6 void tracking( char **nseq, char **aseq, char *seq, int **ijp, int icyc )
9 int iin, ifi, jin, jfi;
10 int lgth = strlen( aseq[0] );
11 int lgth1 = strlen( seq );
14 for( i=0; i<=icyc+1; i++ )
16 nseq[i] += lgth+lgth1;
20 iin = lgth; jin = lgth1;
21 for( k=0; k<=lgth+lgth1; )
23 if( ijp[iin][jin] < 0 )
25 ifi = iin-1; jfi = jin+ijp[iin][jin];
27 else if( ijp[iin][jin] > 0 )
29 ifi = iin-ijp[iin][jin]; jfi = jin-1;
33 ifi = iin-1; jfi = jin-1;
35 if( ifi < 0 ) jfi = -1;
36 if( jfi < 0 ) ifi = -1;
39 for( l=1; l<iin-ifi; l++ )
41 for( i=0; i<icyc+1; i++ )
42 *--nseq[i] = aseq[i][iin-l];
43 *--nseq[icyc+1] = *gap;
46 for( l=1; l<jin-jfi; l++ )
48 for( i=0; i<icyc+1; i++ )
51 *nseq[icyc+1] = seq[jin-l];
54 if( iin <= 0 || jin <= 0 ) break;
55 for( i=0; i<icyc+1; i++ )
56 *--nseq[i] = aseq[i][ifi];
57 *--nseq[icyc+1] = seq[jfi];
64 char **Calignm1( float *wm, char **aseq, char *seq, double *effarr, int icyc, int ex )
66 register int i, j, k, l;
74 static float **v, **w;
81 static int orlgth = 0, orlgth1 = 0;
84 float *AllocateFloatVec( int );
88 float **AllocateFloatTri( int );
89 void FreeFloatTri( float ** );
97 for( i=0; i<icyc+1; i++ ) totaleff += effarr[i];
98 lgth = strlen( aseq[0] );
99 lgth1 = strlen( seq );
103 fprintf( stdout, "effarr in Calignm1 s = \n", ex );
104 for( i=0; i<icyc+1; i++ )
105 fprintf( stdout, " %f", effarr[i] );
109 if( orlgth > 0 && orlgth1 > 0 ) for( i=0; i<njob+1; i++ ) nseq[i] = mseq[i];
114 if( lgth > orlgth || lgth1 > orlgth1 )
118 if( orlgth > 0 && orlgth1 > 0 )
134 FreeFloatVec( gvsa );
135 FreeFloatMtx( scmx );
138 ll1 = MAX( (int)(1.1*lgth), orlgth );
139 ll2 = MAX( (int)(1.1*lgth1), orlgth1 );
140 ll1 = MAX( ll1, ll2 );
143 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
146 v = AllocateFloatMtx( ll1+2, ll2+2 );
148 g = AllocateFloatCub( ll1+2, 3, 3 );
150 gl = AllocateFloatTri( MAX( ll1, ll2 ) + 3 );
151 gs = AllocateFloatTri( MAX( ll1, ll2 ) + 3 );
153 m = AllocateFloatVec( ll1+2 );
154 mp = AllocateIntVec( ll1+2 );
156 mseq = AllocateCharMtx( njob+1, 1 );
157 nseq = AllocateCharMtx( njob+1, ll1+ll2 );
158 for( i=0; i<njob+1; i++ ) mseq[i] = nseq[i];
160 gvsa = AllocateFloatVec( ll1+2 );
162 scmx = AllocateFloatMtx( 26, ll1+2 );
165 fprintf( stderr, "succeeded\n\n" );
171 fprintf( stderr, "orlgth == %d\n", orlgth );
175 if( orlgth > commonAlloc1 || orlgth1 > commonAlloc2 )
179 if( commonAlloc1 && commonAlloc2 )
181 FreeIntMtx( commonIP );
184 ll1 = MAX( commonAlloc1, orlgth );
185 ll2 = MAX( commonAlloc2, orlgth1 );
188 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
191 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
194 fprintf( stderr, "succeeded\n\n" );
203 scmx_calc( icyc, aseq, effarr, scmx );
205 for( i=0; i<lgth; i++ )
208 for( l=0; l<26; l++ )
211 for( k=0; k<26; k++ )
212 scarr[l] += n_dis[k][l] * scmx[k][i];
214 for( j=0; j<lgth1; j++ )
216 v[i][j] = scarr[(int)amino_n[(int)seq[j]]];
220 for( i=0; i<lgth+1; i++ ) v[i][lgth1] = 0.0;
221 for( j=0; j<lgth1+1; j++ ) v[lgth][j] = 0.0;
224 for( j=0; j<lgth+1; j++ ) for( k=0; k<3; k++ ) for( l=0; l<3; l++ )
228 for( i=0; i<icyc+1; i++ )
230 efficient = (float)effarr[i];
233 for( j=0; j<lgth+1; j++ )
237 gc1 = ( aseq[i][j] == '-' );
238 if( j == lgth ) gc1 = 0; /* 0? */
241 g[j][0][0] += ( !gb1 * gc1 ) * efficient * penalty;
242 g[j][1][0] += ( 0 ) * efficient * penalty;
243 g[j][2][0] += ( 0 ) * efficient * penalty;
244 g[j][0][1] += ( gb1 * !gc1 + !gb1 * !gc1 ) * efficient * penalty;
245 g[j][1][1] += ( 0 ) * efficient * penalty;
246 g[j][0][2] += ( !gb1 ) * efficient * penalty; /* ??? */
247 g[j][2][2] += ( 0 ) * efficient * penalty;
251 for( i=0; i<lgth+2; i++ ) for( j=0; j<=i+1; j++ )
253 gs[i][j] = gl[i][j] = 0.0;
256 for( i=0; i<icyc+1; i++ )
258 efficient = (float)effarr[i];
260 for( j=0; j<lgth+1; j++ )
262 if( aseq[i][j] == '-' )
265 gs[j][inttmp] += (float)efficient * penalty;
267 if( aseq[i][j+1] != '-' )
268 gl[j][inttmp] += (float)efficient * penalty;
276 for( i=0; i<lgth+1; i++ )
278 fprintf( stderr, "gl[%d][] = ", i );
279 fprintf( stderr, "\n" );
280 for( j=0; j<=i+1; j++ )
281 fprintf( stderr, " %.0f", i, j, gl[i][j] );
282 fprintf( stderr, "\n" );
286 for( i=0; i<lgth+1; i++ )
288 for( j=1; j<=i; j++ )
289 gs[i][j+1] += gs[i][j];
291 gl[i][j] += gl[i][j+1];
295 for( i=0; i<lgth+1; i++ )
298 fprintf( stderr, "gs[%d][%d] = %f, gl[%d][%d] = %f\n", i, j, gs[i][j], i, j, gl[i][j] );
305 printf( "%d\n", lgth );
306 for( i=0; i<lgth; i++ )
308 printf( "%f %f %f %f %f %f %f\n", g[i][0][0], g[i][1][0], g[i][2][0], g[i][0][1], g[i][1][1], g[i][0][2], g[i][2][2] );
310 printf( "\n\n\n\n\n" );
313 for( i=0; i<lgth; i++ ) for( j=0; j<3; j++ ) for( k=0; k<3; k++ )
315 g[i][j][k] /= ( (float)icyc + 1.0 );
320 fp = fopen( "gapp", "a" );
321 fprintf( fp, "gapmatrix g[i][][] except for %d\n", ex+1 );
322 for( i=0; i<100; i++ )
324 fprintf( fp, "site No.%#3d ", i+1 );
325 fprintf( fp, "00:%#5d ", -(int)g[i][0][0] );
326 fprintf( fp, "10:%#5d ", -(int)g[i][1][0] );
327 fprintf( fp, "20:%#5d ", -(int)g[i][2][0] );
328 fprintf( fp, "01:%#5d ", -(int)g[i][0][1] );
329 fprintf( fp, "11:%#5d ", -(int)g[i][1][1] );
330 fprintf( fp, "02:%#5d ", -(int)g[i][0][2] );
331 fprintf( fp, "22:%#5d ", -(int)g[i][2][2] );
340 w[0][0] += /* scmx[24][0] * penalty */ g[0][0][0]; /* [0][0] ? */
341 w[1][0] += g[0][0][1] + g[1][1][0] + gs[1][2] + gvsa[0]; /* ??? */
343 for( i=2; i<lgth+2; i++ )
345 tmpfloat += g[i-1][1][1] + gl[i-2][i-1] + gvsa[i-1];
346 w[i][0] += g[0][0][1] + gvsa[0] + tmpfloat + g[i][1][0] + gs[i][i+1];
349 w[0][1] += penalty * totaleff + n_dis[24][0] * totaleff;
351 for( j=2; j<lgth1+1; j++ )
353 tmpfloat += n_dis[24][0] * totaleff;
354 w[0][j] += penalty * totaleff + tmpfloat;
357 for( j=0; j<lgth1+1; ++j )
361 for( i=1; i<lgth+1; i++ )
364 for( j=1; j<lgth1+1; j++ )
368 x = w[i-1][j-2] + g[i-0][0][2] + n_dis[24][0] * totaleff;
369 mi += g[i-1][2][2] + n_dis[24][0] * totaleff;
378 mi = w[i-1][0] + g[i-0][0][2]/* + n_dis[24][0] * totaleff */; /* 0.0? */
380 fprintf( stderr, " i == %d j == 1, w[i][0] == %f, mi == %f, g[i][0][2] == %f\n", i, w[i][0], mi, g[i][0][2] );
383 mi = 0.0 + g[i-0][0][2]; * 0.0? *
390 x = w[i-2][j-1] + g[i-1][0][1] + gvsa[i-1];
391 m[j] += g[i-1][1][1] + gl[i-2][i-mp[j]-2] + gvsa[i-1]; /* ??? */
401 m[j] = w[0][j-1]+ g[i][0][1]/* + gvsa[1] */; /* 0.0? */
405 wmax = w[i-1][j-1] + g[i][0][0];
420 ijp[i][j] = -( j - mpi ); /* ijp[][] < 0 -> j ni gap */
423 x = m[j] + g[i][1][0] + gs[i][i-mp[j]];
431 ijp[i][j] = +( i - mp[j] );
438 w[lgth][lgth1] = w[lgth-1][lgth1-1] + g[lgth][0][0];
440 ip[lgth][lgth1] = lgth-1;
441 jp[lgth][lgth1] = lgth1-1;
443 ijp[lgth][lgth1] = 0;
446 *wm = w[lgth][lgth1];
448 for( i=0; i<lgth+1; i++ )
451 ip[i][0] = -1; jp[i][0] = -1;
455 for( j=0; j<lgth1+1; j++ )
458 ip[0][j] = -1; jp[0][j] = -1;
460 ijp[0][j] = -( j + 1 );
463 ip[lgth][lgth1]=iin; jp[lgth][lgth1]=jin;
464 ip[lgth+1][lgth1+1]=lgth; jp[lgth+1][lgth1+1]=lgth1;
467 tracking( nseq, aseq, seq, ip, jp, icyc );
469 tracking( nseq, aseq, seq, ijp, icyc );
471 fp = fopen( "matrx", "a" );
472 fprintf( fp, "matrix v\n" );
473 fprintf( fp, "%#5d", 0 );
474 for( j=0; j<20; j++ ) fprintf( fp, "%#5d", j );
476 for( i=0; i<lgth+1; i++ )
478 fprintf( fp, "%#5d", i );
479 for( j=0; j<20; j++ )
481 fprintf( fp, "%#5d", (int)( w[i][j]/(icyc+1) ) );
488 printf( "seq1[0] = %s\n", nseq[0] );
489 printf( "seq2[0] = %s\n", nseq[icyc+1] );
491 *wm = w[lgth][lgth1];
495 double Cscore_m_1( char **seq, int locnjob, int ex, double **eff )
498 int len = strlen( seq[0] );
499 int gb1, gb2, gc1, gc2;
510 glen1 = AllocateIntVec( locnjob );
511 glen2 = AllocateIntVec( locnjob );
515 printf( "in Cscore_m_1\n" );
516 for( i=0; i<locnjob; i++ )
518 if( i == ex ) continue;
519 fprintf( stdout, "%d-%d:%f\n", ex, i, eff[ex][i] );
522 for( k=0; k<len; k++ )
526 for( i=0; i<locnjob; i++ )
528 double efficient = eff[ex][i];
529 if( i == ex ) continue;
532 gb1 = ( seq[i][k-1] == '-' );
533 gb2 = ( seq[ex][k-1] == '-' );
541 if( gb1 ) glen1[i]++; else glen1[i] = 0;
542 if( gb2 ) glen2[i]++; else glen2[i] = 0;
544 gc1 = ( seq[i][k] == '-' );
545 gc2 = ( seq[ex][k] == '-' );
547 if( glen1[i] >= glen2[i] ) tmp1 = 1; else tmp1 = 0;
548 if( glen1[i] <= glen2[i] ) tmp2 = 1; else tmp2 = 0;
569 pen += (double)cob * penalty * efficient;
570 tmpscore += (double)amino_dis[(int)seq[i][k]][(int)seq[ex][k]] * efficient;
572 nglen += ( !gc1 * !gc2 );
577 printf( "%c<->%c * %f = %f\n", seq[ex][k], seq[i][k], efficient, amino_dis[seq[i][k]][seq[ex][k]] * efficient );
584 if( 1 ) fprintf( stdout, "%c %f ", seq[ex][k], score / (locnjob-1) );
585 if( 1 ) fprintf( stdout, "penalty %f\n", (double)pen / (locnjob - 1 ) );
589 fprintf( stdout, "in Cscore_m_1 %f\n", score );
592 fprintf( stdout, "%s\n", seq[ex] );
596 return( (double)score /*/ nglen + 400.0*/ );