binaries/src/mafft/core/nj.c

   1 #include "mltaln.h"
   2 #define DEBUG 0
   3
   4
   5 void topolcpy( int s1[], int s2[], int *mpt1, int *mpt2 )
   6 {
   7     int i;
   8
   9     *mpt1 = *mpt2;
  10     for( i=0; i<*mpt2; i++ )
  11     {
  12         s1[i] = s2[i];
  13     }
  14 }
  15
  16 void topolcat( int s1[], int s2[], int *mpt1, int *mpt2 )
  17 {
  18     int i;
  19
  20     for( i=*mpt1; i<*mpt1+*mpt2; i++ )
  21     {
  22         s1[i] = s2[i-*mpt1];
  23     }
  24     *mpt1 += *mpt2;
  25 }
  26
  27 void topolsort( int m, int s[] )
  28 {
  29     int i, j, im;
  30     int sm;
  31
  32     for( j=0; j<m-1; j++ )
  33     {
  34         sm = s[j]; im = j;
  35         for( i=j+1; i<m; i++ )
  36         {
  37             if( s[i] < sm )
  38             {
  39                 sm = s[i];
  40                 im = i;
  41             }
  42         }
  43         s[im] = s[j]; s[j] = sm;
  44     }
  45 }
  46
  47 void topolswap( int s1[], int s2[], int *mpt1, int *mpt2 )
  48 {
  49     int i;
  50     int im;
  51     int b;
  52     b = *mpt1; *mpt1 = *mpt2; *mpt2 = b;
  53     im = MAX(*mpt1,*mpt2);
  54     for( i=0; i<im; i++ )
  55     {
  56         b = s1[i]; s1[i] = s2[i]; s2[i] = b;
  57     /*
  58     printf( "s1[%d]=%d\ns2[%d]=%d\n", i, s1[i], i, s2[i] );
  59     */
  60     }
  61 }
  62
  63 void reduc( double **mtx, int nseq, int im, int jm )
  64 {
  65     int i;
  66     for( i=0; i<nseq; i++ )
  67     {
  68         if(    i==im || i==jm
  69             || mtx[MIN(i,im)][MAX(i,im)] == 9999.9
  70             || mtx[MIN(i,jm)][MAX(i,jm)] == 9999.9
  71           ) continue;
  72         mtx[MIN(i,im)][MAX(i,im)]
  73         = 0.5 * ( mtx[MIN(i,im)][MAX(i,im)] + mtx[MIN(i,jm)][MAX(i,jm)]
  74                   - mtx[MIN(im,jm)][MAX(im,jm)] );
  75         mtx[MIN(i,jm)][MAX(i,jm)] = 9999.9;
  76     }
  77     mtx[MIN(im,jm)][MAX(im,jm)] = 9999.9;
  78 }
  79
  80
  81 void  nj( int nseq, double **omtx, int ***topol, double **dis )
  82 {
  83     int i, j, l, n, m;
  84     int count;
  85     double r[M];
  86     double t;
  87     double s, sm;
  88     double totallen = 0.0;
  89     int im=0, jm=0;
  90     double len1, len2;
  91 #if 1
  92         static char **par = NULL;
  93         static double **mtx = NULL;
  94         static int **mem = NULL;
  95         if( par == NULL )
  96         {
  97                 par = AllocateCharMtx( njob, njob );
  98                 mtx = AllocateDoubleMtx( njob, njob );
  99                 mem = AllocateIntMtx( njob, 2 );
 100         }
 101 #else
 102     char par[njob][njob];
 103         double mtx[njob][njob];
 104         int mem[njob][2];
 105 #endif
 106         for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ ) mtx[i][j] = omtx[i][j];
 107     for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ ) par[i][j] = 0;
 108     for( i=0; i<nseq; i++ ) par[i][i] = 1;
 109 //      for( i=0; i<nseq; i++ ) for( j=0; j<2; j++ ) for( l=0; l<nseq+1; l++ ) topol[i][j][l] = -1;
 110         for( i=0; i<nseq; i++ ) for( j=0; j<2; j++ ) for( l=0; l<nseq; l++ ) topol[i][j][l] = -1;
 111     for( n=nseq, m=0; n>2; n--, m=nseq-n )
 112     {
 113         t = 0.0;
 114         for( i=0; i<nseq-1; i++ ) for( j=0; j<nseq; j++ ) if( mtx[i][j] < 9999.9 )
 115             t += mtx[i][j];
 116         for( i=0; i<nseq; i++ )
 117         {
 118             r[i] = 0.0;
 119             for( l=0; l<nseq; l++ )
 120                 if( ( l != i ) && ( mtx[MIN(i,l)][MAX(i,l)] < 9999.9 ) )
 121                     r[i] += mtx[MIN(i,l)][MAX(i,l)];
 122         }
 123         sm = 9999.9;
 124         for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ ) if( mtx[i][j] < 9999.9)
 125         {
 126             s = ( ( 2.0 * t - r[i] - r[j] + (n-2.0)*mtx[i][j] ) ) / ( 2.0*(n-2.0) );
 127             if ( s < sm )
 128             {
 129                 sm = s;
 130                 im = i; jm = j;
 131             }
 132         }
 133         len1 = ( (n-2)*mtx[im][jm] + r[im] - r[jm] ) / (2*(n-2));
 134         len2 = ( (n-2)*mtx[im][jm] - r[im] + r[jm] ) / (2*(n-2));
 135
 136 #if DEBUG
 137         fprintf( stderr, "STEP-%3d  %3d: L = %5.5f\n", m+1, im+1, len1 );
 138         fprintf( stderr, "          %3d: L = %5.5f\n",      jm+1, len2 );
 139 #endif
 140
 141         totallen += len1;
 142         totallen += len2;
 143
 144         dis[m][0] = len1;
 145         dis[m][1] = len2;
 146
 147         for( l=0, count=0; l<nseq; l++ )
 148             if( par[im][l] > 0 )
 149             {
 150                 topol[m][0][count] = l;
 151                 count++;
 152             }
 153         mem[m][0] = count;
 154         for( l=0, count=0; l<nseq; l++ )
 155             if( par[jm][l] > 0 )
 156             {
 157                 topol[m][1][count] = l;
 158                 count++;
 159             }
 160         mem[m][1] = count;
 161         for( l=0; l<nseq; l++ )
 162             par[im][l] += ( par[jm][l] > 0 );
 163         if( n > 3 ) reduc( mtx, nseq, im, jm );
 164     }
 165     for( i=0; i<nseq; i++ )
 166         if( i!=im && i!=jm && mtx[MIN(i,im)][MAX(i,im)]<9999.9 )
 167             break;
 168     len2 = ( mtx[MIN(i,im)][MAX(i,im)] - r[im] + r[i] ) / 2;
 169
 170 /*
 171     printf("          %3d: L = %5.5f\n", i+1, len2 );
 172 */
 173     totallen += len2;
 174
 175     dis[m][0] = len2;
 176     dis[m][1] = 0.0;
 177     for( l=0, count=0; l<nseq; l++ )
 178         if( par[i][l] > 0 )
 179         {
 180             topol[m][0][count] = l;
 181             count++;
 182         }
 183     mem[m][0] = count;
 184     /*
 185     printf( " total length == %f\n", totallen );
 186     */
 187
 188     topolcpy( topol[nseq-2][1], topol[nseq-3][0], mem[nseq-2]+1, mem[nseq-3] );
 189     topolcat( topol[nseq-2][1], topol[nseq-3][1], mem[nseq-2]+1, mem[nseq-3]+1 );
 190     topolsort( mem[nseq-2][1], topol[nseq-2][1] );
 191
 192         if( topol[nseq-2][0][0] > topol[nseq-2][1][0] )
 193                 topolswap( topol[nseq-2][0], topol[nseq-2][1], mem[nseq-2], mem[nseq-2]+1 );
 194
 195 }