7 #define USE_PENALTY_EX 0
9 static void OpeningGapCount( double *ogcp, int clus, char **seq, double *eff )
12 int len = strlen( seq[0] );
13 double totaleff = 0.0;
15 for( i=0; i<len; i++ ) ogcp[i] = 0.0;
16 for( j=0; j<clus; j++ )
19 for( i=0; i<len; i++ )
22 gc = ( seq[j][i] == '-' );
24 if( !gb * gc ) ogcp[i] += eff[j];
29 for( i=0; i<len; i++ )
33 static void FinalGapCount( double *fgcp, int clus, char **seq, double *eff )
36 int len = strlen( seq[0] );
37 double totaleff = 0.0;
39 for( i=0; i<len; i++ ) fgcp[i] = 0.0;
40 for( j=0; j<clus; j++ )
42 gc = ( seq[j][0] == '-' );
43 for( i=1; i<len+1; i++ )
46 gc = ( seq[j][i] == '-' );
48 if( gb * !gc ) fgcp[i-1] += eff[j];
53 for( i=0; i<len; i++ )
56 static float **impmtx = NULL;
57 float imp_match_out_sc( int i1, int j1 )
59 // fprintf( stderr, "imp+match = %f\n", impmtx[i1][j1] * fastathreshold );
60 // fprintf( stderr, "val = %f\n", impmtx[i1][j1] );
61 return( impmtx[i1][j1] );
63 static void imp_match_out_vead_gapmap( float *imp, int i1, int lgth2, int *gapmap2 )
66 float *pt = impmtx[i1];
67 for( j=0; j<lgth2; j++ )
68 *imp++ += pt[gapmap2[j]];
70 static void imp_match_out_vead( float *imp, int i1, int lgth2 )
73 float *pt = impmtx[i1];
74 for( j=0; j<lgth2; j++ )
78 void imp_match_init_strict( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom, int forscore )
80 int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2;
81 static int impalloclen = 0;
85 static char *nocount1 = NULL;
86 static char *nocount2 = NULL;
89 if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
91 if( impmtx ) FreeFloatMtx( impmtx );
92 if( nocount1 ) free( nocount1 );
93 if( nocount2 ) free( nocount2 );
94 impalloclen = MAX( lgth1, lgth2 ) + 2;
95 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
96 nocount1 = AllocateCharVec( impalloclen );
97 nocount2 = AllocateCharVec( impalloclen );
100 for( i=0; i<lgth1; i++ )
102 for( j=0; j<clus1; j++ )
103 if( seq1[j][i] == '-' ) break;
104 if( j != clus1 ) nocount1[i] = 1;
105 else nocount1[i] = 0;
107 for( i=0; i<lgth2; i++ )
109 for( j=0; j<clus2; j++ )
110 if( seq2[j][i] == '-' ) break;
111 if( j != clus2 ) nocount2[i] = 1;
112 else nocount2[i] = 0;
116 fprintf( stderr, "nocount2 =\n" );
117 for( i = 0; i<impalloclen; i++ )
119 fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
126 fprintf( stderr, "eff1 in _init_strict = \n" );
127 for( i=0; i<clus1; i++ )
128 fprintf( stderr, "eff1[] = %f\n", eff1[i] );
129 for( i=0; i<clus2; i++ )
130 fprintf( stderr, "eff2[] = %f\n", eff2[i] );
133 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
135 effijx = fastathreshold;
136 for( i=0; i<clus1; i++ )
138 for( j=0; j<clus2; j++ )
140 effij = eff1[i] * eff2[j] * effijx;
141 tmpptr = localhom[i][j];
144 // fprintf( stderr, "start1 = %d\n", tmpptr->start1 );
145 // fprintf( stderr, "end1 = %d\n", tmpptr->end1 );
146 // fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] );
147 // fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] );
152 if( *pt++ != '-' ) tmpint++;
153 if( tmpint == tmpptr->start1 ) break;
155 start1 = pt - seq1[i] - 1;
159 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
160 if( *pt++ != '-' ) tmpint++;
161 if( tmpint == tmpptr->end1 ) break;
163 end1 = pt - seq1[i] - 1;
169 if( *pt++ != '-' ) tmpint++;
170 if( tmpint == tmpptr->start2 ) break;
172 start2 = pt - seq2[j] - 1;
175 if( *pt++ != '-' ) tmpint++;
176 if( tmpint == tmpptr->end2 ) break;
178 end2 = pt - seq2[j] - 1;
179 // 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] );
180 // fprintf( stderr, "step 0\n" );
181 if( end1 - start1 != end2 - start2 )
183 // fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
187 k1 = start1; k2 = start2;
190 while( *pt1 && *pt2 )
192 if( *pt1 != '-' && *pt2 != '-' )
194 // ½Å¤ß¤òÆó½Å¤Ë¤«¤±¤Ê¤¤¤è¤¦¤ËÃí°Õ¤·¤Æ²¼¤µ¤¤¡£
195 // impmtx[k1][k2] += tmpptr->wimportance * fastathreshold;
196 impmtx[k1][k2] += tmpptr->importance * effij;
197 // fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
201 else if( *pt1 != '-' && *pt2 == '-' )
203 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
206 else if( *pt1 == '-' && *pt2 != '-' )
208 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
211 else if( *pt1 == '-' && *pt2 == '-' )
213 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
217 if( k1 > end1 || k2 > end2 ) break;
220 while( k1 <= end1 && k2 <= end2 )
222 fprintf( stderr, "k1,k2=%d,%d - ", k1, k2 );
223 if( !nocount1[k1] && !nocount2[k2] )
225 impmtx[k1][k2] += tmpptr->wimportance * eff1[i] * eff2[j] * fastathreshold;
226 fprintf( stderr, "marked\n" );
229 fprintf( stderr, "no count\n" );
233 tmpptr = tmpptr->next;
238 fprintf( stderr, "impmtx = \n" );
239 for( k2=0; k2<lgth2; k2++ )
240 fprintf( stderr, "%6.3f ", (double)k2 );
241 fprintf( stderr, "\n" );
242 for( k1=0; k1<lgth1; k1++ )
244 fprintf( stderr, "%d", k1 );
245 for( k2=0; k2<lgth2; k2++ )
246 fprintf( stderr, "%2.1f ", impmtx[k1][k2] );
247 fprintf( stderr, "\n" );
252 void imp_match_init( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom )
254 int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2;
255 static int impalloclen = 0;
258 static char *nocount1 = NULL;
259 static char *nocount2 = NULL;
261 if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
263 if( impmtx ) FreeFloatMtx( impmtx );
264 if( nocount1 ) free( nocount1 );
265 if( nocount2 ) free( nocount2 );
266 impalloclen = MAX( lgth1, lgth2 ) + 2;
267 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
268 nocount1 = AllocateCharVec( impalloclen );
269 nocount2 = AllocateCharVec( impalloclen );
272 for( i=0; i<lgth1; i++ )
274 for( j=0; j<clus1; j++ )
275 if( seq1[j][i] == '-' ) break;
276 if( j != clus1 ) nocount1[i] = 1;
277 else nocount1[i] = 0;
279 for( i=0; i<lgth2; i++ )
281 for( j=0; j<clus2; j++ )
282 if( seq2[j][i] == '-' ) break;
283 if( j != clus2 ) nocount2[i] = 1;
284 else nocount2[i] = 0;
288 fprintf( stderr, "nocount2 =\n" );
289 for( i = 0; i<impalloclen; i++ )
291 fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
295 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
297 for( i=0; i<clus1; i++ )
299 fprintf( stderr, "i = %d, seq1 = %s\n", i, seq1[i] );
300 for( j=0; j<clus2; j++ )
302 fprintf( stderr, "start1 = %d\n", localhom[i][j]->start1 );
303 fprintf( stderr, "end1 = %d\n", localhom[i][j]->end1 );
304 fprintf( stderr, "j = %d, seq2 = %s\n", j, seq2[j] );
309 if( *pt++ != '-' ) tmpint++;
310 if( tmpint == localhom[i][j]->start1 ) break;
312 start1 = pt - seq1[i] - 1;
316 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, localhom[i][j].end1, pt-seq1[i] );
317 if( *pt++ != '-' ) tmpint++;
318 if( tmpint == localhom[i][j]->end1 ) break;
320 end1 = pt - seq1[i] - 1;
326 if( *pt++ != '-' ) tmpint++;
327 if( tmpint == localhom[i][j]->start2 ) break;
329 start2 = pt - seq2[j] - 1;
332 if( *pt++ != '-' ) tmpint++;
333 if( tmpint == localhom[i][j]->end2 ) break;
335 end2 = pt - seq2[j] - 1;
336 // fprintf( stderr, "start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
339 fprintf( stderr, "step 0\n" );
340 while( k1 <= end1 && k2 <= end2 )
343 if( !nocount1[k1] && !nocount2[k2] )
344 impmtx[k1][k2] += localhom[i][j].wimportance * eff1[i] * eff2[j];
347 if( !nocount1[k1] && !nocount2[k2] )
348 impmtx[k1][k2] += localhom[i][j]->wimportance * eff1[i] * eff2[j];
353 dif = ( end1 - start1 ) - ( end2 - start2 );
354 fprintf( stderr, "dif = %d\n", dif );
359 fprintf( stderr, "dif = %d\n", dif );
362 while( k1 <= end1 && k2 <= end2 )
364 if( 0 <= k2 && start2 <= k2 && !nocount1[k1] && !nocount2[k2] )
365 impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
379 if( k1 >= 0 && k1 >= start1 && !nocount1[k1] && !nocount2[k2] )
380 impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
389 fprintf( stderr, "impmtx = \n" );
390 for( k2=0; k2<lgth2; k2++ )
391 fprintf( stderr, "%6.3f ", (double)k2 );
392 fprintf( stderr, "\n" );
393 for( k1=0; k1<lgth1; k1++ )
395 fprintf( stderr, "%d", k1 );
396 for( k2=0; k2<lgth2; k2++ )
397 fprintf( stderr, "%6.3f ", impmtx[k1][k2] );
398 fprintf( stderr, "\n" );
403 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
407 float **cpmxpd = floatwork;
408 int **cpmxpdn = intwork;
413 for( j=0; j<lgth2; j++ )
416 for( l=0; l<26; l++ )
420 cpmxpd[count][j] = cpmx2[l][j];
421 cpmxpdn[count][j] = l;
425 cpmxpdn[count][j] = -1;
429 for( l=0; l<26; l++ )
432 for( k=0; k<26; k++ )
433 scarr[l] += n_dis[k][l] * cpmx1[k][i1];
435 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤
\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
437 float *fpt, **fptpt, *fpt2;
445 ipt=*iptpt,fpt=*fptpt;
447 *fpt2 += scarr[*ipt++] * *fpt++;
448 fpt2++,iptpt++,fptpt++;
452 for( j=0; j<lgth2; j++ )
455 for( k=0; cpmxpdn[k][j]>-1; k++ )
456 match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
461 static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
462 char **seq1, char **seq2,
463 char **mseq1, char **mseq2,
464 float **cpmx1, float **cpmx2,
465 short **ijp, int icyc, int jcyc )
467 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
470 lgth1 = strlen( seq1[0] );
471 lgth2 = strlen( seq2[0] );
474 for( i=0; i<lgth1; i++ )
476 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
484 wm = lastverticalw[0];
485 for( i=0; i<lgth1; i++ )
487 if( lastverticalw[i] >= wm )
489 wm = lastverticalw[i];
490 iin = i; jin = lgth2-1;
491 ijp[lgth1][lgth2] = +( lgth1 - i );
494 for( j=0; j<lgth2; j++ )
496 if( lasthorizontalw[j] >= wm )
498 wm = lasthorizontalw[j];
499 iin = lgth1-1; jin = j;
500 ijp[lgth1][lgth2] = -( lgth2 - j );
505 for( i=0; i<lgth1+1; i++ )
509 for( j=0; j<lgth2+1; j++ )
511 ijp[0][j] = -( j + 1 );
514 for( i=0; i<icyc; i++ )
516 mseq1[i] += lgth1+lgth2;
519 for( j=0; j<jcyc; j++ )
521 mseq2[j] += lgth1+lgth2;
524 iin = lgth1; jin = lgth2;
526 for( k=0; k<=lgth1+lgth2; k++ )
528 if( ijp[iin][jin] < 0 )
530 ifi = iin-1; jfi = jin+ijp[iin][jin];
532 else if( ijp[iin][jin] > 0 )
534 ifi = iin-ijp[iin][jin]; jfi = jin-1;
538 ifi = iin-1; jfi = jin-1;
543 for( i=0; i<icyc; i++ )
544 *--mseq1[i] = seq1[i][ifi+l];
545 for( j=0; j<jcyc; j++ )
552 for( i=0; i<icyc; i++ )
554 for( j=0; j<jcyc; j++ )
555 *--mseq2[j] = seq2[j][jfi+l];
558 if( iin == lgth1 || jin == lgth2 )
562 *impwmpt += imp_match_out_sc( iin, jin );
564 // fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
566 if( iin <= 0 || jin <= 0 ) break;
567 for( i=0; i<icyc; i++ )
568 *--mseq1[i] = seq1[i][ifi];
569 for( j=0; j<jcyc; j++ )
570 *--mseq2[j] = seq2[j][jfi];
572 iin = ifi; jin = jfi;
575 static void Atracking_localhom_gapmap( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
576 char **seq1, char **seq2,
577 char **mseq1, char **mseq2,
578 float **cpmx1, float **cpmx2,
579 short **ijp, int icyc, int jcyc,
580 int *gapmap1, int *gapmap2 )
582 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
585 lgth1 = strlen( seq1[0] );
586 lgth2 = strlen( seq2[0] );
589 for( i=0; i<lgth1; i++ )
591 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
599 wm = lastverticalw[0];
600 for( i=0; i<lgth1; i++ )
602 if( lastverticalw[i] >= wm )
604 wm = lastverticalw[i];
605 iin = i; jin = lgth2-1;
606 ijp[lgth1][lgth2] = +( lgth1 - i );
609 for( j=0; j<lgth2; j++ )
611 if( lasthorizontalw[j] >= wm )
613 wm = lasthorizontalw[j];
614 iin = lgth1-1; jin = j;
615 ijp[lgth1][lgth2] = -( lgth2 - j );
620 for( i=0; i<lgth1+1; i++ )
624 for( j=0; j<lgth2+1; j++ )
626 ijp[0][j] = -( j + 1 );
629 for( i=0; i<icyc; i++ )
631 mseq1[i] += lgth1+lgth2;
634 for( j=0; j<jcyc; j++ )
636 mseq2[j] += lgth1+lgth2;
639 iin = lgth1; jin = lgth2;
641 for( k=0; k<=lgth1+lgth2; k++ )
643 if( ijp[iin][jin] < 0 )
645 ifi = iin-1; jfi = jin+ijp[iin][jin];
647 else if( ijp[iin][jin] > 0 )
649 ifi = iin-ijp[iin][jin]; jfi = jin-1;
653 ifi = iin-1; jfi = jin-1;
658 for( i=0; i<icyc; i++ )
659 *--mseq1[i] = seq1[i][ifi+l];
660 for( j=0; j<jcyc; j++ )
667 for( i=0; i<icyc; i++ )
669 for( j=0; j<jcyc; j++ )
670 *--mseq2[j] = seq2[j][jfi+l];
673 if( iin == lgth1 || jin == lgth2 )
677 *impwmpt += imp_match_out_sc( gapmap1[iin], gapmap2[jin] );
679 // fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
681 if( iin <= 0 || jin <= 0 ) break;
682 for( i=0; i<icyc; i++ )
683 *--mseq1[i] = seq1[i][ifi];
684 for( j=0; j<jcyc; j++ )
685 *--mseq2[j] = seq2[j][jfi];
687 iin = ifi; jin = jfi;
690 static float Atracking( float *lasthorizontalw, float *lastverticalw,
691 char **seq1, char **seq2,
692 char **mseq1, char **mseq2,
693 float **cpmx1, float **cpmx2,
694 short **ijp, int icyc, int jcyc )
696 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
699 lgth1 = strlen( seq1[0] );
700 lgth2 = strlen( seq2[0] );
703 for( i=0; i<lgth1; i++ )
705 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
713 wm = lastverticalw[0];
714 for( i=0; i<lgth1; i++ )
716 if( lastverticalw[i] >= wm )
718 wm = lastverticalw[i];
719 iin = i; jin = lgth2-1;
720 ijp[lgth1][lgth2] = +( lgth1 - i );
723 for( j=0; j<lgth2; j++ )
725 if( lasthorizontalw[j] >= wm )
727 wm = lasthorizontalw[j];
728 iin = lgth1-1; jin = j;
729 ijp[lgth1][lgth2] = -( lgth2 - j );
734 for( i=0; i<lgth1+1; i++ )
738 for( j=0; j<lgth2+1; j++ )
740 ijp[0][j] = -( j + 1 );
743 for( i=0; i<icyc; i++ )
745 mseq1[i] += lgth1+lgth2;
748 for( j=0; j<jcyc; j++ )
750 mseq2[j] += lgth1+lgth2;
753 iin = lgth1; jin = lgth2;
754 for( k=0; k<=lgth1+lgth2; k++ )
756 if( ijp[iin][jin] < 0 )
758 ifi = iin-1; jfi = jin+ijp[iin][jin];
760 else if( ijp[iin][jin] > 0 )
762 ifi = iin-ijp[iin][jin]; jfi = jin-1;
766 ifi = iin-1; jfi = jin-1;
771 for( i=0; i<icyc; i++ )
772 *--mseq1[i] = seq1[i][ifi+l];
773 for( j=0; j<jcyc; j++ )
780 for( i=0; i<icyc; i++ )
782 for( j=0; j<jcyc; j++ )
783 *--mseq2[j] = seq2[j][jfi+l];
786 if( iin <= 0 || jin <= 0 ) break;
787 for( i=0; i<icyc; i++ )
788 *--mseq1[i] = seq1[i][ifi];
789 for( j=0; j<jcyc; j++ )
790 *--mseq2[j] = seq2[j][jfi];
792 iin = ifi; jin = jfi;
797 float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch )
798 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
802 int lasti; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
805 float wm; /* int ?????? */
807 float *currentw, *previousw;
808 float fpenalty = (float)penalty;
809 float fpenalty_ex = (float)penalty_ex;
813 float *mjpt, *prept, *curpt;
819 static float *w1, *w2;
821 static float *initverticalw; /* kufuu sureba iranai */
822 static float *lastverticalw; /* kufuu sureba iranai */
826 static double *ogcp1;
827 static double *ogcp2;
828 static double *fgcp1;
829 static double *fgcp2;
830 static float **cpmx1;
831 static float **cpmx2;
832 static int **intwork;
833 static float **floatwork;
834 static int orlgth1 = 0, orlgth2 = 0;
838 fprintf( stderr, "eff in SA+++align\n" );
839 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
843 mseq1 = AllocateCharMtx( njob, 0 );
844 mseq2 = AllocateCharMtx( njob, 0 );
848 lgth1 = strlen( seq1[0] );
849 lgth2 = strlen( seq2[0] );
851 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
855 if( orlgth1 > 0 && orlgth2 > 0 )
859 FreeFloatVec( match );
860 FreeFloatVec( initverticalw );
861 FreeFloatVec( lastverticalw );
868 FreeDoubleVec( ogcp1 );
869 FreeDoubleVec( ogcp2 );
870 FreeDoubleVec( fgcp1 );
871 FreeDoubleVec( fgcp2 );
874 FreeFloatMtx( cpmx1 );
875 FreeFloatMtx( cpmx2 );
877 FreeFloatMtx( floatwork );
878 FreeIntMtx( intwork );
881 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
882 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
885 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
888 w1 = AllocateFloatVec( ll2+2 );
889 w2 = AllocateFloatVec( ll2+2 );
890 match = AllocateFloatVec( ll2+2 );
892 initverticalw = AllocateFloatVec( ll1+2 );
893 lastverticalw = AllocateFloatVec( ll1+2 );
895 m = AllocateFloatVec( ll2+2 );
896 mp = AllocateIntVec( ll2+2 );
898 mseq = AllocateCharMtx( njob, ll1+ll2 );
900 ogcp1 = AllocateDoubleVec( ll1+2 );
901 ogcp2 = AllocateDoubleVec( ll2+2 );
902 fgcp1 = AllocateDoubleVec( ll1+2 );
903 fgcp2 = AllocateDoubleVec( ll2+2 );
905 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
906 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
908 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
909 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
912 fprintf( stderr, "succeeded\n" );
920 for( i=0; i<icyc; i++ ) mseq1[i] = mseq[i];
921 for( j=0; j<jcyc; j++ ) mseq2[j] = mseq[icyc+j];
924 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
928 if( commonAlloc1 && commonAlloc2 )
930 FreeShortMtx( commonIP );
933 ll1 = MAX( orlgth1, commonAlloc1 );
934 ll2 = MAX( orlgth2, commonAlloc2 );
937 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
940 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
943 fprintf( stderr, "succeeded\n\n" );
954 for( i=0; i<icyc; i++ )
956 fprintf( stderr, "## totaleff = %f\n", t );
959 cpmx_calc_new( seq1, cpmx1, eff1, strlen( seq1[0] ), icyc );
960 cpmx_calc_new( seq2, cpmx2, eff2, strlen( seq2[0] ), jcyc );
962 OpeningGapCount( ogcp1, icyc, seq1, eff1 );
963 OpeningGapCount( ogcp2, jcyc, seq2, eff2 );
964 FinalGapCount( fgcp1, icyc, seq1, eff1 );
965 FinalGapCount( fgcp2, jcyc, seq2, eff2 );
967 for( i=0; i<lgth1; i++ )
969 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] );
970 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] );
972 for( i=0; i<lgth2; i++ )
974 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] );
975 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] );
978 for( i=0; i<lgth1; i++ )
979 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
985 match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
987 match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
990 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
996 for( i=1; i<lgth1+1; i++ )
998 initverticalw[i] += fpenalty * ( ogcp1[0] + fgcp1[i-1] ) ;
1000 for( j=1; j<lgth2+1; j++ )
1002 currentw[j] += fpenalty * ( ogcp2[0] + fgcp2[j-1] ) ;
1008 for( j=1; j<lgth2+1; j++ )
1009 currentw[j] -= offset * j / 2.0;
1010 for( i=1; i<lgth1+1; i++ )
1011 initverticalw[i] -= offset * i / 2.0;
1015 for( j=1; j<lgth2+1; ++j )
1017 m[j] = currentw[j-1] + fpenalty * ogcp1[1]; mp[j] = 0;
1020 lastverticalw[0] = currentw[lgth2-1];
1022 if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1025 fprintf( stderr, "currentw = \n" );
1026 for( i=0; i<lgth1+1; i++ )
1028 fprintf( stderr, "%5.2f ", currentw[i] );
1030 fprintf( stderr, "\n" );
1031 fprintf( stderr, "initverticalw = \n" );
1032 for( i=0; i<lgth2+1; i++ )
1034 fprintf( stderr, "%5.2f ", initverticalw[i] );
1036 fprintf( stderr, "\n" );
1037 fprintf( stderr, "fcgp\n" );
1038 for( i=0; i<lgth1; i++ )
1039 fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1040 for( i=0; i<lgth2; i++ )
1041 fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1044 for( i=1; i<lasti; i++ )
1047 previousw = currentw;
1050 previousw[0] = initverticalw[i-1];
1052 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1054 fprintf( stderr, "\n" );
1055 fprintf( stderr, "i=%d\n", i );
1056 fprintf( stderr, "currentw = \n" );
1057 for( j=0; j<lgth2; j++ )
1059 fprintf( stderr, "%5.2f ", currentw[j] );
1061 fprintf( stderr, "\n" );
1065 // fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1067 imp_match_out_vead( currentw, i, lgth2 );
1069 imp_match_out_vead( currentw, i, lgth2 );
1073 fprintf( stderr, "\n" );
1074 fprintf( stderr, "i=%d\n", i );
1075 fprintf( stderr, "currentw = \n" );
1076 for( j=0; j<lgth2; j++ )
1078 fprintf( stderr, "%5.2f ", currentw[j] );
1080 fprintf( stderr, "\n" );
1082 currentw[0] = initverticalw[i];
1086 mi = previousw[0] + fpenalty * ogcp2[1]; mpi = 0;
1091 curpt = currentw + 1;
1093 for( j=1; j<lgth2+1; j++ )
1099 fprintf( stderr, "%5.0f->", wm );
1101 g = mi + fpenalty * fgcp2[j-1];
1103 fprintf( stderr, "%5.0f?", g );
1108 *ijppt = -( j - mpi );
1110 g = *prept + fpenalty * ogcp2[j];
1120 g = *mjpt + fpenalty * fgcp1[i-1];
1122 fprintf( stderr, "%5.0f?", g );
1127 *ijppt = +( i - *mpjpt );
1129 g = *prept + fpenalty * ogcp1[i];
1136 m[j] += fpenalty_ex;
1140 fprintf( stderr, "%5.0f ", wm );
1149 lastverticalw[i] = currentw[lgth2-1];
1151 floatncpy( previousw, currentw, lgth2+1 );
1152 previousw[0] = initverticalw[i-1];
1155 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1156 currentw[0] = initverticalw[i];
1158 mi = previousw[0] + fpenalty * ogcp2[1]; mpi = 0;
1159 for( j=1; j<lgth2+1; j++ )
1161 wm = previousw[j-1];
1164 g = fpenalty * fgcp2[j-1];
1168 ijp[i][j] = -( j - mpi );
1170 g = fpenalty * ogcp2[j];
1171 if( mi <= previousw[j-1] + g )
1173 mi = previousw[j-1] + g;
1177 g = fpenalty * fgcp1[i-1];
1181 ijp[i][j] = +( i - mp[j] );
1183 g = fpenalty * ogcp1[i];
1184 if( m[j] <= previousw[j-1] + g )
1186 m[j] = previousw[j-1] + g ;
1191 lastverticalw[i] = currentw[lgth2-1];
1199 for( j=1; j<lgth2+1; j++ )
1200 currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1201 for( i=1; i<lgth1+1; i++ )
1202 lastverticalw[i] -= offset * ( lgth1 - i / 2.0);
1207 fprintf( stderr, "\n" );
1208 for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1209 fprintf( stderr, "#####\n" );
1210 for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1211 fprintf( stderr, "====>" );
1212 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1213 for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1217 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1220 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1222 // fprintf( stderr, "### impmatch = %f\n", *impmatch );
1224 resultlen = strlen( mseq1[0] );
1225 if( alloclen < resultlen || resultlen > N )
1227 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1228 ErrorExit( "LENGTH OVER!\n" );
1232 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1233 for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1235 fprintf( stderr, "\n" );
1236 for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1237 fprintf( stderr, "#####\n" );
1238 for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1245 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 )
1246 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1250 int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
1253 float wm; /* int ?????? */
1255 float *currentw, *previousw;
1256 float fpenalty = (float)penalty;
1257 float fpenalty_ex = (float)penalty_ex;
1261 float *mjpt, *prept, *curpt;
1264 static float mi, *m;
1266 static int mpi, *mp;
1267 static float *w1, *w2;
1268 static float *match;
1269 static float *initverticalw; /* kufuu sureba iranai */
1270 static float *lastverticalw; /* kufuu sureba iranai */
1271 static char **mseq1;
1272 static char **mseq2;
1274 static double *ogcp1;
1275 static double *ogcp2;
1276 static double *fgcp1;
1277 static double *fgcp2;
1278 static float **cpmx1;
1279 static float **cpmx2;
1280 static int **intwork;
1281 static float **floatwork;
1282 static int orlgth1 = 0, orlgth2 = 0;
1286 fprintf( stderr, "eff in SA+++align\n" );
1287 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1291 mseq1 = AllocateCharMtx( njob, 0 );
1292 mseq2 = AllocateCharMtx( njob, 0 );
1296 lgth1 = strlen( seq1[0] );
1297 lgth2 = strlen( seq2[0] );
1299 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
1303 if( orlgth1 > 0 && orlgth2 > 0 )
1307 FreeFloatVec( match );
1308 FreeFloatVec( initverticalw );
1309 FreeFloatVec( lastverticalw );
1314 FreeCharMtx( mseq );
1316 FreeDoubleVec( ogcp1 );
1317 FreeDoubleVec( ogcp2 );
1318 FreeDoubleVec( fgcp1 );
1319 FreeDoubleVec( fgcp2 );
1322 FreeFloatMtx( cpmx1 );
1323 FreeFloatMtx( cpmx2 );
1325 FreeFloatMtx( floatwork );
1326 FreeIntMtx( intwork );
1329 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
1330 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
1333 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1336 w1 = AllocateFloatVec( ll2+2 );
1337 w2 = AllocateFloatVec( ll2+2 );
1338 match = AllocateFloatVec( ll2+2 );
1340 initverticalw = AllocateFloatVec( ll1+2 );
1341 lastverticalw = AllocateFloatVec( ll1+2 );
1343 m = AllocateFloatVec( ll2+2 );
1344 mp = AllocateIntVec( ll2+2 );
1346 mseq = AllocateCharMtx( njob, ll1+ll2 );
1348 ogcp1 = AllocateDoubleVec( ll1+2 );
1349 ogcp2 = AllocateDoubleVec( ll2+2 );
1350 fgcp1 = AllocateDoubleVec( ll1+2 );
1351 fgcp2 = AllocateDoubleVec( ll2+2 );
1353 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1354 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1356 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
1357 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
1360 fprintf( stderr, "succeeded\n" );
1363 orlgth1 = ll1 - 100;
1364 orlgth2 = ll2 - 100;
1368 for( i=0; i<icyc; i++ ) mseq1[i] = mseq[i];
1369 for( j=0; j<jcyc; j++ ) mseq2[j] = mseq[icyc+j];
1372 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1376 if( commonAlloc1 && commonAlloc2 )
1378 FreeShortMtx( commonIP );
1381 ll1 = MAX( orlgth1, commonAlloc1 );
1382 ll2 = MAX( orlgth2, commonAlloc2 );
1385 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1388 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
1391 fprintf( stderr, "succeeded\n\n" );
1399 cpmx_calc_new( seq1, cpmx1, eff1, strlen( seq1[0] ), icyc );
1400 cpmx_calc_new( seq2, cpmx2, eff2, strlen( seq2[0] ), jcyc );
1402 OpeningGapCount( ogcp1, icyc, seq1, eff1 );
1403 OpeningGapCount( ogcp2, jcyc, seq2, eff2 );
1404 FinalGapCount( fgcp1, icyc, seq1, eff1 );
1405 FinalGapCount( fgcp2, jcyc, seq2, eff2 );
1407 for( i=0; i<lgth1; i++ )
1409 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] );
1410 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] );
1412 for( i=0; i<lgth2; i++ )
1414 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] );
1415 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] );
1418 for( i=0; i<lgth1; i++ )
1419 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1426 match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1429 match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1430 #if 0 // -> tbfast.c
1432 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1438 for( i=1; i<lgth1+1; i++ )
1440 initverticalw[i] += fpenalty * ( ogcp1[0] + fgcp1[i-1] ) ;
1442 for( j=1; j<lgth2+1; j++ )
1444 currentw[j] += fpenalty * ( ogcp2[0] + fgcp2[j-1] ) ;
1450 for( j=1; j<lgth2+1; j++ )
1451 currentw[j] -= offset * j / 2.0;
1452 for( i=1; i<lgth1+1; i++ )
1453 initverticalw[i] -= offset * i / 2.0;
1457 for( j=1; j<lgth2+1; ++j )
1459 m[j] = currentw[j-1] + fpenalty * ogcp1[1]; mp[j] = 0;
1462 lastverticalw[0] = currentw[lgth2-1];
1464 if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1467 fprintf( stderr, "currentw = \n" );
1468 for( i=0; i<lgth1+1; i++ )
1470 fprintf( stderr, "%5.2f ", currentw[i] );
1472 fprintf( stderr, "\n" );
1473 fprintf( stderr, "initverticalw = \n" );
1474 for( i=0; i<lgth2+1; i++ )
1476 fprintf( stderr, "%5.2f ", initverticalw[i] );
1478 fprintf( stderr, "\n" );
1479 fprintf( stderr, "fcgp\n" );
1480 for( i=0; i<lgth1; i++ )
1481 fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1482 for( i=0; i<lgth2; i++ )
1483 fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1486 for( i=1; i<lasti; i++ )
1489 previousw = currentw;
1492 previousw[0] = initverticalw[i-1];
1494 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1496 fprintf( stderr, "\n" );
1497 fprintf( stderr, "i=%d\n", i );
1498 fprintf( stderr, "currentw = \n" );
1499 for( j=0; j<lgth2; j++ )
1501 fprintf( stderr, "%5.2f ", currentw[j] );
1503 fprintf( stderr, "\n" );
1507 // fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1509 imp_match_out_vead( currentw, i, lgth2 );
1511 imp_match_out_vead_gapmap( currentw, gapmap1[i], lgth2, gapmap2 );
1515 fprintf( stderr, "\n" );
1516 fprintf( stderr, "i=%d\n", i );
1517 fprintf( stderr, "currentw = \n" );
1518 for( j=0; j<lgth2; j++ )
1520 fprintf( stderr, "%5.2f ", currentw[j] );
1522 fprintf( stderr, "\n" );
1524 currentw[0] = initverticalw[i];
1528 mi = previousw[0] + fpenalty * ogcp2[1]; mpi = 0;
1533 curpt = currentw + 1;
1535 for( j=1; j<lgth2+1; j++ )
1541 fprintf( stderr, "%5.0f->", wm );
1543 g = mi + fpenalty * fgcp2[j-1];
1545 fprintf( stderr, "%5.0f?", g );
1550 *ijppt = -( j - mpi );
1552 g = *prept + fpenalty * ogcp2[j];
1562 g = *mjpt + fpenalty * fgcp1[i-1];
1564 fprintf( stderr, "%5.0f?", g );
1569 *ijppt = +( i - *mpjpt );
1571 g = *prept + fpenalty * ogcp1[i];
1578 m[j] += fpenalty_ex;
1582 fprintf( stderr, "%5.0f ", wm );
1591 lastverticalw[i] = currentw[lgth2-1];
1593 floatncpy( previousw, currentw, lgth2+1 );
1594 previousw[0] = initverticalw[i-1];
1597 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1598 currentw[0] = initverticalw[i];
1600 mi = previousw[0] + fpenalty * ogcp2[1]; mpi = 0;
1602 for( j=1; j<lastj; j++ )
1604 wm = previousw[j-1];
1607 g = fpenalty * fgcp2[j-1];
1611 ijp[i][j] = -( j - mpi );
1613 g = fpenalty * ogcp2[j];
1614 if( mi <= previousw[j-1] + g )
1616 mi = previousw[j-1] + g;
1620 g = fpenalty * fgcp1[i-1];
1624 ijp[i][j] = +( i - mp[j] );
1626 g = fpenalty * ogcp1[i];
1627 if( m[j] <= previousw[j-1] + g )
1629 m[j] = previousw[j-1] + g ;
1634 lastverticalw[i] = currentw[lgth2-1];
1642 for( j=1; j<lgth2+1; j++ )
1643 currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1644 for( i=1; i<lgth1+1; i++ )
1645 lastverticalw[i] -= offset * ( lgth1 - i / 2.0);
1650 fprintf( stderr, "\n" );
1651 for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1652 fprintf( stderr, "#####\n" );
1653 for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1654 fprintf( stderr, "====>" );
1655 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1656 for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1660 Atracking_localhom_gapmap( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, gapmap1, gapmap2 );
1663 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1665 // fprintf( stderr, "### impmatch = %f\n", *impmatch );
1667 resultlen = strlen( mseq1[0] );
1668 if( alloclen < resultlen || resultlen > N )
1670 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1671 ErrorExit( "LENGTH OVER!\n" );
1675 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1676 for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1678 fprintf( stderr, "\n" );
1679 for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1680 fprintf( stderr, "#####\n" );
1681 for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1688 float translate_and_Calign( char **mseq1, char **mseq2, double *effarr1, double *effarr2, int clus1, int clus2, int alloclen )
1698 if ( clus1 == 1 && clus2 != 1 )
1700 seq = mseq1[0]; aseq = mseq2; effarr = effarr2; nseq = clus2+1;
1702 printf( "effarr in transl... = \n" );
1703 for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
1706 else if( clus1 != 1 && clus2 == 1 )
1708 seq = mseq2[0]; aseq = mseq1; effarr = effarr1; nseq = clus1+1;
1710 else ErrorExit( "ERROR in translate_and_Calign" );
1712 result = Calignm1( &wm, aseq, seq, effarr, nseq-2, 0 );
1714 resultlen = strlen( result[0] );
1715 if( alloclen < resultlen || resultlen > N )
1717 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1718 ErrorExit( "LENGTH OVER!\n" );
1720 for( i=0; i<nseq-1; i++ ) strcpy( aseq[i], result[i] );
1721 strcpy( seq, result[nseq-1] );