5 void topolcpy( int s1[], int s2[], int *mpt1, int *mpt2 )
10 for( i=0; i<*mpt2; i++ )
16 void topolcat( int s1[], int s2[], int *mpt1, int *mpt2 )
20 for( i=*mpt1; i<*mpt1+*mpt2; i++ )
27 void topolsort( int m, int s[] )
32 for( j=0; j<m-1; j++ )
35 for( i=j+1; i<m; i++ )
43 s[im] = s[j]; s[j] = sm;
47 void topolswap( int s1[], int s2[], int *mpt1, int *mpt2 )
52 b = *mpt1; *mpt1 = *mpt2; *mpt2 = b;
53 im = MAX(*mpt1,*mpt2);
56 b = s1[i]; s1[i] = s2[i]; s2[i] = b;
58 printf( "s1[%d]=%d\ns2[%d]=%d\n", i, s1[i], i, s2[i] );
63 void reduc( double **mtx, int nseq, int im, int jm )
66 for( i=0; i<nseq; i++ )
69 || mtx[MIN(i,im)][MAX(i,im)] == 9999.9
70 || mtx[MIN(i,jm)][MAX(i,jm)] == 9999.9
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;
77 mtx[MIN(im,jm)][MAX(im,jm)] = 9999.9;
81 void nj( int nseq, double **omtx, int ***topol, double **dis )
88 double totallen = 0.0;
92 static char **par = NULL;
93 static double **mtx = NULL;
94 static int **mem = NULL;
97 par = AllocateCharMtx( njob, njob );
98 mtx = AllocateDoubleMtx( njob, njob );
99 mem = AllocateIntMtx( njob, 2 );
102 char par[njob][njob];
103 double mtx[njob][njob];
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 )
114 for( i=0; i<nseq-1; i++ ) for( j=0; j<nseq; j++ ) if( mtx[i][j] < 9999.9 )
116 for( i=0; i<nseq; i++ )
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)];
124 for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ ) if( mtx[i][j] < 9999.9)
126 s = ( ( 2.0 * t - r[i] - r[j] + (n-2.0)*mtx[i][j] ) ) / ( 2.0*(n-2.0) );
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));
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 );
147 for( l=0, count=0; l<nseq; l++ )
150 topol[m][0][count] = l;
154 for( l=0, count=0; l<nseq; l++ )
157 topol[m][1][count] = l;
161 for( l=0; l<nseq; l++ )
162 par[im][l] += ( par[jm][l] > 0 );
163 if( n > 3 ) reduc( mtx, nseq, im, jm );
165 for( i=0; i<nseq; i++ )
166 if( i!=im && i!=jm && mtx[MIN(i,im)][MAX(i,im)]<9999.9 )
168 len2 = ( mtx[MIN(i,im)][MAX(i,im)] - r[im] + r[i] ) / 2;
171 printf(" %3d: L = %5.5f\n", i+1, len2 );
177 for( l=0, count=0; l<nseq; l++ )
180 topol[m][0][count] = l;
185 printf( " total length == %f\n", totallen );
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] );
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 );