binaries/src/mafft/core/Lalignmm.c

   1 #include "mltaln.h"
   2 #include "dp.h"
   3
   4 #define MEMSAVE 1
   5
   6 #define DEBUG 0
   7 #define USE_PENALTY_EX  0
   8 #define STOREWM 1
   9
  10 #define DPTANNI 10
  11
  12 #define LOCAL 0
  13
  14 static int reccycle = 0;
  15
  16 static float localthr;
  17
  18 static void match_ribosum( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
  19 {
  20         int j, k, l;
  21         float scarr[38];
  22         float **cpmxpd = floatwork;
  23         int **cpmxpdn = intwork;
  24         int count = 0;
  25         float *matchpt;
  26         float **cpmxpdpt;
  27         int **cpmxpdnpt;
  28         int cpkd;
  29
  30         if( initialize )
  31         {
  32                 for( j=0; j<lgth2; j++ )
  33                 {
  34                         count = 0;
  35                         for( l=0; l<37; l++ )
  36                         {
  37                                 if( cpmx2[j][l] )
  38                                 {
  39                                         cpmxpd[j][count] = cpmx2[j][l];
  40                                         cpmxpdn[j][count] = l;
  41                                         count++;
  42                                 }
  43                         }
  44                         cpmxpdn[j][count] = -1;
  45                 }
  46         }
  47
  48         for( l=0; l<37; l++ )
  49         {
  50                 scarr[l] = 0.0;
  51                 for( k=0; k<37; k++ )
  52                 {
  53                         scarr[l] += ribosumdis[k][l] * cpmx1[i1][k];
  54                 }
  55         }
  56 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
  57         {
  58                 float *fpt, **fptpt, *fpt2;
  59                 int *ipt, **iptpt;
  60                 fpt2 = match;
  61                 iptpt = cpmxpdn;
  62                 fptpt = cpmxpd;
  63                 while( lgth2-- )
  64                 {
  65                         *fpt2 = 0.0;
  66                         ipt=*iptpt,fpt=*fptpt;
  67                         while( *ipt > -1 )
  68                                 *fpt2 += scarr[*ipt++] * *fpt++;
  69                         fpt2++,iptpt++,fptpt++;
  70                 }
  71         }
  72         for( j=0; j<lgth2; j++ )
  73         {
  74                 match[j] = 0.0;
  75                 for( k=0; cpmxpdn[j][k]>-1; k++ )
  76                         match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
  77         }
  78 #else
  79         matchpt = match;
  80         cpmxpdnpt = cpmxpdn;
  81         cpmxpdpt = cpmxpd;
  82         while( lgth2-- )
  83         {
  84                 *matchpt = 0.0;
  85                 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
  86                         *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
  87                 matchpt++;
  88                 cpmxpdnpt++;
  89                 cpmxpdpt++;
  90         }
  91 #endif
  92 }
  93
  94 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
  95 {
  96         int j, k, l;
  97         float scarr[26];
  98         float **cpmxpd = floatwork;
  99         int **cpmxpdn = intwork;
 100         int count = 0;
 101         float *matchpt;
 102         float **cpmxpdpt;
 103         int **cpmxpdnpt;
 104         int cpkd;
 105
 106         if( initialize )
 107         {
 108                 for( j=0; j<lgth2; j++ )
 109                 {
 110                         count = 0;
 111                         for( l=0; l<26; l++ )
 112                         {
 113                                 if( cpmx2[j][l] )
 114                                 {
 115                                         cpmxpd[j][count] = cpmx2[j][l];
 116                                         cpmxpdn[j][count] = l;
 117                                         count++;
 118                                 }
 119                         }
 120                         cpmxpdn[j][count] = -1;
 121                 }
 122         }
 123
 124         for( l=0; l<26; l++ )
 125         {
 126                 scarr[l] = 0.0;
 127                 for( k=0; k<26; k++ )
 128                 {
 129                         scarr[l] += (n_dis[k][l]-RNAthr) * cpmx1[i1][k];
 130                 }
 131         }
 132 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
 133         {
 134                 float *fpt, **fptpt, *fpt2;
 135                 int *ipt, **iptpt;
 136                 fpt2 = match;
 137                 iptpt = cpmxpdn;
 138                 fptpt = cpmxpd;
 139                 while( lgth2-- )
 140                 {
 141                         *fpt2 = 0.0;
 142                         ipt=*iptpt,fpt=*fptpt;
 143                         while( *ipt > -1 )
 144                                 *fpt2 += scarr[*ipt++] * *fpt++;
 145                         fpt2++,iptpt++,fptpt++;
 146                 }
 147         }
 148         for( j=0; j<lgth2; j++ )
 149         {
 150                 match[j] = 0.0;
 151                 for( k=0; cpmxpdn[j][k]>-1; k++ )
 152                         match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
 153         }
 154 #else
 155         matchpt = match;
 156         cpmxpdnpt = cpmxpdn;
 157         cpmxpdpt = cpmxpd;
 158         while( lgth2-- )
 159         {
 160                 *matchpt = 0.0;
 161                 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
 162                         *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
 163                 matchpt++;
 164                 cpmxpdnpt++;
 165                 cpmxpdpt++;
 166         }
 167 #endif
 168 }
 169
 170 #if 0
 171 static void match_add( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
 172 {
 173         int j, k, l;
 174         float scarr[26];
 175         float **cpmxpd = floatwork;
 176         int **cpmxpdn = intwork;
 177         int count = 0;
 178         float *matchpt;
 179         float **cpmxpdpt;
 180         int **cpmxpdnpt;
 181         int cpkd;
 182
 183
 184         if( initialize )
 185         {
 186                 for( j=0; j<lgth2; j++ )
 187                 {
 188                         count = 0;
 189                         for( l=0; l<26; l++ )
 190                         {
 191                                 if( cpmx2[j][l] )
 192                                 {
 193                                         cpmxpd[j][count] = cpmx2[j][l];
 194                                         cpmxpdn[j][count] = l;
 195                                         count++;
 196                                 }
 197                         }
 198                         cpmxpdn[j][count] = -1;
 199                 }
 200         }
 201
 202         for( l=0; l<26; l++ )
 203         {
 204                 scarr[l] = 0.0;
 205                 for( k=0; k<26; k++ )
 206                 {
 207                         scarr[l] += n_dis[k][l] * cpmx1[i1][k];
 208                 }
 209         }
 210 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
 211         {
 212                 float *fpt, **fptpt, *fpt2;
 213                 int *ipt, **iptpt;
 214                 fpt2 = match;
 215                 iptpt = cpmxpdn;
 216                 fptpt = cpmxpd;
 217                 while( lgth2-- )
 218                 {
 219                         *fpt2 = 0.0;
 220                         ipt=*iptpt,fpt=*fptpt;
 221                         while( *ipt > -1 )
 222                                 *fpt2 += scarr[*ipt++] * *fpt++;
 223                         fpt2++,iptpt++,fptpt++;
 224                 }
 225         }
 226         for( j=0; j<lgth2; j++ )
 227         {
 228                 match[j] = 0.0;
 229                 for( k=0; cpmxpdn[j][k]>-1; k++ )
 230                         match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
 231         }
 232 #else
 233         matchpt = match;
 234         cpmxpdnpt = cpmxpdn;
 235         cpmxpdpt = cpmxpd;
 236         while( lgth2-- )
 237         {
 238 //              *matchpt = 0.0; // add dakara
 239                 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
 240                         *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
 241                 matchpt++;
 242                 cpmxpdnpt++;
 243                 cpmxpdpt++;
 244         }
 245 #endif
 246 }
 247 #endif
 248
 249 #if 0
 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 )
 255 {
 256         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, klim;
 257         char *gaptable1, *gt1bk;
 258         char *gaptable2, *gt2bk;
 259         lgth1 = ien-ist+1;
 260         lgth2 = jen-jst+1;
 261
 262         gt1bk = AllocateCharVec( lgth1+lgth2+1 );
 263         gt2bk = AllocateCharVec( lgth1+lgth2+1 );
 264
 265 #if 0
 266         for( i=0; i<lgth1; i++ )
 267         {
 268                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 269         }
 270 #endif
 271
 272
 273 //      fprintf( stderr, "in Atracking, lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
 274
 275     for( i=0; i<lgth1+1; i++ )
 276     {
 277         ijp[i][0] = i + 1;
 278     }
 279     for( j=0; j<lgth2+1; j++ )
 280     {
 281         ijp[0][j] = -( j + 1 );
 282     }
 283
 284
 285         gaptable1 = gt1bk + lgth1+lgth2;
 286         *gaptable1 = 0;
 287         gaptable2 = gt2bk + lgth1+lgth2;
 288         *gaptable2 = 0;
 289
 290 //      if( lgth2 == 1 ) fprintf( stderr, "in Atracking, mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
 291
 292         iin = lgth1; jin = lgth2;
 293         klim = lgth1+lgth2;
 294         for( k=0; k<=klim; k++ )
 295         {
 296                 if( ijp[iin][jin] < 0 )
 297                 {
 298                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 299                 }
 300                 else if( ijp[iin][jin] > 0 )
 301                 {
 302                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 303                 }
 304                 else
 305                 {
 306                         ifi = iin-1; jfi = jin-1;
 307                 }
 308                 l = iin - ifi;
 309                 while( --l )
 310                 {
 311                         *--gaptable1 = 'o';
 312                         *--gaptable2 = '-';
 313                         k++;
 314                 }
 315                 l= jin - jfi;
 316                 while( --l )
 317                 {
 318                         *--gaptable1 = '-';
 319                         *--gaptable2 = 'o';
 320                         k++;
 321                 }
 322                 if( iin <= 0 || jin <= 0 ) break;
 323                 *--gaptable1 = 'o';
 324                 *--gaptable2 = 'o';
 325                 k++;
 326                 iin = ifi; jin = jfi;
 327
 328         }
 329
 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 );
 332
 333         free( gt1bk );
 334         free( gt2bk );
 335
 336 //      fprintf( stderr, "in Atracking (owari), mseq1 = %s\n", mseq1[0] );
 337 //      fprintf( stderr, "in Atracking (owari), mseq2 = %s\n", mseq2[0] );
 338         return( 0.0 );
 339 }
 340 #endif
 341
 342
 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 */
 345 {
 346         float value = 0.0;
 347         register int i, j;
 348         char **aseq1, **aseq2;
 349         int ll1, ll2;
 350         int lasti, lastj, imid, jmid = 0;
 351         float wm = 0.0;   /* int ?????? */
 352         float g;
 353         float *currentw, *previousw;
 354 #if USE_PENALTY_EX
 355         float fpenalty_ex = (float)penalty_ex;
 356 #endif
 357 //      float fpenalty = (float)penalty;
 358         float *wtmp;
 359 //      short *ijppt;
 360         int *mpjpt;
 361 //      short **ijp;
 362         int *mp;
 363         int mpi;
 364         float *mjpt, *prept, *curpt;
 365         float mi;
 366         float *m;
 367         float *w1, *w2;
 368 //      float *match;
 369         float *initverticalw;    /* kufuu sureba iranai */
 370         float *lastverticalw;    /* kufuu sureba iranai */
 371         int **intwork;
 372         float **floatwork;
 373 //      short **shortmtx;
 374 #if STOREWM
 375         float **WMMTX;
 376         float **WMMTX2;
 377 #endif
 378         float *midw;
 379         float *midm;
 380         float *midn;
 381         int lgth1, lgth2;
 382         float maxwm = 0.0;
 383         int *jumpforwi;
 384         int *jumpforwj;
 385         int *jumpbacki;
 386         int *jumpbackj;
 387         int *jumpdummi; //muda
 388         int *jumpdummj; //muda
 389         int jumpi, jumpj = 0;
 390         char *gaps;
 391         int ijpi, ijpj;
 392         float *ogcp1;
 393         float *fgcp1;
 394         float *ogcp2;
 395         float *fgcp2;
 396         float firstm;
 397         int firstmp;
 398 #if 0
 399         static char ttt1[50000];
 400         static char ttt2[50000];
 401 #endif
 402
 403         localthr = -offset + 500; // 0?
 404
 405         ogcp1 = gapinfo[0] + ist;
 406         fgcp1 = gapinfo[1] + ist;
 407         ogcp2 = gapinfo[2] + jst;
 408         fgcp2 = gapinfo[3] + jst;
 409
 410         depth++;
 411         reccycle++;
 412
 413         lgth1 = ien-ist+1;
 414         lgth2 = jen-jst+1;
 415
 416 //      if( lgth1 < 5 )
 417 //              fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
 418 //      if( lgth2 < 5 )
 419 //              fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
 420 //
 421
 422 #if 0
 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 );
 426         ttt1[lgth1] = 0;
 427         ttt2[lgth2] = 0;
 428         fprintf( stderr, "seq1 = %s\n", ttt1 );
 429         fprintf( stderr, "seq2 = %s\n", ttt2 );
 430 #endif
 431         if( lgth2 <= 0 ) // lgth1 <= 0 ha?
 432         {
 433 //              fprintf( stderr, "\n\n==== jimei\n\n" );
 434 //              exit( 1 );
 435                 for( i=0; i<icyc; i++ )
 436                 {
 437                         strncpy( mseq1[i], seq1[i]+ist, lgth1 );
 438                         mseq1[i][lgth1] = 0;
 439                 }
 440                 for( i=0; i<jcyc; i++ )
 441                 {
 442                         mseq2[i][0] = 0;
 443                         for( j=0; j<lgth1; j++ )
 444                                 strcat( mseq2[i], "-" );
 445                 }
 446
 447 //              fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
 448 //              fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
 449
 450                 return( 0.0 );
 451         }
 452
 453 #if MEMSAVE
 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];
 458 #else
 459         aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
 460         aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
 461 #endif
 462
 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?
 466         {
 467 //              fprintf( stderr, "==== Going to _tanni\n" );
 468
 469 //              value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2, gapinfo );
 470
 471
 472 #if MEMSAVE
 473                 free( aseq1 );
 474                 free( aseq2 );
 475 #else
 476                 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
 477                 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
 478
 479                 FreeCharMtx( aseq1 );
 480                 FreeCharMtx( aseq2 );
 481 #endif
 482
 483 //              fprintf( stderr, "value = %f\n", value );
 484
 485                 return( value );
 486         }
 487 //      fprintf( stderr, "Trying to divide the mtx\n" );
 488
 489         ll1 = ( (int)(lgth1) ) + 100;
 490         ll2 = ( (int)(lgth2) ) + 100;
 491
 492 //      fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
 493
 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 );
 506
 507         initverticalw = AllocateFloatVec( ll1+2 );
 508         lastverticalw = AllocateFloatVec( ll1+2 );
 509
 510         m = AllocateFloatVec( ll2+2 );
 511         mp = AllocateIntVec( ll2+2 );
 512         gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
 513
 514         floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
 515         intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
 516
 517 #if DEBUG
 518         fprintf( stderr, "succeeded\n" );
 519 #endif
 520
 521 #if STOREWM
 522         WMMTX = AllocateFloatMtx( ll1, ll2 );
 523         WMMTX2 = AllocateFloatMtx( ll1, ll2 );
 524 #endif
 525 #if 0
 526         shortmtx = AllocateShortMtx( ll1, ll2 );
 527
 528 #if DEBUG
 529         fprintf( stderr, "succeeded\n\n" );
 530 #endif
 531
 532         ijp = shortmtx;
 533 #endif
 534
 535         currentw = w1;
 536         previousw = w2;
 537
 538         match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
 539
 540         match_ribosum( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
 541
 542         for( i=1; i<lgth1+1; i++ )
 543         {
 544                 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] );
 545         }
 546         for( j=1; j<lgth2+1; j++ )
 547         {
 548                 currentw[j] += ( ogcp2[0] + fgcp2[j-1] );
 549         }
 550
 551 #if STOREWM
 552         WMMTX[0][0] = initverticalw[0];
 553         for( i=1; i<lgth1+1; i++ )
 554         {
 555                 WMMTX[i][0] = initverticalw[i];
 556         }
 557         for( j=1; j<lgth2+1; j++ )
 558         {
 559                 WMMTX[0][j] = currentw[j];
 560         }
 561 #endif
 562
 563
 564         for( j=1; j<lgth2+1; ++j )
 565         {
 566                 m[j] = currentw[j-1] + ogcp1[1];
 567 //              m[j] = currentw[j-1];
 568                 mp[j] = 0;
 569         }
 570
 571         lastverticalw[0] = currentw[lgth2-1];
 572
 573         imid = lgth1 * 0.5;
 574
 575         jumpi = 0; // atode kawaru.
 576         lasti = lgth1+1;
 577 #if STOREWM
 578         for( i=1; i<lasti; i++ )
 579 #else
 580         for( i=1; i<=imid; i++ )
 581 #endif
 582         {
 583                 wtmp = previousw;
 584                 previousw = currentw;
 585                 currentw = wtmp;
 586
 587                 previousw[0] = initverticalw[i-1];
 588
 589                 match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
 590                 currentw[0] = initverticalw[i];
 591
 592                 m[0] = ogcp1[i];
 593 #if STOREM
 594                 WMMTX2[i][0] = m[0];
 595 #endif
 596                 if( i == imid ) midm[0] = m[0];
 597
 598                 mi = previousw[0] + ogcp2[1];
 599 //              mi = previousw[0];
 600                 mpi = 0;
 601
 602
 603 //              ijppt = ijp[i] + 1;
 604                 mjpt = m + 1;
 605                 prept = previousw;
 606                 curpt = currentw + 1;
 607                 mpjpt = mp + 1;
 608
 609
 610                 lastj = lgth2+1;
 611                 for( j=1; j<lastj; j++ )
 612                 {
 613
 614                         wm = *prept;
 615
 616 #if 0
 617                         fprintf( stderr, "%5.0f->", wm );
 618 #endif
 619                         g = mi + fgcp2[j-1];
 620 //                      g = mi + fpenalty;
 621 #if 0
 622                         fprintf( stderr, "%5.0f?", g );
 623 #endif
 624                         if( g > wm )
 625                         {
 626                                 wm = g;
 627 //                              *ijppt = -( j - mpi );
 628                         }
 629                         g = *prept + ogcp2[j];
 630 //                      g = *prept;
 631                         if( g >= mi )
 632                         {
 633                                 mi = g;
 634                                 mpi = j-1;
 635                         }
 636 #if USE_PENALTY_EX
 637                         mi += fpenalty_ex;
 638 #endif
 639
 640                         g = *mjpt + fgcp1[i-1];
 641 //                      g = *mjpt + fpenalty;
 642 #if 0
 643                         fprintf( stderr, "%5.0f?", g );
 644 #endif
 645                         if( g > wm )
 646                         {
 647                                 wm = g;
 648 //                              *ijppt = +( i - *mpjpt );
 649                         }
 650
 651
 652                         g = *prept + ogcp1[i];
 653 //                      g = *prept;
 654                         if( g >= *mjpt )
 655                         {
 656                                 *mjpt = g;
 657                                 *mpjpt = i-1;
 658                         }
 659 #if USE_PENALTY_EX
 660                         m[j] += fpenalty_ex;
 661 #endif
 662 #if LOCAL
 663             if( wm < localthr )
 664             {
 665 //                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
 666                                 wm = 0;
 667             }
 668 #endif
 669
 670 #if 0
 671                         fprintf( stderr, "%5.0f ", wm );
 672 #endif
 673                         *curpt += wm;
 674
 675
 676 #if STOREWM
 677                         WMMTX[i][j] = *curpt;
 678                         WMMTX2[i][j] = *mjpt;
 679 #endif
 680
 681                         if( i == imid ) //muda
 682                         {
 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 );
 687                                 midw[j] = *curpt;
 688                                 midm[j] = *mjpt;
 689                                 midn[j] = mi;
 690                         }
 691
 692 //                      fprintf( stderr, "m[%d] = %f\n", j, m[j] );
 693                         mjpt++;
 694                         prept++;
 695                         mpjpt++;
 696                         curpt++;
 697
 698                 }
 699                 lastverticalw[i] = currentw[lgth2-1];
 700
 701 #if STOREWM
 702                 WMMTX2[i][lgth2] = m[lgth2-1];
 703 #endif
 704
 705 #if 0  // ue
 706                 if( i == imid )
 707                 {
 708                         for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
 709                         for( j=0; j<lgth2; j++ ) midm[j] = m[j];
 710                 }
 711 #endif
 712         }
 713 //      for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
 714 //      for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
 715
 716 #if 0
 717     for( i=0; i<lgth1; i++ )
 718     {
 719         for( j=0; j<lgth2; j++ )
 720         {
 721             fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
 722         }
 723         fprintf( stderr, "\n" );
 724     }
 725         fprintf( stderr, "\n" );
 726         fprintf( stderr, "WMMTX2 = \n" );
 727     for( i=0; i<lgth1; i++ )
 728     {
 729         for( j=0; j<lgth2; j++ )
 730         {
 731             fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
 732         }
 733         fprintf( stderr, "\n" );
 734     }
 735         fprintf( stderr, "\n" );
 736 #endif
 737
 738 // gyakudp
 739
 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 );
 742
 743         for( i=0; i<lgth1-1; i++ )
 744         {
 745                 initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
 746 //              initverticalw[i] += fpenalty;
 747         }
 748         for( j=0; j<lgth2-1; j++ )
 749         {
 750                 currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
 751 //              currentw[j] += fpenalty;
 752         }
 753
 754 #if STOREWM
 755         for( i=0; i<lgth1-1; i++ )
 756         {
 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] );
 759         }
 760         for( j=0; j<lgth2-1; j++ )
 761         {
 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] );
 764         }
 765 #endif
 766
 767
 768
 769
 770
 771
 772         for( j=lgth2-1; j>0; --j )
 773         {
 774                 m[j-1] = currentw[j] + fgcp2[lgth2-2];
 775 //              m[j-1] = currentw[j];
 776                 mp[j] = lgth1-1;
 777         }
 778
 779 //      for( j=0; j<lgth2; j++ ) m[j] = 0.0;
 780         // m[lgth2-1] ha irunoka?
 781
 782
 783 //      for( i=lgth1-2; i>=imid; i-- )
 784         firstm = -9999999.9;
 785         firstmp = lgth1-1;
 786         for( i=lgth1-2; i>-1; i-- )
 787         {
 788                 wtmp = previousw;
 789                 previousw = currentw;
 790                 currentw = wtmp;
 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 );
 794
 795                 currentw[lgth2-1] = initverticalw[i];
 796
 797 //              m[lgth2] = fgcp1[i];
 798 //              WMMTX2[i][lgth2] += m[lgth2];
 799 //              fprintf( stderr, "m[] = %f\n", m[lgth2] );
 800
 801                 mi = previousw[lgth2-1] + fgcp2[lgth2-2];
 802 //              mi = previousw[lgth2-1];
 803                 mpi = lgth2 - 1;
 804
 805                 mjpt = m + lgth2 - 2;
 806                 prept = previousw + lgth2 - 1;
 807                 curpt = currentw + lgth2 - 2;
 808                 mpjpt = mp + lgth2 - 2;
 809
 810
 811                 for( j=lgth2-2; j>-1; j-- )
 812                 {
 813                         wm = *prept;
 814                         ijpi = i+1;
 815                         ijpj = j+1;
 816
 817                         g = mi + ogcp2[j+1];
 818 //                      g = mi + fpenalty;
 819                         if( g > wm )
 820                         {
 821                                 wm = g;
 822                                 ijpj = mpi;
 823                                 ijpi = i+1;
 824                         }
 825
 826                         g = *prept + fgcp2[j];
 827 //                      g = *prept;
 828                         if( g >= mi )
 829                         {
 830 //                              fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
 831                                 mi = g;
 832                                 mpi = j + 1;
 833                         }
 834
 835 #if USE_PENALTY_EX
 836                         mi += fpenalty_ex;
 837 #endif
 838
 839 //                      fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
 840                         g = *mjpt + ogcp1[i+1];
 841 //                      g = *mjpt + fpenalty;
 842                         if( g > wm )
 843                         {
 844                                 wm = g;
 845                                 ijpi = *mpjpt;
 846                                 ijpj = j+1;
 847                         }
 848
 849 //                      if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
 850                         g = *prept + fgcp1[i];
 851 //                      g = *prept;
 852                         if( g >= *mjpt )
 853                         {
 854                                 *mjpt = g;
 855                                 *mpjpt = i + 1;
 856                         }
 857
 858 #if USE_PENALTY_EX
 859                         m[j] += fpenalty_ex;
 860 #endif
 861
 862                         if( i == jumpi || i == imid - 1 )
 863                         {
 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 );
 868                         }
 869                         if( i == imid ) // muda
 870                         {
 871                                 midw[j] += wm;
 872 //                              midm[j+1] += *mjpt + fpenalty; //??????
 873                                 midm[j+1] += *mjpt; //??????
 874                         }
 875                         if( i == imid - 1 )
 876                         {
 877 //                              midn[j] += mi + fpenalty;  //????
 878                                 midn[j] += mi;  //????
 879                         }
 880 #if LOCAL
 881             if( wm < localthr )
 882             {
 883 //                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
 884                                 wm = 0;
 885             }
 886 #endif
 887
 888 #if STOREWM
 889                         WMMTX[i][j] += wm;
 890 //                      WMMTX2[i][j+1] += *mjpt + fpenalty;
 891                         WMMTX2[i][j] += *curpt;
 892 #endif
 893                         *curpt += wm;
 894
 895                         mjpt--;
 896                         prept--;
 897                         mpjpt--;
 898                         curpt--;
 899                 }
 900 //              fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
 901                 g = *prept + fgcp1[i];
 902                 if( firstm < g )
 903                 {
 904                         firstm = g;
 905                         firstmp = i + 1;
 906                 }
 907 #if STOREWM
 908 //              WMMTX2[i][j+1] += firstm;
 909 #endif
 910                 if( i == imid ) midm[j+1] += firstm;
 911
 912                 if( i == imid - 1 )
 913                 {
 914                         maxwm = midw[1];
 915                         jmid = 0;
 916 //                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
 917                         for( j=2; j<lgth2-1; j++ )
 918                         {
 919                                 wm = midw[j];
 920                                 if( wm > maxwm )
 921                                 {
 922                                         jmid = j;
 923                                         maxwm = wm;
 924                                 }
 925 //                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
 926                         }
 927                         for( j=0; j<lgth2+1; j++ )
 928                         {
 929                                 wm = midm[j];
 930                                 if( wm > maxwm )
 931                                 {
 932                                         jmid = j;
 933                                         maxwm = wm;
 934                                 }
 935 //                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
 936                         }
 937
 938 //                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
 939
 940
 941 //                      fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
 942                         wm = midw[jmid];
 943                         jumpi = imid-1;
 944                         jumpj = jmid-1;
 945                         if( jmid > 0 && midn[jmid-1] > wm ) //060413
 946                         {
 947                                 jumpi = imid-1;
 948                                 jumpj = jumpbacki[jmid];
 949                                 wm = midn[jmid-1];
 950 //                              fprintf( stderr, "rejump (n)\n" );
 951                         }
 952                         if( midm[jmid] > wm )
 953                         {
 954                                 jumpi = jumpbackj[jmid];
 955                                 jumpj = jmid-1;
 956                                 wm = midm[jmid];
 957 //                              fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
 958                         }
 959
 960
 961 //                      fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
 962 //                      fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
 963 #if 0
 964                         fprintf( stderr, "imid = %d\n", imid );
 965                         fprintf( stderr, "midn = \n" );
 966                         for( j=0; j<lgth2; j++ )
 967                         {
 968                                 fprintf( stderr, "% 7.1f ", midn[j] );
 969                         }
 970                         fprintf( stderr, "\n" );
 971                         fprintf( stderr, "midw = \n" );
 972                         for( j=0; j<lgth2; j++ )
 973                         {
 974                                 fprintf( stderr, "% 7.1f ", midw[j] );
 975                         }
 976                         fprintf( stderr, "\n" );
 977                         fprintf( stderr, "midm = \n" );
 978                         for( j=0; j<lgth2; j++ )
 979                         {
 980                                 fprintf( stderr, "% 7.1f ", midm[j] );
 981                         }
 982                         fprintf( stderr, "\n" );
 983 #endif
 984 //                      fprintf( stderr, "maxwm = %f\n", maxwm );
 985                 }
 986                 if( i == jumpi ) //saki?
 987                 {
 988 //                      fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
 989 //                      fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
 990                         if( jmid == 0 )
 991                         {
 992 //                              fprintf( stderr, "CHUI2!\n" );
 993                                 jumpj = 0; jmid = 1;
 994                                 jumpi = firstmp - 1;
 995                                 imid = firstmp;
 996                         }
 997
 998 #if 0
 999                         else if( jmid == lgth2 )
1000                         {
1001                                 fprintf( stderr, "CHUI1!\n" );
1002                                 jumpi=0; jumpj=0;
1003                                 imid=jumpforwi[0]; jmid=lgth2-1;
1004                         }
1005 #else // 060414
1006                         else if( jmid >= lgth2 )
1007                         {
1008 //                              fprintf( stderr, "CHUI1!\n" );
1009                                 jumpi=imid-1; jmid=lgth2;
1010                                 jumpj = lgth2-1;
1011                         }
1012 #endif
1013                         else
1014                         {
1015                                 imid = jumpforwi[jumpj];
1016                                 jmid = jumpforwj[jumpj];
1017                         }
1018 #if 0
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 );
1023 #endif
1024
1025 #if STOREWM
1026 //                      break;
1027 #else
1028                         break;
1029 #endif
1030                 }
1031         }
1032 #if 0
1033                 jumpi=0; jumpj=0;
1034                 imid=lgth1-1; jmid=lgth2-1;
1035         }
1036 #endif
1037
1038 //      fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1039 //      fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1040
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];
1045
1046
1047
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;
1051
1052 #if STOREWM
1053
1054 #if 0
1055         for( i=0; i<lgth1; i++ )
1056         {
1057                 float maxpairscore = -9999.9;
1058                 float tmpscore;
1059
1060                 for( j=0; j<lgth2; j++ )
1061                 {
1062                         if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1063                         {
1064                                 map12[i].pos = j;
1065                                 map12[i].score = tmpscore;
1066                                 maxpairscore = tmpscore;
1067                         }
1068                 }
1069
1070                 for( k=0; k<lgth1; k++ )
1071                 {
1072                         if( i == k ) continue;
1073                         if( map12[i].score <= WMMTX[k][map12[i].pos] )
1074                                 break;
1075                 }
1076                 if( k != lgth1 )
1077                 {
1078                         map12[i].pos = -1;
1079                         map12[i].score = -1.0;
1080                 }
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] );
1082         }
1083         for( j=0; j<lgth2; j++ )
1084         {
1085                 float maxpairscore = -9999.9;
1086                 float tmpscore;
1087
1088                 for( i=0; i<lgth1; i++ )
1089                 {
1090                         if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1091                         {
1092                                 map21[j].pos = i;
1093                                 map21[j].score = tmpscore;
1094                                 maxpairscore = tmpscore;
1095                         }
1096                 }
1097
1098                 for( k=0; k<lgth2; k++ )
1099                 {
1100                         if( j == k ) continue;
1101                         if( map21[j].score <= WMMTX[map21[j].pos][k] )
1102                                 break;
1103                 }
1104                 if( k != lgth2 )
1105                 {
1106                         map21[j].pos = -1;
1107                         map21[j].score = -1.0;
1108                 }
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] );
1110         }
1111
1112         for( i=0; i<lgth1; i++ )
1113         {
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 );
1116         }
1117 #endif
1118
1119 #if 0
1120         fprintf( stderr, "WMMTX = \n" );
1121     for( i=0; i<lgth1; i++ )
1122     {
1123         fprintf( stderr, "%d ", i );
1124         for( j=0; j<lgth2; j++ )
1125         {
1126             fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
1127         }
1128         fprintf( stderr, "\n" );
1129     }
1130 //      fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
1131     for( i=0; i<lgth1; i++ )
1132     {
1133         fprintf( stderr, "%d ", i );
1134         for( j=0; j<lgth2+1; j++ )
1135         {
1136             fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
1137         }
1138         fprintf( stderr, "\n" );
1139     }
1140 #endif
1141
1142 #if 0
1143     fprintf( stdout, "#WMMTX = \n" );
1144     for( i=0; i<lgth1; i++ )
1145     {
1146 //        fprintf( stdout, "%d ", i );
1147         for( j=0; j<lgth2; j++ )
1148         {
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" );
1153                         else
1154                 fprintf( stdout, "\n" );
1155         }
1156         fprintf( stdout, "\n" );
1157     }
1158 #endif
1159
1160 #if 0
1161
1162         fprintf( stderr, "jumpbacki = \n" );
1163         for( j=0; j<lgth2; j++ )
1164         {
1165                 fprintf( stderr, "% 7d ", jumpbacki[j] );
1166         }
1167         fprintf( stderr, "\n" );
1168         fprintf( stderr, "jumpbackj = \n" );
1169         for( j=0; j<lgth2; j++ )
1170         {
1171                 fprintf( stderr, "% 7d ", jumpbackj[j] );
1172         }
1173         fprintf( stderr, "\n" );
1174         fprintf( stderr, "jumpforwi = \n" );
1175         for( j=0; j<lgth2; j++ )
1176         {
1177                 fprintf( stderr, "% 7d ", jumpforwi[j] );
1178         }
1179         fprintf( stderr, "\n" );
1180         fprintf( stderr, "jumpforwj = \n" );
1181         for( j=0; j<lgth2; j++ )
1182         {
1183                 fprintf( stderr, "% 7d ", jumpforwj[j] );
1184         }
1185         fprintf( stderr, "\n" );
1186 #endif
1187
1188
1189 #endif
1190
1191
1192 //      Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1193
1194 #if 0 // irukamo
1195         resultlen = strlen( mseq1[0] );
1196         if( alloclen < resultlen || resultlen > N )
1197         {
1198                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1199                 ErrorExit( "LENGTH OVER!\n" );
1200         }
1201 #endif
1202
1203
1204
1205 #if 0
1206         fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
1207         fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
1208
1209         fprintf( stderr, "imid = %d\n", imid );
1210         fprintf( stderr, "jmid = %d\n", jmid );
1211 #endif
1212
1213
1214         FreeFloatVec( w1 );
1215         FreeFloatVec( w2 );
1216         FreeFloatVec( initverticalw );
1217         FreeFloatVec( lastverticalw );
1218         FreeFloatVec( midw );
1219         FreeFloatVec( midm );
1220         FreeFloatVec( midn );
1221
1222         FreeIntVec( jumpbacki );
1223         FreeIntVec( jumpbackj );
1224         FreeIntVec( jumpforwi );
1225         FreeIntVec( jumpforwj );
1226         FreeIntVec( jumpdummi );
1227         FreeIntVec( jumpdummj );
1228
1229         FreeFloatVec( m );
1230         FreeIntVec( mp );
1231
1232         FreeFloatMtx( floatwork );
1233         FreeIntMtx( intwork );
1234
1235 #if STOREWM
1236         FreeFloatMtx( WMMTX );
1237         FreeFloatMtx( WMMTX2 );
1238 #endif
1239
1240         return( value );
1241
1242 }
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 */
1245 {
1246         int alnlen;
1247         float value = 0.0;
1248         register int i, j;
1249         char **aseq1, **aseq2;
1250         int ll1, ll2, l, len;
1251         int lasti, lastj, imid, jmid=0;
1252         float wm = 0.0;   /* int ?????? */
1253         float g;
1254         float *currentw, *previousw;
1255 #if USE_PENALTY_EX
1256         float fpenalty_ex = (float)RNApenalty_ex;
1257 #endif
1258 //      float fpenalty = (float)penalty;
1259         float *wtmp;
1260 //      short *ijppt;
1261         int *mpjpt;
1262 //      short **ijp;
1263         int *mp;
1264         int mpi;
1265         float *mjpt, *prept, *curpt;
1266         float mi;
1267         float *m;
1268         float *w1, *w2;
1269 //      float *match;
1270         float *initverticalw;    /* kufuu sureba iranai */
1271         float *lastverticalw;    /* kufuu sureba iranai */
1272         int **intwork;
1273         float **floatwork;
1274 //      short **shortmtx;
1275 #if STOREWM
1276         float **WMMTX;
1277         float **WMMTX2;
1278 #endif
1279         float *midw;
1280         float *midm;
1281         float *midn;
1282         int lgth1, lgth2;
1283         float maxwm = 0.0;
1284         int *jumpforwi;
1285         int *jumpforwj;
1286         int *jumpbacki;
1287         int *jumpbackj;
1288         int *jumpdummi; //muda
1289         int *jumpdummj; //muda
1290         int jumpi, jumpj = 0;
1291         char *gaps;
1292         int ijpi, ijpj;
1293         float *ogcp1;
1294         float *fgcp1;
1295         float *ogcp2;
1296         float *fgcp2;
1297         float firstm;
1298         int firstmp;
1299 #if 0
1300         static char ttt1[50000];
1301         static char ttt2[50000];
1302 #endif
1303
1304         localthr = -offset + 500; // 0?
1305
1306         ogcp1 = gapinfo[0] + ist;
1307         fgcp1 = gapinfo[1] + ist;
1308         ogcp2 = gapinfo[2] + jst;
1309         fgcp2 = gapinfo[3] + jst;
1310
1311         depth++;
1312         reccycle++;
1313
1314         lgth1 = ien-ist+1;
1315         lgth2 = jen-jst+1;
1316
1317 //      if( lgth1 < 5 )
1318 //              fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
1319 //      if( lgth2 < 5 )
1320 //              fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
1321 //
1322
1323 #if 0
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 );
1327         ttt1[lgth1] = 0;
1328         ttt2[lgth2] = 0;
1329         fprintf( stderr, "seq1 = %s\n", ttt1 );
1330         fprintf( stderr, "seq2 = %s\n", ttt2 );
1331 #endif
1332         if( lgth2 <= 0 ) // lgth1 <= 0 ha?
1333         {
1334 //              fprintf( stderr, "\n\n==== jimei\n\n" );
1335 //              exit( 1 );
1336                 for( i=0; i<icyc; i++ )
1337                 {
1338                         strncpy( mseq1[i], seq1[i]+ist, lgth1 );
1339                         mseq1[i][lgth1] = 0;
1340                 }
1341                 for( i=0; i<jcyc; i++ )
1342                 {
1343                         mseq2[i][0] = 0;
1344                         for( j=0; j<lgth1; j++ )
1345                                 strcat( mseq2[i], "-" );
1346                 }
1347
1348 //              fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
1349 //              fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
1350
1351                 return( 0.0 );
1352         }
1353
1354 #if MEMSAVE
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];
1359 #else
1360         aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
1361         aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
1362 #endif
1363
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?
1367         {
1368 //              fprintf( stderr, "==== Going to _tanni\n" );
1369
1370 //              value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2, gapinfo );
1371
1372
1373 #if MEMSAVE
1374                 free( aseq1 );
1375                 free( aseq2 );
1376 #else
1377                 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
1378                 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
1379
1380                 FreeCharMtx( aseq1 );
1381                 FreeCharMtx( aseq2 );
1382 #endif
1383
1384 //              fprintf( stderr, "value = %f\n", value );
1385
1386                 return( value );
1387         }
1388 //      fprintf( stderr, "Trying to divide the mtx\n" );
1389
1390         ll1 = ( (int)(lgth1) ) + 100;
1391         ll2 = ( (int)(lgth2) ) + 100;
1392
1393 //      fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
1394
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 );
1407
1408         initverticalw = AllocateFloatVec( ll1+2 );
1409         lastverticalw = AllocateFloatVec( ll1+2 );
1410
1411         m = AllocateFloatVec( ll2+2 );
1412         mp = AllocateIntVec( ll2+2 );
1413         gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
1414
1415         floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
1416         intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
1417
1418 #if DEBUG
1419         fprintf( stderr, "succeeded\n" );
1420 #endif
1421
1422 #if STOREWM
1423         WMMTX = AllocateFloatMtx( ll1, ll2 );
1424         WMMTX2 = AllocateFloatMtx( ll1, ll2 );
1425 #endif
1426 #if 0
1427         shortmtx = AllocateShortMtx( ll1, ll2 );
1428
1429 #if DEBUG
1430         fprintf( stderr, "succeeded\n\n" );
1431 #endif
1432
1433         ijp = shortmtx;
1434 #endif
1435
1436         currentw = w1;
1437         previousw = w2;
1438
1439         match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
1440
1441         match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
1442
1443         for( i=1; i<lgth1+1; i++ )
1444         {
1445                 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] );
1446         }
1447         for( j=1; j<lgth2+1; j++ )
1448         {
1449                 currentw[j] += ( ogcp2[0] + fgcp2[j-1] );
1450         }
1451
1452 #if STOREWM
1453         WMMTX[0][0] = initverticalw[0];
1454         for( i=1; i<lgth1+1; i++ )
1455         {
1456                 WMMTX[i][0] = initverticalw[i];
1457         }
1458         for( j=1; j<lgth2+1; j++ )
1459         {
1460                 WMMTX[0][j] = currentw[j];
1461         }
1462 #endif
1463
1464
1465         for( j=1; j<lgth2+1; ++j )
1466         {
1467                 m[j] = currentw[j-1] + ogcp1[1];
1468 //              m[j] = currentw[j-1];
1469                 mp[j] = 0;
1470         }
1471
1472         lastverticalw[0] = currentw[lgth2-1];
1473
1474         imid = lgth1 * 0.5;
1475
1476         jumpi = 0; // atode kawaru.
1477         lasti = lgth1+1;
1478 #if STOREWM
1479         for( i=1; i<lasti; i++ )
1480 #else
1481         for( i=1; i<=imid; i++ )
1482 #endif
1483         {
1484                 wtmp = previousw;
1485                 previousw = currentw;
1486                 currentw = wtmp;
1487
1488                 previousw[0] = initverticalw[i-1];
1489
1490                 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
1491                 currentw[0] = initverticalw[i];
1492
1493                 m[0] = ogcp1[i];
1494 #if STOREM
1495                 WMMTX2[i][0] = m[0];
1496 #endif
1497                 if( i == imid ) midm[0] = m[0];
1498
1499                 mi = previousw[0] + ogcp2[1];
1500 //              mi = previousw[0];
1501                 mpi = 0;
1502
1503
1504 //              ijppt = ijp[i] + 1;
1505                 mjpt = m + 1;
1506                 prept = previousw;
1507                 curpt = currentw + 1;
1508                 mpjpt = mp + 1;
1509
1510
1511                 lastj = lgth2+1;
1512                 for( j=1; j<lastj; j++ )
1513                 {
1514
1515                         wm = *prept;
1516
1517 #if 0
1518                         fprintf( stderr, "%5.0f->", wm );
1519 #endif
1520                         g = mi + fgcp2[j-1];
1521 //                      g = mi + fpenalty;
1522 #if 0
1523                         fprintf( stderr, "%5.0f?", g );
1524 #endif
1525                         if( g > wm )
1526                         {
1527                                 wm = g;
1528 //                              *ijppt = -( j - mpi );
1529                         }
1530                         g = *prept + ogcp2[j];
1531 //                      g = *prept;
1532                         if( g >= mi )
1533                         {
1534                                 mi = g;
1535                                 mpi = j-1;
1536                         }
1537 #if USE_PENALTY_EX
1538                         mi += fpenalty_ex;
1539 #endif
1540
1541                         g = *mjpt + fgcp1[i-1];
1542 //                      g = *mjpt + fpenalty;
1543 #if 0
1544                         fprintf( stderr, "%5.0f?", g );
1545 #endif
1546                         if( g > wm )
1547                         {
1548                                 wm = g;
1549 //                              *ijppt = +( i - *mpjpt );
1550                         }
1551
1552
1553                         g = *prept + ogcp1[i];
1554 //                      g = *prept;
1555                         if( g >= *mjpt )
1556                         {
1557                                 *mjpt = g;
1558                                 *mpjpt = i-1;
1559                         }
1560 #if USE_PENALTY_EX
1561                         m[j] += fpenalty_ex;
1562 #endif
1563 #if LOCAL
1564             if( wm < localthr )
1565             {
1566 //                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
1567                                 wm = 0;
1568             }
1569 #endif
1570
1571 #if 0
1572                         fprintf( stderr, "%5.0f ", wm );
1573 #endif
1574                         *curpt += wm;
1575
1576
1577 #if STOREWM
1578                         WMMTX[i][j] = *curpt;
1579                         WMMTX2[i][j] = *mjpt;
1580 #endif
1581
1582                         if( i == imid ) //muda
1583                         {
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 );
1588                                 midw[j] = *curpt;
1589                                 midm[j] = *mjpt;
1590                                 midn[j] = mi;
1591                         }
1592
1593 //                      fprintf( stderr, "m[%d] = %f\n", j, m[j] );
1594                         mjpt++;
1595                         prept++;
1596                         mpjpt++;
1597                         curpt++;
1598
1599                 }
1600                 lastverticalw[i] = currentw[lgth2-1];
1601
1602 #if STOREWM
1603                 WMMTX2[i][lgth2] = m[lgth2-1];
1604 #endif
1605
1606 #if 0  // ue
1607                 if( i == imid )
1608                 {
1609                         for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
1610                         for( j=0; j<lgth2; j++ ) midm[j] = m[j];
1611                 }
1612 #endif
1613         }
1614 //      for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
1615 //      for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
1616
1617 #if 0
1618     for( i=0; i<lgth1; i++ )
1619     {
1620         for( j=0; j<lgth2; j++ )
1621         {
1622             fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
1623         }
1624         fprintf( stderr, "\n" );
1625     }
1626         fprintf( stderr, "\n" );
1627         fprintf( stderr, "WMMTX2 = \n" );
1628     for( i=0; i<lgth1; i++ )
1629     {
1630         for( j=0; j<lgth2; j++ )
1631         {
1632             fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
1633         }
1634         fprintf( stderr, "\n" );
1635     }
1636         fprintf( stderr, "\n" );
1637 #endif
1638
1639 // gyakudp
1640
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 );
1643
1644         for( i=0; i<lgth1-1; i++ )
1645         {
1646                 initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
1647 //              initverticalw[i] += fpenalty;
1648         }
1649         for( j=0; j<lgth2-1; j++ )
1650         {
1651                 currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
1652 //              currentw[j] += fpenalty;
1653         }
1654
1655 #if STOREWM
1656         for( i=0; i<lgth1-1; i++ )
1657         {
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] );
1660         }
1661         for( j=0; j<lgth2-1; j++ )
1662         {
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] );
1665         }
1666 #endif
1667
1668
1669
1670
1671
1672
1673         for( j=lgth2-1; j>0; --j )
1674         {
1675                 m[j-1] = currentw[j] + fgcp2[lgth2-2];
1676 //              m[j-1] = currentw[j];
1677                 mp[j] = lgth1-1;
1678         }
1679
1680 //      for( j=0; j<lgth2; j++ ) m[j] = 0.0;
1681         // m[lgth2-1] ha irunoka?
1682
1683
1684 //      for( i=lgth1-2; i>=imid; i-- )
1685         firstm = -9999999.9;
1686         firstmp = lgth1-1;
1687         for( i=lgth1-2; i>-1; i-- )
1688         {
1689                 wtmp = previousw;
1690                 previousw = currentw;
1691                 currentw = wtmp;
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 );
1695
1696                 currentw[lgth2-1] = initverticalw[i];
1697
1698 //              m[lgth2] = fgcp1[i];
1699 //              WMMTX2[i][lgth2] += m[lgth2];
1700 //              fprintf( stderr, "m[] = %f\n", m[lgth2] );
1701
1702                 mi = previousw[lgth2-1] + fgcp2[lgth2-2];
1703 //              mi = previousw[lgth2-1];
1704                 mpi = lgth2 - 1;
1705
1706                 mjpt = m + lgth2 - 2;
1707                 prept = previousw + lgth2 - 1;
1708                 curpt = currentw + lgth2 - 2;
1709                 mpjpt = mp + lgth2 - 2;
1710
1711
1712                 for( j=lgth2-2; j>-1; j-- )
1713                 {
1714                         wm = *prept;
1715                         ijpi = i+1;
1716                         ijpj = j+1;
1717
1718                         g = mi + ogcp2[j+1];
1719 //                      g = mi + fpenalty;
1720                         if( g > wm )
1721                         {
1722                                 wm = g;
1723                                 ijpj = mpi;
1724                                 ijpi = i+1;
1725                         }
1726
1727                         g = *prept + fgcp2[j];
1728 //                      g = *prept;
1729                         if( g >= mi )
1730                         {
1731 //                              fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
1732                                 mi = g;
1733                                 mpi = j + 1;
1734                         }
1735
1736 #if USE_PENALTY_EX
1737                         mi += fpenalty_ex;
1738 #endif
1739
1740 //                      fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
1741                         g = *mjpt + ogcp1[i+1];
1742 //                      g = *mjpt + fpenalty;
1743                         if( g > wm )
1744                         {
1745                                 wm = g;
1746                                 ijpi = *mpjpt;
1747                                 ijpj = j+1;
1748                         }
1749
1750 //                      if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
1751                         g = *prept + fgcp1[i];
1752 //                      g = *prept;
1753                         if( g >= *mjpt )
1754                         {
1755                                 *mjpt = g;
1756                                 *mpjpt = i + 1;
1757                         }
1758
1759 #if USE_PENALTY_EX
1760                         m[j] += fpenalty_ex;
1761 #endif
1762
1763                         if( i == jumpi || i == imid - 1 )
1764                         {
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 );
1769                         }
1770                         if( i == imid ) // muda
1771                         {
1772                                 midw[j] += wm;
1773 //                              midm[j+1] += *mjpt + fpenalty; //??????
1774                                 midm[j+1] += *mjpt; //??????
1775                         }
1776                         if( i == imid - 1 )
1777                         {
1778 //                              midn[j] += mi + fpenalty;  //????
1779                                 midn[j] += mi;  //????
1780                         }
1781 #if LOCAL
1782             if( wm < localthr )
1783             {
1784 //                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
1785                                 wm = 0;
1786             }
1787 #endif
1788
1789 #if STOREWM
1790                         WMMTX[i][j] += wm;
1791 //                      WMMTX2[i][j+1] += *mjpt + fpenalty;
1792                         WMMTX2[i][j] += *curpt;
1793 #endif
1794                         *curpt += wm;
1795
1796                         mjpt--;
1797                         prept--;
1798                         mpjpt--;
1799                         curpt--;
1800                 }
1801 //              fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
1802                 g = *prept + fgcp1[i];
1803                 if( firstm < g )
1804                 {
1805                         firstm = g;
1806                         firstmp = i + 1;
1807                 }
1808 #if STOREWM
1809 //              WMMTX2[i][j+1] += firstm;
1810 #endif
1811                 if( i == imid ) midm[j+1] += firstm;
1812
1813                 if( i == imid - 1 )
1814                 {
1815                         maxwm = midw[1];
1816                         jmid = 0;
1817 //                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1818                         for( j=2; j<lgth2-1; j++ )
1819                         {
1820                                 wm = midw[j];
1821                                 if( wm > maxwm )
1822                                 {
1823                                         jmid = j;
1824                                         maxwm = wm;
1825                                 }
1826 //                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1827                         }
1828                         for( j=0; j<lgth2+1; j++ )
1829                         {
1830                                 wm = midm[j];
1831                                 if( wm > maxwm )
1832                                 {
1833                                         jmid = j;
1834                                         maxwm = wm;
1835                                 }
1836 //                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1837                         }
1838
1839 //                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1840
1841
1842 //                      fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
1843                         wm = midw[jmid];
1844                         jumpi = imid-1;
1845                         jumpj = jmid-1;
1846                         if( jmid > 0 && midn[jmid-1] > wm ) //060413
1847                         {
1848                                 jumpi = imid-1;
1849                                 jumpj = jumpbacki[jmid];
1850                                 wm = midn[jmid-1];
1851 //                              fprintf( stderr, "rejump (n)\n" );
1852                         }
1853                         if( midm[jmid] > wm )
1854                         {
1855                                 jumpi = jumpbackj[jmid];
1856                                 jumpj = jmid-1;
1857                                 wm = midm[jmid];
1858 //                              fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
1859                         }
1860
1861
1862 //                      fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
1863 //                      fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
1864 #if 0
1865                         fprintf( stderr, "imid = %d\n", imid );
1866                         fprintf( stderr, "midn = \n" );
1867                         for( j=0; j<lgth2; j++ )
1868                         {
1869                                 fprintf( stderr, "% 7.1f ", midn[j] );
1870                         }
1871                         fprintf( stderr, "\n" );
1872                         fprintf( stderr, "midw = \n" );
1873                         for( j=0; j<lgth2; j++ )
1874                         {
1875                                 fprintf( stderr, "% 7.1f ", midw[j] );
1876                         }
1877                         fprintf( stderr, "\n" );
1878                         fprintf( stderr, "midm = \n" );
1879                         for( j=0; j<lgth2; j++ )
1880                         {
1881                                 fprintf( stderr, "% 7.1f ", midm[j] );
1882                         }
1883                         fprintf( stderr, "\n" );
1884 #endif
1885 //                      fprintf( stderr, "maxwm = %f\n", maxwm );
1886                 }
1887                 if( i == jumpi ) //saki?
1888                 {
1889 //                      fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
1890 //                      fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
1891                         if( jmid == 0 )
1892                         {
1893 //                              fprintf( stderr, "CHUI2!\n" );
1894                                 jumpj = 0; jmid = 1;
1895                                 jumpi = firstmp - 1;
1896                                 imid = firstmp;
1897                         }
1898
1899 #if 0
1900                         else if( jmid == lgth2 )
1901                         {
1902                                 fprintf( stderr, "CHUI1!\n" );
1903                                 jumpi=0; jumpj=0;
1904                                 imid=jumpforwi[0]; jmid=lgth2-1;
1905                         }
1906 #else // 060414
1907                         else if( jmid >= lgth2 )
1908                         {
1909 //                              fprintf( stderr, "CHUI1!\n" );
1910                                 jumpi=imid-1; jmid=lgth2;
1911                                 jumpj = lgth2-1;
1912                         }
1913 #endif
1914                         else
1915                         {
1916                                 imid = jumpforwi[jumpj];
1917                                 jmid = jumpforwj[jumpj];
1918                         }
1919 #if 0
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 );
1924 #endif
1925
1926 #if STOREWM
1927 //                      break;
1928 #else
1929                         break;
1930 #endif
1931                 }
1932         }
1933 #if 0
1934                 jumpi=0; jumpj=0;
1935                 imid=lgth1-1; jmid=lgth2-1;
1936         }
1937 #endif
1938
1939 //      fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1940 //      fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1941
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];
1946
1947
1948
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;
1952
1953 #if STOREWM
1954
1955 #if 0
1956         for( i=0; i<lgth1; i++ )
1957         {
1958                 float maxpairscore = -9999.9;
1959                 float tmpscore;
1960
1961                 for( j=0; j<lgth2; j++ )
1962                 {
1963                         if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1964                         {
1965                                 map12[i].pos = j;
1966                                 map12[i].score = tmpscore;
1967                                 maxpairscore = tmpscore;
1968                         }
1969                 }
1970
1971                 for( k=0; k<lgth1; k++ )
1972                 {
1973                         if( i == k ) continue;
1974                         if( map12[i].score <= WMMTX[k][map12[i].pos] )
1975                                 break;
1976                 }
1977                 if( k != lgth1 )
1978                 {
1979                         map12[i].pos = -1;
1980                         map12[i].score = -1.0;
1981                 }
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] );
1983         }
1984         for( j=0; j<lgth2; j++ )
1985         {
1986                 float maxpairscore = -9999.9;
1987                 float tmpscore;
1988
1989                 for( i=0; i<lgth1; i++ )
1990                 {
1991                         if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1992                         {
1993                                 map21[j].pos = i;
1994                                 map21[j].score = tmpscore;
1995                                 maxpairscore = tmpscore;
1996                         }
1997                 }
1998
1999                 for( k=0; k<lgth2; k++ )
2000                 {
2001                         if( j == k ) continue;
2002                         if( map21[j].score <= WMMTX[map21[j].pos][k] )
2003                                 break;
2004                 }
2005                 if( k != lgth2 )
2006                 {
2007                         map21[j].pos = -1;
2008                         map21[j].score = -1.0;
2009                 }
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] );
2011         }
2012
2013         for( i=0; i<lgth1; i++ )
2014         {
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 );
2017         }
2018 #endif
2019
2020 #if 0
2021         fprintf( stderr, "WMMTX = \n" );
2022     for( i=0; i<lgth1; i++ )
2023     {
2024         fprintf( stderr, "%d ", i );
2025         for( j=0; j<lgth2; j++ )
2026         {
2027             fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
2028         }
2029         fprintf( stderr, "\n" );
2030     }
2031 //      fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
2032     for( i=0; i<lgth1; i++ )
2033     {
2034         fprintf( stderr, "%d ", i );
2035         for( j=0; j<lgth2+1; j++ )
2036         {
2037             fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
2038         }
2039         fprintf( stderr, "\n" );
2040     }
2041         exit( 1 );
2042 #endif
2043
2044 #if 0
2045     fprintf( stdout, "#WMMTX = \n" );
2046     for( i=0; i<lgth1; i++ )
2047     {
2048 //        fprintf( stdout, "%d ", i );
2049         for( j=0; j<lgth2; j++ )
2050         {
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" );
2055                         else
2056                 fprintf( stdout, "\n" );
2057         }
2058         fprintf( stdout, "\n" );
2059     }
2060         exit( 1 );
2061 #endif
2062
2063 #if 0
2064
2065         fprintf( stderr, "jumpbacki = \n" );
2066         for( j=0; j<lgth2; j++ )
2067         {
2068                 fprintf( stderr, "% 7d ", jumpbacki[j] );
2069         }
2070         fprintf( stderr, "\n" );
2071         fprintf( stderr, "jumpbackj = \n" );
2072         for( j=0; j<lgth2; j++ )
2073         {
2074                 fprintf( stderr, "% 7d ", jumpbackj[j] );
2075         }
2076         fprintf( stderr, "\n" );
2077         fprintf( stderr, "jumpforwi = \n" );
2078         for( j=0; j<lgth2; j++ )
2079         {
2080                 fprintf( stderr, "% 7d ", jumpforwi[j] );
2081         }
2082         fprintf( stderr, "\n" );
2083         fprintf( stderr, "jumpforwj = \n" );
2084         for( j=0; j<lgth2; j++ )
2085         {
2086                 fprintf( stderr, "% 7d ", jumpforwj[j] );
2087         }
2088         fprintf( stderr, "\n" );
2089 #endif
2090
2091
2092 #endif
2093
2094
2095 //      Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
2096
2097 #if 0 // irukamo
2098         resultlen = strlen( mseq1[0] );
2099         if( alloclen < resultlen || resultlen > N )
2100         {
2101                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
2102                 ErrorExit( "LENGTH OVER!\n" );
2103         }
2104 #endif
2105
2106
2107
2108 #if 0
2109         fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
2110         fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
2111
2112         fprintf( stderr, "imid = %d\n", imid );
2113         fprintf( stderr, "jmid = %d\n", jmid );
2114 #endif
2115
2116
2117         FreeFloatVec( w1 );
2118         FreeFloatVec( w2 );
2119         FreeFloatVec( initverticalw );
2120         FreeFloatVec( lastverticalw );
2121         FreeFloatVec( midw );
2122         FreeFloatVec( midm );
2123         FreeFloatVec( midn );
2124
2125         FreeIntVec( jumpbacki );
2126         FreeIntVec( jumpbackj );
2127         FreeIntVec( jumpforwi );
2128         FreeIntVec( jumpforwj );
2129         FreeIntVec( jumpdummi );
2130         FreeIntVec( jumpdummj );
2131
2132         FreeFloatVec( m );
2133         FreeIntVec( mp );
2134
2135         FreeFloatMtx( floatwork );
2136         FreeIntMtx( intwork );
2137
2138 #if STOREWM
2139         FreeFloatMtx( WMMTX );
2140         FreeFloatMtx( WMMTX2 );
2141 #endif
2142
2143         return( value );
2144
2145 //      fprintf( stderr, "==== calling myself (first)\n" );
2146
2147 #if 0
2148                 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2149                 fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
2150 #endif
2151         value = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ist+jumpi, jst, jst+jumpj, alloclen, aseq1, aseq2, depth, gapinfo, map );
2152 #if 0
2153                 fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
2154                 fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
2155 #endif
2156 #if MEMSAVE
2157 #else
2158         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
2159         for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
2160 #endif
2161
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 );
2164
2165         len = strlen( mseq1[0] );
2166 //      fprintf( stderr, "len = %d\n", len );
2167         l = jmid - jumpj - 1;
2168 //      fprintf( stderr, "l=%d\n", l );
2169         if( l > 0 )
2170         {
2171                 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
2172                 for( i=0; i<icyc; i++ )
2173                 {
2174                         strcat( mseq1[i], gaps );
2175                         mseq1[i][len+l] = 0;
2176                 }
2177                 for( j=0; j<jcyc; j++ )
2178                 {
2179                         strncat( mseq2[j], seq2[j]+jst+jumpj+1, l );
2180                         mseq2[j][len+l] = 0;
2181                 }
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;
2185         }
2186         len = strlen( mseq1[0] );
2187         l = imid - jumpi - 1;
2188 //      fprintf( stderr, "l=%d\n", l );
2189         if( l > 0 )
2190         {
2191                 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
2192                 for( i=0; i<icyc; i++ )
2193                 {
2194                         strncat( mseq1[i], seq1[i]+ist+jumpi+1, l );
2195                         mseq1[i][len+l] = 0;
2196                 }
2197                 for( j=0; j<jcyc; j++ )
2198                 {
2199                         strcat( mseq2[j], gaps );
2200                         mseq2[j][len+l] = 0;
2201                 }
2202
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] );
2205
2206
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;
2210         }
2211 #if 0
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] );
2214 #endif
2215
2216 //      fprintf( stderr, "==== calling myself (second)\n" );
2217
2218 #if MEMSAVE
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;
2222 #endif
2223
2224 #if 0
2225                 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2226                 fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
2227 #endif
2228         value += MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist+imid, ien, jst+jmid, jen, alloclen, aseq1, aseq2, depth, gapinfo, map );
2229 #if 0
2230                 fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
2231                 fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
2232 #endif
2233
2234
2235
2236 #if DEBUG
2237         if( value - maxwm > 1 || maxwm - value > 1 )
2238         {
2239                 fprintf( stderr, "WARNING value  = %f, but maxwm = %f\n", value, maxwm );
2240                 for( i=0; i<icyc; i++ )
2241                 {
2242                         fprintf( stderr, ">1-%d\n%s\n", i, mseq1[i] );
2243                         fprintf( stderr, "%s\n", aseq1[i] );
2244                 }
2245                 for( i=0; i<jcyc; i++ )
2246                 {
2247                         fprintf( stderr, ">2-%d\n%s\n", i, mseq2[i] );
2248                         fprintf( stderr, "%s\n", aseq2[i] );
2249                 }
2250
2251 //              exit( 1 );
2252         }
2253         else
2254         {
2255                 fprintf( stderr, "value = %.0f, maxwm = %.0f -> ok\n", value, maxwm );
2256         }
2257 #endif
2258
2259 #if MEMSAVE
2260 #else
2261         for( i=0; i<icyc; i++ ) strcat( mseq1[i], aseq1[i] );
2262         for( i=0; i<jcyc; i++ ) strcat( mseq2[i], aseq2[i] );
2263 #endif
2264
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 );
2267
2268         free( gaps );
2269 #if MEMSAVE
2270         free( aseq1 );
2271         free( aseq2 );
2272 #else
2273         FreeCharMtx( aseq1 );
2274         FreeCharMtx( aseq2 );
2275 #endif
2276
2277         return( value );
2278 }
2279
2280
2281
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 */
2284 {
2285 //      int k;
2286         int i, j;
2287         int ll1, ll2;
2288         int lgth1, lgth2;
2289         float wm = 0.0;   /* int ?????? */
2290         char **mseq1;
2291         char **mseq2;
2292 //      char **mseq;
2293         float *ogcp1;
2294         float *ogcp2;
2295         float *fgcp1;
2296         float *fgcp2;
2297         float **cpmx1;
2298         float **cpmx2;
2299         float **gapinfo;
2300 //      float fpenalty;
2301         float fpenalty = (float)RNApenalty;
2302         int nglen1, nglen2;
2303
2304
2305
2306
2307
2308 #if 0
2309         fprintf( stderr, "eff in SA+++align\n" );
2310         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
2311 #endif
2312
2313         nglen1 = seqlen( seq1[0] );
2314         nglen2 = seqlen( seq2[0] );
2315
2316         lgth1 = strlen( seq1[0] );
2317         lgth2 = strlen( seq2[0] );
2318
2319         ll1 = ( (int)(lgth1) ) + 100;
2320         ll2 = ( (int)(lgth2) ) + 100;
2321
2322         mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
2323         mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
2324
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 );
2330
2331
2332         cpmx1 = AllocateFloatMtx( ll1+2, 27 );
2333         cpmx2 = AllocateFloatMtx( ll2+2, 27 );
2334
2335         for( i=0; i<icyc; i++ )
2336         {
2337                 if( strlen( seq1[i] ) != lgth1 )
2338                 {
2339                         fprintf( stderr, "i = %d / %d\n", i, icyc );
2340                         fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2341                         exit( 1 );
2342                 }
2343         }
2344         for( j=0; j<jcyc; j++ )
2345         {
2346                 if( strlen( seq2[j] ) != lgth2 )
2347                 {
2348                         fprintf( stderr, "j = %d / %d\n", j, jcyc );
2349                         fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2350                         exit( 1 );
2351                 }
2352         }
2353
2354         MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
2355         MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
2356
2357
2358 #if 1
2359
2360         if( sgap1 )
2361         {
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 );
2366         }
2367         else
2368         {
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 );
2373         }
2374
2375 #if 1
2376         for( i=0; i<lgth1; i++ )
2377         {
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] );
2381         }
2382         for( i=0; i<lgth2; i++ )
2383         {
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] );
2387         }
2388 #else
2389         for( i=0; i<lgth1; i++ )
2390         {
2391                 ogcp1[i] = 0.5 * fpenalty;
2392                 fgcp1[i] = 0.5 * fpenalty;
2393         }
2394         for( i=0; i<lgth2; i++ )
2395         {
2396                 ogcp2[i] = 0.5 * fpenalty;
2397                 fgcp2[i] = 0.5 * fpenalty;
2398         }
2399 #endif
2400
2401         gapinfo[0] = ogcp1;
2402         gapinfo[1] = fgcp1;
2403         gapinfo[2] = ogcp2;
2404         gapinfo[3] = fgcp2;
2405 #endif
2406
2407 #if 0
2408         fprintf( stdout, "in MSalignmm.c\n" );
2409         for( i=0; i<icyc; i++ )
2410         {
2411                 fprintf( stdout, ">%d of GROUP1\n", i );
2412                 fprintf( stdout, "%s\n", seq1[i] );
2413         }
2414         for( i=0; i<jcyc; i++ )
2415         {
2416                 fprintf( stdout, ">%d of GROUP2\n", i );
2417                 fprintf( stdout, "%s\n", seq2[i] );
2418         }
2419         fflush( stdout );
2420 #endif
2421
2422         wm = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo, map );
2423 #if DEBUG
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] );
2428 #endif
2429
2430 //      fprintf( stderr, "wm = %f\n", wm );
2431
2432 #if 0
2433
2434         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2435         for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
2436
2437         if( seqlen( seq1[0] ) != nglen1 )
2438         {
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] );
2441                 exit( 1 );
2442         }
2443         if( seqlen( seq2[0] ) != nglen2 )
2444         {
2445                 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
2446                 exit( 1 );
2447         }
2448 #endif
2449
2450         FreeFloatVec( ogcp1 );
2451         FreeFloatVec( ogcp2 );
2452         FreeFloatVec( fgcp1 );
2453         FreeFloatVec( fgcp2 );
2454         FreeFloatMtx( cpmx1 );
2455         FreeFloatMtx( cpmx2 );
2456         free( (void *)gapinfo );
2457
2458         FreeCharMtx( mseq1 );
2459         FreeCharMtx( mseq2 );
2460
2461         lgth1 = strlen( seq1[0] );
2462         lgth2 = strlen( seq2[0] );
2463         for( i=0; i<icyc; i++ )
2464         {
2465                 if( strlen( seq1[i] ) != lgth1 )
2466                 {
2467                         fprintf( stderr, "i = %d / %d\n", i, icyc );
2468                         fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2469                         exit( 1 );
2470                 }
2471         }
2472         for( j=0; j<jcyc; j++ )
2473         {
2474                 if( strlen( seq2[j] ) != lgth2 )
2475                 {
2476                         fprintf( stderr, "j = %d / %d\n", j, jcyc );
2477                         fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2478                         exit( 1 );
2479                 }
2480         }
2481
2482         return( wm );
2483 }
2484
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 */
2487 {
2488 //      int k;
2489         int i, j;
2490         int ll1, ll2;
2491         int lgth1, lgth2;
2492         float wm = 0.0;   /* int ?????? */
2493         char **mseq1;
2494         char **mseq2;
2495         float *ogcp1;
2496         float *ogcp2;
2497         float *fgcp1;
2498         float *fgcp2;
2499         float **cpmx1;
2500         float **cpmx2;
2501         float **gapinfo;
2502         float fpenalty = (float)penalty;
2503         int nglen1, nglen2;
2504
2505 #if 0
2506         fprintf( stderr, "eff in SA+++align\n" );
2507         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
2508 #endif
2509
2510         nglen1 = seqlen( seq1[0] );
2511         nglen2 = seqlen( seq2[0] );
2512
2513         lgth1 = strlen( seq1[0] );
2514         lgth2 = strlen( seq2[0] );
2515
2516         ll1 = ( (int)(lgth1) ) + 100;
2517         ll2 = ( (int)(lgth2) ) + 100;
2518
2519         mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
2520         mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
2521
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 );
2527
2528
2529         cpmx1 = AllocateFloatMtx( ll1+2, 39 );
2530         cpmx2 = AllocateFloatMtx( ll2+2, 39 );
2531
2532         for( i=0; i<icyc; i++ )
2533         {
2534                 if( strlen( seq1[i] ) != lgth1 )
2535                 {
2536                         fprintf( stderr, "i = %d / %d\n", i, icyc );
2537                         fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2538                         exit( 1 );
2539                 }
2540         }
2541         for( j=0; j<jcyc; j++ )
2542         {
2543                 if( strlen( seq2[j] ) != lgth2 )
2544                 {
2545                         fprintf( stderr, "j = %d / %d\n", j, jcyc );
2546                         fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2547                         exit( 1 );
2548                 }
2549         }
2550
2551 #if 0
2552         MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
2553         MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
2554 #else
2555         cpmx_ribosum( seq1, seq1r, dir1, cpmx1, eff1, lgth1, icyc );
2556         cpmx_ribosum( seq2, seq2r, dir2, cpmx2, eff2, lgth2, jcyc );
2557 #endif
2558
2559
2560 #if 1
2561
2562         if( sgap1 )
2563         {
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 );
2568         }
2569         else
2570         {
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 );
2575         }
2576
2577 #if 1
2578         for( i=0; i<lgth1; i++ )
2579         {
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] );
2583         }
2584         for( i=0; i<lgth2; i++ )
2585         {
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] );
2589         }
2590 #else
2591         for( i=0; i<lgth1; i++ )
2592         {
2593                 ogcp1[i] = 0.5 * fpenalty;
2594                 fgcp1[i] = 0.5 * fpenalty;
2595         }
2596         for( i=0; i<lgth2; i++ )
2597         {
2598                 ogcp2[i] = 0.5 * fpenalty;
2599                 fgcp2[i] = 0.5 * fpenalty;
2600         }
2601 #endif
2602
2603         gapinfo[0] = ogcp1;
2604         gapinfo[1] = fgcp1;
2605         gapinfo[2] = ogcp2;
2606         gapinfo[3] = fgcp2;
2607 #endif
2608
2609 #if 0
2610         fprintf( stdout, "in MSalignmm.c\n" );
2611         for( i=0; i<icyc; i++ )
2612         {
2613                 fprintf( stdout, ">%d of GROUP1\n", i );
2614                 fprintf( stdout, "%s\n", seq1[i] );
2615         }
2616         for( i=0; i<jcyc; i++ )
2617         {
2618                 fprintf( stdout, ">%d of GROUP2\n", i );
2619                 fprintf( stdout, "%s\n", seq2[i] );
2620         }
2621         fflush( stdout );
2622 #endif
2623
2624         wm = MSalign2m2m_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo, map );
2625 #if DEBUG
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] );
2630 #endif
2631
2632 //      fprintf( stderr, "wm = %f\n", wm );
2633
2634 #if 0
2635
2636         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2637         for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
2638
2639         if( seqlen( seq1[0] ) != nglen1 )
2640         {
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] );
2643                 exit( 1 );
2644         }
2645         if( seqlen( seq2[0] ) != nglen2 )
2646         {
2647                 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
2648                 exit( 1 );
2649         }
2650 #endif
2651
2652         FreeFloatVec( ogcp1 );
2653         FreeFloatVec( ogcp2 );
2654         FreeFloatVec( fgcp1 );
2655         FreeFloatVec( fgcp2 );
2656         FreeFloatMtx( cpmx1 );
2657         FreeFloatMtx( cpmx2 );
2658         free( (void *)gapinfo );
2659
2660         FreeCharMtx( mseq1 );
2661         FreeCharMtx( mseq2 );
2662
2663         lgth1 = strlen( seq1[0] );
2664         lgth2 = strlen( seq2[0] );
2665         for( i=0; i<icyc; i++ )
2666         {
2667                 if( strlen( seq1[i] ) != lgth1 )
2668                 {
2669                         fprintf( stderr, "i = %d / %d\n", i, icyc );
2670                         fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2671                         exit( 1 );
2672                 }
2673         }
2674         for( j=0; j<jcyc; j++ )
2675         {
2676                 if( strlen( seq2[j] ) != lgth2 )
2677                 {
2678                         fprintf( stderr, "j = %d / %d\n", j, jcyc );
2679                         fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2680                         exit( 1 );
2681                 }
2682         }
2683
2684         return( wm );
2685 }