Mac binaries

[jabaws.git] / website / archive / binaries / mac / src / mafft / core / nj.c

diff --git a/website/archive/binaries/mac/src/mafft/core/nj.c b/website/archive/binaries/mac/src/mafft/core/nj.c
new file mode 100644 (file)
index 0000000..0adf4e5
--- /dev/null
+++ b/website/archive/binaries/mac/src/mafft/core/nj.c
@@ -0,0 +1,195 @@
+#include "mltaln.h"
+#define DEBUG 0
+
+
+void topolcpy( int s1[], int s2[], int *mpt1, int *mpt2 )
+{
+    int i;
+
+    *mpt1 = *mpt2;
+    for( i=0; i<*mpt2; i++ )
+    {
+        s1[i] = s2[i];
+    }
+}
+
+void topolcat( int s1[], int s2[], int *mpt1, int *mpt2 )
+{
+    int i;
+
+    for( i=*mpt1; i<*mpt1+*mpt2; i++ )
+    {
+        s1[i] = s2[i-*mpt1];
+    }
+    *mpt1 += *mpt2;
+}
+   
+void topolsort( int m, int s[] )
+{
+    int i, j, im;
+    int sm;
+
+    for( j=0; j<m-1; j++ )
+    {
+        sm = s[j]; im = j;
+        for( i=j+1; i<m; i++ )
+        {
+            if( s[i] < sm )
+            {
+                sm = s[i];
+                im = i;
+            }
+        }
+        s[im] = s[j]; s[j] = sm;
+    }
+}
+
+void topolswap( int s1[], int s2[], int *mpt1, int *mpt2 )
+{
+    int i;
+    int im;
+    int b;
+    b = *mpt1; *mpt1 = *mpt2; *mpt2 = b;
+    im = MAX(*mpt1,*mpt2);
+    for( i=0; i<im; i++ )
+    {
+        b = s1[i]; s1[i] = s2[i]; s2[i] = b;
+    /*
+    printf( "s1[%d]=%d\ns2[%d]=%d\n", i, s1[i], i, s2[i] );
+    */
+    }
+}
+
+void reduc( double **mtx, int nseq, int im, int jm )
+{
+    int i;
+    for( i=0; i<nseq; i++ )
+    {
+        if(    i==im || i==jm
+            || mtx[MIN(i,im)][MAX(i,im)] == 9999.9
+            || mtx[MIN(i,jm)][MAX(i,jm)] == 9999.9
+          ) continue;
+        mtx[MIN(i,im)][MAX(i,im)]
+        = 0.5 * ( mtx[MIN(i,im)][MAX(i,im)] + mtx[MIN(i,jm)][MAX(i,jm)]
+                  - mtx[MIN(im,jm)][MAX(im,jm)] );
+        mtx[MIN(i,jm)][MAX(i,jm)] = 9999.9;
+    }
+    mtx[MIN(im,jm)][MAX(im,jm)] = 9999.9;
+}
+
+
+void  nj( int nseq, double **omtx, int ***topol, double **dis )
+{
+    int i, j, l, n, m;
+    int count;
+    double r[M];
+    double t;
+    double s, sm;
+    double totallen = 0.0;
+    int im=0, jm=0;
+    double len1, len2;
+#if 1
+       static char **par = NULL;
+       static double **mtx = NULL;
+       static int **mem = NULL;
+       if( par == NULL )
+       {
+               par = AllocateCharMtx( njob, njob );
+               mtx = AllocateDoubleMtx( njob, njob );
+               mem = AllocateIntMtx( njob, 2 );
+       }
+#else
+    char par[njob][njob];
+       double mtx[njob][njob];
+       int mem[njob][2];
+#endif
+       for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ ) mtx[i][j] = omtx[i][j];
+    for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ ) par[i][j] = 0;
+    for( i=0; i<nseq; i++ ) par[i][i] = 1;
+//     for( i=0; i<nseq; i++ ) for( j=0; j<2; j++ ) for( l=0; l<nseq+1; l++ ) topol[i][j][l] = -1;
+       for( i=0; i<nseq; i++ ) for( j=0; j<2; j++ ) for( l=0; l<nseq; l++ ) topol[i][j][l] = -1;
+    for( n=nseq, m=0; n>2; n--, m=nseq-n )
+    {
+        t = 0.0;
+        for( i=0; i<nseq-1; i++ ) for( j=0; j<nseq; j++ ) if( mtx[i][j] < 9999.9 )
+            t += mtx[i][j];
+        for( i=0; i<nseq; i++ )
+        {
+            r[i] = 0.0;
+            for( l=0; l<nseq; l++ ) 
+                if( ( l != i ) && ( mtx[MIN(i,l)][MAX(i,l)] < 9999.9 ) )
+                    r[i] += mtx[MIN(i,l)][MAX(i,l)];
+        }
+        sm = 9999.9;
+        for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ ) if( mtx[i][j] < 9999.9)
+        {
+            s = ( ( 2.0 * t - r[i] - r[j] + (n-2.0)*mtx[i][j] ) ) / ( 2.0*(n-2.0) );
+            if ( s < sm )
+            {
+                sm = s;
+                im = i; jm = j;
+            }
+        }
+        len1 = ( (n-2)*mtx[im][jm] + r[im] - r[jm] ) / (2*(n-2));
+        len2 = ( (n-2)*mtx[im][jm] - r[im] + r[jm] ) / (2*(n-2));
+
+#if DEBUG
+        fprintf( stderr, "STEP-%3d  %3d: L = %5.5f\n", m+1, im+1, len1 );
+        fprintf( stderr, "          %3d: L = %5.5f\n",      jm+1, len2 );
+#endif
+
+        totallen += len1;
+        totallen += len2;
+
+        dis[m][0] = len1;
+        dis[m][1] = len2;
+
+        for( l=0, count=0; l<nseq; l++ )
+            if( par[im][l] > 0 )
+            {
+                topol[m][0][count] = l;
+                count++;
+            }
+        mem[m][0] = count;
+        for( l=0, count=0; l<nseq; l++ )
+            if( par[jm][l] > 0 )
+            {
+                topol[m][1][count] = l;
+                count++;
+            }
+        mem[m][1] = count;
+        for( l=0; l<nseq; l++ )
+            par[im][l] += ( par[jm][l] > 0 );
+        if( n > 3 ) reduc( mtx, nseq, im, jm );
+    }
+    for( i=0; i<nseq; i++ )
+        if( i!=im && i!=jm && mtx[MIN(i,im)][MAX(i,im)]<9999.9 )
+            break;
+    len2 = ( mtx[MIN(i,im)][MAX(i,im)] - r[im] + r[i] ) / 2;
+
+/*
+    printf("          %3d: L = %5.5f\n", i+1, len2 );
+*/
+    totallen += len2;
+
+    dis[m][0] = len2;
+    dis[m][1] = 0.0;
+    for( l=0, count=0; l<nseq; l++ )
+        if( par[i][l] > 0 )
+        {
+            topol[m][0][count] = l;
+            count++;
+        }
+    mem[m][0] = count;
+    /*
+    printf( " total length == %f\n", totallen );
+    */
+
+    topolcpy( topol[nseq-2][1], topol[nseq-3][0], mem[nseq-2]+1, mem[nseq-3] );
+    topolcat( topol[nseq-2][1], topol[nseq-3][1], mem[nseq-2]+1, mem[nseq-3]+1 );
+    topolsort( mem[nseq-2][1], topol[nseq-2][1] );
+       
+       if( topol[nseq-2][0][0] > topol[nseq-2][1][0] )
+               topolswap( topol[nseq-2][0], topol[nseq-2][1], mem[nseq-2], mem[nseq-2]+1 );
+
+}