forester/java/src/org/forester/evoinference/distance/NeighborJoining.java

   1 // $Id:
   2 // FORESTER -- software libraries and applications
   3 // for evolutionary biology research and applications.
   4 //
   5 // Copyright (C) 2008-2009 Christian M. Zmasek
   6 // Copyright (C) 2008-2009 Burnham Institute for Medical Research
   7 // All rights reserved
   8 //
   9 // This library is free software; you can redistribute it and/or
  10 // modify it under the terms of the GNU Lesser General Public
  11 // License as published by the Free Software Foundation; either
  12 // version 2.1 of the License, or (at your option) any later version.
  13 //
  14 // This library is distributed in the hope that it will be useful,
  15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  17 // Lesser General Public License for more details.
  18 //
  19 // You should have received a copy of the GNU Lesser General Public
  20 // License along with this library; if not, write to the Free Software
  21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  22 //
  23 // Contact: phylosoft @ gmail . com
  24 // WWW: www.phylosoft.org/forester
  25
  26 package org.forester.evoinference.distance;
  27
  28 import java.util.ArrayList;
  29 import java.util.List;
  30
  31 import org.forester.evoinference.matrix.distance.BasicSymmetricalDistanceMatrix;
  32 import org.forester.evoinference.matrix.distance.DistanceMatrix;
  33 import org.forester.phylogeny.Phylogeny;
  34 import org.forester.phylogeny.PhylogenyNode;
  35 import org.forester.phylogeny.PhylogenyNodeI;
  36 import org.forester.util.ForesterUtil;
  37
  38 public class NeighborJoining {
  39
  40     private final static boolean VERBOSE_DEFAULT = false;
  41     private DistanceMatrix       _d;
  42     private DistanceMatrix       _m;
  43     private double[]             _r;
  44     private int                  _n;
  45     private PhylogenyNodeI[]     _external_nodes;
  46     private int[]                _mappings;
  47     private boolean              _verbose;
  48
  49     public NeighborJoining() {
  50         init();
  51     }
  52
  53     private void calculateDistancesFromNewNode( final int otu1, final int otu2, final double d ) {
  54         for( int i = 0; i < _n; ++i ) {
  55             if ( ( i == otu1 ) || ( i == otu2 ) ) {
  56                 continue;
  57             }
  58             final double nd = ( getValueFromD( otu1, i ) + getValueFromD( i, otu2 ) - d ) / 2;
  59             setValueInD( nd, otu1, i );
  60         }
  61     }
  62
  63     private double calculateM( final int i, final int j ) {
  64         return getValueFromD( i, j ) - ( _r[ i ] + _r[ j ] ) / ( _n - 2 );
  65     }
  66
  67     private void calculateNetDivergences() {
  68         for( int i = 0; i < _n; ++i ) {
  69             double d = 0.0;
  70             for( int n = 0; n < _n; ++n ) {
  71                 d += getValueFromD( i, n );
  72             }
  73             _r[ i ] = d;
  74         }
  75     }
  76
  77     public Phylogeny execute( final DistanceMatrix distance ) {
  78         reset( distance );
  79         final Phylogeny phylogeny = new Phylogeny();
  80         while ( _n > 2 ) {
  81             updateM();
  82             final int[] s = findMinimalDistance();
  83             final int otu1 = s[ 0 ];
  84             final int otu2 = s[ 1 ];
  85             // It is a condition that otu1 < otu2.
  86             if ( otu1 > otu2 ) {
  87                 throw new AssertionError( "NJ code is faulty: otu1 > otu2" );
  88             }
  89             final PhylogenyNodeI node = new PhylogenyNode();
  90             final double d = getValueFromD( otu1, otu2 );
  91             final double d1 = ( d / 2 ) + ( ( _r[ otu1 ] - _r[ otu2 ] ) / ( 2 * ( _n - 2 ) ) );
  92             final double d2 = d - d1;
  93             getExternalPhylogenyNode( otu1 ).setDistanceToParent( d1 );
  94             getExternalPhylogenyNode( otu2 ).setDistanceToParent( d2 );
  95             node.addAsChild( getExternalPhylogenyNode( otu1 ) );
  96             node.addAsChild( getExternalPhylogenyNode( otu2 ) );
  97             if ( isVerbose() ) {
  98                 printProgress( otu1, otu2 );
  99             }
 100             calculateDistancesFromNewNode( otu1, otu2, d );
 101             setExternalPhylogenyNode( node, otu1 );
 102             updateMappings( otu2 );
 103             --_n;
 104         }
 105         final double d = getValueFromD( 0, 1 ) / 2;
 106         getExternalPhylogenyNode( 0 ).setDistanceToParent( d );
 107         getExternalPhylogenyNode( 1 ).setDistanceToParent( d );
 108         final PhylogenyNodeI root = new PhylogenyNode();
 109         root.addAsChild( getExternalPhylogenyNode( 0 ) );
 110         root.addAsChild( getExternalPhylogenyNode( 1 ) );
 111         if ( isVerbose() ) {
 112             printProgress( 0, 1 );
 113         }
 114         phylogeny.setRoot( ( PhylogenyNode ) root );
 115         phylogeny.setRooted( false );
 116         return phylogeny;
 117     }
 118
 119     public List<Phylogeny> execute( final List<DistanceMatrix> distances_list ) {
 120         final List<Phylogeny> pl = new ArrayList<Phylogeny>();
 121         for( final DistanceMatrix distances : distances_list ) {
 122             pl.add( execute( distances ) );
 123         }
 124         return pl;
 125     }
 126
 127     private int[] findMinimalDistance() {
 128         // if more than one minimal distances, always the first found is
 129         // returned
 130         // i could randomize this, so that any would be returned in a randomized
 131         // fashion...
 132         double minimum = Double.MAX_VALUE;
 133         int otu_1 = -1;
 134         int otu_2 = -1;
 135         for( int j = 1; j < _n; ++j ) {
 136             for( int i = 0; i < j; ++i ) {
 137                 if ( _m.getValue( i, j ) < minimum ) {
 138                     minimum = _m.getValue( i, j );
 139                     otu_1 = i;
 140                     otu_2 = j;
 141                 }
 142             }
 143         }
 144         return new int[] { otu_1, otu_2 };
 145     }
 146
 147     private PhylogenyNodeI getExternalPhylogenyNode( final int i ) {
 148         return _external_nodes[ _mappings[ i ] ];
 149     }
 150
 151     private double getValueFromD( final int otu1, final int otu2 ) {
 152         return _d.getValue( _mappings[ otu1 ], _mappings[ otu2 ] );
 153     }
 154
 155     private void init() {
 156         setVerbose( VERBOSE_DEFAULT );
 157     }
 158
 159     private void initExternalNodes() {
 160         _external_nodes = new PhylogenyNodeI[ _n ];
 161         for( int i = 0; i < _n; ++i ) {
 162             _external_nodes[ i ] = new PhylogenyNode();
 163             final String id = _d.getIdentifier( i );
 164             if ( id != null ) {
 165                 _external_nodes[ i ].setName( id );
 166             }
 167             else {
 168                 _external_nodes[ i ].setName( "" + i );
 169             }
 170             _mappings[ i ] = i;
 171         }
 172     }
 173
 174     private boolean isVerbose() {
 175         return _verbose;
 176     }
 177
 178     private void printProgress( final int otu1, final int otu2 ) {
 179         final PhylogenyNodeI n1 = getExternalPhylogenyNode( otu1 );
 180         final PhylogenyNodeI n2 = getExternalPhylogenyNode( otu2 );
 181         System.out.println( "Node " + ( ForesterUtil.isEmpty( n1.getName() ) ? n1.getId() : n1.getName() ) + " joins "
 182                 + ( ForesterUtil.isEmpty( n2.getName() ) ? n2.getId() : n2.getName() ) );
 183     }
 184
 185     // only the values in the lower triangle are used.
 186     // !matrix values will be changed!
 187     private void reset( final DistanceMatrix distances ) {
 188         _n = distances.getSize();
 189         _d = distances;
 190         _m = new BasicSymmetricalDistanceMatrix( _n );
 191         _r = new double[ _n ];
 192         _mappings = new int[ _n ];
 193         initExternalNodes();
 194     }
 195
 196     private void setExternalPhylogenyNode( final PhylogenyNodeI node, final int i ) {
 197         _external_nodes[ _mappings[ i ] ] = node;
 198     }
 199
 200     private void setValueInD( final double d, final int otu1, final int otu2 ) {
 201         _d.setValue( _mappings[ otu1 ], _mappings[ otu2 ], d );
 202     }
 203
 204     public void setVerbose( final boolean verbose ) {
 205         _verbose = verbose;
 206     }
 207
 208     private void updateM() {
 209         calculateNetDivergences();
 210         for( int j = 1; j < _n; ++j ) {
 211             for( int i = 0; i < j; ++i ) {
 212                 _m.setValue( i, j, calculateM( i, j ) );
 213             }
 214         }
 215     }
 216
 217     // otu2 will, in effect, be "deleted" from the matrix.
 218     private void updateMappings( final int otu2 ) {
 219         for( int i = otu2; i < _mappings.length - 1; ++i ) {
 220             _mappings[ i ] = _mappings[ i + 1 ];
 221         }
 222     }
 223
 224     public static NeighborJoining createInstance() {
 225         return new NeighborJoining();
 226     }
 227 }