6 #define USE_PENALTY_EX 1
10 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
13 int *intptr = amino_dis[(int)s1[0][i1]];
16 *match++ = intptr[(int)*seq2++];
18 static void match_calc_mtx( int mtx[0x80][0x80], float *match, char **s1, char **s2, int i1, int lgth2 )
21 int *intptr = mtx[(int)s1[0][i1]];
24 *match++ = intptr[(int)*seq2++];
27 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
31 for( j=0; j<lgth2; j++ )
32 match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
36 static float Atracking( float *lasthorizontalw, float *lastverticalw,
37 char **seq1, char **seq2,
38 char **mseq1, char **mseq2,
42 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
46 lgth1 = strlen( seq1[0] );
47 lgth2 = strlen( seq2[0] );
52 for( i=0; i<lgth1; i++ )
54 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
58 for( i=0; i<lgth1+1; i++ )
62 for( j=0; j<lgth2+1; j++ )
64 ijp[0][j] = -( j + 1 );
71 wm = lastverticalw[0];
72 for( i=0; i<lgth1; i++ )
74 if( lastverticalw[i] >= wm )
76 wm = lastverticalw[i];
77 iin = i; jin = lgth2-1;
78 ijp[lgth1][lgth2] = +( lgth1 - i );
81 for( j=0; j<lgth2; j++ )
83 if( lasthorizontalw[j] >= wm )
85 wm = lasthorizontalw[j];
86 iin = lgth1-1; jin = j;
87 ijp[lgth1][lgth2] = -( lgth2 - j );
94 mseq1[0] += lgth1+lgth2;
96 mseq2[0] += lgth1+lgth2;
100 iin = lgth1; jin = lgth2;
101 limk = lgth1+lgth2 + 1;
102 for( k=0; k<limk; k++ )
104 if( ijp[iin][jin] < 0 )
106 ifi = iin-1; jfi = jin+ijp[iin][jin];
108 else if( ijp[iin][jin] > 0 )
110 ifi = iin-ijp[iin][jin]; jfi = jin-1;
114 ifi = iin-1; jfi = jin-1;
119 *--mseq1[0] = seq1[0][ifi+l];
127 *--mseq2[0] = seq2[0][jfi+l];
130 if( iin <= 0 || jin <= 0 ) break;
131 *--mseq1[0] = seq1[0][ifi];
132 *--mseq2[0] = seq2[0][jfi];
134 iin = ifi; jin = jfi;
140 float G__align11( char **seq1, char **seq2, int alloclen, int headgp, int tailgp )
141 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
145 int lasti; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
148 float wm; /* int ?????? */
150 float *currentw, *previousw;
151 float fpenalty = (float)penalty;
153 float fpenalty_ex = (float)penalty_ex;
158 float *mjpt, *prept, *curpt;
161 static TLS float mi, *m;
162 static TLS int **ijp;
163 static TLS int mpi, *mp;
164 static TLS float *w1, *w2;
165 static TLS float *match;
166 static TLS float *initverticalw; /* kufuu sureba iranai */
167 static TLS float *lastverticalw; /* kufuu sureba iranai */
168 static TLS char **mseq1;
169 static TLS char **mseq2;
170 static TLS char **mseq;
171 static TLS int **intwork;
172 static TLS float **floatwork;
173 static TLS int orlgth1 = 0, orlgth2 = 0;
177 if( orlgth1 > 0 && orlgth2 > 0 )
185 FreeFloatVec( match );
186 FreeFloatVec( initverticalw );
187 FreeFloatVec( lastverticalw );
196 FreeFloatMtx( floatwork );
197 FreeIntMtx( intwork );
202 lgth1 = strlen( seq1[0] );
203 lgth2 = strlen( seq2[0] );
205 if( lgth1 <= 0 || lgth2 <= 0 )
207 fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
211 if( lgth1 == 0 && lgth2 == 0 )
217 while( lgth2 ) seq1[0][--lgth2] = '-';
218 // fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
225 while( lgth1 ) seq2[0][--lgth1] = '-';
226 // fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
236 mseq1 = AllocateCharMtx( njob, 0 );
237 mseq2 = AllocateCharMtx( njob, 0 );
242 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
246 if( orlgth1 > 0 && orlgth2 > 0 )
250 FreeFloatVec( match );
251 FreeFloatVec( initverticalw );
252 FreeFloatVec( lastverticalw );
261 FreeFloatMtx( floatwork );
262 FreeIntMtx( intwork );
265 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
266 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
269 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
272 w1 = AllocateFloatVec( ll2+2 );
273 w2 = AllocateFloatVec( ll2+2 );
274 match = AllocateFloatVec( ll2+2 );
276 initverticalw = AllocateFloatVec( ll1+2 );
277 lastverticalw = AllocateFloatVec( ll1+2 );
279 m = AllocateFloatVec( ll2+2 );
280 mp = AllocateIntVec( ll2+2 );
282 mseq = AllocateCharMtx( njob, ll1+ll2 );
285 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
286 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
289 fprintf( stderr, "succeeded\n" );
301 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
305 if( commonAlloc1 && commonAlloc2 )
307 FreeIntMtx( commonIP );
310 ll1 = MAX( orlgth1, commonAlloc1 );
311 ll2 = MAX( orlgth2, commonAlloc2 );
314 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
317 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
320 fprintf( stderr, "succeeded\n\n" );
330 for( i=0; i<lgth1; i++ )
331 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
338 match_calc( initverticalw, seq2, seq1, 0, lgth1 );
341 match_calc( currentw, seq1, seq2, 0, lgth2 );
345 for( i=1; i<lgth1+1; i++ )
347 initverticalw[i] += fpenalty;
349 for( j=1; j<lgth2+1; j++ )
351 currentw[j] += fpenalty;
355 for( j=1; j<lgth2+1; ++j )
357 m[j] = currentw[j-1]; mp[j] = 0;
361 lastverticalw[0] = 0.0; // lgth2==0 no toki error
363 lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
365 if( tailgp ) lasti = lgth1+1; else lasti = lgth1;
368 fprintf( stderr, "currentw = \n" );
369 for( i=0; i<lgth1+1; i++ )
371 fprintf( stderr, "%5.2f ", currentw[i] );
373 fprintf( stderr, "\n" );
374 fprintf( stderr, "initverticalw = \n" );
375 for( i=0; i<lgth2+1; i++ )
377 fprintf( stderr, "%5.2f ", initverticalw[i] );
379 fprintf( stderr, "\n" );
382 for( i=1; i<lasti; i++ )
385 previousw = currentw;
388 previousw[0] = initverticalw[i-1];
390 match_calc( currentw, seq1, seq2, i, lgth2 );
392 fprintf( stderr, "\n" );
393 fprintf( stderr, "i=%d\n", i );
394 fprintf( stderr, "currentw = \n" );
395 for( j=0; j<lgth2; j++ )
397 fprintf( stderr, "%5.2f ", currentw[j] );
399 fprintf( stderr, "\n" );
402 fprintf( stderr, "\n" );
403 fprintf( stderr, "i=%d\n", i );
404 fprintf( stderr, "currentw = \n" );
405 for( j=0; j<lgth2; j++ )
407 fprintf( stderr, "%5.2f ", currentw[j] );
409 fprintf( stderr, "\n" );
411 currentw[0] = initverticalw[i];
413 mi = previousw[0]; mpi = 0;
418 curpt = currentw + 1;
420 for( j=1; j<lgth2+1; j++ )
426 fprintf( stderr, "%5.0f->", wm );
429 fprintf( stderr, "%5.0f?", g );
431 if( (g=mi+fpenalty) > wm )
434 *ijppt = -( j - mpi );
436 if( (g=*prept) >= mi )
446 fprintf( stderr, "%5.0f?", g );
448 if( (g=*mjpt + fpenalty) > wm )
451 *ijppt = +( i - *mpjpt );
453 if( (g=*prept) >= *mjpt )
463 fprintf( stderr, "%5.0f ", wm );
471 lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
474 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, tailgp );
476 resultlen = strlen( mseq1[0] );
477 if( alloclen < resultlen || resultlen > N )
479 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
480 ErrorExit( "LENGTH OVER!\n" );
484 strcpy( seq1[0], mseq1[0] );
485 strcpy( seq2[0], mseq2[0] );
487 fprintf( stderr, "\n" );
488 fprintf( stderr, ">\n%s\n", mseq1[0] );
489 fprintf( stderr, ">\n%s\n", mseq2[0] );
490 fprintf( stderr, "wm = %f\n", wm );
496 float G__align11_noalign( int scoremtx[0x80][0x80], int penal, int penal_ex, char **seq1, char **seq2, int alloclen )
497 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
501 int lasti; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
504 float wm; /* int ?????? */
506 float *currentw, *previousw;
507 float fpenalty = (float)penal;
509 float fpenalty_ex = (float)penal_ex;
513 float *mjpt, *prept, *curpt;
516 static TLS float mi, *m;
517 static TLS float *w1, *w2;
518 static TLS float *match;
519 static TLS float *initverticalw; /* kufuu sureba iranai */
520 static TLS float *lastverticalw; /* kufuu sureba iranai */
521 static TLS int **intwork;
522 static TLS float **floatwork;
523 static TLS int orlgth1 = 0, orlgth2 = 0;
527 if( orlgth1 > 0 && orlgth2 > 0 )
533 FreeFloatVec( match );
534 FreeFloatVec( initverticalw );
535 FreeFloatVec( lastverticalw );
539 FreeFloatMtx( floatwork );
540 FreeIntMtx( intwork );
549 lgth1 = strlen( seq1[0] );
550 lgth2 = strlen( seq2[0] );
554 if( lgth1 <= 0 || lgth2 <= 0 )
556 fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
559 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
563 if( orlgth1 > 0 && orlgth2 > 0 )
567 FreeFloatVec( match );
568 FreeFloatVec( initverticalw );
569 FreeFloatVec( lastverticalw );
576 FreeFloatMtx( floatwork );
577 FreeIntMtx( intwork );
580 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
581 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
584 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
587 w1 = AllocateFloatVec( ll2+2 );
588 w2 = AllocateFloatVec( ll2+2 );
589 match = AllocateFloatVec( ll2+2 );
591 initverticalw = AllocateFloatVec( ll1+2 );
592 lastverticalw = AllocateFloatVec( ll1+2 );
594 m = AllocateFloatVec( ll2+2 );
598 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
599 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
602 fprintf( stderr, "succeeded\n" );
615 for( i=0; i<lgth1; i++ )
616 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
623 match_calc_mtx( scoremtx, initverticalw, seq2, seq1, 0, lgth1 );
626 match_calc_mtx( scoremtx, currentw, seq1, seq2, 0, lgth2 );
628 if( 1 ) // tsuneni outgap-1
630 for( i=1; i<lgth1+1; i++ )
632 initverticalw[i] += fpenalty;
634 for( j=1; j<lgth2+1; j++ )
636 currentw[j] += fpenalty;
640 for( j=1; j<lgth2+1; ++j )
642 m[j] = currentw[j-1];
646 lastverticalw[0] = 0.0; // lgth2==0 no toki error
648 lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
650 if( 1 ) lasti = lgth1+1; else lasti = lgth1; // tsuneni outgap-1
653 fprintf( stderr, "currentw = \n" );
654 for( i=0; i<lgth1+1; i++ )
656 fprintf( stderr, "%5.2f ", currentw[i] );
658 fprintf( stderr, "\n" );
659 fprintf( stderr, "initverticalw = \n" );
660 for( i=0; i<lgth2+1; i++ )
662 fprintf( stderr, "%5.2f ", initverticalw[i] );
664 fprintf( stderr, "\n" );
667 for( i=1; i<lasti; i++ )
670 previousw = currentw;
673 previousw[0] = initverticalw[i-1];
675 match_calc_mtx( scoremtx, currentw, seq1, seq2, i, lgth2 );
677 fprintf( stderr, "\n" );
678 fprintf( stderr, "i=%d\n", i );
679 fprintf( stderr, "currentw = \n" );
680 for( j=0; j<lgth2; j++ )
682 fprintf( stderr, "%5.2f ", currentw[j] );
684 fprintf( stderr, "\n" );
687 fprintf( stderr, "\n" );
688 fprintf( stderr, "i=%d\n", i );
689 fprintf( stderr, "currentw = \n" );
690 for( j=0; j<lgth2; j++ )
692 fprintf( stderr, "%5.2f ", currentw[j] );
694 fprintf( stderr, "\n" );
696 currentw[0] = initverticalw[i];
702 curpt = currentw + 1;
703 for( j=1; j<lgth2+1; j++ )
708 fprintf( stderr, "%5.0f->", wm );
711 fprintf( stderr, "%5.0f?", g );
713 if( (g=mi+fpenalty) > wm )
717 if( (g=*prept) >= mi )
726 fprintf( stderr, "%5.0f?", g );
728 if( (g=*mjpt + fpenalty) > wm )
732 if( (g=*prept) >= *mjpt )
741 fprintf( stderr, "%5.0f ", wm );
747 lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
751 fprintf( stderr, "\n" );
752 fprintf( stderr, ">\n%s\n", mseq1[0] );
753 fprintf( stderr, ">\n%s\n", mseq2[0] );
754 fprintf( stderr, "wm = %f\n", wm );