7 #define USE_PENALTY_EX 1
9 static short localstop;
12 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
17 intptr = amino_dis[s1[0][i1]];
19 *match++ = intptr[*seq2++];
22 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
26 for( j=0; j<lgth2; j++ )
27 match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
31 static void match_calc_bk( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
35 float **cpmxpd = floatwork;
36 int **cpmxpdn = intwork;
41 for( j=0; j<lgth2; j++ )
48 cpmxpd[count][j] = cpmx2[l][j];
49 cpmxpdn[count][j] = l;
53 cpmxpdn[count][j] = -1;
61 scarr[l] += n_dis[k][l] * cpmx1[k][i1];
63 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤
\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
65 float *fpt, **fptpt, *fpt2;
73 ipt=*iptpt,fpt=*fptpt;
75 *fpt2 += scarr[*ipt++] * *fpt++;
76 fpt2++,iptpt++,fptpt++;
80 for( j=0; j<lgth2; j++ )
83 for( k=0; cpmxpdn[k][j]>-1; k++ )
84 match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
89 static float Ltracking( float *lasthorizontalw, float *lastverticalw,
90 char **seq1, char **seq2,
91 char **mseq1, char **mseq2,
92 float **cpmx1, float **cpmx2,
93 short **ijp, int *off1pt, int *off2pt, int endi, int endj )
95 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
97 lgth1 = strlen( seq1[0] );
98 lgth2 = strlen( seq2[0] );
101 for( i=0; i<lgth1; i++ )
103 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
107 for( i=0; i<lgth1+1; i++ )
109 ijp[i][0] = localstop;
111 for( j=0; j<lgth2+1; j++ )
113 ijp[0][j] = localstop;
116 mseq1[0] += lgth1+lgth2;
118 mseq2[0] += lgth1+lgth2;
120 iin = endi; jin = endj;
122 for( k=0; k<=limk; k++ )
124 if( ijp[iin][jin] < 0 )
126 ifi = iin-1; jfi = jin+ijp[iin][jin];
128 else if( ijp[iin][jin] > 0 )
130 ifi = iin-ijp[iin][jin]; jfi = jin-1;
134 ifi = iin-1; jfi = jin-1;
139 *--mseq1[0] = seq1[0][ifi+l];
147 *--mseq2[0] = seq2[0][jfi+l];
151 if( iin <= 0 || jin <= 0 ) break;
152 *--mseq1[0] = seq1[0][ifi];
153 *--mseq2[0] = seq2[0][jfi];
154 if( ijp[ifi][jfi] == localstop ) break;
156 iin = ifi; jin = jfi;
158 if( ifi == -1 ) *off1pt = 0; else *off1pt = ifi;
159 if( jfi == -1 ) *off2pt = 0; else *off2pt = jfi;
161 // fprintf( stderr, "ifn = %d, jfn = %d\n", ifi, jfi );
168 float L__align11( char **seq1, char **seq2, int alloclen, int *off1pt, int *off2pt )
169 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
173 int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
176 float wm = 0.0; /* int ?????? */
178 float *currentw, *previousw;
182 float *mjpt, *prept, *curpt;
188 static float *w1, *w2;
190 static float *initverticalw; /* kufuu sureba iranai */
191 static float *lastverticalw; /* kufuu sureba iranai */
195 static float **cpmx1;
196 static float **cpmx2;
197 static int **intwork;
198 static float **floatwork;
199 static int orlgth1 = 0, orlgth2 = 0;
202 float localthr = -offset;
203 float localthr2 = -offset;
204 // float localthr = 100;
205 // float localthr2 = 100;
206 float fpenalty = (float)penalty;
207 float fpenalty_ex = (float)penalty_ex;
212 mseq1 = AllocateCharMtx( njob, 0 );
213 mseq2 = AllocateCharMtx( njob, 0 );
217 lgth1 = strlen( seq1[0] );
218 lgth2 = strlen( seq2[0] );
220 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
224 if( orlgth1 > 0 && orlgth2 > 0 )
228 FreeFloatVec( match );
229 FreeFloatVec( initverticalw );
230 FreeFloatVec( lastverticalw );
238 FreeFloatMtx( cpmx1 );
239 FreeFloatMtx( cpmx2 );
241 FreeFloatMtx( floatwork );
242 FreeIntMtx( intwork );
245 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
246 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
249 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
252 w1 = AllocateFloatVec( ll2+2 );
253 w2 = AllocateFloatVec( ll2+2 );
254 match = AllocateFloatVec( ll2+2 );
256 initverticalw = AllocateFloatVec( ll1+2 );
257 lastverticalw = AllocateFloatVec( ll1+2 );
259 m = AllocateFloatVec( ll2+2 );
260 mp = AllocateIntVec( ll2+2 );
262 mseq = AllocateCharMtx( njob, ll1+ll2 );
264 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
265 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
267 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
268 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
271 fprintf( stderr, "succeeded\n" );
283 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
287 if( commonAlloc1 && commonAlloc2 )
289 FreeShortMtx( commonIP );
292 ll1 = MAX( orlgth1, commonAlloc1 );
293 ll2 = MAX( orlgth2, commonAlloc2 );
296 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
299 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
302 fprintf( stderr, "succeeded\n\n" );
312 for( i=0; i<lgth1; i++ )
313 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
319 match_calc( initverticalw, seq2, seq1, 0, lgth1 );
321 match_calc( currentw, seq1, seq2, 0, lgth2 );
325 for( j=1; j<lasti; ++j )
327 m[j] = currentw[j-1]; mp[j] = 0;
329 if( m[j] < localthr ) m[j] = localthr2;
333 lastverticalw[0] = currentw[lgth2-1];
338 fprintf( stderr, "currentw = \n" );
339 for( i=0; i<lgth1+1; i++ )
341 fprintf( stderr, "%5.2f ", currentw[i] );
343 fprintf( stderr, "\n" );
344 fprintf( stderr, "initverticalw = \n" );
345 for( i=0; i<lgth2+1; i++ )
347 fprintf( stderr, "%5.2f ", initverticalw[i] );
349 fprintf( stderr, "\n" );
352 fprintf( stderr, "\n" );
353 fprintf( stderr, " " );
354 for( j=0; j<lgth2; j++ )
355 fprintf( stderr, "%c ", seq2[0][j] );
356 fprintf( stderr, "\n" );
359 localstop = lgth1+lgth2+1;
362 fprintf( stderr, "\n" );
363 fprintf( stderr, "%c ", seq1[0][0] );
365 for( j=0; j<lgth2+1; j++ )
366 fprintf( stderr, "%5.0f ", currentw[j] );
367 fprintf( stderr, "\n" );
370 for( i=1; i<lasti; i++ )
373 previousw = currentw;
376 previousw[0] = initverticalw[i-1];
378 match_calc( currentw, seq1, seq2, i, lgth2 );
380 fprintf( stderr, "%c ", seq1[0][i] );
381 fprintf( stderr, "%5.0f ", currentw[0] );
385 fprintf( stderr, "\n" );
386 fprintf( stderr, "i=%d\n", i );
387 fprintf( stderr, "currentw = \n" );
388 for( j=0; j<lgth2; j++ )
390 fprintf( stderr, "%5.2f ", currentw[j] );
392 fprintf( stderr, "\n" );
395 fprintf( stderr, "\n" );
396 fprintf( stderr, "i=%d\n", i );
397 fprintf( stderr, "currentw = \n" );
398 for( j=0; j<lgth2; j++ )
400 fprintf( stderr, "%5.2f ", currentw[j] );
402 fprintf( stderr, "\n" );
404 currentw[0] = initverticalw[i];
406 mi = previousw[0]; mpi = 0;
409 if( mi < localthr ) mi = localthr2;
415 curpt = currentw + 1;
418 for( j=1; j<lastj; j++ )
424 fprintf( stderr, "%5.0f->", wm );
427 fprintf( stderr, "%5.0f?", g );
429 if( (g=mi+fpenalty) > wm )
432 *ijppt = -( j - mpi );
445 fprintf( stderr, "%5.0f?", g );
447 if( (g=*mjpt+fpenalty) > wm )
450 *ijppt = +( i - *mpjpt );
458 *mjpt += fpenalty_ex;
470 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
476 fprintf( stderr, "%5.0f ", *curpt );
479 fprintf( stderr, "%5.0f ", wm );
480 // fprintf( stderr, "%c-%c *ijppt = %d, localstop = %d\n", seq1[0][i], seq2[0][j], *ijppt, localstop );
490 fprintf( stderr, "\n" );
493 lastverticalw[i] = currentw[lgth2-1];
498 fprintf( stderr, "maxwm = %f\n", maxwm );
499 fprintf( stderr, "endali = %d\n", endali );
500 fprintf( stderr, "endalj = %d\n", endalj );
503 Ltracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, off1pt, off2pt, endali, endalj );
506 resultlen = strlen( mseq1[0] );
507 if( alloclen < resultlen || resultlen > N )
509 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
510 ErrorExit( "LENGTH OVER!\n" );
514 strcpy( seq1[0], mseq1[0] );
515 strcpy( seq2[0], mseq2[0] );
518 fprintf( stderr, "wm=%f\n", wm );
519 fprintf( stderr, ">\n%s\n", mseq1[0] );
520 fprintf( stderr, ">\n%s\n", mseq2[0] );
527 float L__align11_noalign( char **seq1, char **seq2, int alloclen, int *off1pt, int *off2pt )
528 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
532 int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
535 float wm = 0.0; /* int ?????? */
537 float *currentw, *previousw;
541 float *mjpt, *prept, *curpt;
547 static float *w1, *w2;
549 static float *initverticalw; /* kufuu sureba iranai */
550 static float *lastverticalw; /* kufuu sureba iranai */
554 static float **cpmx1;
555 static float **cpmx2;
556 static int **intwork;
557 static float **floatwork;
558 static int orlgth1 = 0, orlgth2 = 0;
561 float localthr = -offset;
562 float localthr2 = -offset;
563 // float localthr = 100;
564 // float localthr2 = 100;
565 float fpenalty = (float)penalty;
566 float fpenalty_ex = (float)penalty_ex;
571 mseq1 = AllocateCharMtx( njob, 0 );
572 mseq2 = AllocateCharMtx( njob, 0 );
576 lgth1 = strlen( seq1[0] );
577 lgth2 = strlen( seq2[0] );
579 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
583 if( orlgth1 > 0 && orlgth2 > 0 )
587 FreeFloatVec( match );
588 FreeFloatVec( initverticalw );
589 FreeFloatVec( lastverticalw );
597 FreeFloatMtx( cpmx1 );
598 FreeFloatMtx( cpmx2 );
600 FreeFloatMtx( floatwork );
601 FreeIntMtx( intwork );
604 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
605 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
608 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
611 w1 = AllocateFloatVec( ll2+2 );
612 w2 = AllocateFloatVec( ll2+2 );
613 match = AllocateFloatVec( ll2+2 );
615 initverticalw = AllocateFloatVec( ll1+2 );
616 lastverticalw = AllocateFloatVec( ll1+2 );
618 m = AllocateFloatVec( ll2+2 );
619 mp = AllocateIntVec( ll2+2 );
621 mseq = AllocateCharMtx( njob, ll1+ll2 );
623 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
624 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
626 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
627 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
630 fprintf( stderr, "succeeded\n" );
643 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
647 if( commonAlloc1 && commonAlloc2 )
649 FreeShortMtx( commonIP );
652 ll1 = MAX( orlgth1, commonAlloc1 );
653 ll2 = MAX( orlgth2, commonAlloc2 );
656 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
659 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
662 fprintf( stderr, "succeeded\n\n" );
673 for( i=0; i<lgth1; i++ )
674 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
680 match_calc( initverticalw, seq2, seq1, 0, lgth1 );
682 match_calc( currentw, seq1, seq2, 0, lgth2 );
686 for( j=1; j<lasti; ++j )
688 m[j] = currentw[j-1]; mp[j] = 0;
690 if( m[j] < localthr ) m[j] = localthr2;
694 lastverticalw[0] = currentw[lgth2-1];
699 fprintf( stderr, "currentw = \n" );
700 for( i=0; i<lgth1+1; i++ )
702 fprintf( stderr, "%5.2f ", currentw[i] );
704 fprintf( stderr, "\n" );
705 fprintf( stderr, "initverticalw = \n" );
706 for( i=0; i<lgth2+1; i++ )
708 fprintf( stderr, "%5.2f ", initverticalw[i] );
710 fprintf( stderr, "\n" );
713 fprintf( stderr, "\n" );
714 fprintf( stderr, " " );
715 for( j=0; j<lgth2; j++ )
716 fprintf( stderr, "%c ", seq2[0][j] );
717 fprintf( stderr, "\n" );
720 localstop = lgth1+lgth2+1;
723 fprintf( stderr, "\n" );
724 fprintf( stderr, "%c ", seq1[0][0] );
726 for( j=0; j<lgth2+1; j++ )
727 fprintf( stderr, "%5.0f ", currentw[j] );
728 fprintf( stderr, "\n" );
731 for( i=1; i<lasti; i++ )
734 previousw = currentw;
737 previousw[0] = initverticalw[i-1];
739 match_calc( currentw, seq1, seq2, i, lgth2 );
741 fprintf( stderr, "%c ", seq1[0][i] );
742 fprintf( stderr, "%5.0f ", currentw[0] );
746 fprintf( stderr, "\n" );
747 fprintf( stderr, "i=%d\n", i );
748 fprintf( stderr, "currentw = \n" );
749 for( j=0; j<lgth2; j++ )
751 fprintf( stderr, "%5.2f ", currentw[j] );
753 fprintf( stderr, "\n" );
756 fprintf( stderr, "\n" );
757 fprintf( stderr, "i=%d\n", i );
758 fprintf( stderr, "currentw = \n" );
759 for( j=0; j<lgth2; j++ )
761 fprintf( stderr, "%5.2f ", currentw[j] );
763 fprintf( stderr, "\n" );
765 currentw[0] = initverticalw[i];
767 mi = previousw[0]; mpi = 0;
770 if( mi < localthr ) mi = localthr2;
773 // ijppt = ijp[i] + 1;
776 curpt = currentw + 1;
779 for( j=1; j<lastj; j++ )
785 fprintf( stderr, "%5.0f->", wm );
788 fprintf( stderr, "%5.0f?", g );
790 if( (g=mi+fpenalty) > wm )
793 // *ijppt = -( j - mpi );
806 fprintf( stderr, "%5.0f?", g );
808 if( (g=*mjpt+fpenalty) > wm )
811 // *ijppt = +( i - *mpjpt );
819 *mjpt += fpenalty_ex;
831 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
832 // *ijppt = localstop;
837 fprintf( stderr, "%5.0f ", *curpt );
840 fprintf( stderr, "%5.0f ", wm );
841 // fprintf( stderr, "%c-%c *ijppt = %d, localstop = %d\n", seq1[0][i], seq2[0][j], *ijppt, localstop );
851 fprintf( stderr, "\n" );
854 lastverticalw[i] = currentw[lgth2-1];
859 fprintf( stderr, "maxwm = %f\n", maxwm );
860 fprintf( stderr, "endali = %d\n", endali );
861 fprintf( stderr, "endalj = %d\n", endalj );
866 Ltracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, off1pt, off2pt, endali, endalj );
869 resultlen = strlen( mseq1[0] );
870 if( alloclen < resultlen || resultlen > N )
872 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
873 ErrorExit( "LENGTH OVER!\n" );
877 strcpy( seq1[0], mseq1[0] );
878 strcpy( seq2[0], mseq2[0] );
880 fprintf( stderr, "wm=%f\n", wm );
881 fprintf( stderr, ">\n%s\n", mseq1[0] );
882 fprintf( stderr, ">\n%s\n", mseq2[0] );