7 #define USE_PENALTY_EX 0
14 static int reccycle = 0;
16 static float localthr;
18 static void match_ribosum( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
22 float **cpmxpd = floatwork;
23 int **cpmxpdn = intwork;
32 for( j=0; j<lgth2; j++ )
39 cpmxpd[j][count] = cpmx2[j][l];
40 cpmxpdn[j][count] = l;
44 cpmxpdn[j][count] = -1;
53 scarr[l] += ribosumdis[k][l] * cpmx1[i1][k];
56 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤
\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
58 float *fpt, **fptpt, *fpt2;
66 ipt=*iptpt,fpt=*fptpt;
68 *fpt2 += scarr[*ipt++] * *fpt++;
69 fpt2++,iptpt++,fptpt++;
72 for( j=0; j<lgth2; j++ )
75 for( k=0; cpmxpdn[j][k]>-1; k++ )
76 match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
85 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
86 *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
94 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
98 float **cpmxpd = floatwork;
99 int **cpmxpdn = intwork;
108 for( j=0; j<lgth2; j++ )
111 for( l=0; l<26; l++ )
115 cpmxpd[j][count] = cpmx2[j][l];
116 cpmxpdn[j][count] = l;
120 cpmxpdn[j][count] = -1;
124 for( l=0; l<26; l++ )
127 for( k=0; k<26; k++ )
129 scarr[l] += (n_dis[k][l]-RNAthr) * cpmx1[i1][k];
132 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤
\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
134 float *fpt, **fptpt, *fpt2;
142 ipt=*iptpt,fpt=*fptpt;
144 *fpt2 += scarr[*ipt++] * *fpt++;
145 fpt2++,iptpt++,fptpt++;
148 for( j=0; j<lgth2; j++ )
151 for( k=0; cpmxpdn[j][k]>-1; k++ )
152 match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
161 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
162 *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
171 static void match_add( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
175 float **cpmxpd = floatwork;
176 int **cpmxpdn = intwork;
186 for( j=0; j<lgth2; j++ )
189 for( l=0; l<26; l++ )
193 cpmxpd[j][count] = cpmx2[j][l];
194 cpmxpdn[j][count] = l;
198 cpmxpdn[j][count] = -1;
202 for( l=0; l<26; l++ )
205 for( k=0; k<26; k++ )
207 scarr[l] += n_dis[k][l] * cpmx1[i1][k];
210 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤
\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
212 float *fpt, **fptpt, *fpt2;
220 ipt=*iptpt,fpt=*fptpt;
222 *fpt2 += scarr[*ipt++] * *fpt++;
223 fpt2++,iptpt++,fptpt++;
226 for( j=0; j<lgth2; j++ )
229 for( k=0; cpmxpdn[j][k]>-1; k++ )
230 match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
238 // *matchpt = 0.0; // add dakara
239 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
240 *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
250 static float Atracking(
251 char **seq1, char **seq2,
252 char **mseq1, char **mseq2,
253 int **ijp, int icyc, int jcyc,
254 int ist, int ien, int jst, int jen )
256 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, klim;
257 char *gaptable1, *gt1bk;
258 char *gaptable2, *gt2bk;
262 gt1bk = AllocateCharVec( lgth1+lgth2+1 );
263 gt2bk = AllocateCharVec( lgth1+lgth2+1 );
266 for( i=0; i<lgth1; i++ )
268 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
273 // fprintf( stderr, "in Atracking, lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
275 for( i=0; i<lgth1+1; i++ )
279 for( j=0; j<lgth2+1; j++ )
281 ijp[0][j] = -( j + 1 );
285 gaptable1 = gt1bk + lgth1+lgth2;
287 gaptable2 = gt2bk + lgth1+lgth2;
290 // if( lgth2 == 1 ) fprintf( stderr, "in Atracking, mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
292 iin = lgth1; jin = lgth2;
294 for( k=0; k<=klim; k++ )
296 if( ijp[iin][jin] < 0 )
298 ifi = iin-1; jfi = jin+ijp[iin][jin];
300 else if( ijp[iin][jin] > 0 )
302 ifi = iin-ijp[iin][jin]; jfi = jin-1;
306 ifi = iin-1; jfi = jin-1;
322 if( iin <= 0 || jin <= 0 ) break;
326 iin = ifi; jin = jfi;
330 for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i]+ist, gaptable1 );
331 for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j]+jst, gaptable2 );
336 // fprintf( stderr, "in Atracking (owari), mseq1 = %s\n", mseq1[0] );
337 // fprintf( stderr, "in Atracking (owari), mseq2 = %s\n", mseq2[0] );
343 static float MSalign2m2m_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, float **map )
344 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
348 char **aseq1, **aseq2;
350 int lasti, lastj, imid, jmid = 0;
351 float wm = 0.0; /* int ?????? */
353 float *currentw, *previousw;
355 float fpenalty_ex = (float)penalty_ex;
357 // float fpenalty = (float)penalty;
364 float *mjpt, *prept, *curpt;
369 float *initverticalw; /* kufuu sureba iranai */
370 float *lastverticalw; /* kufuu sureba iranai */
387 int *jumpdummi; //muda
388 int *jumpdummj; //muda
389 int jumpi, jumpj = 0;
399 static char ttt1[50000];
400 static char ttt2[50000];
403 localthr = -offset + 500; // 0?
405 ogcp1 = gapinfo[0] + ist;
406 fgcp1 = gapinfo[1] + ist;
407 ogcp2 = gapinfo[2] + jst;
408 fgcp2 = gapinfo[3] + jst;
417 // fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
419 // fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
423 fprintf( stderr, "==== MSalign (depth=%d, reccycle=%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, reccycle, ist, ien, jst, jen );
424 strncpy( ttt1, seq1[0]+ist, lgth1 );
425 strncpy( ttt2, seq2[0]+jst, lgth2 );
428 fprintf( stderr, "seq1 = %s\n", ttt1 );
429 fprintf( stderr, "seq2 = %s\n", ttt2 );
431 if( lgth2 <= 0 ) // lgth1 <= 0 ha?
433 // fprintf( stderr, "\n\n==== jimei\n\n" );
435 for( i=0; i<icyc; i++ )
437 strncpy( mseq1[i], seq1[i]+ist, lgth1 );
440 for( i=0; i<jcyc; i++ )
443 for( j=0; j<lgth1; j++ )
444 strcat( mseq2[i], "-" );
447 // fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
448 // fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
454 aseq1 = AllocateCharMtx( icyc, 0 );
455 aseq2 = AllocateCharMtx( jcyc, 0 );
456 for( i=0; i<icyc; i++ ) aseq1[i] = mseq1[i];
457 for( i=0; i<jcyc; i++ ) aseq2[i] = mseq2[i];
459 aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
460 aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
463 // if( lgth1 < DPTANNI && lgth2 < DPTANNI ) // & dato lgth ==1 no kanousei ga arunode yokunai
464 // if( lgth1 < DPTANNI ) // kore mo lgth2 ga mijikasugiru kanousei ari
465 if( lgth1 < DPTANNI || lgth2 < DPTANNI ) // zettai ni anzen ka?
467 // fprintf( stderr, "==== Going to _tanni\n" );
469 // value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2, gapinfo );
476 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
477 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
479 FreeCharMtx( aseq1 );
480 FreeCharMtx( aseq2 );
483 // fprintf( stderr, "value = %f\n", value );
487 // fprintf( stderr, "Trying to divide the mtx\n" );
489 ll1 = ( (int)(lgth1) ) + 100;
490 ll2 = ( (int)(lgth2) ) + 100;
492 // fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
494 w1 = AllocateFloatVec( ll2+2 );
495 w2 = AllocateFloatVec( ll2+2 );
496 // match = AllocateFloatVec( ll2+2 );
497 midw = AllocateFloatVec( ll2+2 );
498 midn = AllocateFloatVec( ll2+2 );
499 midm = AllocateFloatVec( ll2+2 );
500 jumpbacki = AllocateIntVec( ll2+2 );
501 jumpbackj = AllocateIntVec( ll2+2 );
502 jumpforwi = AllocateIntVec( ll2+2 );
503 jumpforwj = AllocateIntVec( ll2+2 );
504 jumpdummi = AllocateIntVec( ll2+2 );
505 jumpdummj = AllocateIntVec( ll2+2 );
507 initverticalw = AllocateFloatVec( ll1+2 );
508 lastverticalw = AllocateFloatVec( ll1+2 );
510 m = AllocateFloatVec( ll2+2 );
511 mp = AllocateIntVec( ll2+2 );
512 gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
514 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
515 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
518 fprintf( stderr, "succeeded\n" );
522 WMMTX = AllocateFloatMtx( ll1, ll2 );
523 WMMTX2 = AllocateFloatMtx( ll1, ll2 );
526 shortmtx = AllocateShortMtx( ll1, ll2 );
529 fprintf( stderr, "succeeded\n\n" );
538 match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
540 match_ribosum( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
542 for( i=1; i<lgth1+1; i++ )
544 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] );
546 for( j=1; j<lgth2+1; j++ )
548 currentw[j] += ( ogcp2[0] + fgcp2[j-1] );
552 WMMTX[0][0] = initverticalw[0];
553 for( i=1; i<lgth1+1; i++ )
555 WMMTX[i][0] = initverticalw[i];
557 for( j=1; j<lgth2+1; j++ )
559 WMMTX[0][j] = currentw[j];
564 for( j=1; j<lgth2+1; ++j )
566 m[j] = currentw[j-1] + ogcp1[1];
567 // m[j] = currentw[j-1];
571 lastverticalw[0] = currentw[lgth2-1];
575 jumpi = 0; // atode kawaru.
578 for( i=1; i<lasti; i++ )
580 for( i=1; i<=imid; i++ )
584 previousw = currentw;
587 previousw[0] = initverticalw[i-1];
589 match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
590 currentw[0] = initverticalw[i];
596 if( i == imid ) midm[0] = m[0];
598 mi = previousw[0] + ogcp2[1];
599 // mi = previousw[0];
603 // ijppt = ijp[i] + 1;
606 curpt = currentw + 1;
611 for( j=1; j<lastj; j++ )
617 fprintf( stderr, "%5.0f->", wm );
620 // g = mi + fpenalty;
622 fprintf( stderr, "%5.0f?", g );
627 // *ijppt = -( j - mpi );
629 g = *prept + ogcp2[j];
640 g = *mjpt + fgcp1[i-1];
641 // g = *mjpt + fpenalty;
643 fprintf( stderr, "%5.0f?", g );
648 // *ijppt = +( i - *mpjpt );
652 g = *prept + ogcp1[i];
665 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
671 fprintf( stderr, "%5.0f ", wm );
677 WMMTX[i][j] = *curpt;
678 WMMTX2[i][j] = *mjpt;
681 if( i == imid ) //muda
683 jumpbackj[j] = *mpjpt; // muda atode matomeru
684 jumpbacki[j] = mpi; // muda atode matomeru
685 // fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
686 // fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
692 // fprintf( stderr, "m[%d] = %f\n", j, m[j] );
699 lastverticalw[i] = currentw[lgth2-1];
702 WMMTX2[i][lgth2] = m[lgth2-1];
708 for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
709 for( j=0; j<lgth2; j++ ) midm[j] = m[j];
713 // for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
714 // for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
717 for( i=0; i<lgth1; i++ )
719 for( j=0; j<lgth2; j++ )
721 fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
723 fprintf( stderr, "\n" );
725 fprintf( stderr, "\n" );
726 fprintf( stderr, "WMMTX2 = \n" );
727 for( i=0; i<lgth1; i++ )
729 for( j=0; j<lgth2; j++ )
731 fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
733 fprintf( stderr, "\n" );
735 fprintf( stderr, "\n" );
740 match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
741 match_ribosum( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
743 for( i=0; i<lgth1-1; i++ )
745 initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
746 // initverticalw[i] += fpenalty;
748 for( j=0; j<lgth2-1; j++ )
750 currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
751 // currentw[j] += fpenalty;
755 for( i=0; i<lgth1-1; i++ )
757 WMMTX[i][lgth2-1] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
758 // fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
760 for( j=0; j<lgth2-1; j++ )
762 WMMTX[lgth1-1][j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
763 // fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
772 for( j=lgth2-1; j>0; --j )
774 m[j-1] = currentw[j] + fgcp2[lgth2-2];
775 // m[j-1] = currentw[j];
779 // for( j=0; j<lgth2; j++ ) m[j] = 0.0;
780 // m[lgth2-1] ha irunoka?
783 // for( i=lgth1-2; i>=imid; i-- )
786 for( i=lgth1-2; i>-1; i-- )
789 previousw = currentw;
791 previousw[lgth2-1] = initverticalw[i+1];
792 // match_calc( currentw, seq1, seq2, i, lgth2 );
793 match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
795 currentw[lgth2-1] = initverticalw[i];
797 // m[lgth2] = fgcp1[i];
798 // WMMTX2[i][lgth2] += m[lgth2];
799 // fprintf( stderr, "m[] = %f\n", m[lgth2] );
801 mi = previousw[lgth2-1] + fgcp2[lgth2-2];
802 // mi = previousw[lgth2-1];
805 mjpt = m + lgth2 - 2;
806 prept = previousw + lgth2 - 1;
807 curpt = currentw + lgth2 - 2;
808 mpjpt = mp + lgth2 - 2;
811 for( j=lgth2-2; j>-1; j-- )
818 // g = mi + fpenalty;
826 g = *prept + fgcp2[j];
830 // fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
839 // fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
840 g = *mjpt + ogcp1[i+1];
841 // g = *mjpt + fpenalty;
849 // if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
850 g = *prept + fgcp1[i];
862 if( i == jumpi || i == imid - 1 )
864 jumpforwi[j] = ijpi; //muda
865 jumpforwj[j] = ijpj; //muda
866 // fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
867 // fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
869 if( i == imid ) // muda
872 // midm[j+1] += *mjpt + fpenalty; //??????
873 midm[j+1] += *mjpt; //??????
877 // midn[j] += mi + fpenalty; //????
878 midn[j] += mi; //????
883 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
890 // WMMTX2[i][j+1] += *mjpt + fpenalty;
891 WMMTX2[i][j] += *curpt;
900 // fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
901 g = *prept + fgcp1[i];
908 // WMMTX2[i][j+1] += firstm;
910 if( i == imid ) midm[j+1] += firstm;
916 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
917 for( j=2; j<lgth2-1; j++ )
925 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
927 for( j=0; j<lgth2+1; j++ )
935 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
938 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
941 // fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
945 if( jmid > 0 && midn[jmid-1] > wm ) //060413
948 jumpj = jumpbacki[jmid];
950 // fprintf( stderr, "rejump (n)\n" );
952 if( midm[jmid] > wm )
954 jumpi = jumpbackj[jmid];
957 // fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
961 // fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
962 // fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
964 fprintf( stderr, "imid = %d\n", imid );
965 fprintf( stderr, "midn = \n" );
966 for( j=0; j<lgth2; j++ )
968 fprintf( stderr, "% 7.1f ", midn[j] );
970 fprintf( stderr, "\n" );
971 fprintf( stderr, "midw = \n" );
972 for( j=0; j<lgth2; j++ )
974 fprintf( stderr, "% 7.1f ", midw[j] );
976 fprintf( stderr, "\n" );
977 fprintf( stderr, "midm = \n" );
978 for( j=0; j<lgth2; j++ )
980 fprintf( stderr, "% 7.1f ", midm[j] );
982 fprintf( stderr, "\n" );
984 // fprintf( stderr, "maxwm = %f\n", maxwm );
986 if( i == jumpi ) //saki?
988 // fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
989 // fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
992 // fprintf( stderr, "CHUI2!\n" );
999 else if( jmid == lgth2 )
1001 fprintf( stderr, "CHUI1!\n" );
1003 imid=jumpforwi[0]; jmid=lgth2-1;
1006 else if( jmid >= lgth2 )
1008 // fprintf( stderr, "CHUI1!\n" );
1009 jumpi=imid-1; jmid=lgth2;
1015 imid = jumpforwi[jumpj];
1016 jmid = jumpforwj[jumpj];
1019 fprintf( stderr, "jumpi -> %d\n", jumpi );
1020 fprintf( stderr, "jumpj -> %d\n", jumpj );
1021 fprintf( stderr, "imid -> %d\n", imid );
1022 fprintf( stderr, "jmid -> %d\n", jmid );
1034 imid=lgth1-1; jmid=lgth2-1;
1038 // fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1039 // fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1041 // for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1042 // for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1043 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1044 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1048 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
1049 map[i][j] = WMMTX[i][j] / maxwm;
1050 // map[i][j] = WMMTX2[i][j] / maxwm;
1055 for( i=0; i<lgth1; i++ )
1057 float maxpairscore = -9999.9;
1060 for( j=0; j<lgth2; j++ )
1062 if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1065 map12[i].score = tmpscore;
1066 maxpairscore = tmpscore;
1070 for( k=0; k<lgth1; k++ )
1072 if( i == k ) continue;
1073 if( map12[i].score <= WMMTX[k][map12[i].pos] )
1079 map12[i].score = -1.0;
1081 fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, map12[i].pos, map12[i].score, seq1[0][i], seq2[0][map12[i].pos] );
1083 for( j=0; j<lgth2; j++ )
1085 float maxpairscore = -9999.9;
1088 for( i=0; i<lgth1; i++ )
1090 if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1093 map21[j].score = tmpscore;
1094 maxpairscore = tmpscore;
1098 for( k=0; k<lgth2; k++ )
1100 if( j == k ) continue;
1101 if( map21[j].score <= WMMTX[map21[j].pos][k] )
1107 map21[j].score = -1.0;
1109 fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", j, map21[j].pos, map21[j].score, seq2[0][j], seq1[0][map21[j].pos] );
1112 for( i=0; i<lgth1; i++ )
1114 if( map12[i].pos != -1 ) if( map21[map12[i].pos].pos != i )
1115 fprintf( stderr, "ERROR i=%d, but map12[i].pos=%d and map21[map12[i].pos]=%d\n", i, map12[i].pos, map21[map12[i].pos].pos );
1120 fprintf( stderr, "WMMTX = \n" );
1121 for( i=0; i<lgth1; i++ )
1123 fprintf( stderr, "%d ", i );
1124 for( j=0; j<lgth2; j++ )
1126 fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
1128 fprintf( stderr, "\n" );
1130 // fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
1131 for( i=0; i<lgth1; i++ )
1133 fprintf( stderr, "%d ", i );
1134 for( j=0; j<lgth2+1; j++ )
1136 fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
1138 fprintf( stderr, "\n" );
1143 fprintf( stdout, "#WMMTX = \n" );
1144 for( i=0; i<lgth1; i++ )
1146 // fprintf( stdout, "%d ", i );
1147 for( j=0; j<lgth2; j++ )
1149 // if( WMMTX[i][j] > amino_dis['a']['g'] -1 )
1150 fprintf( stdout, "%d %d %8.1f", i, j, WMMTX[i][j] );
1151 if( WMMTX[i][j] == maxwm )
1152 fprintf( stdout, "selected \n" );
1154 fprintf( stdout, "\n" );
1156 fprintf( stdout, "\n" );
1162 fprintf( stderr, "jumpbacki = \n" );
1163 for( j=0; j<lgth2; j++ )
1165 fprintf( stderr, "% 7d ", jumpbacki[j] );
1167 fprintf( stderr, "\n" );
1168 fprintf( stderr, "jumpbackj = \n" );
1169 for( j=0; j<lgth2; j++ )
1171 fprintf( stderr, "% 7d ", jumpbackj[j] );
1173 fprintf( stderr, "\n" );
1174 fprintf( stderr, "jumpforwi = \n" );
1175 for( j=0; j<lgth2; j++ )
1177 fprintf( stderr, "% 7d ", jumpforwi[j] );
1179 fprintf( stderr, "\n" );
1180 fprintf( stderr, "jumpforwj = \n" );
1181 for( j=0; j<lgth2; j++ )
1183 fprintf( stderr, "% 7d ", jumpforwj[j] );
1185 fprintf( stderr, "\n" );
1192 // Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1195 resultlen = strlen( mseq1[0] );
1196 if( alloclen < resultlen || resultlen > N )
1198 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1199 ErrorExit( "LENGTH OVER!\n" );
1206 fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
1207 fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
1209 fprintf( stderr, "imid = %d\n", imid );
1210 fprintf( stderr, "jmid = %d\n", jmid );
1216 FreeFloatVec( initverticalw );
1217 FreeFloatVec( lastverticalw );
1218 FreeFloatVec( midw );
1219 FreeFloatVec( midm );
1220 FreeFloatVec( midn );
1222 FreeIntVec( jumpbacki );
1223 FreeIntVec( jumpbackj );
1224 FreeIntVec( jumpforwi );
1225 FreeIntVec( jumpforwj );
1226 FreeIntVec( jumpdummi );
1227 FreeIntVec( jumpdummj );
1232 FreeFloatMtx( floatwork );
1233 FreeIntMtx( intwork );
1236 FreeFloatMtx( WMMTX );
1237 FreeFloatMtx( WMMTX2 );
1243 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, float **map )
1244 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1249 char **aseq1, **aseq2;
1250 int ll1, ll2, l, len;
1251 int lasti, lastj, imid, jmid=0;
1252 float wm = 0.0; /* int ?????? */
1254 float *currentw, *previousw;
1256 float fpenalty_ex = (float)RNApenalty_ex;
1258 // float fpenalty = (float)penalty;
1265 float *mjpt, *prept, *curpt;
1270 float *initverticalw; /* kufuu sureba iranai */
1271 float *lastverticalw; /* kufuu sureba iranai */
1274 // short **shortmtx;
1288 int *jumpdummi; //muda
1289 int *jumpdummj; //muda
1290 int jumpi, jumpj = 0;
1300 static char ttt1[50000];
1301 static char ttt2[50000];
1304 localthr = -offset + 500; // 0?
1306 ogcp1 = gapinfo[0] + ist;
1307 fgcp1 = gapinfo[1] + ist;
1308 ogcp2 = gapinfo[2] + jst;
1309 fgcp2 = gapinfo[3] + jst;
1318 // fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
1320 // fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
1324 fprintf( stderr, "==== MSalign (depth=%d, reccycle=%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, reccycle, ist, ien, jst, jen );
1325 strncpy( ttt1, seq1[0]+ist, lgth1 );
1326 strncpy( ttt2, seq2[0]+jst, lgth2 );
1329 fprintf( stderr, "seq1 = %s\n", ttt1 );
1330 fprintf( stderr, "seq2 = %s\n", ttt2 );
1332 if( lgth2 <= 0 ) // lgth1 <= 0 ha?
1334 // fprintf( stderr, "\n\n==== jimei\n\n" );
1336 for( i=0; i<icyc; i++ )
1338 strncpy( mseq1[i], seq1[i]+ist, lgth1 );
1339 mseq1[i][lgth1] = 0;
1341 for( i=0; i<jcyc; i++ )
1344 for( j=0; j<lgth1; j++ )
1345 strcat( mseq2[i], "-" );
1348 // fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
1349 // fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
1355 aseq1 = AllocateCharMtx( icyc, 0 );
1356 aseq2 = AllocateCharMtx( jcyc, 0 );
1357 for( i=0; i<icyc; i++ ) aseq1[i] = mseq1[i];
1358 for( i=0; i<jcyc; i++ ) aseq2[i] = mseq2[i];
1360 aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
1361 aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
1364 // if( lgth1 < DPTANNI && lgth2 < DPTANNI ) // & dato lgth ==1 no kanousei ga arunode yokunai
1365 // if( lgth1 < DPTANNI ) // kore mo lgth2 ga mijikasugiru kanousei ari
1366 if( lgth1 < DPTANNI || lgth2 < DPTANNI ) // zettai ni anzen ka?
1368 // fprintf( stderr, "==== Going to _tanni\n" );
1370 // value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2, gapinfo );
1377 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
1378 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
1380 FreeCharMtx( aseq1 );
1381 FreeCharMtx( aseq2 );
1384 // fprintf( stderr, "value = %f\n", value );
1388 // fprintf( stderr, "Trying to divide the mtx\n" );
1390 ll1 = ( (int)(lgth1) ) + 100;
1391 ll2 = ( (int)(lgth2) ) + 100;
1393 // fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
1395 w1 = AllocateFloatVec( ll2+2 );
1396 w2 = AllocateFloatVec( ll2+2 );
1397 // match = AllocateFloatVec( ll2+2 );
1398 midw = AllocateFloatVec( ll2+2 );
1399 midn = AllocateFloatVec( ll2+2 );
1400 midm = AllocateFloatVec( ll2+2 );
1401 jumpbacki = AllocateIntVec( ll2+2 );
1402 jumpbackj = AllocateIntVec( ll2+2 );
1403 jumpforwi = AllocateIntVec( ll2+2 );
1404 jumpforwj = AllocateIntVec( ll2+2 );
1405 jumpdummi = AllocateIntVec( ll2+2 );
1406 jumpdummj = AllocateIntVec( ll2+2 );
1408 initverticalw = AllocateFloatVec( ll1+2 );
1409 lastverticalw = AllocateFloatVec( ll1+2 );
1411 m = AllocateFloatVec( ll2+2 );
1412 mp = AllocateIntVec( ll2+2 );
1413 gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
1415 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
1416 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
1419 fprintf( stderr, "succeeded\n" );
1423 WMMTX = AllocateFloatMtx( ll1, ll2 );
1424 WMMTX2 = AllocateFloatMtx( ll1, ll2 );
1427 shortmtx = AllocateShortMtx( ll1, ll2 );
1430 fprintf( stderr, "succeeded\n\n" );
1439 match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
1441 match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
1443 for( i=1; i<lgth1+1; i++ )
1445 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] );
1447 for( j=1; j<lgth2+1; j++ )
1449 currentw[j] += ( ogcp2[0] + fgcp2[j-1] );
1453 WMMTX[0][0] = initverticalw[0];
1454 for( i=1; i<lgth1+1; i++ )
1456 WMMTX[i][0] = initverticalw[i];
1458 for( j=1; j<lgth2+1; j++ )
1460 WMMTX[0][j] = currentw[j];
1465 for( j=1; j<lgth2+1; ++j )
1467 m[j] = currentw[j-1] + ogcp1[1];
1468 // m[j] = currentw[j-1];
1472 lastverticalw[0] = currentw[lgth2-1];
1476 jumpi = 0; // atode kawaru.
1479 for( i=1; i<lasti; i++ )
1481 for( i=1; i<=imid; i++ )
1485 previousw = currentw;
1488 previousw[0] = initverticalw[i-1];
1490 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
1491 currentw[0] = initverticalw[i];
1495 WMMTX2[i][0] = m[0];
1497 if( i == imid ) midm[0] = m[0];
1499 mi = previousw[0] + ogcp2[1];
1500 // mi = previousw[0];
1504 // ijppt = ijp[i] + 1;
1507 curpt = currentw + 1;
1512 for( j=1; j<lastj; j++ )
1518 fprintf( stderr, "%5.0f->", wm );
1520 g = mi + fgcp2[j-1];
1521 // g = mi + fpenalty;
1523 fprintf( stderr, "%5.0f?", g );
1528 // *ijppt = -( j - mpi );
1530 g = *prept + ogcp2[j];
1541 g = *mjpt + fgcp1[i-1];
1542 // g = *mjpt + fpenalty;
1544 fprintf( stderr, "%5.0f?", g );
1549 // *ijppt = +( i - *mpjpt );
1553 g = *prept + ogcp1[i];
1561 m[j] += fpenalty_ex;
1566 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
1572 fprintf( stderr, "%5.0f ", wm );
1578 WMMTX[i][j] = *curpt;
1579 WMMTX2[i][j] = *mjpt;
1582 if( i == imid ) //muda
1584 jumpbackj[j] = *mpjpt; // muda atode matomeru
1585 jumpbacki[j] = mpi; // muda atode matomeru
1586 // fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
1587 // fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
1593 // fprintf( stderr, "m[%d] = %f\n", j, m[j] );
1600 lastverticalw[i] = currentw[lgth2-1];
1603 WMMTX2[i][lgth2] = m[lgth2-1];
1609 for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
1610 for( j=0; j<lgth2; j++ ) midm[j] = m[j];
1614 // for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
1615 // for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
1618 for( i=0; i<lgth1; i++ )
1620 for( j=0; j<lgth2; j++ )
1622 fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
1624 fprintf( stderr, "\n" );
1626 fprintf( stderr, "\n" );
1627 fprintf( stderr, "WMMTX2 = \n" );
1628 for( i=0; i<lgth1; i++ )
1630 for( j=0; j<lgth2; j++ )
1632 fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
1634 fprintf( stderr, "\n" );
1636 fprintf( stderr, "\n" );
1641 match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
1642 match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
1644 for( i=0; i<lgth1-1; i++ )
1646 initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
1647 // initverticalw[i] += fpenalty;
1649 for( j=0; j<lgth2-1; j++ )
1651 currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
1652 // currentw[j] += fpenalty;
1656 for( i=0; i<lgth1-1; i++ )
1658 WMMTX[i][lgth2-1] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
1659 // fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
1661 for( j=0; j<lgth2-1; j++ )
1663 WMMTX[lgth1-1][j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
1664 // fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
1673 for( j=lgth2-1; j>0; --j )
1675 m[j-1] = currentw[j] + fgcp2[lgth2-2];
1676 // m[j-1] = currentw[j];
1680 // for( j=0; j<lgth2; j++ ) m[j] = 0.0;
1681 // m[lgth2-1] ha irunoka?
1684 // for( i=lgth1-2; i>=imid; i-- )
1685 firstm = -9999999.9;
1687 for( i=lgth1-2; i>-1; i-- )
1690 previousw = currentw;
1692 previousw[lgth2-1] = initverticalw[i+1];
1693 // match_calc( currentw, seq1, seq2, i, lgth2 );
1694 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
1696 currentw[lgth2-1] = initverticalw[i];
1698 // m[lgth2] = fgcp1[i];
1699 // WMMTX2[i][lgth2] += m[lgth2];
1700 // fprintf( stderr, "m[] = %f\n", m[lgth2] );
1702 mi = previousw[lgth2-1] + fgcp2[lgth2-2];
1703 // mi = previousw[lgth2-1];
1706 mjpt = m + lgth2 - 2;
1707 prept = previousw + lgth2 - 1;
1708 curpt = currentw + lgth2 - 2;
1709 mpjpt = mp + lgth2 - 2;
1712 for( j=lgth2-2; j>-1; j-- )
1718 g = mi + ogcp2[j+1];
1719 // g = mi + fpenalty;
1727 g = *prept + fgcp2[j];
1731 // fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
1740 // fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
1741 g = *mjpt + ogcp1[i+1];
1742 // g = *mjpt + fpenalty;
1750 // if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
1751 g = *prept + fgcp1[i];
1760 m[j] += fpenalty_ex;
1763 if( i == jumpi || i == imid - 1 )
1765 jumpforwi[j] = ijpi; //muda
1766 jumpforwj[j] = ijpj; //muda
1767 // fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
1768 // fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
1770 if( i == imid ) // muda
1773 // midm[j+1] += *mjpt + fpenalty; //??????
1774 midm[j+1] += *mjpt; //??????
1778 // midn[j] += mi + fpenalty; //????
1779 midn[j] += mi; //????
1784 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
1791 // WMMTX2[i][j+1] += *mjpt + fpenalty;
1792 WMMTX2[i][j] += *curpt;
1801 // fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
1802 g = *prept + fgcp1[i];
1809 // WMMTX2[i][j+1] += firstm;
1811 if( i == imid ) midm[j+1] += firstm;
1817 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1818 for( j=2; j<lgth2-1; j++ )
1826 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1828 for( j=0; j<lgth2+1; j++ )
1836 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1839 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1842 // fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
1846 if( jmid > 0 && midn[jmid-1] > wm ) //060413
1849 jumpj = jumpbacki[jmid];
1851 // fprintf( stderr, "rejump (n)\n" );
1853 if( midm[jmid] > wm )
1855 jumpi = jumpbackj[jmid];
1858 // fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
1862 // fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
1863 // fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
1865 fprintf( stderr, "imid = %d\n", imid );
1866 fprintf( stderr, "midn = \n" );
1867 for( j=0; j<lgth2; j++ )
1869 fprintf( stderr, "% 7.1f ", midn[j] );
1871 fprintf( stderr, "\n" );
1872 fprintf( stderr, "midw = \n" );
1873 for( j=0; j<lgth2; j++ )
1875 fprintf( stderr, "% 7.1f ", midw[j] );
1877 fprintf( stderr, "\n" );
1878 fprintf( stderr, "midm = \n" );
1879 for( j=0; j<lgth2; j++ )
1881 fprintf( stderr, "% 7.1f ", midm[j] );
1883 fprintf( stderr, "\n" );
1885 // fprintf( stderr, "maxwm = %f\n", maxwm );
1887 if( i == jumpi ) //saki?
1889 // fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
1890 // fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
1893 // fprintf( stderr, "CHUI2!\n" );
1894 jumpj = 0; jmid = 1;
1895 jumpi = firstmp - 1;
1900 else if( jmid == lgth2 )
1902 fprintf( stderr, "CHUI1!\n" );
1904 imid=jumpforwi[0]; jmid=lgth2-1;
1907 else if( jmid >= lgth2 )
1909 // fprintf( stderr, "CHUI1!\n" );
1910 jumpi=imid-1; jmid=lgth2;
1916 imid = jumpforwi[jumpj];
1917 jmid = jumpforwj[jumpj];
1920 fprintf( stderr, "jumpi -> %d\n", jumpi );
1921 fprintf( stderr, "jumpj -> %d\n", jumpj );
1922 fprintf( stderr, "imid -> %d\n", imid );
1923 fprintf( stderr, "jmid -> %d\n", jmid );
1935 imid=lgth1-1; jmid=lgth2-1;
1939 // fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1940 // fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1942 // for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1943 // for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1944 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1945 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1949 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
1950 map[i][j] = WMMTX[i][j] / maxwm;
1951 // map[i][j] = WMMTX2[i][j] / maxwm;
1956 for( i=0; i<lgth1; i++ )
1958 float maxpairscore = -9999.9;
1961 for( j=0; j<lgth2; j++ )
1963 if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1966 map12[i].score = tmpscore;
1967 maxpairscore = tmpscore;
1971 for( k=0; k<lgth1; k++ )
1973 if( i == k ) continue;
1974 if( map12[i].score <= WMMTX[k][map12[i].pos] )
1980 map12[i].score = -1.0;
1982 fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, map12[i].pos, map12[i].score, seq1[0][i], seq2[0][map12[i].pos] );
1984 for( j=0; j<lgth2; j++ )
1986 float maxpairscore = -9999.9;
1989 for( i=0; i<lgth1; i++ )
1991 if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1994 map21[j].score = tmpscore;
1995 maxpairscore = tmpscore;
1999 for( k=0; k<lgth2; k++ )
2001 if( j == k ) continue;
2002 if( map21[j].score <= WMMTX[map21[j].pos][k] )
2008 map21[j].score = -1.0;
2010 fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", j, map21[j].pos, map21[j].score, seq2[0][j], seq1[0][map21[j].pos] );
2013 for( i=0; i<lgth1; i++ )
2015 if( map12[i].pos != -1 ) if( map21[map12[i].pos].pos != i )
2016 fprintf( stderr, "ERROR i=%d, but map12[i].pos=%d and map21[map12[i].pos]=%d\n", i, map12[i].pos, map21[map12[i].pos].pos );
2021 fprintf( stderr, "WMMTX = \n" );
2022 for( i=0; i<lgth1; i++ )
2024 fprintf( stderr, "%d ", i );
2025 for( j=0; j<lgth2; j++ )
2027 fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
2029 fprintf( stderr, "\n" );
2031 // fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
2032 for( i=0; i<lgth1; i++ )
2034 fprintf( stderr, "%d ", i );
2035 for( j=0; j<lgth2+1; j++ )
2037 fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
2039 fprintf( stderr, "\n" );
2045 fprintf( stdout, "#WMMTX = \n" );
2046 for( i=0; i<lgth1; i++ )
2048 // fprintf( stdout, "%d ", i );
2049 for( j=0; j<lgth2; j++ )
2051 // if( WMMTX[i][j] > amino_dis['a']['g'] -1 )
2052 fprintf( stdout, "%d %d %8.1f", i, j, WMMTX[i][j] );
2053 if( WMMTX[i][j] == maxwm )
2054 fprintf( stdout, "selected \n" );
2056 fprintf( stdout, "\n" );
2058 fprintf( stdout, "\n" );
2065 fprintf( stderr, "jumpbacki = \n" );
2066 for( j=0; j<lgth2; j++ )
2068 fprintf( stderr, "% 7d ", jumpbacki[j] );
2070 fprintf( stderr, "\n" );
2071 fprintf( stderr, "jumpbackj = \n" );
2072 for( j=0; j<lgth2; j++ )
2074 fprintf( stderr, "% 7d ", jumpbackj[j] );
2076 fprintf( stderr, "\n" );
2077 fprintf( stderr, "jumpforwi = \n" );
2078 for( j=0; j<lgth2; j++ )
2080 fprintf( stderr, "% 7d ", jumpforwi[j] );
2082 fprintf( stderr, "\n" );
2083 fprintf( stderr, "jumpforwj = \n" );
2084 for( j=0; j<lgth2; j++ )
2086 fprintf( stderr, "% 7d ", jumpforwj[j] );
2088 fprintf( stderr, "\n" );
2095 // Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
2098 resultlen = strlen( mseq1[0] );
2099 if( alloclen < resultlen || resultlen > N )
2101 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
2102 ErrorExit( "LENGTH OVER!\n" );
2109 fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
2110 fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
2112 fprintf( stderr, "imid = %d\n", imid );
2113 fprintf( stderr, "jmid = %d\n", jmid );
2119 FreeFloatVec( initverticalw );
2120 FreeFloatVec( lastverticalw );
2121 FreeFloatVec( midw );
2122 FreeFloatVec( midm );
2123 FreeFloatVec( midn );
2125 FreeIntVec( jumpbacki );
2126 FreeIntVec( jumpbackj );
2127 FreeIntVec( jumpforwi );
2128 FreeIntVec( jumpforwj );
2129 FreeIntVec( jumpdummi );
2130 FreeIntVec( jumpdummj );
2135 FreeFloatMtx( floatwork );
2136 FreeIntMtx( intwork );
2139 FreeFloatMtx( WMMTX );
2140 FreeFloatMtx( WMMTX2 );
2145 // fprintf( stderr, "==== calling myself (first)\n" );
2148 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2149 fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
2151 value = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ist+jumpi, jst, jst+jumpj, alloclen, aseq1, aseq2, depth, gapinfo, map );
2153 fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
2154 fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
2158 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
2159 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
2162 // fprintf( stderr, "====(f) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
2163 // fprintf( stderr, "====(f) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
2165 len = strlen( mseq1[0] );
2166 // fprintf( stderr, "len = %d\n", len );
2167 l = jmid - jumpj - 1;
2168 // fprintf( stderr, "l=%d\n", l );
2171 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
2172 for( i=0; i<icyc; i++ )
2174 strcat( mseq1[i], gaps );
2175 mseq1[i][len+l] = 0;
2177 for( j=0; j<jcyc; j++ )
2179 strncat( mseq2[j], seq2[j]+jst+jumpj+1, l );
2180 mseq2[j][len+l] = 0;
2182 // fprintf( stderr, "penalizing (2) .. %f(%d), %f(%d)\n", ogcp2[jumpj+1], jumpj+1, fgcp2[jmid-1], jmid-1 );
2183 value += ( ogcp2[jumpj+1] + fgcp2[jmid-1] );
2184 // value += fpenalty;
2186 len = strlen( mseq1[0] );
2187 l = imid - jumpi - 1;
2188 // fprintf( stderr, "l=%d\n", l );
2191 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
2192 for( i=0; i<icyc; i++ )
2194 strncat( mseq1[i], seq1[i]+ist+jumpi+1, l );
2195 mseq1[i][len+l] = 0;
2197 for( j=0; j<jcyc; j++ )
2199 strcat( mseq2[j], gaps );
2200 mseq2[j][len+l] = 0;
2203 // for( i=0; i<lgth1; i++ ) fprintf( stderr, "ogcp1[%d] = %f\n", i, ogcp1[i] );
2204 // for( i=0; i<lgth1; i++ ) fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
2207 // fprintf( stderr, "penalizing (1) .. ogcp1[%d] = %f, fgcp1[%d] = %f\n", jumpi+1, ogcp1[jumpi+1], imid-1, fgcp1[imid-1] );
2208 value += ( ogcp1[jumpi+1] + fgcp1[imid-1] );
2209 // value += fpenalty;
2212 for( i=0; i<icyc; i++ ) fprintf( stderr, "after gapfill mseq1[%d]=%s\n", i, mseq1[i] );
2213 for( i=0; i<jcyc; i++ ) fprintf( stderr, "after gapfill mseq2[%d]=%s\n", i, mseq2[i] );
2216 // fprintf( stderr, "==== calling myself (second)\n" );
2219 alnlen = strlen( aseq1[0] );
2220 for( i=0; i<icyc; i++ ) aseq1[i] += alnlen;
2221 for( i=0; i<jcyc; i++ ) aseq2[i] += alnlen;
2225 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2226 fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
2228 value += MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist+imid, ien, jst+jmid, jen, alloclen, aseq1, aseq2, depth, gapinfo, map );
2230 fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
2231 fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
2237 if( value - maxwm > 1 || maxwm - value > 1 )
2239 fprintf( stderr, "WARNING value = %f, but maxwm = %f\n", value, maxwm );
2240 for( i=0; i<icyc; i++ )
2242 fprintf( stderr, ">1-%d\n%s\n", i, mseq1[i] );
2243 fprintf( stderr, "%s\n", aseq1[i] );
2245 for( i=0; i<jcyc; i++ )
2247 fprintf( stderr, ">2-%d\n%s\n", i, mseq2[i] );
2248 fprintf( stderr, "%s\n", aseq2[i] );
2255 fprintf( stderr, "value = %.0f, maxwm = %.0f -> ok\n", value, maxwm );
2261 for( i=0; i<icyc; i++ ) strcat( mseq1[i], aseq1[i] );
2262 for( i=0; i<jcyc; i++ ) strcat( mseq2[i], aseq2[i] );
2265 // fprintf( stderr, "====(s) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
2266 // fprintf( stderr, "====(s) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
2273 FreeCharMtx( aseq1 );
2274 FreeCharMtx( aseq2 );
2282 float Lalignmm_hmout( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, float **map )
2283 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
2289 float wm = 0.0; /* int ?????? */
2301 float fpenalty = (float)RNApenalty;
2309 fprintf( stderr, "eff in SA+++align\n" );
2310 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
2313 nglen1 = seqlen( seq1[0] );
2314 nglen2 = seqlen( seq2[0] );
2316 lgth1 = strlen( seq1[0] );
2317 lgth2 = strlen( seq2[0] );
2319 ll1 = ( (int)(lgth1) ) + 100;
2320 ll2 = ( (int)(lgth2) ) + 100;
2322 mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
2323 mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
2325 gapinfo = AllocateFloatMtx( 4, 0 );
2326 ogcp1 = AllocateFloatVec( ll1+2 );
2327 ogcp2 = AllocateFloatVec( ll2+2 );
2328 fgcp1 = AllocateFloatVec( ll1+2 );
2329 fgcp2 = AllocateFloatVec( ll2+2 );
2332 cpmx1 = AllocateFloatMtx( ll1+2, 27 );
2333 cpmx2 = AllocateFloatMtx( ll2+2, 27 );
2335 for( i=0; i<icyc; i++ )
2337 if( strlen( seq1[i] ) != lgth1 )
2339 fprintf( stderr, "i = %d / %d\n", i, icyc );
2340 fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2344 for( j=0; j<jcyc; j++ )
2346 if( strlen( seq2[j] ) != lgth2 )
2348 fprintf( stderr, "j = %d / %d\n", j, jcyc );
2349 fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2354 MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
2355 MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
2362 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
2363 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
2364 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap2 );
2365 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
2369 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
2370 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
2371 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
2372 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
2376 for( i=0; i<lgth1; i++ )
2378 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
2379 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
2380 // fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
2382 for( i=0; i<lgth2; i++ )
2384 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
2385 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
2386 // fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
2389 for( i=0; i<lgth1; i++ )
2391 ogcp1[i] = 0.5 * fpenalty;
2392 fgcp1[i] = 0.5 * fpenalty;
2394 for( i=0; i<lgth2; i++ )
2396 ogcp2[i] = 0.5 * fpenalty;
2397 fgcp2[i] = 0.5 * fpenalty;
2408 fprintf( stdout, "in MSalignmm.c\n" );
2409 for( i=0; i<icyc; i++ )
2411 fprintf( stdout, ">%d of GROUP1\n", i );
2412 fprintf( stdout, "%s\n", seq1[i] );
2414 for( i=0; i<jcyc; i++ )
2416 fprintf( stdout, ">%d of GROUP2\n", i );
2417 fprintf( stdout, "%s\n", seq2[i] );
2422 wm = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo, map );
2424 fprintf( stderr, " seq1[0] = %s\n", seq1[0] );
2425 fprintf( stderr, " seq2[0] = %s\n", seq2[0] );
2426 fprintf( stderr, "mseq1[0] = %s\n", mseq1[0] );
2427 fprintf( stderr, "mseq2[0] = %s\n", mseq2[0] );
2430 // fprintf( stderr, "wm = %f\n", wm );
2434 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2435 for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
2437 if( seqlen( seq1[0] ) != nglen1 )
2439 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
2440 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2443 if( seqlen( seq2[0] ) != nglen2 )
2445 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
2450 FreeFloatVec( ogcp1 );
2451 FreeFloatVec( ogcp2 );
2452 FreeFloatVec( fgcp1 );
2453 FreeFloatVec( fgcp2 );
2454 FreeFloatMtx( cpmx1 );
2455 FreeFloatMtx( cpmx2 );
2456 free( (void *)gapinfo );
2458 FreeCharMtx( mseq1 );
2459 FreeCharMtx( mseq2 );
2461 lgth1 = strlen( seq1[0] );
2462 lgth2 = strlen( seq2[0] );
2463 for( i=0; i<icyc; i++ )
2465 if( strlen( seq1[i] ) != lgth1 )
2467 fprintf( stderr, "i = %d / %d\n", i, icyc );
2468 fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2472 for( j=0; j<jcyc; j++ )
2474 if( strlen( seq2[j] ) != lgth2 )
2476 fprintf( stderr, "j = %d / %d\n", j, jcyc );
2477 fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2485 float Lalign2m2m_hmout( char **seq1, char **seq2, char **seq1r, char **seq2r, char *dir1, char *dir2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, float **map )
2486 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
2492 float wm = 0.0; /* int ?????? */
2502 float fpenalty = (float)penalty;
2506 fprintf( stderr, "eff in SA+++align\n" );
2507 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
2510 nglen1 = seqlen( seq1[0] );
2511 nglen2 = seqlen( seq2[0] );
2513 lgth1 = strlen( seq1[0] );
2514 lgth2 = strlen( seq2[0] );
2516 ll1 = ( (int)(lgth1) ) + 100;
2517 ll2 = ( (int)(lgth2) ) + 100;
2519 mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
2520 mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
2522 gapinfo = AllocateFloatMtx( 4, 0 );
2523 ogcp1 = AllocateFloatVec( ll1+2 );
2524 ogcp2 = AllocateFloatVec( ll2+2 );
2525 fgcp1 = AllocateFloatVec( ll1+2 );
2526 fgcp2 = AllocateFloatVec( ll2+2 );
2529 cpmx1 = AllocateFloatMtx( ll1+2, 39 );
2530 cpmx2 = AllocateFloatMtx( ll2+2, 39 );
2532 for( i=0; i<icyc; i++ )
2534 if( strlen( seq1[i] ) != lgth1 )
2536 fprintf( stderr, "i = %d / %d\n", i, icyc );
2537 fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2541 for( j=0; j<jcyc; j++ )
2543 if( strlen( seq2[j] ) != lgth2 )
2545 fprintf( stderr, "j = %d / %d\n", j, jcyc );
2546 fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2552 MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
2553 MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
2555 cpmx_ribosum( seq1, seq1r, dir1, cpmx1, eff1, lgth1, icyc );
2556 cpmx_ribosum( seq2, seq2r, dir2, cpmx2, eff2, lgth2, jcyc );
2564 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
2565 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
2566 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap2 );
2567 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
2571 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
2572 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
2573 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
2574 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
2578 for( i=0; i<lgth1; i++ )
2580 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
2581 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
2582 // fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
2584 for( i=0; i<lgth2; i++ )
2586 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
2587 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
2588 // fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
2591 for( i=0; i<lgth1; i++ )
2593 ogcp1[i] = 0.5 * fpenalty;
2594 fgcp1[i] = 0.5 * fpenalty;
2596 for( i=0; i<lgth2; i++ )
2598 ogcp2[i] = 0.5 * fpenalty;
2599 fgcp2[i] = 0.5 * fpenalty;
2610 fprintf( stdout, "in MSalignmm.c\n" );
2611 for( i=0; i<icyc; i++ )
2613 fprintf( stdout, ">%d of GROUP1\n", i );
2614 fprintf( stdout, "%s\n", seq1[i] );
2616 for( i=0; i<jcyc; i++ )
2618 fprintf( stdout, ">%d of GROUP2\n", i );
2619 fprintf( stdout, "%s\n", seq2[i] );
2624 wm = MSalign2m2m_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo, map );
2626 fprintf( stderr, " seq1[0] = %s\n", seq1[0] );
2627 fprintf( stderr, " seq2[0] = %s\n", seq2[0] );
2628 fprintf( stderr, "mseq1[0] = %s\n", mseq1[0] );
2629 fprintf( stderr, "mseq2[0] = %s\n", mseq2[0] );
2632 // fprintf( stderr, "wm = %f\n", wm );
2636 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2637 for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
2639 if( seqlen( seq1[0] ) != nglen1 )
2641 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
2642 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2645 if( seqlen( seq2[0] ) != nglen2 )
2647 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
2652 FreeFloatVec( ogcp1 );
2653 FreeFloatVec( ogcp2 );
2654 FreeFloatVec( fgcp1 );
2655 FreeFloatVec( fgcp2 );
2656 FreeFloatMtx( cpmx1 );
2657 FreeFloatMtx( cpmx2 );
2658 free( (void *)gapinfo );
2660 FreeCharMtx( mseq1 );
2661 FreeCharMtx( mseq2 );
2663 lgth1 = strlen( seq1[0] );
2664 lgth2 = strlen( seq2[0] );
2665 for( i=0; i<icyc; i++ )
2667 if( strlen( seq1[i] ) != lgth1 )
2669 fprintf( stderr, "i = %d / %d\n", i, icyc );
2670 fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2674 for( j=0; j<jcyc; j++ )
2676 if( strlen( seq2[j] ) != lgth2 )
2678 fprintf( stderr, "j = %d / %d\n", j, jcyc );
2679 fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );