8 #define USE_PENALTY_EX 0
9 #define FASTMATCHCALC 1
12 static TLS float **impmtx = NULL;
13 static TLS int impalloclen = 0;
14 float imp_match_out_sc( int i1, int j1 )
16 // fprintf( stderr, "imp+match = %f\n", impmtx[i1][j1] * fastathreshold );
17 // fprintf( stderr, "val = %f\n", impmtx[i1][j1] );
18 return( impmtx[i1][j1] );
20 static void imp_match_out_vead_gapmap( float *imp, int i1, int lgth2, int *gapmap2 )
23 float *pt = impmtx[i1];
24 int *gapmappt = gapmap2;
26 *imp++ += pt[*gapmappt++];
29 float *pt = impmtx[i1];
30 for( j=0; j<lgth2; j++ )
31 *imp++ += pt[gapmap2[j]];
36 static void imp_match_out_vead_tate_gapmap( float *imp, int j1, int lgth1, int *gapmap1 )
39 int *gapmappt = gapmap1;
41 *imp++ += impmtx[*gapmappt++][j1];
44 for( i=0; i<lgth1; i++ )
45 *imp++ += impmtx[gapmap1[i]][j1];
49 static void imp_match_out_vead( float *imp, int i1, int lgth2 )
52 float *pt = impmtx[i1];
57 float *pt = impmtx[i1];
58 for( j=0; j<lgth2; j++ )
62 static void imp_match_out_vead_tate( float *imp, int j1, int lgth1 )
65 for( i=0; i<lgth1; i++ )
66 *imp++ += impmtx[i][j1];
69 void imp_rna( int nseq1, int nseq2, char **seq1, char **seq2, double *eff1, double *eff2, RNApair ***grouprna1, RNApair ***grouprna2, int *gapmap1, int *gapmap2, RNApair *pair )
71 foldrna( nseq1, nseq2, seq1, seq2, eff1, eff2, grouprna1, grouprna2, impmtx, gapmap1, gapmap2, pair );
75 void imp_match_init_strict( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, double *eff1_kozo, double *eff2_kozo, LocalHom ***localhom, int forscore )
77 int i, j, k1, k2, tmpint, start1, start2, end1, end2;
82 static TLS char *nocount1 = NULL;
83 static TLS char *nocount2 = NULL;
88 if( impmtx ) FreeFloatMtx( impmtx );
90 if( nocount1 ) free( nocount1 );
92 if( nocount2 ) free( nocount2 );
98 if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
100 if( impmtx ) FreeFloatMtx( impmtx );
101 if( nocount1 ) free( nocount1 );
102 if( nocount2 ) free( nocount2 );
103 impalloclen = MAX( lgth1, lgth2 ) + 2;
104 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
105 nocount1 = AllocateCharVec( impalloclen );
106 nocount2 = AllocateCharVec( impalloclen );
109 for( i=0; i<lgth1; i++ )
111 for( j=0; j<clus1; j++ )
112 if( seq1[j][i] == '-' ) break;
113 if( j != clus1 ) nocount1[i] = 1;
114 else nocount1[i] = 0;
116 for( i=0; i<lgth2; i++ )
118 for( j=0; j<clus2; j++ )
119 if( seq2[j][i] == '-' ) break;
120 if( j != clus2 ) nocount2[i] = 1;
121 else nocount2[i] = 0;
125 fprintf( stderr, "nocount2 =\n" );
126 for( i = 0; i<impalloclen; i++ )
128 fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
135 fprintf( stderr, "eff1 in _init_strict = \n" );
136 for( i=0; i<clus1; i++ )
137 fprintf( stderr, "eff1[] = %f\n", eff1[i] );
138 for( i=0; i<clus2; i++ )
139 fprintf( stderr, "eff2[] = %f\n", eff2[i] );
142 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
144 effijx = fastathreshold;
145 for( i=0; i<clus1; i++ )
147 for( j=0; j<clus2; j++ )
149 effij = (float)( eff1[i] * eff2[j] * effijx );
150 effij_kozo = (float)( eff1_kozo[i] * eff2_kozo[j] * effijx );
151 tmpptr = localhom[i][j];
154 // fprintf( stderr, "start1 = %d\n", tmpptr->start1 );
155 // fprintf( stderr, "end1 = %d\n", tmpptr->end1 );
156 // fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] );
157 // fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] );
162 if( *pt++ != '-' ) tmpint++;
163 if( tmpint == tmpptr->start1 ) break;
165 start1 = pt - seq1[i] - 1;
167 if( tmpptr->start1 == tmpptr->end1 ) end1 = start1;
173 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
174 if( tmpint == tmpptr->end1 ) break;
175 if( *pt++ != '-' ) tmpint++;
177 end1 = pt - seq1[i] - 0;
181 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
182 if( *pt++ != '-' ) tmpint++;
183 if( tmpint == tmpptr->end1 ) break;
185 end1 = pt - seq1[i] - 1;
193 if( *pt++ != '-' ) tmpint++;
194 if( tmpint == tmpptr->start2 ) break;
196 start2 = pt - seq2[j] - 1;
197 if( tmpptr->start2 == tmpptr->end2 ) end2 = start2;
203 if( tmpint == tmpptr->end2 ) break;
204 if( *pt++ != '-' ) tmpint++;
206 end2 = pt - seq2[j] - 0;
210 if( *pt++ != '-' ) tmpint++;
211 if( tmpint == tmpptr->end2 ) break;
213 end2 = pt - seq2[j] - 1;
216 // fprintf( stderr, "start1 = %d (%c), end1 = %d (%c), start2 = %d (%c), end2 = %d (%c)\n", start1, seq1[i][start1], end1, seq1[i][end1], start2, seq2[j][start2], end2, seq2[j][end2] );
217 // fprintf( stderr, "step 0\n" );
218 if( end1 - start1 != end2 - start2 )
220 // fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
224 k1 = start1; k2 = start2;
227 while( *pt1 && *pt2 )
229 if( *pt1 != '-' && *pt2 != '-' )
231 // ½Å¤ß¤òÆó½Å¤Ë¤«¤±¤Ê¤¤¤è¤¦¤ËÃðդ·¤Æ²¼¤µ¤¤¡£
232 // impmtx[k1][k2] += tmpptr->wimportance * fastathreshold;
233 // impmtx[k1][k2] += tmpptr->importance * effij;
234 // impmtx[k1][k2] += tmpptr->fimportance * effij;
235 if( tmpptr->korh == 'k' )
236 impmtx[k1][k2] += tmpptr->fimportance * effij_kozo;
238 impmtx[k1][k2] += tmpptr->fimportance * effij;
240 // fprintf( stderr, "#### impmtx[k1][k2] = %f, tmpptr->fimportance=%f, effij=%f\n", impmtx[k1][k2], tmpptr->fimportance, effij );
241 // fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
242 // fprintf( stderr, "%d (%c) - %d (%c) - %f\n", k1, *pt1, k2, *pt2, tmpptr->fimportance * effij );
246 else if( *pt1 != '-' && *pt2 == '-' )
248 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
251 else if( *pt1 == '-' && *pt2 != '-' )
253 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
256 else if( *pt1 == '-' && *pt2 == '-' )
258 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
262 if( k1 > end1 || k2 > end2 ) break;
265 while( k1 <= end1 && k2 <= end2 )
267 fprintf( stderr, "k1,k2=%d,%d - ", k1, k2 );
268 if( !nocount1[k1] && !nocount2[k2] )
270 impmtx[k1][k2] += tmpptr->wimportance * eff1[i] * eff2[j] * fastathreshold;
271 fprintf( stderr, "marked\n" );
274 fprintf( stderr, "no count\n" );
278 tmpptr = tmpptr->next;
284 if( clus1 == 1 && clus2 == 1 )
286 fprintf( stderr, "writing impmtx\n" );
287 fprintf( stderr, "\n" );
288 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
289 fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
290 fprintf( stderr, "impmtx = \n" );
291 for( k2=0; k2<lgth2; k2++ )
292 fprintf( stderr, "%6.3f ", (double)k2 );
293 fprintf( stderr, "\n" );
294 for( k1=0; k1<lgth1; k1++ )
296 fprintf( stderr, "%d ", k1 );
297 for( k2=0; k2<30; k2++ )
298 fprintf( stderr, "%2.1f ", impmtx[k1][k2] );
299 fprintf( stderr, "\n" );
307 void imp_match_init( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom )
309 int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2;
310 static TLS int impalloclen = 0;
313 static TLS char *nocount1 = NULL;
314 static TLS char *nocount2 = NULL;
316 if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
318 if( impmtx ) FreeFloatMtx( impmtx );
319 if( nocount1 ) free( nocount1 );
320 if( nocount2 ) free( nocount2 );
321 impalloclen = MAX( lgth1, lgth2 ) + 2;
322 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
323 nocount1 = AllocateCharVec( impalloclen );
324 nocount2 = AllocateCharVec( impalloclen );
327 for( i=0; i<lgth1; i++ )
329 for( j=0; j<clus1; j++ )
330 if( seq1[j][i] == '-' ) break;
331 if( j != clus1 ) nocount1[i] = 1;
332 else nocount1[i] = 0;
334 for( i=0; i<lgth2; i++ )
336 for( j=0; j<clus2; j++ )
337 if( seq2[j][i] == '-' ) break;
338 if( j != clus2 ) nocount2[i] = 1;
339 else nocount2[i] = 0;
343 fprintf( stderr, "nocount2 =\n" );
344 for( i = 0; i<impalloclen; i++ )
346 fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
350 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
352 for( i=0; i<clus1; i++ )
354 fprintf( stderr, "i = %d, seq1 = %s\n", i, seq1[i] );
355 for( j=0; j<clus2; j++ )
357 fprintf( stderr, "start1 = %d\n", localhom[i][j]->start1 );
358 fprintf( stderr, "end1 = %d\n", localhom[i][j]->end1 );
359 fprintf( stderr, "j = %d, seq2 = %s\n", j, seq2[j] );
364 if( *pt++ != '-' ) tmpint++;
365 if( tmpint == localhom[i][j]->start1 ) break;
367 start1 = pt - seq1[i] - 1;
371 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, localhom[i][j].end1, pt-seq1[i] );
372 if( *pt++ != '-' ) tmpint++;
373 if( tmpint == localhom[i][j]->end1 ) break;
375 end1 = pt - seq1[i] - 1;
381 if( *pt++ != '-' ) tmpint++;
382 if( tmpint == localhom[i][j]->start2 ) break;
384 start2 = pt - seq2[j] - 1;
387 if( *pt++ != '-' ) tmpint++;
388 if( tmpint == localhom[i][j]->end2 ) break;
390 end2 = pt - seq2[j] - 1;
391 // fprintf( stderr, "start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
394 fprintf( stderr, "step 0\n" );
395 while( k1 <= end1 && k2 <= end2 )
398 if( !nocount1[k1] && !nocount2[k2] )
399 impmtx[k1][k2] += localhom[i][j].wimportance * eff1[i] * eff2[j];
402 if( !nocount1[k1] && !nocount2[k2] )
403 impmtx[k1][k2] += localhom[i][j]->wimportance * eff1[i] * eff2[j];
408 dif = ( end1 - start1 ) - ( end2 - start2 );
409 fprintf( stderr, "dif = %d\n", dif );
414 fprintf( stderr, "dif = %d\n", dif );
417 while( k1 <= end1 && k2 <= end2 )
419 if( 0 <= k2 && start2 <= k2 && !nocount1[k1] && !nocount2[k2] )
420 impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
434 if( k1 >= 0 && k1 >= start1 && !nocount1[k1] && !nocount2[k2] )
435 impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
444 fprintf( stderr, "impmtx = \n" );
445 for( k2=0; k2<lgth2; k2++ )
446 fprintf( stderr, "%6.3f ", (double)k2 );
447 fprintf( stderr, "\n" );
448 for( k1=0; k1<lgth1; k1++ )
450 fprintf( stderr, "%d", k1 );
451 for( k2=0; k2<lgth2; k2++ )
452 fprintf( stderr, "%6.3f ", impmtx[k1][k2] );
453 fprintf( stderr, "\n" );
459 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
464 float **cpmxpd = floatwork;
465 int **cpmxpdn = intwork;
466 float *matchpt, *cpmxpdpt, **cpmxpdptpt;
467 int *cpmxpdnpt, **cpmxpdnptpt;
471 for( j=0; j<lgth2; j++ )
474 for( l=0; l<26; l++ )
478 cpmxpd[j][count] = cpmx2[l][j];
479 cpmxpdn[j][count] = l;
483 cpmxpdn[j][count] = -1;
488 for( l=0; l<26; l++ )
491 for( j=0; j<26; j++ )
492 // scarr[l] += n_dis[j][l] * cpmx1[j][i1];
493 scarr[l] += n_dis_consweight_multi[j][l] * cpmx1[j][i1];
496 cpmxpdnptpt = cpmxpdn;
501 cpmxpdnpt = *cpmxpdnptpt++;
502 cpmxpdpt = *cpmxpdptpt++;
503 while( *cpmxpdnpt>-1 )
504 *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;
511 float **cpmxpd = floatwork;
512 int **cpmxpdn = intwork;
517 for( j=0; j<lgth2; j++ )
520 for( l=0; l<26; l++ )
524 cpmxpd[count][j] = cpmx2[l][j];
525 cpmxpdn[count][j] = l;
529 cpmxpdn[count][j] = -1;
532 for( l=0; l<26; l++ )
535 for( k=0; k<26; k++ )
536 scarr[l] += n_dis_consweight_multi[k][l] * cpmx1[k][i1];
537 // scarr[l] += n_dis[k][l] * cpmx1[k][i1];
539 for( j=0; j<lgth2; j++ )
542 for( k=0; cpmxpdn[k][j]>-1; k++ )
543 match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
548 static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
549 char **seq1, char **seq2,
550 char **mseq1, char **mseq2,
551 float **cpmx1, float **cpmx2,
552 int **ijp, int icyc, int jcyc )
554 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
556 char *gaptable1, *gt1bk;
557 char *gaptable2, *gt2bk;
558 lgth1 = strlen( seq1[0] );
559 lgth2 = strlen( seq2[0] );
560 gt1bk = AllocateCharVec( lgth1+lgth2+1 );
561 gt2bk = AllocateCharVec( lgth1+lgth2+1 );
564 for( i=0; i<lgth1; i++ )
566 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
574 wm = lastverticalw[0];
575 for( i=0; i<lgth1; i++ )
577 if( lastverticalw[i] >= wm )
579 wm = lastverticalw[i];
580 iin = i; jin = lgth2-1;
581 ijp[lgth1][lgth2] = +( lgth1 - i );
584 for( j=0; j<lgth2; j++ )
586 if( lasthorizontalw[j] >= wm )
588 wm = lasthorizontalw[j];
589 iin = lgth1-1; jin = j;
590 ijp[lgth1][lgth2] = -( lgth2 - j );
595 for( i=0; i<lgth1+1; i++ )
599 for( j=0; j<lgth2+1; j++ )
601 ijp[0][j] = -( j + 1 );
604 gaptable1 = gt1bk + lgth1+lgth2;
606 gaptable2 = gt2bk + lgth1+lgth2;
609 iin = lgth1; jin = lgth2;
611 for( k=0; k<=lgth1+lgth2; k++ )
613 if( ijp[iin][jin] < 0 )
615 ifi = iin-1; jfi = jin+ijp[iin][jin];
617 else if( ijp[iin][jin] > 0 )
619 ifi = iin-ijp[iin][jin]; jfi = jin-1;
623 ifi = iin-1; jfi = jin-1;
639 if( iin == lgth1 || jin == lgth2 )
643 *impwmpt += imp_match_out_sc( iin, jin );
645 // fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
647 if( iin <= 0 || jin <= 0 ) break;
651 iin = ifi; jin = jfi;
654 for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );
655 for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );
660 static void Atracking_localhom_gapmap( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
661 char **seq1, char **seq2,
662 char **mseq1, char **mseq2,
663 float **cpmx1, float **cpmx2,
664 int **ijp, int icyc, int jcyc,
665 int *gapmap1, int *gapmap2 )
667 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
669 char *gaptable1, *gt1bk;
670 char *gaptable2, *gt2bk;
671 lgth1 = strlen( seq1[0] );
672 lgth2 = strlen( seq2[0] );
673 gt1bk = AllocateCharVec( lgth1+lgth2+1 );
674 gt2bk = AllocateCharVec( lgth1+lgth2+1 );
677 for( i=0; i<lgth1; i++ )
679 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
687 wm = lastverticalw[0];
688 for( i=0; i<lgth1; i++ )
690 if( lastverticalw[i] >= wm )
692 wm = lastverticalw[i];
693 iin = i; jin = lgth2-1;
694 ijp[lgth1][lgth2] = +( lgth1 - i );
697 for( j=0; j<lgth2; j++ )
699 if( lasthorizontalw[j] >= wm )
701 wm = lasthorizontalw[j];
702 iin = lgth1-1; jin = j;
703 ijp[lgth1][lgth2] = -( lgth2 - j );
708 for( i=0; i<lgth1+1; i++ )
712 for( j=0; j<lgth2+1; j++ )
714 ijp[0][j] = -( j + 1 );
717 gaptable1 = gt1bk + lgth1+lgth2;
719 gaptable2 = gt2bk + lgth1+lgth2;
722 iin = lgth1; jin = lgth2;
724 for( k=0; k<=lgth1+lgth2; k++ )
726 if( ijp[iin][jin] < 0 )
728 ifi = iin-1; jfi = jin+ijp[iin][jin];
730 else if( ijp[iin][jin] > 0 )
732 ifi = iin-ijp[iin][jin]; jfi = jin-1;
736 ifi = iin-1; jfi = jin-1;
752 if( iin == lgth1 || jin == lgth2 )
756 *impwmpt += imp_match_out_sc( gapmap1[iin], gapmap2[jin] );
758 // fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
760 if( iin <= 0 || jin <= 0 ) break;
764 iin = ifi; jin = jfi;
766 for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );
767 for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );
772 static float Atracking( float *lasthorizontalw, float *lastverticalw,
773 char **seq1, char **seq2,
774 char **mseq1, char **mseq2,
775 float **cpmx1, float **cpmx2,
776 int **ijp, int icyc, int jcyc,
779 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
781 char *gaptable1, *gt1bk;
782 char *gaptable2, *gt2bk;
783 lgth1 = strlen( seq1[0] );
784 lgth2 = strlen( seq2[0] );
786 gt1bk = AllocateCharVec( lgth1+lgth2+1 );
787 gt2bk = AllocateCharVec( lgth1+lgth2+1 );
790 for( i=0; i<lgth1; i++ )
792 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
800 wm = lastverticalw[0];
801 for( i=0; i<lgth1; i++ )
803 if( lastverticalw[i] >= wm )
805 wm = lastverticalw[i];
806 iin = i; jin = lgth2-1;
807 ijp[lgth1][lgth2] = +( lgth1 - i );
810 for( j=0; j<lgth2; j++ )
812 if( lasthorizontalw[j] >= wm )
814 wm = lasthorizontalw[j];
815 iin = lgth1-1; jin = j;
816 ijp[lgth1][lgth2] = -( lgth2 - j );
821 for( i=0; i<lgth1+1; i++ )
825 for( j=0; j<lgth2+1; j++ )
827 ijp[0][j] = -( j + 1 );
830 gaptable1 = gt1bk + lgth1+lgth2;
832 gaptable2 = gt2bk + lgth1+lgth2;
835 iin = lgth1; jin = lgth2;
836 for( k=0; k<=lgth1+lgth2; k++ )
838 if( ijp[iin][jin] < 0 )
840 ifi = iin-1; jfi = jin+ijp[iin][jin];
842 else if( ijp[iin][jin] > 0 )
844 ifi = iin-ijp[iin][jin]; jfi = jin-1;
848 ifi = iin-1; jfi = jin-1;
864 if( iin <= 0 || jin <= 0 ) break;
868 iin = ifi; jin = jfi;
871 for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );
872 for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );
880 float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, char *sgap1, char *sgap2, char *egap1, char *egap2, int *chudanpt, int chudanref, int *chudanres, int headgp, int tailgp )
881 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
885 int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
888 float wm = 0.0; /* int ?????? */
890 float *currentw, *previousw;
891 // float fpenalty = (float)penalty;
893 float fpenalty_ex = (float)penalty_ex;
898 float *mjpt, *prept, *curpt;
901 static TLS float mi, *m;
902 static TLS int **ijp;
903 static TLS int mpi, *mp;
904 static TLS float *w1, *w2;
905 static TLS float *match;
906 static TLS float *initverticalw; /* kufuu sureba iranai */
907 static TLS float *lastverticalw; /* kufuu sureba iranai */
908 static TLS char **mseq1;
909 static TLS char **mseq2;
910 static TLS char **mseq;
911 static TLS float *ogcp1;
912 static TLS float *ogcp2;
913 static TLS float *fgcp1;
914 static TLS float *fgcp2;
915 static TLS float **cpmx1;
916 static TLS float **cpmx2;
917 static TLS int **intwork;
918 static TLS float **floatwork;
919 static TLS int orlgth1 = 0, orlgth2 = 0;
920 float fpenalty = (float)penalty;
931 // fprintf( stderr, "## Freeing local arrays in A__align\n" );
935 imp_match_init_strict( NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0 );
941 FreeFloatVec( match );
942 FreeFloatVec( initverticalw );
943 FreeFloatVec( lastverticalw );
950 FreeFloatVec( ogcp1 );
951 FreeFloatVec( ogcp2 );
952 FreeFloatVec( fgcp1 );
953 FreeFloatVec( fgcp2 );
956 FreeFloatMtx( cpmx1 );
957 FreeFloatMtx( cpmx2 );
959 FreeFloatMtx( floatwork );
960 FreeIntMtx( intwork );
965 // fprintf( stderr, "## Not allocated\n" );
970 lgth1 = strlen( seq1[0] );
971 lgth2 = strlen( seq2[0] );
973 if( lgth1 == 0 || lgth2 == 0 )
975 fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
977 if( lgth1 == 0 && lgth2 == 0 )
982 for( i=0; i<icyc; i++ )
986 while( j ) seq1[i][--j] = '-';
987 // fprintf( stderr, "seq1[i] = %s\n", seq1[i] );
994 for( i=0; i<jcyc; i++ )
998 while( j ) seq2[i][--j] = '-';
999 // fprintf( stderr, "seq2[i] = %s\n", seq2[i] );
1007 fprintf( stderr, "#### eff in SA+++align\n" );
1008 fprintf( stderr, "#### seq1[0] = %s\n", seq1[0] );
1009 fprintf( stderr, "#### strlen( seq1[0] ) = %d\n", strlen( seq1[0] ) );
1010 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1011 fprintf( stderr, "#### seq2[0] = %s\n", seq2[0] );
1012 fprintf( stderr, "#### strlen( seq2[0] ) = %d\n", strlen( seq2[0] ) );
1013 for( i=0; i<jcyc; i++ ) fprintf( stderr, "eff2[%d] = %f\n", i, eff2[i] );
1017 mseq1 = AllocateCharMtx( njob, 0 );
1018 mseq2 = AllocateCharMtx( njob, 0 );
1023 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
1027 if( orlgth1 > 0 && orlgth2 > 0 )
1031 FreeFloatVec( match );
1032 FreeFloatVec( initverticalw );
1033 FreeFloatVec( lastverticalw );
1038 FreeCharMtx( mseq );
1040 FreeFloatVec( ogcp1 );
1041 FreeFloatVec( ogcp2 );
1042 FreeFloatVec( fgcp1 );
1043 FreeFloatVec( fgcp2 );
1046 FreeFloatMtx( cpmx1 );
1047 FreeFloatMtx( cpmx2 );
1049 FreeFloatMtx( floatwork );
1050 FreeIntMtx( intwork );
1053 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
1054 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
1057 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1060 w1 = AllocateFloatVec( ll2+2 );
1061 w2 = AllocateFloatVec( ll2+2 );
1062 match = AllocateFloatVec( ll2+2 );
1064 initverticalw = AllocateFloatVec( ll1+2 );
1065 lastverticalw = AllocateFloatVec( ll1+2 );
1067 m = AllocateFloatVec( ll2+2 );
1068 mp = AllocateIntVec( ll2+2 );
1070 mseq = AllocateCharMtx( njob, ll1+ll2 );
1072 ogcp1 = AllocateFloatVec( ll1+2 );
1073 ogcp2 = AllocateFloatVec( ll2+2 );
1074 fgcp1 = AllocateFloatVec( ll1+2 );
1075 fgcp2 = AllocateFloatVec( ll2+2 );
1077 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1078 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1081 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
1082 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 );
1084 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
1085 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
1089 fprintf( stderr, "succeeded\n" );
1092 orlgth1 = ll1 - 100;
1093 orlgth2 = ll2 - 100;
1097 for( i=0; i<icyc; i++ )
1102 for( j=0; j<jcyc; j++ )
1104 mseq2[j] = mseq[icyc+j];
1109 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1113 if( commonAlloc1 && commonAlloc2 )
1115 FreeIntMtx( commonIP );
1118 ll1 = MAX( orlgth1, commonAlloc1 );
1119 ll2 = MAX( orlgth2, commonAlloc2 );
1122 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1125 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
1128 fprintf( stderr, "succeeded\n\n" );
1139 for( i=0; i<icyc; i++ )
1141 fprintf( stderr, "## totaleff = %f\n", t );
1145 cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1146 cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1150 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
1151 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
1152 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap1 );
1153 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
1157 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1158 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1159 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1160 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1163 for( i=0; i<lgth1; i++ )
1165 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
1166 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
1168 for( i=0; i<lgth2; i++ )
1170 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
1171 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
1174 for( i=0; i<lgth1; i++ )
1175 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1181 match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1183 imp_match_out_vead_tate( initverticalw, 0, lgth1 ); // 060306
1185 match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1187 imp_match_out_vead( currentw, 0, lgth2 ); // 060306
1188 #if 0 // -> tbfast.c
1190 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1196 for( i=1; i<lgth1+1; i++ )
1198 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1200 for( j=1; j<lgth2+1; j++ )
1202 currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1208 for( j=1; j<lgth2+1; j++ )
1209 currentw[j] -= offset * j / 2.0;
1210 for( i=1; i<lgth1+1; i++ )
1211 initverticalw[i] -= offset * i / 2.0;
1215 for( j=1; j<lgth2+1; ++j )
1217 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1220 lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari
1222 lastverticalw[0] = currentw[lgth2-1];
1224 if( tailgp ) lasti = lgth1+1; else lasti = lgth1;
1227 fprintf( stderr, "currentw = \n" );
1228 for( i=0; i<lgth1+1; i++ )
1230 fprintf( stderr, "%5.2f ", currentw[i] );
1232 fprintf( stderr, "\n" );
1233 fprintf( stderr, "initverticalw = \n" );
1234 for( i=0; i<lgth2+1; i++ )
1236 fprintf( stderr, "%5.2f ", initverticalw[i] );
1238 fprintf( stderr, "\n" );
1239 fprintf( stderr, "fcgp\n" );
1240 for( i=0; i<lgth1; i++ )
1241 fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1242 for( i=0; i<lgth2; i++ )
1243 fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1246 for( i=1; i<lasti; i++ )
1249 #ifdef enablemultithread
1250 // fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
1251 if( chudanpt && *chudanpt != chudanref )
1253 // fprintf( stderr, "\n\n## CHUUDAN!!! S\n" );
1259 previousw = currentw;
1262 previousw[0] = initverticalw[i-1];
1264 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1266 fprintf( stderr, "\n" );
1267 fprintf( stderr, "i=%d\n", i );
1268 fprintf( stderr, "currentw = \n" );
1269 for( j=0; j<lgth2; j++ )
1271 fprintf( stderr, "%5.2f ", currentw[j] );
1273 fprintf( stderr, "\n" );
1277 // fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1279 imp_match_out_vead( currentw, i, lgth2 );
1281 imp_match_out_vead( currentw, i, lgth2 );
1285 fprintf( stderr, "\n" );
1286 fprintf( stderr, "i=%d\n", i );
1287 fprintf( stderr, "currentw = \n" );
1288 for( j=0; j<lgth2; j++ )
1290 fprintf( stderr, "%5.2f ", currentw[j] );
1292 fprintf( stderr, "\n" );
1294 currentw[0] = initverticalw[i];
1297 mi = previousw[0] + ogcp2[1]; mpi = 0;
1301 curpt = currentw + 1;
1304 ogcp2pt = ogcp2 + 1;
1305 fgcp1va = fgcp1[i-1];
1308 for( j=1; j<lastj; j++ )
1310 #ifdef xxxenablemultithread
1311 // fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
1312 if( chudanpt && *chudanpt != chudanref )
1314 // fprintf( stderr, "\n\n## CHUUDAN!!! S\n" );
1323 fprintf( stderr, "%5.0f->", wm );
1326 fprintf( stderr, "%5.0f?", g );
1328 if( (g=mi+*fgcp2pt) > wm )
1331 *ijppt = -( j - mpi );
1333 if( (g=*prept+*ogcp2pt) >= mi )
1343 fprintf( stderr, "%5.0f?", g );
1345 if( (g=*mjpt+fgcp1va) > wm )
1348 *ijppt = +( i - *mpjpt );
1350 if( (g=*prept+ogcp1va) >= *mjpt )
1356 m[j] += fpenalty_ex;
1360 fprintf( stderr, "%5.0f ", wm );
1370 lastverticalw[i] = currentw[lgth2-1];
1373 // fprintf( stderr, "wm = %f\n", wm );
1378 for( j=1; j<lgth2+1; j++ )
1379 currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1380 for( i=1; i<lgth1+1; i++ )
1381 lastverticalw[i] -= offset * ( lgth1 - i / 2.0);
1386 fprintf( stderr, "\n" );
1387 for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1388 fprintf( stderr, "#####\n" );
1389 for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1390 fprintf( stderr, "====>" );
1391 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1392 for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1396 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1399 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, tailgp );
1401 // fprintf( stderr, "### impmatch = %f\n", *impmatch );
1403 resultlen = strlen( mseq1[0] );
1404 if( alloclen < resultlen || resultlen > N )
1406 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1407 ErrorExit( "LENGTH OVER!\n" );
1411 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1412 for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1414 fprintf( stderr, "\n" );
1415 for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1416 fprintf( stderr, "#####\n" );
1417 for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1420 // fprintf( stderr, "wm = %f\n", wm );
1425 float A__align_gapmap( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, int *gapmap1, int *gapmap2 )
1426 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1430 int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
1433 float wm = 0.0; /* int ?????? */
1435 float *currentw, *previousw;
1436 // float fpenalty = (float)penalty;
1438 float fpenalty_ex = (float)penalty_ex;
1443 float *mjpt, *prept, *curpt;
1446 static TLS float mi, *m;
1447 static TLS int **ijp;
1448 static TLS int mpi, *mp;
1449 static TLS float *w1, *w2;
1450 static TLS float *match;
1451 static TLS float *initverticalw; /* kufuu sureba iranai */
1452 static TLS float *lastverticalw; /* kufuu sureba iranai */
1453 static TLS char **mseq1;
1454 static TLS char **mseq2;
1455 static TLS char **mseq;
1456 static TLS float *ogcp1;
1457 static TLS float *ogcp2;
1458 static TLS float *fgcp1;
1459 static TLS float *fgcp2;
1460 static TLS float **cpmx1;
1461 static TLS float **cpmx2;
1462 static TLS int **intwork;
1463 static TLS float **floatwork;
1464 static TLS int orlgth1 = 0, orlgth2 = 0;
1472 fprintf( stderr, "eff in SA+++align\n" );
1473 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1477 mseq1 = AllocateCharMtx( njob, 0 );
1478 mseq2 = AllocateCharMtx( njob, 0 );
1482 lgth1 = strlen( seq1[0] );
1483 lgth2 = strlen( seq2[0] );
1485 if( lgth1 == 0 || lgth2 == 0 )
1487 fprintf( stderr, "WARNING (Aalign_gapmap): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
1491 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
1495 if( orlgth1 > 0 && orlgth2 > 0 )
1499 FreeFloatVec( match );
1500 FreeFloatVec( initverticalw );
1501 FreeFloatVec( lastverticalw );
1506 FreeCharMtx( mseq );
1508 FreeFloatVec( ogcp1 );
1509 FreeFloatVec( ogcp2 );
1510 FreeFloatVec( fgcp1 );
1511 FreeFloatVec( fgcp2 );
1514 FreeFloatMtx( cpmx1 );
1515 FreeFloatMtx( cpmx2 );
1517 FreeFloatMtx( floatwork );
1518 FreeIntMtx( intwork );
1521 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
1522 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
1525 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1528 w1 = AllocateFloatVec( ll2+2 );
1529 w2 = AllocateFloatVec( ll2+2 );
1530 match = AllocateFloatVec( ll2+2 );
1532 initverticalw = AllocateFloatVec( ll1+2 );
1533 lastverticalw = AllocateFloatVec( ll1+2 );
1535 m = AllocateFloatVec( ll2+2 );
1536 mp = AllocateIntVec( ll2+2 );
1538 mseq = AllocateCharMtx( njob, ll1+ll2 );
1540 ogcp1 = AllocateFloatVec( ll1+2 );
1541 ogcp2 = AllocateFloatVec( ll2+2 );
1542 fgcp1 = AllocateFloatVec( ll1+2 );
1543 fgcp2 = AllocateFloatVec( ll2+2 );
1545 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1546 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1549 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
1550 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
1552 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
1553 intwork = AllocateIntMtx( 27, MAX( ll1, ll2 )+2 );
1557 fprintf( stderr, "succeeded\n" );
1560 orlgth1 = ll1 - 100;
1561 orlgth2 = ll2 - 100;
1565 for( i=0; i<icyc; i++ )
1570 for( j=0; j<jcyc; j++ )
1572 mseq2[j] = mseq[icyc+j];
1577 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1581 if( commonAlloc1 && commonAlloc2 )
1583 FreeIntMtx( commonIP );
1586 ll1 = MAX( orlgth1, commonAlloc1 );
1587 ll2 = MAX( orlgth2, commonAlloc2 );
1590 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1593 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
1596 fprintf( stderr, "succeeded\n\n" );
1604 cpmx_calc_new( seq1, cpmx1, eff1, strlen( seq1[0] ), icyc );
1605 cpmx_calc_new( seq2, cpmx2, eff2, strlen( seq2[0] ), jcyc );
1607 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1608 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1609 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1610 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1612 for( i=0; i<lgth1; i++ )
1614 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] );
1615 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] );
1617 for( i=0; i<lgth2; i++ )
1619 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] );
1620 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] );
1623 for( i=0; i<lgth1; i++ )
1624 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1631 match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1633 imp_match_out_vead_tate_gapmap( initverticalw, gapmap2[0], lgth1, gapmap1 ); // 060306
1636 match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1638 imp_match_out_vead_gapmap( currentw, gapmap1[0], lgth2, gapmap2 ); // 060306
1639 #if 0 // -> tbfast.c
1641 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1645 if( 1 ) // tsuneni outgap=1
1647 for( i=1; i<lgth1+1; i++ )
1649 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1651 for( j=1; j<lgth2+1; j++ )
1653 currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1659 for( j=1; j<lgth2+1; j++ )
1660 currentw[j] -= offset * j / 2.0;
1661 for( i=1; i<lgth1+1; i++ )
1662 initverticalw[i] -= offset * i / 2.0;
1666 for( j=1; j<lgth2+1; ++j )
1668 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1672 lastverticalw[0] = 0.0;
1674 lastverticalw[0] = currentw[lgth2-1];
1676 if( 1 ) lasti = lgth1+1; else lasti = lgth1; // tsuneni outgap=1
1679 fprintf( stderr, "currentw = \n" );
1680 for( i=0; i<lgth1+1; i++ )
1682 fprintf( stderr, "%5.2f ", currentw[i] );
1684 fprintf( stderr, "\n" );
1685 fprintf( stderr, "initverticalw = \n" );
1686 for( i=0; i<lgth2+1; i++ )
1688 fprintf( stderr, "%5.2f ", initverticalw[i] );
1690 fprintf( stderr, "\n" );
1691 fprintf( stderr, "fcgp\n" );
1692 for( i=0; i<lgth1; i++ )
1693 fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1694 for( i=0; i<lgth2; i++ )
1695 fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1698 for( i=1; i<lasti; i++ )
1701 previousw = currentw;
1704 previousw[0] = initverticalw[i-1];
1706 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1708 fprintf( stderr, "\n" );
1709 fprintf( stderr, "i=%d\n", i );
1710 fprintf( stderr, "currentw = \n" );
1711 for( j=0; j<lgth2; j++ )
1713 fprintf( stderr, "%5.2f ", currentw[j] );
1715 fprintf( stderr, "\n" );
1719 // fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1721 imp_match_out_vead( currentw, i, lgth2 );
1723 imp_match_out_vead_gapmap( currentw, gapmap1[i], lgth2, gapmap2 );
1727 fprintf( stderr, "\n" );
1728 fprintf( stderr, "i=%d\n", i );
1729 fprintf( stderr, "currentw = \n" );
1730 for( j=0; j<lgth2; j++ )
1732 fprintf( stderr, "%5.2f ", currentw[j] );
1734 fprintf( stderr, "\n" );
1736 currentw[0] = initverticalw[i];
1739 mi = previousw[0] + ogcp2[1]; mpi = 0;
1744 curpt = currentw + 1;
1747 ogcp2pt = ogcp2 + 1;
1748 fgcp1va = fgcp1[i-1];
1751 for( j=1; j<lastj; j++ )
1757 fprintf( stderr, "%5.0f->", wm );
1761 fprintf( stderr, "%5.0f?", g );
1766 *ijppt = -( j - mpi );
1768 g = *prept + *ogcp2pt;
1778 g = *mjpt + fgcp1va;
1780 fprintf( stderr, "%5.0f?", g );
1785 *ijppt = +( i - *mpjpt );
1787 g = *prept + ogcp1va;
1794 m[j] += fpenalty_ex;
1798 fprintf( stderr, "%5.0f ", wm );
1808 lastverticalw[i] = currentw[lgth2-1];
1814 for( j=1; j<lgth2+1; j++ )
1815 currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1816 for( i=1; i<lgth1+1; i++ )
1817 lastverticalw[i] -= offset * ( lgth1 - i / 2.0);
1822 fprintf( stderr, "\n" );
1823 for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1824 fprintf( stderr, "#####\n" );
1825 for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1826 fprintf( stderr, "====>" );
1827 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1828 for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1832 Atracking_localhom_gapmap( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, gapmap1, gapmap2 );
1835 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, 1 );
1837 // fprintf( stderr, "### impmatch = %f\n", *impmatch );
1839 resultlen = strlen( mseq1[0] );
1840 if( alloclen < resultlen || resultlen > N )
1842 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1843 ErrorExit( "LENGTH OVER!\n" );
1847 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1848 for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1850 fprintf( stderr, "\n" );
1851 for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1852 fprintf( stderr, "#####\n" );
1853 for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1856 // fprintf( stderr, "wm = %f\n", wm );