inprogress

diff --git a/forester/java/src/org/forester/evoinference/TestPhylogenyReconstruction.java b/forester/java/src/org/forester/evoinference/TestPhylogenyReconstruction.java
index 6073985..1963491 100644 (file)
--- a/forester/java/src/org/forester/evoinference/TestPhylogenyReconstruction.java
+++ b/forester/java/src/org/forester/evoinference/TestPhylogenyReconstruction.java
@@ -63,7 +63,8 @@ public class TestPhylogenyReconstruction {
     }
 
     public static void main( final String[] args ) {
-        timeNeighborJoining();
+        testNeighborJoining();
+        // timeNeighborJoining();
     }
 
     public static boolean test( final File test_dir ) {
@@ -1913,6 +1914,19 @@ public class TestPhylogenyReconstruction {
 
     private static boolean testNeighborJoining() {
         try {
+            final NeighborJoining nj = NeighborJoining.createInstance();
+            //            BasicSymmetricalDistanceMatrix m0 = new BasicSymmetricalDistanceMatrix( 4 );
+            //            m0.setIdentifier( 0, "A" );
+            //            m0.setIdentifier( 1, "B" );
+            //            m0.setIdentifier( 2, "C" );
+            //            m0.setIdentifier( 3, "D" );
+            //            m0.setRow( "5 ", 1 );
+            //            m0.setRow( "3 6 ", 2 );
+            //            m0.setRow( "7.5 10.5 5.5", 3 );
+            //            System.out.println( m0.toString() );
+            //            final Phylogeny p0 = nj.execute( m0 );
+            //            Archaeopteryx.createApplication( p0 );
+            //
             BasicSymmetricalDistanceMatrix m = new BasicSymmetricalDistanceMatrix( 6 );
             m.setRow( "5", 1 );
             m.setRow( "4 7", 2 );
@@ -1925,64 +1939,66 @@ public class TestPhylogenyReconstruction {
             m.setIdentifier( 3, "D" );
             m.setIdentifier( 4, "E" );
             m.setIdentifier( 5, "F" );
-            final NeighborJoining nj = NeighborJoining.createInstance();
-            nj.execute( m );
-            m = new BasicSymmetricalDistanceMatrix( 7 );
-            m.setIdentifier( 0, "Bovine" );
-            m.setIdentifier( 1, "Mouse" );
-            m.setIdentifier( 2, "Gibbon" );
-            m.setIdentifier( 3, "Orang" );
-            m.setIdentifier( 4, "Gorilla" );
-            m.setIdentifier( 5, "Chimp" );
-            m.setIdentifier( 6, "Human" );
-            m.setRow( "0.00000 1.68660 1.71980 1.66060 1.52430 1.60430 1.59050", 0 );
-            m.setRow( "1.68660 0.00000 1.52320 1.48410 1.44650 1.43890 1.46290", 1 );
-            m.setRow( "1.71980 1.52320 0.00000 0.71150 0.59580 0.61790 0.55830", 2 );
-            m.setRow( "1.66060 1.48410 0.71150 0.00000 0.46310 0.50610 0.47100", 3 );
-            m.setRow( "1.52430 1.44650 0.59580 0.46310 0.00000 0.34840 0.30830", 4 );
-            m.setRow( "1.60430 1.43890 0.61790 0.50610 0.34840 0.00000 0.26920", 5 );
-            m.setRow( "1.59050 1.46290 0.55830 0.47100 0.30830 0.26920 0.00000", 6 );
             System.out.println( m.toString() );
-            final Phylogeny p2 = nj.execute( m );
-            p2.reRoot( p2.getNode( "Bovine" ) );
-            System.out.println( p2.toString() );
-            // from phylip Neighbor-Joining/UPGMA method version 3.69:
-            // ((((((Chimp:0.15167,Human:0.11753):0.03982,Gorilla:0.15393):0.02696,Orang:0.28469):0.04648,Gibbon:0.35793):0.42027,Mouse:0.76891):0.458845,Bovine:0.458845);
-            Archaeopteryx.createApplication( p2 );
-            m = new BasicSymmetricalDistanceMatrix( 4 );
-            m.setIdentifier( 0, "A" );
-            m.setIdentifier( 1, "B" );
-            m.setIdentifier( 2, "C" );
-            m.setIdentifier( 3, "D" );
-            m.setRow( "0.00 0.95 0.17 0.98", 0 );
-            m.setRow( "0.95 0.00 1.02 1.83", 1 );
-            m.setRow( "0.17 1.02 0.00 1.01", 2 );
-            m.setRow( "0.98 1.83 1.01 0.00", 3 );
-            final Phylogeny p3 = nj.execute( m );
-            //
-            // -- A 0.05
-            // - |0.01
-            // ----------------------- B 0.90
-            //
-            // --- C 0.10
-            // - |0.01
-            // ------------------------- D 0.91
-            p3.reRoot( p3.getNode( "C" ).getParent() );
-            if ( !isEqual( p3.getNode( "A" ).getDistanceToParent(), 0.05 ) ) {
-                return false;
-            }
-            if ( !isEqual( p3.getNode( "B" ).getDistanceToParent(), 0.90 ) ) {
-                return false;
-            }
-            if ( !isEqual( p3.getNode( "C" ).getDistanceToParent(), 0.10 ) ) {
-                return false;
-            }
-            if ( !isEqual( p3.getNode( "D" ).getDistanceToParent(), 0.91 ) ) {
-                return false;
-            }
-            if ( TIME ) {
-                timeNeighborJoining();
-            }
+            final Phylogeny p1 = nj.execute( m );
+            Archaeopteryx.createApplication( p1 );
+            //            m = new BasicSymmetricalDistanceMatrix( 7 );
+            //            m.setIdentifier( 0, "Bovine" );
+            //            m.setIdentifier( 1, "Mouse" );
+            //            m.setIdentifier( 2, "Gibbon" );
+            //            m.setIdentifier( 3, "Orang" );
+            //            m.setIdentifier( 4, "Gorilla" );
+            //            m.setIdentifier( 5, "Chimp" );
+            //            m.setIdentifier( 6, "Human" );
+            //            m.setRow( "0.00000 1.68660 1.71980 1.66060 1.52430 1.60430 1.59050", 0 );
+            //            m.setRow( "1.68660 0.00000 1.52320 1.48410 1.44650 1.43890 1.46290", 1 );
+            //            m.setRow( "1.71980 1.52320 0.00000 0.71150 0.59580 0.61790 0.55830", 2 );
+            //            m.setRow( "1.66060 1.48410 0.71150 0.00000 0.46310 0.50610 0.47100", 3 );
+            //            m.setRow( "1.52430 1.44650 0.59580 0.46310 0.00000 0.34840 0.30830", 4 );
+            //            m.setRow( "1.60430 1.43890 0.61790 0.50610 0.34840 0.00000 0.26920", 5 );
+            //            m.setRow( "1.59050 1.46290 0.55830 0.47100 0.30830 0.26920 0.00000", 6 );
+            //            System.out.println( m.toString() );
+            //            final Phylogeny p2 = nj.execute( m );
+            //            p2.reRoot( p2.getNode( "Bovine" ) );
+            //            System.out.println( p2.toString() );
+            //            Archaeopteryx.createApplication( p2 );
+            //            // from phylip Neighbor-Joining/UPGMA method version 3.69:
+            //            // ((((((Chimp:0.15167,Human:0.11753):0.03982,Gorilla:0.15393):0.02696,Orang:0.28469):0.04648,Gibbon:0.35793):0.42027,Mouse:0.76891):0.458845,Bovine:0.458845);
+            //            Archaeopteryx.createApplication( p2 );
+            //            m = new BasicSymmetricalDistanceMatrix( 4 );
+            //            m.setIdentifier( 0, "A" );
+            //            m.setIdentifier( 1, "B" );
+            //            m.setIdentifier( 2, "C" );
+            //            m.setIdentifier( 3, "D" );
+            //            m.setRow( "0.00 0.95 0.17 0.98", 0 );
+            //            m.setRow( "0.95 0.00 1.02 1.83", 1 );
+            //            m.setRow( "0.17 1.02 0.00 1.01", 2 );
+            //            m.setRow( "0.98 1.83 1.01 0.00", 3 );
+            //            final Phylogeny p3 = nj.execute( m );
+            //            //
+            //            // -- A 0.05
+            //            // - |0.01
+            //            // ----------------------- B 0.90
+            //            //
+            //            // --- C 0.10
+            //            // - |0.01
+            //            // ------------------------- D 0.91
+            //            p3.reRoot( p3.getNode( "C" ).getParent() );
+            //            if ( !isEqual( p3.getNode( "A" ).getDistanceToParent(), 0.05 ) ) {
+            //                return false;
+            //            }
+            //            if ( !isEqual( p3.getNode( "B" ).getDistanceToParent(), 0.90 ) ) {
+            //                return false;
+            //            }
+            //            if ( !isEqual( p3.getNode( "C" ).getDistanceToParent(), 0.10 ) ) {
+            //                return false;
+            //            }
+            //            if ( !isEqual( p3.getNode( "D" ).getDistanceToParent(), 0.91 ) ) {
+            //                return false;
+            //            }
+            //            if ( TIME ) {
+            //                timeNeighborJoining();
+            //            }
         }
         catch ( final Exception e ) {
             e.printStackTrace( System.out );
@@ -2353,7 +2369,7 @@ public class TestPhylogenyReconstruction {
 
     private static void timeNeighborJoining() {
         final NeighborJoining nj = NeighborJoining.createInstance();
-        for( int n = 3; n <= 13; ++n ) {
+        for( int n = 3; n <= 6; ++n ) {
             final int x = ( int ) Math.pow( 2, n );
             final BasicSymmetricalDistanceMatrix mt = new BasicSymmetricalDistanceMatrix( x );
             mt.randomize( new Date().getTime() );

diff --git a/forester/java/src/org/forester/evoinference/distance/NeighborJoining.java b/forester/java/src/org/forester/evoinference/distance/NeighborJoining.java
index bcb6df1..2dc0ffc 100644 (file)
--- a/forester/java/src/org/forester/evoinference/distance/NeighborJoining.java
+++ b/forester/java/src/org/forester/evoinference/distance/NeighborJoining.java
@@ -44,12 +44,34 @@ public final class NeighborJoining {
     private PhylogenyNode[]                _external_nodes;
     private int[]                          _mappings;
     private final boolean                  _verbose;
-    private final static boolean           DEBUG = false;
+    private final static boolean           DEBUG = true;
 
     private NeighborJoining( final boolean verbose ) {
         _verbose = verbose;
     }
 
+    private final void printM() {
+        for( int i = 0; i < _m_values.length; i++ ) {
+            for( int j = 0; j < _m_values.length; j++ ) {
+                System.out.print( _m_values[ i ][ j ] );
+                System.out.print( " " );
+            }
+            System.out.println();
+        }
+        System.out.println();
+    }
+
+    private final void printD() {
+        for( int i = 0; i < _d_values.length; i++ ) {
+            for( int j = 0; j < _d_values.length; j++ ) {
+                System.out.print( _d_values[ i ][ j ] );
+                System.out.print( " " );
+            }
+            System.out.println();
+        }
+        System.out.println();
+    }
+
     private final void calculateDistancesFromNewNode( final int otu1, final int otu2, final double d ) {
         final int otu1_m = _mappings[ otu1 ];
         final int otu2_m = _mappings[ otu2 ];
@@ -66,15 +88,35 @@ public final class NeighborJoining {
 
     private final void calculateNetDivergences() {
         double d;
-        int i_m;
+        //  int i_m;
         for( int i = 0; i < _n; ++i ) {
             d = 0;
-            i_m = _mappings[ i ];
+            //i_m = _mappings[ i ];
             for( int n = 0; n < _n; ++n ) {
                 //d += _d_values[ i_m ][ _mappings[ n ] ];
-                d += getValueFromD( i, n );
+                if ( i != n ) {
+                    if ( i < n ) {
+                        d += getValueFromD( i, n );
+                        System.out.print( "+" );
+                        System.out.print( getValueFromD( i, n ) );
+                    }
+                    else {
+                        d += getValueFromD( n, i );
+                        System.out.print( "+" );
+                        System.out.print( getValueFromD( n, i ) );
+                    }
+                }
+                else {
+                    if ( DEBUG ) {
+                        if ( getValueFromD( i, n ) != 0 ) {
+                            throw new RuntimeException( "faulty NJ code" );
+                        }
+                    }
+                }
             }
             _r[ i ] = d;
+            System.out.print( "=" );
+            System.out.println( d );
         }
     }
 
@@ -147,6 +189,10 @@ public final class NeighborJoining {
     }
 
     private final double getValueFromD( final int otu1, final int otu2 ) {
+        if ( otu1 > otu2 ) {
+            //throw new IllegalStateException();
+            return _d_values[ _mappings[ otu2 ] ][ _mappings[ otu1 ] ];
+        }
         return _d_values[ _mappings[ otu1 ] ][ _mappings[ otu2 ] ];
     }
 
@@ -191,13 +237,15 @@ public final class NeighborJoining {
         int j_m;
         final int _n_2 = _n - 2;
         for( int j = 1; j < _n; ++j ) {
-            r_j = _r[ j ];
-            j_m = _mappings[ j ];
+            //r_j = _r[ j ];
+            // j_m = _mappings[ j ];
             for( int i = 0; i < j; ++i ) {
-                _m_values[ i ][ j ] = getValueFromD( i, j ) - ( _r[ i ] + r_j ) / ( _n - 2 );
+                _m_values[ i ][ j ] = getValueFromD( i, j ) - ( ( _r[ i ] + _r[ j ] ) / ( _n - 2 ) );
                 //_m_values[ i ][ j ] = _d_values[ _mappings[ i ] ][ j_m ] - ( ( _r[ i ] + r_j ) / ( _n_2 ) );
             }
         }
+        printM();
+        printD();
     }
 
     //  private final double getValueFromD( final int otu1, final int otu2 ) {