8 #define USE_PENALTY_EX 0
16 static void OpeningGapCount( float *ogcp, int clus, char **seq, double *eff, int len )
21 for( i=0; i<len; i++ ) ogcp[i] = 0.0;
22 for( j=0; j<clus; j++ )
26 // gc = ( gappat[j] == '-' );
27 for( i=0; i<len; i++ )
30 gc = ( seq[j][i] == '-' );
32 if( !gb * gc ) ogcp[i] += feff;
36 // for( i=0; i<len; i++ )
38 // fprintf( stderr, "ogcp[%d] = %30.20f\n", i, ogcp[i] );
42 static void FinalGapCount( float *fgcp, int clus, char **seq, double *eff, int len )
47 for( i=0; i<len; i++ ) fgcp[i] = 0.0;
48 for( j=0; j<clus; j++ )
51 // gc = ( seq[j][0] == '-' );
53 for( i=1; i<len; i++ )
56 gc = ( seq[j][i] == '-' );
68 // fprintf( stderr, "changed gc=%d, gb=%d\n", gc, gb );
69 if( gb * !gc ) fgcp[i-1] += feff;
72 // fprintf( stderr, "eff[%d] = %30.20f\n", j, eff[j] );
74 // for( i=0; i<len; i++ )
76 // fprintf( stderr, "fgcp[%d] = %30.20f\n", i, fgcp[i] );
80 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
84 float **cpmxpd = floatwork;
85 int **cpmxpdn = intwork;
90 for( j=0; j<lgth2; j++ )
97 cpmxpd[j][count] = cpmx2[j][l];
98 cpmxpdn[j][count] = l;
102 cpmxpdn[j][count] = -1;
106 for( l=0; l<26; l++ )
109 for( k=0; k<26; k++ )
111 scarr[l] += n_dis[k][l] * cpmx1[i1][k];
114 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤
\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
116 float *fpt, **fptpt, *fpt2;
124 ipt=*iptpt,fpt=*fptpt;
126 *fpt2 += scarr[*ipt++] * *fpt++;
127 fpt2++,iptpt++,fptpt++;
131 for( j=0; j<lgth2; j++ )
134 for( k=0; cpmxpdn[j][k]>-1; k++ )
135 match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
140 static float Atracking(
141 char **seq1, char **seq2,
142 char **mseq1, char **mseq2,
143 short **ijp, int icyc, int jcyc,
144 int ist, int ien, int jst, int jen )
146 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, klim;
152 for( i=0; i<lgth1; i++ )
154 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
159 // fprintf( stderr, "in Atracking, lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
161 for( i=0; i<lgth1+1; i++ )
165 for( j=0; j<lgth2+1; j++ )
167 ijp[0][j] = -( j + 1 );
171 for( i=0; i<icyc; i++ )
173 mseq1[i] += lgth1+lgth2;
176 for( j=0; j<jcyc; j++ )
178 mseq2[j] += lgth1+lgth2;
182 // if( lgth2 == 1 ) fprintf( stderr, "in Atracking, mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
184 iin = lgth1; jin = lgth2;
186 for( k=0; k<=klim; k++ )
188 if( ijp[iin][jin] < 0 )
190 ifi = iin-1; jfi = jin+ijp[iin][jin];
192 else if( ijp[iin][jin] > 0 )
194 ifi = iin-ijp[iin][jin]; jfi = jin-1;
198 ifi = iin-1; jfi = jin-1;
203 for( i=0; i<icyc; i++ )
204 *--mseq1[i] = seq1[i][ist+ifi+l];
205 for( j=0; j<jcyc; j++ )
212 for( i=0; i<icyc; i++ )
214 for( j=0; j<jcyc; j++ )
215 *--mseq2[j] = seq2[j][jst+jfi+l];
218 if( iin <= 0 || jin <= 0 ) break;
219 for( i=0; i<icyc; i++ )
220 *--mseq1[i] = seq1[i][ist+ifi];
221 for( j=0; j<jcyc; j++ )
222 *--mseq2[j] = seq2[j][jst+jfi];
224 iin = ifi; jin = jfi;
226 // if( lgth2 == 1 ) fprintf( stderr, "in Atracking, mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
228 // fprintf( stderr, "in Atracking (owari), mseq1 = %s\n", mseq1[0] );
229 // fprintf( stderr, "in Atracking (owari), mseq2 = %s\n", mseq2[0] );
233 static float MSalignmm_tanni( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, float **cpmx1, float **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, float **gapinfo )
234 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
241 float wm = 0.0; /* int ?????? */
243 float *currentw, *previousw;
244 float fpenalty = (float)penalty;
245 float fpenalty_ex = (float)penalty_ex;
249 float *mjpt, *prept, *curpt;
256 float *initverticalw; /* kufuu sureba iranai */
257 float *lastverticalw; /* kufuu sureba iranai */
269 char **aseq1bk, **aseq2bk;
271 ogcp1 = gapinfo[0] + ist;
272 fgcp1 = gapinfo[1] + ist;
273 ogcp2 = gapinfo[2] + jst;
274 fgcp2 = gapinfo[3] + jst;
277 char ttt1[10000], ttt2[10000];
285 strncpy( ttt1, seq1[0]+ist, lgth1 ); ttt1[lgth1] = 0;
286 strncpy( ttt2, seq2[0]+jst, lgth2 ); ttt2[lgth2] = 0;
288 fprintf( stderr, "in _tanni ist,ien = %d,%d, lgth1=%d\n", ist, ien, lgth1 );
289 fprintf( stderr, "in _tanni jst,jen = %d,%d, lgth2=%d\n", jst, jen, lgth2 );
290 fprintf( stderr, "ttt1 = %s\n", ttt1 );
291 fprintf( stderr, "ttt2 = %s\n", ttt2 );
295 ll1 = ( (int)(1.3*lgth1) ) + 100;
296 ll2 = ( (int)(1.3*lgth2) ) + 100;
298 aseq1 = AllocateCharMtx( icyc, 0 );
299 aseq2 = AllocateCharMtx( jcyc, 0 );
300 aseq1bk = AllocateCharMtx( icyc, lgth1+lgth2+100 );
301 aseq2bk = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
302 for( i=0; i<icyc; i++ ) aseq1[i] = aseq1bk[i];
303 for( i=0; i<jcyc; i++ ) aseq2[i] = aseq2bk[i];
305 w1 = AllocateFloatVec( ll2+2 );
306 w2 = AllocateFloatVec( ll2+2 );
308 initverticalw = AllocateFloatVec( ll1+2 );
309 lastverticalw = AllocateFloatVec( ll1+2 );
311 m = AllocateFloatVec( ll2+2 );
312 mp = AllocateIntVec( ll2+2 );
314 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 27 );
315 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 );
318 shortmtx = AllocateShortMtx( ll1+1, ll2+1 );
325 match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
327 match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
329 for( i=1; i<lgth1+1; i++ )
331 initverticalw[i] += fpenalty * ( ogcp1[0] + fgcp1[i-1] );
333 for( j=1; j<lgth2+1; j++ )
335 currentw[j] += fpenalty * ( ogcp2[0] + fgcp2[j-1] );
338 for( j=1; j<lgth2+1; ++j )
340 m[j] = currentw[j-1] + fpenalty * ogcp1[1]; mp[j] = 0;;
343 lastverticalw[0] = currentw[lgth2-1];
348 for( i=1; i<lasti; i++ )
352 previousw = currentw;
355 previousw[0] = initverticalw[i-1];
357 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
358 currentw[0] = initverticalw[i];
360 mi = previousw[0] + fpenalty * ogcp2[1];
366 curpt = currentw + 1;
369 for( j=1; j<lastj; j++ )
375 fprintf( stderr, "%5.0f->", wm );
377 g = mi + fpenalty * fgcp2[j-1];
379 fprintf( stderr, "%5.0f?", g );
384 *ijppt = -( j - mpi );
386 g = *prept + fpenalty * ogcp2[j];
396 g = *mjpt + fpenalty * fgcp1[i-1];
398 fprintf( stderr, "%5.0f?", g );
403 *ijppt = +( i - *mpjpt );
405 g = *prept + fpenalty * ogcp1[i];
416 fprintf( stderr, "%5.0f ", wm );
427 lastverticalw[i] = currentw[lgth2-1];
430 // fprintf( stderr, "wm = %f\n", wm );
432 Atracking( seq1, seq2, aseq1, aseq2, ijp, icyc, jcyc, ist, ien, jst, jen );
434 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
435 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
437 // fprintf( stderr, "in _tanni, aseq1 = %s\n", aseq1[0] );
438 // fprintf( stderr, "in _tanni, mseq1 = %s\n", mseq1[0] );
442 FreeFloatVec( initverticalw );
443 FreeFloatVec( lastverticalw );
449 FreeFloatMtx( floatwork );
450 FreeIntMtx( intwork );
452 FreeShortMtx( shortmtx );
454 FreeCharMtx( aseq1bk );
455 FreeCharMtx( aseq2bk );
464 static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, float **cpmx1, float **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, int depth, float **gapinfo )
465 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
471 char **aseq1, **aseq2;
472 int ll1, ll2, l, len;
473 int lasti, lastj, imid, jmid;
475 float wm = 0.0; /* int ?????? */
477 float *currentw, *previousw;
478 float fpenalty = (float)penalty;
479 float fpenalty_ex = (float)penalty_ex;
486 float *mjpt, *prept, *curpt;
491 float *initverticalw; /* kufuu sureba iranai */
492 float *lastverticalw; /* kufuu sureba iranai */
510 int *jumpdummi; //muda
511 int *jumpdummj; //muda
522 static char ttt1[50000];
523 static char ttt2[50000];
526 ogcp1 = gapinfo[0] + ist;
527 fgcp1 = gapinfo[1] + ist;
528 ogcp2 = gapinfo[2] + jst;
529 fgcp2 = gapinfo[3] + jst;
539 fprintf( stderr, "==== MSalign (depth=%d, reccycle=%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, reccycle, ist, ien, jst, jen );
540 strncpy( ttt1, seq1[0]+ist, lgth1 );
541 strncpy( ttt2, seq2[0]+jst, lgth2 );
544 fprintf( stderr, "seq1 = %s\n", ttt1 );
545 fprintf( stderr, "seq2 = %s\n", ttt2 );
547 if( lgth2 <= 0 ) // lgth1 <= 0 ha?
549 // fprintf( stderr, "\n\n==== jimei\n\n" );
551 for( i=0; i<icyc; i++ )
553 strncpy( mseq1[i], seq1[i]+ist, lgth1 );
556 for( i=0; i<jcyc; i++ )
559 for( j=0; j<lgth1; j++ )
560 strcat( mseq2[i], "-" );
563 // fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
564 // fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
570 aseq1 = AllocateCharMtx( icyc, 0 );
571 aseq2 = AllocateCharMtx( jcyc, 0 );
572 for( i=0; i<icyc; i++ ) aseq1[i] = mseq1[i];
573 for( i=0; i<jcyc; i++ ) aseq2[i] = mseq2[i];
575 aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
576 aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
579 if( lgth1 < DPTANNI || lgth2 < DPTANNI )
581 // fprintf( stderr, "==== Going to _tanni\n" );
583 value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2, gapinfo );
588 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
589 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
591 FreeCharMtx( aseq1 );
592 FreeCharMtx( aseq2 );
595 // fprintf( stderr, "value = %f\n", value );
599 // fprintf( stderr, "Trying to divide the mtx\n" );
601 ll1 = ( (int)(1.3*lgth1) ) + 100;
602 ll2 = ( (int)(1.3*lgth2) ) + 100;
604 // fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
606 w1 = AllocateFloatVec( ll2+2 );
607 w2 = AllocateFloatVec( ll2+2 );
608 // match = AllocateFloatVec( ll2+2 );
609 midw = AllocateFloatVec( ll2+2 );
610 midn = AllocateFloatVec( ll2+2 );
611 midm = AllocateFloatVec( ll2+2 );
612 jumpbacki = AllocateIntVec( ll2+2 );
613 jumpbackj = AllocateIntVec( ll2+2 );
614 jumpforwi = AllocateIntVec( ll2+2 );
615 jumpforwj = AllocateIntVec( ll2+2 );
616 jumpdummi = AllocateIntVec( ll2+2 );
617 jumpdummj = AllocateIntVec( ll2+2 );
619 initverticalw = AllocateFloatVec( ll1+2 );
620 lastverticalw = AllocateFloatVec( ll1+2 );
622 m = AllocateFloatVec( ll2+2 );
623 mp = AllocateIntVec( ll2+2 );
624 gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
626 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
627 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
630 fprintf( stderr, "succeeded\n" );
634 WMMTX = AllocateFloatMtx( ll1, ll2 );
635 WMMTX2 = AllocateFloatMtx( ll1, ll2 );
638 shortmtx = AllocateShortMtx( ll1, ll2 );
641 fprintf( stderr, "succeeded\n\n" );
650 match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
652 match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
654 for( i=1; i<lgth1+1; i++ )
656 initverticalw[i] += fpenalty * ( ogcp1[0] + fgcp1[i-1] );
657 // initverticalw[i] += fpenalty;
659 for( j=1; j<lgth2+1; j++ )
661 currentw[j] += fpenalty * ( ogcp2[0] + fgcp2[j-1] );
662 // currentw[j] += fpenalty;
666 WMMTX[0][0] = initverticalw[0];
667 for( i=1; i<lgth1+1; i++ )
669 WMMTX[i][0] = initverticalw[i];
671 for( j=1; j<lgth2+1; j++ )
673 WMMTX[0][j] = currentw[j];
678 for( j=1; j<lgth2+1; ++j )
680 m[j] = currentw[j-1] + fpenalty * ogcp1[1];
681 // m[j] = currentw[j-1];
685 lastverticalw[0] = currentw[lgth2-1];
689 jumpi = 0; // atode kawaru.
692 for( i=1; i<lasti; i++ )
694 for( i=1; i<=imid; i++ )
698 previousw = currentw;
701 previousw[0] = initverticalw[i-1];
703 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
704 currentw[0] = initverticalw[i];
706 m[0] = fpenalty * ogcp1[i];
710 if( i == imid ) midm[0] = m[0];
712 mi = previousw[0] + fpenalty * ogcp2[1];
713 // mi = previousw[0];
717 // ijppt = ijp[i] + 1;
720 curpt = currentw + 1;
725 for( j=1; j<lastj; j++ )
731 fprintf( stderr, "%5.0f->", wm );
733 g = mi + fpenalty * fgcp2[j-1];
734 // g = mi + fpenalty;
736 fprintf( stderr, "%5.0f?", g );
741 // *ijppt = -( j - mpi );
743 g = *prept + fpenalty * ogcp2[j];
754 g = *mjpt + fpenalty * fgcp1[i-1];
755 // g = *mjpt + fpenalty;
757 fprintf( stderr, "%5.0f?", g );
762 // *ijppt = +( i - *mpjpt );
766 g = *prept + fpenalty * ogcp1[i];
778 fprintf( stderr, "%5.0f ", wm );
783 WMMTX[i][j] = *curpt;
784 WMMTX2[i][j] = *mjpt;
787 if( i == imid ) //muda
789 jumpbackj[j] = *mpjpt; // muda atode matomeru
790 jumpbacki[j] = mpi; // muda atode matomeru
791 // fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
792 // fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
798 // fprintf( stderr, "m[%d] = %f\n", j, m[j] );
805 lastverticalw[i] = currentw[lgth2-1];
808 WMMTX2[i][lgth2] = m[lgth2-1];
814 for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
815 for( j=0; j<lgth2; j++ ) midm[j] = m[j];
819 // for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
820 // for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
823 for( i=0; i<lgth1; i++ )
825 for( j=0; j<lgth2; j++ )
827 fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
829 fprintf( stderr, "\n" );
831 fprintf( stderr, "\n" );
832 fprintf( stderr, "WMMTX2 = \n" );
833 for( i=0; i<lgth1; i++ )
835 for( j=0; j<lgth2; j++ )
837 fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
839 fprintf( stderr, "\n" );
841 fprintf( stderr, "\n" );
846 match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
847 match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
849 for( i=0; i<lgth1-1; i++ )
851 initverticalw[i] += fpenalty * ( fgcp1[lgth1-1] + ogcp1[i+1] );
852 // initverticalw[i] += fpenalty;
854 for( j=0; j<lgth2-1; j++ )
856 currentw[j] += fpenalty * ( fgcp2[lgth2-1] + ogcp2[j+1] );
857 // currentw[j] += fpenalty;
861 for( i=0; i<lgth1-1; i++ )
863 WMMTX[i][lgth2-1] += fpenalty * ( fgcp1[lgth1-1] + ogcp1[i+1] );
864 fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
866 for( j=0; j<lgth2-1; j++ )
868 WMMTX[lgth1-1][j] += fpenalty * ( fgcp2[lgth2-1] + ogcp2[j+1] );
869 fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
878 for( j=lgth2-1; j>0; --j )
880 m[j-1] = currentw[j] + fpenalty * fgcp2[lgth2-2];
881 // m[j-1] = currentw[j];
885 // for( j=0; j<lgth2; j++ ) m[j] = 0.0;
886 // m[lgth2-1] ha irunoka?
889 // for( i=lgth1-2; i>=imid; i-- )
892 for( i=lgth1-2; i>-1; i-- )
895 previousw = currentw;
897 previousw[lgth2-1] = initverticalw[i+1];
898 // match_calc( currentw, seq1, seq2, i, lgth2 );
899 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
901 currentw[lgth2-1] = initverticalw[i];
903 // m[lgth2] = fpenalty * fgcp1[i];
904 // WMMTX2[i][lgth2] += m[lgth2];
905 // fprintf( stderr, "m[] = %f\n", m[lgth2] );
907 mi = previousw[lgth2-1] + fpenalty * fgcp2[lgth2-2];
908 // mi = previousw[lgth2-1];
911 mjpt = m + lgth2 - 2;
912 prept = previousw + lgth2 - 1;
913 curpt = currentw + lgth2 - 2;
914 mpjpt = mp + lgth2 - 2;
917 for( j=lgth2-2; j>-1; j-- )
923 g = mi + fpenalty * ogcp2[j+1];
924 // g = mi + fpenalty;
932 g = *prept + fpenalty * fgcp2[j];
936 // fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
945 // fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
946 g = *mjpt + fpenalty * ogcp1[i+1];
947 // g = *mjpt + fpenalty;
955 // if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
956 g = *prept + fpenalty * fgcp1[i];
968 if( i == jumpi || i == imid - 1 )
970 jumpforwi[j] = ijpi; //muda
971 jumpforwj[j] = ijpj; //muda
972 // fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
973 // fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
975 if( i == imid ) // muda
978 // midm[j+1] += *mjpt + fpenalty; //??????
979 midm[j+1] += *mjpt; //??????
983 // midn[j] += mi + fpenalty; //????
984 midn[j] += mi; //????
988 // WMMTX2[i][j+1] += *mjpt + fpenalty;
989 WMMTX2[i][j+1] += *mjpt;
998 // fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
999 g = *prept + fpenalty * fgcp1[i];
1006 WMMTX2[i][j+1] += firstm;
1008 if( i == imid ) midm[j+1] += firstm;
1014 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1015 for( j=2; j<lgth2-1; j++ )
1023 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1025 for( j=0; j<lgth2+1; j++ )
1033 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1036 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1039 // fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
1043 if( jmid > 0 && midn[jmid-1] > wm ) //060413
1046 jumpj = jumpbacki[jmid];
1048 // fprintf( stderr, "rejump (n)\n" );
1050 if( midm[jmid] > wm )
1052 jumpi = jumpbackj[jmid];
1055 // fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
1059 // fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
1060 // fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
1062 fprintf( stderr, "imid = %d\n", imid );
1063 fprintf( stderr, "midn = \n" );
1064 for( j=0; j<lgth2; j++ )
1066 fprintf( stderr, "% 7.1f ", midn[j] );
1068 fprintf( stderr, "\n" );
1069 fprintf( stderr, "midw = \n" );
1070 for( j=0; j<lgth2; j++ )
1072 fprintf( stderr, "% 7.1f ", midw[j] );
1074 fprintf( stderr, "\n" );
1075 fprintf( stderr, "midm = \n" );
1076 for( j=0; j<lgth2; j++ )
1078 fprintf( stderr, "% 7.1f ", midm[j] );
1080 fprintf( stderr, "\n" );
1082 // fprintf( stderr, "maxwm = %f\n", maxwm );
1084 if( i == jumpi ) //saki?
1086 // fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
1087 // fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
1090 // fprintf( stderr, "CHUI2!\n" );
1091 jumpj = 0; jmid = 1;
1092 jumpi = firstmp - 1;
1097 else if( jmid == lgth2 )
1099 fprintf( stderr, "CHUI1!\n" );
1101 imid=jumpforwi[0]; jmid=lgth2-1;
1104 else if( jmid >= lgth2 )
1106 fprintf( stderr, "CHUI1!\n" );
1107 jumpi=imid-1; jmid=lgth2;
1113 imid = jumpforwi[jumpj];
1114 jmid = jumpforwj[jumpj];
1117 fprintf( stderr, "jumpi -> %d\n", jumpi );
1118 fprintf( stderr, "jumpj -> %d\n", jumpj );
1119 fprintf( stderr, "imid -> %d\n", imid );
1120 fprintf( stderr, "jmid -> %d\n", jmid );
1132 imid=lgth1-1; jmid=lgth2-1;
1136 // fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1137 // fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1139 // for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1140 // for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1141 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1142 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1146 fprintf( stderr, "WMMTX = \n" );
1147 for( i=0; i<lgth1; i++ )
1149 fprintf( stderr, "%d ", i );
1150 for( j=0; j<lgth2; j++ )
1152 fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
1154 fprintf( stderr, "\n" );
1156 fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
1157 for( i=0; i<lgth1; i++ )
1159 fprintf( stderr, "%d ", i );
1160 for( j=0; j<lgth2+1; j++ )
1162 fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
1164 fprintf( stderr, "\n" );
1167 fprintf( stderr, "jumpbacki = \n" );
1168 for( j=0; j<lgth2; j++ )
1170 fprintf( stderr, "% 7d ", jumpbacki[j] );
1172 fprintf( stderr, "\n" );
1173 fprintf( stderr, "jumpbackj = \n" );
1174 for( j=0; j<lgth2; j++ )
1176 fprintf( stderr, "% 7d ", jumpbackj[j] );
1178 fprintf( stderr, "\n" );
1179 fprintf( stderr, "jumpforwi = \n" );
1180 for( j=0; j<lgth2; j++ )
1182 fprintf( stderr, "% 7d ", jumpforwi[j] );
1184 fprintf( stderr, "\n" );
1185 fprintf( stderr, "jumpforwj = \n" );
1186 for( j=0; j<lgth2; j++ )
1188 fprintf( stderr, "% 7d ", jumpforwj[j] );
1190 fprintf( stderr, "\n" );
1196 // Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1199 resultlen = strlen( mseq1[0] );
1200 if( alloclen < resultlen || resultlen > N )
1202 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1203 ErrorExit( "LENGTH OVER!\n" );
1210 fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
1211 fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
1213 fprintf( stderr, "imid = %d\n", imid );
1214 fprintf( stderr, "jmid = %d\n", jmid );
1220 FreeFloatVec( initverticalw );
1221 FreeFloatVec( lastverticalw );
1222 FreeFloatVec( midw );
1223 FreeFloatVec( midm );
1224 FreeFloatVec( midn );
1226 FreeIntVec( jumpbacki );
1227 FreeIntVec( jumpbackj );
1228 FreeIntVec( jumpforwi );
1229 FreeIntVec( jumpforwj );
1230 FreeIntVec( jumpdummi );
1231 FreeIntVec( jumpdummj );
1236 FreeFloatMtx( floatwork );
1237 FreeIntMtx( intwork );
1240 FreeFloatMtx( WMMTX );
1241 FreeFloatMtx( WMMTX2 );
1244 // fprintf( stderr, "==== calling myself (first)\n" );
1246 // fprintf( stderr, "jumpi = %d\n", jumpi );
1247 value = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ist+jumpi, jst, jst+jumpj, alloclen, aseq1, aseq2, depth, gapinfo );
1249 alnlen = strlen( aseq1[0] );
1250 for( i=0; i<icyc; i++ ) aseq1[i] += alnlen;
1251 for( i=0; i<jcyc; i++ ) aseq2[i] += alnlen;
1253 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
1254 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
1257 // fprintf( stderr, "====(f) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
1258 // fprintf( stderr, "====(f) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
1260 len = strlen( mseq1[0] );
1261 // fprintf( stderr, "len = %d\n", len );
1262 l = jmid - jumpj - 1;
1263 // fprintf( stderr, "l=%d\n", l );
1266 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
1267 for( i=0; i<icyc; i++ )
1269 strcat( mseq1[i], gaps );
1270 mseq1[i][len+l] = 0;
1272 for( j=0; j<jcyc; j++ )
1274 strncat( mseq2[j], seq2[j]+jst+jumpj+1, l );
1275 mseq2[j][len+l] = 0;
1277 // fprintf( stderr, "penalizing (2) .. %f(%d), %f(%d)\n", ogcp2[jumpj+1], jumpj+1, fgcp2[jmid-1], jmid-1 );
1278 value += fpenalty * ( ogcp2[jumpj+1] + fgcp2[jmid-1] );
1279 // value += fpenalty;
1281 len = strlen( mseq1[0] );
1282 l = imid - jumpi - 1;
1283 // fprintf( stderr, "l=%d\n", l );
1286 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
1287 for( i=0; i<icyc; i++ )
1289 strncat( mseq1[i], seq1[i]+ist+jumpi+1, l );
1290 mseq1[i][len+l] = 0;
1292 for( j=0; j<jcyc; j++ )
1294 strcat( mseq2[j], gaps );
1295 mseq2[j][len+l] = 0;
1298 // for( i=0; i<lgth1; i++ ) fprintf( stderr, "ogcp1[%d] = %f\n", i, ogcp1[i] );
1299 // for( i=0; i<lgth1; i++ ) fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
1302 // fprintf( stderr, "penalizing (1) .. ogcp1[%d] = %f, fgcp1[%d] = %f\n", jumpi+1, ogcp1[jumpi+1], imid-1, fgcp1[imid-1] );
1303 value += fpenalty * ( ogcp1[jumpi+1] + fgcp1[imid-1] );
1304 // value += fpenalty;
1307 for( i=0; i<icyc; i++ ) fprintf( stderr, "after gapfill mseq1[%d]=%s\n", i, mseq1[i] );
1308 for( i=0; i<jcyc; i++ ) fprintf( stderr, "after gapfill mseq2[%d]=%s\n", i, mseq2[i] );
1311 // fprintf( stderr, "==== calling myself (second)\n" );
1314 value += MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist+imid, ien, jst+jmid, jen, alloclen, aseq1, aseq2, depth, gapinfo );
1319 if( value - maxwm > 1 || maxwm - value > 1 )
1321 fprintf( stderr, "WARNING value = %f, but maxwm = %f penalty = %f\n", value, maxwm, fpenalty );
1322 for( i=0; i<icyc; i++ )
1324 fprintf( stderr, ">1-%d\n%s\n", i, mseq1[i] );
1325 fprintf( stderr, "%s\n", aseq1[i] );
1327 for( i=0; i<jcyc; i++ )
1329 fprintf( stderr, ">2-%d\n%s\n", i, mseq2[i] );
1330 fprintf( stderr, "%s\n", aseq2[i] );
1337 fprintf( stderr, "value = %.0f, maxwm = %.0f -> ok\n", value, maxwm );
1343 for( i=0; i<icyc; i++ ) strcat( mseq1[i], aseq1[i] );
1344 for( i=0; i<jcyc; i++ ) strcat( mseq2[i], aseq2[i] );
1347 // fprintf( stderr, "====(s) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
1348 // fprintf( stderr, "====(s) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
1353 FreeCharMtx( aseq1 );
1354 FreeCharMtx( aseq2 );
1362 float MSalignmm( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen )
1363 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1369 float wm = 0.0; /* int ?????? */
1382 fprintf( stderr, "eff in SA+++align\n" );
1383 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1386 lgth1 = strlen( seq1[0] );
1387 lgth2 = strlen( seq2[0] );
1389 ll1 = ( (int)(1.3*lgth1) ) + 100;
1390 ll2 = ( (int)(1.3*lgth2) ) + 100;
1392 mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
1393 mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
1395 gapinfo = AllocateFloatMtx( 4, 0 );
1396 ogcp1 = AllocateFloatVec( ll1+2 );
1397 ogcp2 = AllocateFloatVec( ll2+2 );
1398 fgcp1 = AllocateFloatVec( ll1+2 );
1399 fgcp2 = AllocateFloatVec( ll2+2 );
1402 cpmx1 = AllocateFloatMtx( ll1+2, 27 );
1403 cpmx2 = AllocateFloatMtx( ll2+2, 27 );
1405 for( i=0; i<icyc; i++ )
1409 for( j=0; j<jcyc; j++ )
1414 MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1415 MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1420 OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1421 OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1422 FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1423 FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1426 for( i=0; i<lgth1; i++ )
1428 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] );
1429 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] );
1430 // fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
1432 for( i=0; i<lgth2; i++ )
1434 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] );
1435 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] );
1436 // fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
1439 for( i=0; i<lgth1; i++ )
1444 for( i=0; i<lgth2; i++ )
1458 fprintf( stdout, "in MSalignmm.c\n" );
1459 for( i=0; i<icyc; i++ )
1461 fprintf( stdout, ">%d of GROUP1\n", i );
1462 fprintf( stdout, "%s\n", seq1[i] );
1464 for( i=0; i<jcyc; i++ )
1466 fprintf( stdout, ">%d of GROUP2\n", i );
1467 fprintf( stdout, "%s\n", seq2[i] );
1472 wm = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo );
1474 // fprintf( stderr, "wm = %f\n", wm );
1476 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1477 for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
1480 FreeFloatVec( ogcp1 );
1481 FreeFloatVec( ogcp2 );
1482 FreeFloatVec( fgcp1 );
1483 FreeFloatVec( fgcp2 );
1484 FreeFloatMtx( cpmx1 );
1485 FreeFloatMtx( cpmx2 );
1486 free( (void *)gapinfo );
1488 FreeCharMtx( mseq1 );
1489 FreeCharMtx( mseq2 );