binaries/src/mafft/core/MSalignmm.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 0
   9
  10 #define DPTANNI 100
  11
  12
  13 static TLS int reccycle = 0;
  14
  15
  16 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
  17 {
  18         int j, k, l;
  19         float scarr[26];
  20         float **cpmxpd = floatwork;
  21         int **cpmxpdn = intwork;
  22         int count = 0;
  23         float *matchpt;
  24         float **cpmxpdpt;
  25         int **cpmxpdnpt;
  26         int cpkd;
  27
  28         if( initialize )
  29         {
  30                 for( j=0; j<lgth2; j++ )
  31                 {
  32                         count = 0;
  33                         for( l=0; l<26; l++ )
  34                         {
  35                                 if( cpmx2[j][l] )
  36                                 {
  37                                         cpmxpd[j][count] = cpmx2[j][l];
  38                                         cpmxpdn[j][count] = l;
  39                                         count++;
  40                                 }
  41                         }
  42                         cpmxpdn[j][count] = -1;
  43                 }
  44         }
  45
  46         for( l=0; l<26; l++ )
  47         {
  48                 scarr[l] = 0.0;
  49                 for( k=0; k<26; k++ )
  50                 {
  51                         scarr[l] += n_dis[k][l] * cpmx1[i1][k];
  52                 }
  53         }
  54 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
  55         {
  56                 float *fpt, **fptpt, *fpt2;
  57                 int *ipt, **iptpt;
  58                 fpt2 = match;
  59                 iptpt = cpmxpdn;
  60                 fptpt = cpmxpd;
  61                 while( lgth2-- )
  62                 {
  63                         *fpt2 = 0.0;
  64                         ipt=*iptpt,fpt=*fptpt;
  65                         while( *ipt > -1 )
  66                                 *fpt2 += scarr[*ipt++] * *fpt++;
  67                         fpt2++,iptpt++,fptpt++;
  68                 }
  69         }
  70         for( j=0; j<lgth2; j++ )
  71         {
  72                 match[j] = 0.0;
  73                 for( k=0; cpmxpdn[j][k]>-1; k++ )
  74                         match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
  75         }
  76 #else
  77         matchpt = match;
  78         cpmxpdnpt = cpmxpdn;
  79         cpmxpdpt = cpmxpd;
  80         while( lgth2-- )
  81         {
  82                 *matchpt = 0.0;
  83                 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
  84                         *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
  85                 matchpt++;
  86                 cpmxpdnpt++;
  87                 cpmxpdpt++;
  88         }
  89 #endif
  90 }
  91
  92 static float Atracking( float *lasthorizontalw, float *lastverticalw,
  93                                                 char **seq1, char **seq2,
  94                         char **mseq1, char **mseq2,
  95                         int **ijp, int icyc, int jcyc,
  96                                                 int ist, int ien, int jst, int jen,
  97                                                 int fulllen1, int fulllen2, int tailgp )
  98 {
  99         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, klim;
 100         char *gaptable1, *gt1bk;
 101         char *gaptable2, *gt2bk;
 102         float wm;
 103         lgth1 = ien-ist+1;
 104         lgth2 = jen-jst+1;
 105
 106         gt1bk = AllocateCharVec( lgth1+lgth2+1 );
 107         gt2bk = AllocateCharVec( lgth1+lgth2+1 );
 108
 109 #if 0
 110         for( i=0; i<lgth1; i++ )
 111         {
 112                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 113         }
 114 #endif
 115
 116
 117 //      fprintf( stderr, "in Atracking, tailgp=%d, ien=%d, fulllen1=%d, jen=%d, fulllen2=%d\n", tailgp, ien, fulllen1, jen, fulllen2 );
 118
 119         if( tailgp == 1 )
 120                 ;
 121         else if( ien == fulllen1-1 || jen == fulllen2-1 )
 122         {
 123 //              fprintf( stderr, "searching lasthorizontalw and lasthorizontalw\n" );
 124                 wm = lastverticalw[0];
 125                 for( i=0; i<lgth1; i++ )
 126                 {
 127                         if( lastverticalw[i] >= wm )
 128                         {
 129                                 wm = lastverticalw[i];
 130                                 iin = i; jin = lgth2-1;
 131                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 132                         }
 133                 }
 134                 for( j=0; j<lgth2; j++ )
 135                 {
 136                         if( lasthorizontalw[j] >= wm )
 137                         {
 138                                 wm = lasthorizontalw[j];
 139                                 iin = lgth1-1; jin = j;
 140                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 141                         }
 142                 }
 143         }
 144 #if 0
 145         else if( jen == fulllen2-1 )
 146         {
 147                 fprintf( stderr, "searching lastverticalw\n" );
 148                 wm = lastverticalw[0];
 149                 for( i=0; i<lgth1; i++ )
 150                 {
 151                         if( lastverticalw[i] >= wm )
 152                         {
 153                                 wm = lastverticalw[i];
 154                                 iin = i; jin = lgth2-1;
 155                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 156                         }
 157                 }
 158         }
 159         else if( ien == fulllen1-1 )
 160         {
 161                 fprintf( stderr, "searching lasthorizontalw\n" );
 162                 wm = lasthorizontalw[0];
 163                 for( j=0; j<lgth2; j++ )
 164                 {
 165                         if( lasthorizontalw[j] >= wm )
 166                         {
 167                                 wm = lasthorizontalw[j];
 168                                 iin = lgth1-1; jin = j;
 169                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 170                         }
 171                 }
 172         }
 173 #endif
 174
 175     for( i=0; i<lgth1+1; i++ )
 176     {
 177         ijp[i][0] = i + 1;
 178     }
 179     for( j=0; j<lgth2+1; j++ )
 180     {
 181         ijp[0][j] = -( j + 1 );
 182     }
 183
 184
 185         gaptable1 = gt1bk + lgth1+lgth2;
 186         *gaptable1 = 0;
 187         gaptable2 = gt2bk + lgth1+lgth2;
 188         *gaptable2 = 0;
 189
 190 //      if( lgth2 == 1 ) fprintf( stderr, "in Atracking, mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
 191
 192         iin = lgth1; jin = lgth2;
 193         klim = lgth1+lgth2;
 194         for( k=0; k<=klim; k++ )
 195         {
 196                 if( ijp[iin][jin] < 0 )
 197                 {
 198                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 199                 }
 200                 else if( ijp[iin][jin] > 0 )
 201                 {
 202                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 203                 }
 204                 else
 205                 {
 206                         ifi = iin-1; jfi = jin-1;
 207                 }
 208                 l = iin - ifi;
 209                 while( --l )
 210                 {
 211                         *--gaptable1 = 'o';
 212                         *--gaptable2 = '-';
 213                         k++;
 214                 }
 215                 l= jin - jfi;
 216                 while( --l )
 217                 {
 218                         *--gaptable1 = '-';
 219                         *--gaptable2 = 'o';
 220                         k++;
 221                 }
 222                 if( iin <= 0 || jin <= 0 ) break;
 223                 *--gaptable1 = 'o';
 224                 *--gaptable2 = 'o';
 225                 k++;
 226                 iin = ifi; jin = jfi;
 227
 228         }
 229
 230         for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i]+ist, gaptable1 );
 231         for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j]+jst, gaptable2 );
 232
 233         free( gt1bk );
 234         free( gt2bk );
 235
 236 //      fprintf( stderr, "in Atracking (owari), mseq1 = %s\n", mseq1[0] );
 237 //      fprintf( stderr, "in Atracking (owari), mseq2 = %s\n", mseq2[0] );
 238         return( 0.0 );
 239 }
 240
 241 static float MSalignmm_tanni( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, float **cpmx1, float **cpmx2, int ist, int ien, int jst, int jen, int alloclen, int fulllen1, int fulllen2, char **mseq1, char **mseq2, float **gapinfo, int headgp, int tailgp )
 242 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 243 {
 244 //      int k;
 245         register int i, j;
 246         int ll1, ll2;
 247         int lasti, lastj;
 248         float wm = 0.0;   /* int ?????? */
 249         float g;
 250         float *currentw, *previousw;
 251 #if USE_PENALTY_EX
 252         float fpenalty_ex = (float)penalty_ex;
 253 #endif
 254 #if 1
 255         float *wtmp;
 256         int *ijppt;
 257         float *mjpt, *prept, *curpt;
 258         int *mpjpt;
 259 #endif
 260         float mi, *m;
 261         int **ijp;
 262         int mpi, *mp;
 263         float *w1, *w2;
 264         float *initverticalw;    /* kufuu sureba iranai */
 265         float *lastverticalw;    /* kufuu sureba iranai */
 266         int **intwork;
 267         float **floatwork;
 268         int **intmtx;
 269         int lgth1, lgth2;
 270         float *ogcp1;
 271         float *fgcp1;
 272         float *ogcp2;
 273         float *fgcp2;
 274 //      char **aseq1;
 275 //      char **aseq2;
 276 //      char **aseq1bk, **aseq2bk;
 277
 278
 279         ogcp1 = gapinfo[0] + ist;
 280         fgcp1 = gapinfo[1] + ist;
 281         ogcp2 = gapinfo[2] + jst;
 282         fgcp2 = gapinfo[3] + jst;
 283
 284 #if STOREWM
 285         char ttt1[10000], ttt2[10000];
 286 #endif
 287
 288
 289         lgth1 = ien-ist+1;
 290         lgth2 = jen-jst+1;
 291
 292 #if STOREWM
 293         strncpy( ttt1, seq1[0]+ist, lgth1 ); ttt1[lgth1] = 0;
 294         strncpy( ttt2, seq2[0]+jst, lgth2 ); ttt2[lgth2] = 0;
 295
 296         fprintf( stderr, "in _tanni ist,ien = %d,%d, lgth1=%d\n", ist, ien, lgth1 );
 297         fprintf( stderr, "in _tanni jst,jen = %d,%d, lgth2=%d\n", jst, jen, lgth2 );
 298         fprintf( stderr, "ttt1 = %s\n", ttt1 );
 299         fprintf( stderr, "ttt2 = %s\n", ttt2 );
 300 #endif
 301
 302 #if 0
 303         fprintf( stderr, "in _tanni ist,ien = %d,%d, fulllen1=%d\n", ist, ien, fulllen1 );
 304         fprintf( stderr, "in _tanni jst,jen = %d,%d, fulllen2=%d\n", jst, jen, fulllen2 );
 305         fprintf( stderr, "in _tanni seq1[0] = %-*.*s\n", ien-ist+1, ien-ist+1, seq1[0]+ist );
 306         fprintf( stderr, "in _tanni seq2[0] = %-*.*s\n", jen-jst+1, jen-jst+1, seq2[0]+jst );
 307 #endif
 308
 309
 310         ll1 = ( (int)(lgth1) ) + 100;
 311         ll2 = ( (int)(lgth2) ) + 100;
 312
 313 //      aseq1 = AllocateCharMtx( icyc, 0 );
 314 //      aseq2 = AllocateCharMtx( jcyc, 0 );
 315 //      aseq1bk = AllocateCharMtx( icyc, lgth1+lgth2+100 );
 316 //      aseq2bk = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
 317 //      for( i=0; i<icyc; i++ ) aseq1[i] = aseq1bk[i];
 318 //      for( i=0; i<jcyc; i++ ) aseq2[i] = aseq2bk[i];
 319
 320         w1 = AllocateFloatVec( ll2+2 );
 321         w2 = AllocateFloatVec( ll2+2 );
 322
 323         initverticalw = AllocateFloatVec( ll1+2 );
 324         lastverticalw = AllocateFloatVec( ll1+2 );
 325
 326         m = AllocateFloatVec( ll2+2 );
 327         mp = AllocateIntVec( ll2+2 );
 328
 329         floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 27 );
 330         intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 );
 331
 332
 333         intmtx = AllocateIntMtx( ll1+1, ll2+1 );
 334
 335         ijp = intmtx;
 336
 337         currentw = w1;
 338         previousw = w2;
 339
 340         match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
 341
 342         match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
 343
 344         if( headgp || ist != 0 )
 345         {
 346                 for( i=1; i<lgth1+1; i++ )
 347                 {
 348                         initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
 349                 }
 350         }
 351         if( headgp || jst != 0 )
 352         {
 353                 for( j=1; j<lgth2+1; j++ )
 354                 {
 355                         currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
 356                 }
 357         }
 358
 359         for( j=1; j<lgth2+1; ++j )
 360         {
 361                 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;;
 362         }
 363
 364         lastverticalw[0] = currentw[lgth2-1];
 365
 366
 367
 368         if( tailgp || jen != fulllen2-1 ) lasti = lgth1+1; else lasti = lgth1;
 369 //      if( 1 ) lasti = lgth1+1; else lasti = lgth1;
 370         for( i=1; i<lasti; i++ )
 371         {
 372                 wtmp = previousw;
 373                 previousw = currentw;
 374                 currentw = wtmp;
 375
 376                 previousw[0] = initverticalw[i-1];
 377
 378                 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
 379                 currentw[0] = initverticalw[i];
 380
 381                 mi = previousw[0] + ogcp2[1];
 382                 mpi = 0;
 383
 384                 ijppt = ijp[i] + 1;
 385                 mjpt = m + 1;
 386                 prept = previousw;
 387                 curpt = currentw + 1;
 388                 mpjpt = mp + 1;
 389 //              if( tailgp && jen != fulllen2-1 ) lastj = lgth2+1; else lastj = lgth2;
 390                 lastj = lgth2+1;
 391                 for( j=1; j<lastj; j++ )
 392                 {
 393                         wm = *prept;
 394                         *ijppt = 0;
 395
 396 #if 0
 397                         fprintf( stderr, "%5.0f->", wm );
 398 #endif
 399                         g = mi + fgcp2[j-1];
 400 #if 0
 401                         fprintf( stderr, "%5.0f?", g );
 402 #endif
 403                         if( g > wm )
 404                         {
 405                                 wm = g;
 406                                 *ijppt = -( j - mpi );
 407                         }
 408                         g = *prept + ogcp2[j];
 409                         if( g >= mi )
 410                         {
 411                                 mi = g;
 412                                 mpi = j-1;
 413                         }
 414 #if USE_PENALTY_EX
 415                         mi += fpenalty_ex;
 416 #endif
 417
 418                         g = *mjpt + fgcp1[i-1];
 419 #if 0
 420                         fprintf( stderr, "%5.0f?", g );
 421 #endif
 422                         if( g > wm )
 423                         {
 424                                 wm = g;
 425                                 *ijppt = +( i - *mpjpt );
 426                         }
 427                         g = *prept + ogcp1[i];
 428                         if( g >= *mjpt )
 429                         {
 430                                 *mjpt = g;
 431                                 *mpjpt = i-1;
 432                         }
 433 #if USE_PENALTY_EX
 434                         m[j] += fpenalty_ex;
 435 #endif
 436
 437 #if 0
 438                         fprintf( stderr, "%5.0f ", wm );
 439 #endif
 440                         *curpt += wm;
 441
 442
 443                         ijppt++;
 444                         mjpt++;
 445                         prept++;
 446                         mpjpt++;
 447                         curpt++;
 448                 }
 449                 lastverticalw[i] = currentw[lgth2-1];
 450         }
 451
 452 //      fprintf( stderr, "wm = %f\n", wm );
 453
 454         Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, ist, ien, jst, jen, fulllen1, fulllen2, tailgp );
 455 #if 0
 456         fprintf( stderr, "res in _tanni mseq1[0] = %s\n", mseq1[0] );
 457         fprintf( stderr, "res in _tanni mseq2[0] = %s\n", mseq2[0] );
 458 #endif
 459
 460 //      for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
 461 //      for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
 462
 463 //      fprintf( stderr, "in _tanni, aseq1 = %s\n", aseq1[0] );
 464 //      fprintf( stderr, "in _tanni, mseq1 = %s\n", mseq1[0] );
 465
 466         FreeFloatVec( w1 );
 467         FreeFloatVec( w2 );
 468         FreeFloatVec( initverticalw );
 469         FreeFloatVec( lastverticalw );
 470
 471         FreeFloatVec( m );
 472         FreeIntVec( mp );
 473
 474
 475         FreeFloatMtx( floatwork );
 476         FreeIntMtx( intwork );
 477
 478         FreeIntMtx( intmtx );
 479
 480
 481 //      FreeCharMtx( aseq1bk );
 482 //      FreeCharMtx( aseq2bk );
 483
 484 //      free( aseq1 );
 485 //      free( aseq2 );
 486
 487         return( wm );
 488
 489 }
 490
 491 static void freearrays_rec1(
 492         float *w1, float *w2, float *initverticalw, float *lastverticalw,
 493         float *midw, float *midm, float *midn,
 494         int *jumpbacki, int *jumpbackj, int *jumpforwi, int *jumpforwj, int *jumpdummi, int *jumpdummj,
 495         float *m, int *mp,
 496         float **floatwork, int **intwork
 497 #if STOREWM
 498         , float **WMMTX, float **WMMTX2
 499 #endif
 500 )
 501 {
 502
 503         FreeFloatVec( w1 );
 504         FreeFloatVec( w2 );
 505         FreeFloatVec( initverticalw );
 506         FreeFloatVec( lastverticalw );
 507         FreeFloatVec( midw );
 508         FreeFloatVec( midm );
 509         FreeFloatVec( midn );
 510
 511         FreeIntVec( jumpbacki );
 512         FreeIntVec( jumpbackj );
 513         FreeIntVec( jumpforwi );
 514         FreeIntVec( jumpforwj );
 515         FreeIntVec( jumpdummi );
 516         FreeIntVec( jumpdummj );
 517
 518         FreeFloatVec( m );
 519         FreeIntVec( mp );
 520
 521         FreeFloatMtx( floatwork );
 522         FreeIntMtx( intwork );
 523
 524 #if STOREWM
 525         FreeFloatMtx( WMMTX );
 526         FreeFloatMtx( WMMTX2 );
 527 #endif
 528 }
 529
 530 static void freearrays_rec2( char *gaps, char **aseq1, char **aseq2 )
 531 {
 532         free( gaps );
 533 #if MEMSAVE
 534         free( aseq1 );
 535         free( aseq2 );
 536 #else
 537         FreeCharMtx( aseq1 );
 538         FreeCharMtx( aseq2 );
 539 #endif
 540 }
 541
 542 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, int fulllen1, int fulllen2, char **mseq1, char **mseq2, int depth, float **gapinfo, int *chudanpt, int chudanref, int *chudanres, int headgp, int tailgp )
 543 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 544 {
 545 //      int k;
 546         int alnlen;
 547         float value = 0.0;
 548         register int i, j;
 549         char **aseq1, **aseq2;
 550         int ll1, ll2, l, len;
 551         int lasti, lastj, imid;
 552         int jmid = 0; // by D.Mathog, a guess
 553         float wm = 0.0;   /* int ?????? */
 554         float g;
 555         float *currentw, *previousw;
 556 #if USE_PENALTY_EX
 557         float fpenalty_ex = (float)penalty_ex;
 558 #endif
 559         float *wtmp;
 560 //      short *ijppt;
 561         int *mpjpt;
 562 //      short **ijp;
 563         int *mp;
 564         int mpi;
 565         float *mjpt, *prept, *curpt;
 566         float mi;
 567         float *m;
 568         float *w1, *w2;
 569 //      float *match;
 570         float *initverticalw;    /* kufuu sureba iranai */
 571         float *lastverticalw;    /* kufuu sureba iranai */
 572         int **intwork;
 573         float **floatwork;
 574 //      short **shortmtx;
 575 #if STOREWM
 576         float **WMMTX;
 577         float **WMMTX2;
 578 #endif
 579         float *midw;
 580         float *midm;
 581         float *midn;
 582         int lgth1, lgth2;
 583         float maxwm;
 584         int *jumpforwi;
 585         int *jumpforwj;
 586         int *jumpbacki;
 587         int *jumpbackj;
 588         int *jumpdummi; //muda
 589         int *jumpdummj = NULL; // by D.Mathog, a guess
 590         int jumpi, jumpj = 0; // by D.Mathog, a guess
 591         char *gaps;
 592         int ijpi, ijpj;
 593         float *ogcp1;
 594         float *fgcp1;
 595         float *ogcp2;
 596         float *fgcp2;
 597         float firstm;
 598         int firstmp;
 599 #if STOREWM
 600         static TLS char ttt1[50000];
 601         static TLS char ttt2[50000];
 602 #endif
 603
 604 #if 0
 605         int nglen1, nglen2;
 606         nglen1 = seqlen( seq1[0] );
 607         nglen2 = seqlen( seq2[0] );
 608 #endif
 609
 610 //      fprintf( stderr, "fulllen1 = %d, fulllen2 = %d, headgp = %d, tailgp = %d\n", fulllen1, fulllen2, headgp, tailgp );
 611
 612         ogcp1 = gapinfo[0] + ist;
 613         fgcp1 = gapinfo[1] + ist;
 614         ogcp2 = gapinfo[2] + jst;
 615         fgcp2 = gapinfo[3] + jst;
 616
 617         depth++;
 618         reccycle++;
 619
 620         lgth1 = ien-ist+1;
 621         lgth2 = jen-jst+1;
 622
 623 //      if( lgth1 < 5 )
 624 //              fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
 625 //      if( lgth2 < 5 )
 626 //              fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
 627 //
 628
 629
 630 #if STOREWM
 631         fprintf( stderr, "==== MSalign (depth=%d, reccycle=%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, reccycle, ist, ien, jst, jen );
 632         strncpy( ttt1, seq1[0]+ist, lgth1 );
 633         strncpy( ttt2, seq2[0]+jst, lgth2 );
 634         ttt1[lgth1] = 0;
 635         ttt2[lgth2] = 0;
 636         fprintf( stderr, "seq1 = %s\n", ttt1 );
 637         fprintf( stderr, "seq2 = %s\n", ttt2 );
 638 #endif
 639         if( lgth2 <= 0 ) // lgth1 <= 0 ha?
 640         {
 641 //              fprintf( stderr, "\n\n==== jimei\n\n" );
 642 //              exit( 1 );
 643                 for( i=0; i<icyc; i++ )
 644                 {
 645                         strncpy( mseq1[i], seq1[i]+ist, lgth1 );
 646                         mseq1[i][lgth1] = 0;
 647                 }
 648                 for( i=0; i<jcyc; i++ )
 649                 {
 650                         mseq2[i][0] = 0;
 651                         for( j=0; j<lgth1; j++ )
 652                                 strcat( mseq2[i], "-" );
 653                 }
 654
 655 //              fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
 656 //              fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
 657
 658                 return( 0.0 );
 659         }
 660
 661 #if MEMSAVE
 662         aseq1 = AllocateCharMtx( icyc, 0 );
 663         aseq2 = AllocateCharMtx( jcyc, 0 );
 664         for( i=0; i<icyc; i++ ) aseq1[i] = mseq1[i];
 665         for( i=0; i<jcyc; i++ ) aseq2[i] = mseq2[i];
 666 #else
 667         aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
 668         aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
 669 #endif
 670
 671 //      fprintf( stderr, "####(s) seq1[0] (%d) = \n%-*.*s\n (a%d-a%d)\n", depth, ien-ist+1, ien-ist+1, seq1[0]+ist, ist, ien );
 672 //      fprintf( stderr, "####(s) seq2[0] (%d) = \n%-*.*s\n (b%d-b%d)\n", depth, jen-jst+1, jen-jst+1, seq2[0]+jst, jst, jen );
 673
 674 //  if( lgth1 < DPTANNI && lgth2 < DPTANNI ) // & dato lgth ==1 no kanousei ga arunode yokunai
 675 //    if( lgth1 < DPTANNI ) // kore mo lgth2 ga mijikasugiru kanousei ari
 676     if( lgth1 < DPTANNI || lgth2 < DPTANNI ) // zettai ni anzen ka?
 677         {
 678 //              fprintf( stderr, "==== Going to _tanni\n" );
 679
 680                 value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, fulllen1, fulllen2, aseq1, aseq2, gapinfo, headgp, tailgp );
 681
 682
 683 #if MEMSAVE
 684                 free( aseq1 );
 685                 free( aseq2 );
 686 #else
 687                 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
 688                 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
 689
 690                 FreeCharMtx( aseq1 );
 691                 FreeCharMtx( aseq2 );
 692 #endif
 693
 694 //              fprintf( stderr, "value = %f\n", value );
 695
 696                 return( value );
 697         }
 698 //      fprintf( stderr, "Trying to divide the mtx\n" );
 699
 700         ll1 = ( (int)(lgth1) ) + 100;
 701         ll2 = ( (int)(lgth2) ) + 100;
 702
 703 //      fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
 704
 705         w1 = AllocateFloatVec( ll2+2 );
 706         w2 = AllocateFloatVec( ll2+2 );
 707 //      match = AllocateFloatVec( ll2+2 );
 708         midw = AllocateFloatVec( ll2+2 );
 709         midn = AllocateFloatVec( ll2+2 );
 710         midm = AllocateFloatVec( ll2+2 );
 711         jumpbacki = AllocateIntVec( ll2+2 );
 712         jumpbackj = AllocateIntVec( ll2+2 );
 713         jumpforwi = AllocateIntVec( ll2+2 );
 714         jumpforwj = AllocateIntVec( ll2+2 );
 715         jumpdummi = AllocateIntVec( ll2+2 );
 716         jumpdummj = AllocateIntVec( ll2+2 );
 717
 718         initverticalw = AllocateFloatVec( ll1+2 );
 719         lastverticalw = AllocateFloatVec( ll1+2 );
 720
 721         m = AllocateFloatVec( ll2+2 );
 722         mp = AllocateIntVec( ll2+2 );
 723         gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
 724
 725         floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
 726         intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
 727
 728 #if DEBUG
 729         fprintf( stderr, "succeeded\n" );
 730 #endif
 731
 732 #if STOREWM
 733         WMMTX = AllocateFloatMtx( ll1, ll2 );
 734         WMMTX2 = AllocateFloatMtx( ll1, ll2 );
 735 #endif
 736 #if 0
 737         shortmtx = AllocateShortMtx( ll1, ll2 );
 738
 739 #if DEBUG
 740         fprintf( stderr, "succeeded\n\n" );
 741 #endif
 742
 743         ijp = shortmtx;
 744 #endif
 745
 746         currentw = w1;
 747         previousw = w2;
 748
 749         match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
 750
 751         match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
 752
 753         for( i=1; i<lgth1+1; i++ )
 754         {
 755                 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
 756         }
 757         for( j=1; j<lgth2+1; j++ )
 758         {
 759                 currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
 760         }
 761
 762 #if STOREWM
 763         WMMTX[0][0] = initverticalw[0];
 764         for( i=1; i<lgth1+1; i++ )
 765         {
 766                 WMMTX[i][0] = initverticalw[i];
 767         }
 768         for( j=1; j<lgth2+1; j++ )
 769         {
 770                 WMMTX[0][j] = currentw[j];
 771         }
 772 #endif
 773
 774
 775         for( j=1; j<lgth2+1; ++j )
 776         {
 777                 m[j] = currentw[j-1] + ogcp1[1];
 778 //              m[j] = currentw[j-1];
 779                 mp[j] = 0;
 780         }
 781
 782         lastverticalw[0] = currentw[lgth2-1];
 783
 784         imid = lgth1 * 0.5;
 785
 786         jumpi = 0; // atode kawaru.
 787         lasti = lgth1+1;
 788 #if STOREWM
 789         for( i=1; i<lasti; i++ )
 790 #else
 791         for( i=1; i<=imid; i++ )
 792 #endif
 793         {
 794 #ifdef enablemultithread
 795 //              fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
 796                 if( chudanpt && *chudanpt != chudanref )
 797                 {
 798 //                      fprintf( stderr, "\n\n## CHUUDAN!!! zenhan\n" );
 799                         *chudanres = 1;
 800                         freearrays_rec1
 801                         (
 802                                 w1, w2, initverticalw, lastverticalw, midw, midm, midn,
 803                                 jumpbacki, jumpbackj, jumpforwi, jumpforwj, jumpdummi, jumpdummj,
 804                                 m, mp,
 805                                 floatwork, intwork
 806 #if STOREWM
 807                                 , WMMTX, WMMTX2
 808 #endif
 809                         );
 810                         freearrays_rec2( gaps, aseq1, aseq2 );
 811                         return( -1.0 );
 812                 }
 813 #endif
 814                 wtmp = previousw;
 815                 previousw = currentw;
 816                 currentw = wtmp;
 817
 818                 previousw[0] = initverticalw[i-1];
 819
 820                 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
 821                 currentw[0] = initverticalw[i];
 822
 823                 m[0] = ogcp1[i];
 824 #if STOREM
 825                 WMMTX2[i][0] = m[0];
 826 #endif
 827                 if( i == imid ) midm[0] = m[0];
 828
 829                 mi = previousw[0] + ogcp2[1];
 830 //              mi = previousw[0];
 831                 mpi = 0;
 832
 833
 834 //              ijppt = ijp[i] + 1;
 835                 mjpt = m + 1;
 836                 prept = previousw;
 837                 curpt = currentw + 1;
 838                 mpjpt = mp + 1;
 839
 840
 841                 lastj = lgth2+1;
 842                 for( j=1; j<lastj; j++ )
 843                 {
 844
 845                         wm = *prept;
 846
 847 #if 0
 848                         fprintf( stderr, "%5.0f->", wm );
 849 #endif
 850                         g = mi + fgcp2[j-1];
 851 //                      g = mi + fpenalty;
 852 #if 0
 853                         fprintf( stderr, "%5.0f?", g );
 854 #endif
 855                         if( g > wm )
 856                         {
 857                                 wm = g;
 858 //                              *ijppt = -( j - mpi );
 859                         }
 860                         g = *prept + ogcp2[j];
 861 //                      g = *prept;
 862                         if( g >= mi )
 863                         {
 864                                 mi = g;
 865                                 mpi = j-1;
 866                         }
 867 #if USE_PENALTY_EX
 868                         mi += fpenalty_ex;
 869 #endif
 870
 871                         g = *mjpt + fgcp1[i-1];
 872 //                      g = *mjpt + fpenalty;
 873 #if 0
 874                         fprintf( stderr, "%5.0f?", g );
 875 #endif
 876                         if( g > wm )
 877                         {
 878                                 wm = g;
 879 //                              *ijppt = +( i - *mpjpt );
 880                         }
 881
 882
 883                         g = *prept + ogcp1[i];
 884 //                      g = *prept;
 885                         if( g >= *mjpt )
 886                         {
 887                                 *mjpt = g;
 888                                 *mpjpt = i-1;
 889                         }
 890 #if USE_PENALTY_EX
 891                         m[j] += fpenalty_ex;
 892 #endif
 893
 894 #if 0
 895                         fprintf( stderr, "%5.0f ", wm );
 896 #endif
 897                         *curpt += wm;
 898
 899 #if STOREWM
 900                         WMMTX[i][j] = *curpt;
 901                         WMMTX2[i][j] = *mjpt;
 902 #endif
 903
 904                         if( i == imid ) //muda
 905                         {
 906                                 jumpbackj[j] = *mpjpt; // muda atode matomeru
 907                                 jumpbacki[j] = mpi; // muda atode matomeru
 908 //                              fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
 909 //                              fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
 910                                 midw[j] = *curpt;
 911                                 midm[j] = *mjpt;
 912                                 midn[j] = mi;
 913                         }
 914
 915 //                      fprintf( stderr, "m[%d] = %f\n", j, m[j] );
 916                         mjpt++;
 917                         prept++;
 918                         mpjpt++;
 919                         curpt++;
 920
 921                 }
 922                 lastverticalw[i] = currentw[lgth2-1];
 923
 924 #if STOREWM
 925                 WMMTX2[i][lgth2] = m[lgth2-1];
 926 #endif
 927
 928 #if 0  // ue
 929                 if( i == imid )
 930                 {
 931                         for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
 932                         for( j=0; j<lgth2; j++ ) midm[j] = m[j];
 933                 }
 934 #endif
 935         }
 936 //      for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
 937 //      for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
 938
 939 #if 0
 940     for( i=0; i<lgth1; i++ )
 941     {
 942         for( j=0; j<lgth2; j++ )
 943         {
 944             fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
 945         }
 946         fprintf( stderr, "\n" );
 947     }
 948         fprintf( stderr, "\n" );
 949         fprintf( stderr, "WMMTX2 = \n" );
 950     for( i=0; i<lgth1; i++ )
 951     {
 952         for( j=0; j<lgth2; j++ )
 953         {
 954             fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
 955         }
 956         fprintf( stderr, "\n" );
 957     }
 958         fprintf( stderr, "\n" );
 959 #endif
 960
 961 // gyakudp
 962
 963         match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
 964         match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
 965
 966         for( i=0; i<lgth1-1; i++ )
 967         {
 968                 initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
 969 //              initverticalw[i] += fpenalty;
 970         }
 971         for( j=0; j<lgth2-1; j++ )
 972         {
 973                 currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
 974 //              currentw[j] += fpenalty;
 975         }
 976
 977 #if STOREWM
 978         for( i=0; i<lgth1-1; i++ )
 979         {
 980                 WMMTX[i][lgth2-1] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
 981                 fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
 982         }
 983         for( j=0; j<lgth2-1; j++ )
 984         {
 985                 WMMTX[lgth1-1][j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
 986                 fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
 987         }
 988 #endif
 989
 990
 991
 992
 993
 994
 995         for( j=lgth2-1; j>0; --j )
 996         {
 997                 m[j-1] = currentw[j] + fgcp2[lgth2-2];
 998 //              m[j-1] = currentw[j];
 999                 mp[j] = lgth1-1;
1000         }
1001
1002 //      for( j=0; j<lgth2; j++ ) m[j] = 0.0;
1003         // m[lgth2-1] ha irunoka?
1004
1005
1006 //      for( i=lgth1-2; i>=imid; i-- )
1007         firstm = -9999999.9;
1008 //      firstmp = lgth1-1;
1009         firstmp = lgth1;
1010         for( i=lgth1-2; i>-1; i-- )
1011         {
1012 #ifdef enablemultithread
1013 //              fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
1014                 if( chudanpt && *chudanpt != chudanref )
1015                 {
1016 //                      fprintf( stderr, "\n\n## CHUUDAN!!! kouhan\n" );
1017                         *chudanres = 1;
1018                         freearrays_rec1
1019                         (
1020                                 w1, w2, initverticalw, lastverticalw, midw, midm, midn,
1021                                 jumpbacki, jumpbackj, jumpforwi, jumpforwj, jumpdummi, jumpdummj,
1022                                 m, mp,
1023                                 floatwork, intwork
1024 #if STOREWM
1025                                 , WMMTX, WMMTX2
1026 #endif
1027                         );
1028                         freearrays_rec2( gaps, aseq1, aseq2 );
1029                         return( -1.0 );
1030                 }
1031 #endif
1032                 wtmp = previousw;
1033                 previousw = currentw;
1034                 currentw = wtmp;
1035                 previousw[lgth2-1] = initverticalw[i+1];
1036 //              match_calc( currentw, seq1, seq2, i, lgth2 );
1037                 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
1038
1039                 currentw[lgth2-1] = initverticalw[i];
1040
1041 //              m[lgth2] = fgcp1[i];
1042 //              WMMTX2[i][lgth2] += m[lgth2];
1043 //              fprintf( stderr, "m[] = %f\n", m[lgth2] );
1044
1045                 mi = previousw[lgth2-1] + fgcp2[lgth2-2];
1046 //              mi = previousw[lgth2-1];
1047                 mpi = lgth2 - 1;
1048
1049                 mjpt = m + lgth2 - 2;
1050                 prept = previousw + lgth2 - 1;
1051                 curpt = currentw + lgth2 - 2;
1052                 mpjpt = mp + lgth2 - 2;
1053
1054
1055                 for( j=lgth2-2; j>-1; j-- )
1056                 {
1057                         wm = *prept;
1058                         ijpi = i+1;
1059                         ijpj = j+1;
1060
1061                         g = mi + ogcp2[j+1];
1062 //                      g = mi + fpenalty;
1063                         if( g > wm )
1064                         {
1065                                 wm = g;
1066                                 ijpj = mpi;
1067                                 ijpi = i+1;
1068                         }
1069
1070                         g = *prept + fgcp2[j];
1071 //                      g = *prept;
1072                         if( g >= mi )
1073                         {
1074 //                              fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
1075                                 mi = g;
1076                                 mpi = j + 1;
1077                         }
1078
1079 #if USE_PENALTY_EX
1080                         mi += fpenalty_ex;
1081 #endif
1082
1083 //                      fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
1084                         g = *mjpt + ogcp1[i+1];
1085 //                      g = *mjpt + fpenalty;
1086                         if( g > wm )
1087                         {
1088                                 wm = g;
1089                                 ijpi = *mpjpt;
1090                                 ijpj = j+1;
1091                         }
1092
1093 //                      if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
1094                         g = *prept + fgcp1[i];
1095 //                      g = *prept;
1096                         if( g >= *mjpt )
1097                         {
1098                                 *mjpt = g;
1099                                 *mpjpt = i + 1;
1100                         }
1101
1102 #if USE_PENALTY_EX
1103                         m[j] += fpenalty_ex;
1104 #endif
1105
1106                         if( i == jumpi || i == imid - 1 )
1107                         {
1108                                 jumpforwi[j] = ijpi; //muda
1109                                 jumpforwj[j] = ijpj; //muda
1110 //                              fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
1111 //                              fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
1112                         }
1113                         if( i == imid ) // muda
1114                         {
1115                                 midw[j] += wm;
1116 //                              midm[j+1] += *mjpt + fpenalty; //??????
1117                                 midm[j+1] += *mjpt; //??????
1118                         }
1119                         if( i == imid - 1 )
1120                         {
1121 //                              midn[j] += mi + fpenalty;  //????
1122                                 midn[j] += mi;  //????
1123                         }
1124 #if STOREWM
1125                         WMMTX[i][j] += wm;
1126 //                      WMMTX2[i][j+1] += *mjpt + fpenalty;
1127                         WMMTX2[i][j+1] += *mjpt;
1128 #endif
1129                         *curpt += wm;
1130
1131                         mjpt--;
1132                         prept--;
1133                         mpjpt--;
1134                         curpt--;
1135                 }
1136 //              fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
1137                 g = *prept + fgcp1[i];
1138                 if( firstm < g )
1139                 {
1140                         firstm = g;
1141                         firstmp = i + 1;
1142                 }
1143 #if STOREWM
1144                 WMMTX2[i][j+1] += firstm;
1145 #endif
1146                 if( i == imid ) midm[j+1] += firstm;
1147
1148
1149                 if( i == imid - 1 )
1150                 {
1151                         maxwm = midw[1];
1152                         jmid = 0;
1153 //                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1154                         for( j=2; j<lgth2-1; j++ )
1155                         {
1156                                 wm = midw[j];
1157                                 if( wm > maxwm )
1158                                 {
1159                                         jmid = j;
1160                                         maxwm = wm;
1161                                 }
1162 //                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1163                         }
1164                         for( j=0; j<lgth2+1; j++ )
1165                         {
1166                                 wm = midm[j];
1167                                 if( wm > maxwm )
1168                                 {
1169                                         jmid = j;
1170                                         maxwm = wm;
1171                                 }
1172 //                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1173                         }
1174
1175 //                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1176
1177
1178 //                      fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
1179                         wm = midw[jmid];
1180                         jumpi = imid-1;
1181                         jumpj = jmid-1;
1182                         if( jmid > 0 && midn[jmid-1] > wm ) //060413
1183                         {
1184                                 jumpi = imid-1;
1185                                 jumpj = jumpbacki[jmid];
1186                                 wm = midn[jmid-1];
1187 //                              fprintf( stderr, "rejump (n)\n" );
1188                         }
1189                         if( midm[jmid] > wm )
1190                         {
1191                                 jumpi = jumpbackj[jmid];
1192                                 jumpj = jmid-1;
1193                                 wm = midm[jmid];
1194 //                              fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
1195                         }
1196
1197
1198 //                      fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
1199 //                      fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
1200 #if STOREWM
1201                         fprintf( stderr, "imid = %d\n", imid );
1202                         fprintf( stderr, "midn = \n" );
1203                         for( j=0; j<lgth2; j++ )
1204                         {
1205                                 fprintf( stderr, "% 7.1f ", midn[j] );
1206                         }
1207                         fprintf( stderr, "\n" );
1208                         fprintf( stderr, "midw = \n" );
1209                         for( j=0; j<lgth2; j++ )
1210                         {
1211                                 fprintf( stderr, "% 7.1f ", midw[j] );
1212                         }
1213                         fprintf( stderr, "\n" );
1214                         fprintf( stderr, "midm = \n" );
1215                         for( j=0; j<lgth2; j++ )
1216                         {
1217                                 fprintf( stderr, "% 7.1f ", midm[j] );
1218                         }
1219                         fprintf( stderr, "\n" );
1220 #endif
1221 //                      fprintf( stderr, "maxwm = %f\n", maxwm );
1222                 }
1223                 if( i == jumpi ) //saki?
1224                 {
1225 //                      fprintf( stderr, "#### FIRST i=%d, jumpi<imid=%d<%d, ist=%d, ien=%d, firstmp=%d, lgth1=%d\n", i, jumpi, imid, ist, ien, firstmp, lgth1 );
1226 //                      fprintf( stderr, "#### mark 1\n" );
1227 //                      fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
1228 //                      fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
1229
1230                         if( jmid == 0 )
1231                         {
1232 //                              fprintf( stderr, "#### CHUI2!\n" );
1233                                 jumpj = 0; jmid = 1;
1234 #if 0 // v6.823 made
1235                                 jumpi = firstmp-1;
1236                                 imid = firstmp;
1237 #endif
1238 #if 0
1239                                 jumpi = 0;
1240                                 imid = 1;
1241 #else
1242 //                              if( 1 || firstmp > 100 ) // naze 100
1243                                 if( imid < firstmp-1 ) // naze 100
1244                                 {
1245                                         jumpi = firstmp;
1246                                         imid = firstmp+1;
1247                                 }
1248 #if 0
1249                                 else
1250                                 {
1251                                         jumpi = 0;
1252                                         imid = 1;
1253                                 }
1254 #endif
1255 #endif
1256                         }
1257
1258 #if 0
1259                         else if( jmid == lgth2 )
1260                         {
1261                                 fprintf( stderr, "CHUI1!\n" );
1262                                 jumpi=0; jumpj=0;
1263                                 imid=jumpforwi[0]; jmid=lgth2-1;
1264                         }
1265 #else // 060414
1266                         else if( jmid >= lgth2 )
1267                         {
1268 //                              fprintf( stderr, "CHUI1!\n" );
1269                                 jumpi=imid-1; jmid=lgth2;
1270                                 jumpj = lgth2-1;
1271                         }
1272 #endif
1273                         else
1274                         {
1275 //                              fprintf( stderr, "#### CHUI3!\n" );
1276                                 imid = jumpforwi[jumpj];
1277                                 jmid = jumpforwj[jumpj];
1278                                 if( imid == jumpi ) jumpi = imid-1;
1279                         }
1280 #if 0
1281                         fprintf( stderr, "jumpi -> %d\n", jumpi );
1282                         fprintf( stderr, "jumpj -> %d\n", jumpj );
1283                         fprintf( stderr, "imid -> %d\n", imid );
1284                         fprintf( stderr, "jmid -> %d\n", jmid );
1285 #endif
1286 //                      fprintf( stderr, "#### FINAL i=%d, jumpi<imid=%d<%d, ist=%d, ien=%d, firstmp=%d\n", i, jumpi, imid, ist, ien, firstmp );
1287
1288 #if STOREWM
1289                         break;
1290 #else
1291                         break;
1292 #endif
1293                 }
1294         }
1295
1296 //      fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1297 //      fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1298
1299 //      for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1300 //      for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1301 //      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1302 //      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1303
1304
1305 #if STOREWM
1306         fprintf( stderr, "WMMTX = \n" );
1307     for( i=0; i<lgth1; i++ )
1308     {
1309         fprintf( stderr, "%d ", i );
1310         for( j=0; j<lgth2; j++ )
1311         {
1312             fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
1313         }
1314         fprintf( stderr, "\n" );
1315     }
1316 //      fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
1317     for( i=0; i<lgth1; i++ )
1318     {
1319         fprintf( stderr, "%d ", i );
1320         for( j=0; j<lgth2+1; j++ )
1321         {
1322             fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
1323         }
1324         fprintf( stderr, "\n" );
1325     }
1326
1327         fprintf( stderr, "jumpbacki = \n" );
1328         for( j=0; j<lgth2; j++ )
1329         {
1330                 fprintf( stderr, "% 7d ", jumpbacki[j] );
1331         }
1332         fprintf( stderr, "\n" );
1333         fprintf( stderr, "jumpbackj = \n" );
1334         for( j=0; j<lgth2; j++ )
1335         {
1336                 fprintf( stderr, "% 7d ", jumpbackj[j] );
1337         }
1338         fprintf( stderr, "\n" );
1339         fprintf( stderr, "jumpforwi = \n" );
1340         for( j=0; j<lgth2; j++ )
1341         {
1342                 fprintf( stderr, "% 7d ", jumpforwi[j] );
1343         }
1344         fprintf( stderr, "\n" );
1345         fprintf( stderr, "jumpforwj = \n" );
1346         for( j=0; j<lgth2; j++ )
1347         {
1348                 fprintf( stderr, "% 7d ", jumpforwj[j] );
1349         }
1350         fprintf( stderr, "\n" );
1351
1352
1353 #endif
1354
1355
1356 //      Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1357
1358 #if 0 // irukamo
1359         resultlen = strlen( mseq1[0] );
1360         if( alloclen < resultlen || resultlen > N )
1361         {
1362                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1363                 ErrorExit( "LENGTH OVER!\n" );
1364         }
1365 #endif
1366
1367
1368
1369 #if 0
1370         fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
1371         fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
1372
1373         fprintf( stderr, "imid = %d\n", imid );
1374         fprintf( stderr, "jmid = %d\n", jmid );
1375 #endif
1376
1377         freearrays_rec1
1378         (
1379                 w1, w2, initverticalw, lastverticalw, midw, midm, midn,
1380                 jumpbacki, jumpbackj, jumpforwi, jumpforwj, jumpdummi, jumpdummj,
1381                 m, mp,
1382                 floatwork, intwork
1383 #if STOREWM
1384                 , WMMTX, WMMTX2
1385 #endif
1386         );
1387
1388
1389 //      fprintf( stderr, "==== calling myself (first)\n" );
1390
1391         value = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ist+jumpi, jst, jst+jumpj, alloclen, fulllen1, fulllen2, aseq1, aseq2, depth, gapinfo, NULL, 0, NULL, headgp, tailgp ); // chudan mada
1392 #if 0
1393                 fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
1394                 fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
1395 #endif
1396 #if MEMSAVE
1397 #else
1398         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
1399         for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
1400 #endif
1401
1402 //      fprintf( stderr, "====(f) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
1403 //      fprintf( stderr, "====(f) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
1404
1405         len = strlen( mseq1[0] );
1406 //      fprintf( stderr, "len = %d\n", len );
1407         l = jmid - jumpj - 1;
1408 //      fprintf( stderr, "l=%d\n", l );
1409         if( l > 0 )
1410         {
1411                 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
1412                 for( i=0; i<icyc; i++ )
1413                 {
1414                         strcat( mseq1[i], gaps );
1415                         mseq1[i][len+l] = 0;
1416                 }
1417                 for( j=0; j<jcyc; j++ )
1418                 {
1419                         strncat( mseq2[j], seq2[j]+jst+jumpj+1, l );
1420                         mseq2[j][len+l] = 0;
1421                 }
1422 //              fprintf( stderr, "penalizing (2) .. %f(%d), %f(%d)\n", ogcp2[jumpj+1], jumpj+1, fgcp2[jmid-1], jmid-1 );
1423                 value +=  ( ogcp2[jumpj+1] + fgcp2[jmid-1] );
1424 //              value += fpenalty;
1425         }
1426         len = strlen( mseq1[0] );
1427         l = imid - jumpi - 1;
1428 //      fprintf( stderr, "l=%d\n", l );
1429         if( l > 0 )
1430         {
1431                 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
1432                 for( i=0; i<icyc; i++ )
1433                 {
1434                         strncat( mseq1[i], seq1[i]+ist+jumpi+1, l );
1435                         mseq1[i][len+l] = 0;
1436                 }
1437                 for( j=0; j<jcyc; j++ )
1438                 {
1439                         strcat( mseq2[j], gaps );
1440                         mseq2[j][len+l] = 0;
1441                 }
1442
1443 //              for( i=0; i<lgth1; i++ ) fprintf( stderr, "ogcp1[%d] = %f\n", i, ogcp1[i] );
1444 //              for( i=0; i<lgth1; i++ ) fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
1445
1446
1447 //              fprintf( stderr, "penalizing (1) .. ogcp1[%d] = %f, fgcp1[%d] = %f\n", jumpi+1, ogcp1[jumpi+1], imid-1, fgcp1[imid-1] );
1448                 value += ( ogcp1[jumpi+1] + fgcp1[imid-1] );
1449 //              value += fpenalty;
1450         }
1451 #if 0
1452         for( i=0; i<icyc; i++ ) fprintf( stderr, "after gapfill mseq1[%d]=%s\n", i, mseq1[i] );
1453         for( i=0; i<jcyc; i++ ) fprintf( stderr, "after gapfill mseq2[%d]=%s\n", i, mseq2[i] );
1454 #endif
1455
1456 //      fprintf( stderr, "==== calling myself (second)\n" );
1457
1458 #if MEMSAVE
1459         alnlen = strlen( aseq1[0] );
1460         for( i=0; i<icyc; i++ ) aseq1[i] += alnlen;
1461         for( i=0; i<jcyc; i++ ) aseq2[i] += alnlen;
1462 #endif
1463
1464         value += MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist+imid, ien, jst+jmid, jen, alloclen, fulllen1, fulllen2, aseq1, aseq2, depth, gapinfo, NULL, 0, NULL, headgp, tailgp ); // chudan mada
1465 #if 0
1466                 fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
1467                 fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
1468 #endif
1469
1470
1471
1472 #if DEBUG
1473         if( value - maxwm > 1 || maxwm - value > 1 )
1474         {
1475                 fprintf( stderr, "WARNING value  = %f, but maxwm = %f\n", value, maxwm );
1476                 for( i=0; i<icyc; i++ )
1477                 {
1478                         fprintf( stderr, ">1-%d\n%s\n", i, mseq1[i] );
1479                         fprintf( stderr, "%s\n", aseq1[i] );
1480                 }
1481                 for( i=0; i<jcyc; i++ )
1482                 {
1483                         fprintf( stderr, ">2-%d\n%s\n", i, mseq2[i] );
1484                         fprintf( stderr, "%s\n", aseq2[i] );
1485                 }
1486
1487 //              exit( 1 );
1488         }
1489         else
1490         {
1491                 fprintf( stderr, "value = %.0f, maxwm = %.0f -> ok\n", value, maxwm );
1492         }
1493 #endif
1494
1495 #if MEMSAVE
1496 #else
1497         for( i=0; i<icyc; i++ ) strcat( mseq1[i], aseq1[i] );
1498         for( i=0; i<jcyc; i++ ) strcat( mseq2[i], aseq2[i] );
1499 #endif
1500
1501 //      fprintf( stderr, "====(s) aseq1[0] (%d) = \n%s\n (a%d-a%d)\n", depth, mseq1[0], ist, ien );
1502 //      fprintf( stderr, "====(s) aseq2[0] (%d) = \n%s\n (b%d-b%d)\n", depth, mseq2[0], jst, jen );
1503
1504         freearrays_rec2( gaps, aseq1, aseq2 );
1505
1506 #if 0
1507         if( seqlen( seq1[0] ) != nglen1 )
1508         {
1509                 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
1510                 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
1511                 exit( 1 );
1512         }
1513         else
1514                 fprintf( stderr, "nglen1 is ok in _rec\n" );
1515         if( seqlen( seq2[0] ) != nglen2 )
1516         {
1517                 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
1518                 exit( 1 );
1519         }
1520         else
1521                 fprintf( stderr, "nglen2 is ok in _rec\n" );
1522 #endif
1523
1524         return( value );
1525 }
1526
1527 static void freearrays(
1528         float *ogcp1,
1529         float *ogcp2,
1530         float *fgcp1,
1531         float *fgcp2,
1532         float **cpmx1,
1533         float **cpmx2,
1534         float **gapinfo,
1535         char **mseq1,
1536         char **mseq2
1537 )
1538 {
1539         FreeFloatVec( ogcp1 );
1540         FreeFloatVec( ogcp2 );
1541         FreeFloatVec( fgcp1 );
1542         FreeFloatVec( fgcp2 );
1543         FreeFloatMtx( cpmx1 );
1544         FreeFloatMtx( cpmx2 );
1545         free( (void *)gapinfo );
1546
1547         FreeCharMtx( mseq1 );
1548         FreeCharMtx( mseq2 );
1549 }
1550
1551 float MSalignmm( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, int *chudanpt, int chudanref, int *chudanres, int headgp, int tailgp )
1552 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1553 {
1554 //      int k;
1555         int i, j;
1556         int ll1, ll2;
1557         int lgth1, lgth2;
1558         float wm = 0.0;   /* int ?????? */
1559         char **mseq1;
1560         char **mseq2;
1561         float *ogcp1;
1562         float *ogcp2;
1563         float *fgcp1;
1564         float *fgcp2;
1565         float **cpmx1;
1566         float **cpmx2;
1567         float **gapinfo;
1568         float fpenalty = (float)penalty;
1569         int nglen1, nglen2;
1570
1571 #if 0
1572         fprintf( stderr, "eff in SA+++align\n" );
1573         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1574 #endif
1575
1576         nglen1 = seqlen( seq1[0] );
1577         nglen2 = seqlen( seq2[0] );
1578
1579 #if 0
1580         fprintf( stderr, "\n" );
1581         for( i=0; i<icyc; i++ ) fprintf( stderr, "seq1[%d] at root = %s\n", i, seq1[i] );
1582         for( j=0; j<jcyc; j++ ) fprintf( stderr, "seq2[%d] at root = %s\n", j, seq2[j] );
1583         fprintf( stderr, "\n" );
1584 #endif
1585
1586         lgth1 = strlen( seq1[0] );
1587         lgth2 = strlen( seq2[0] );
1588
1589         ll1 = ( (int)(lgth1) ) + 100;
1590         ll2 = ( (int)(lgth2) ) + 100;
1591
1592         mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
1593         mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
1594
1595         gapinfo = AllocateFloatMtx( 4, 0 );
1596         ogcp1 = AllocateFloatVec( ll1+2 );
1597         ogcp2 = AllocateFloatVec( ll2+2 );
1598         fgcp1 = AllocateFloatVec( ll1+2 );
1599         fgcp2 = AllocateFloatVec( ll2+2 );
1600
1601
1602         cpmx1 = AllocateFloatMtx( ll1+2, 27 );
1603         cpmx2 = AllocateFloatMtx( ll2+2, 27 );
1604
1605         for( i=0; i<icyc; i++ )
1606         {
1607                 if( strlen( seq1[i] ) != lgth1 )
1608                 {
1609                         fprintf( stderr, "i = %d / %d\n", i, icyc );
1610                         fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
1611                         exit( 1 );
1612                 }
1613         }
1614         for( j=0; j<jcyc; j++ )
1615         {
1616                 if( strlen( seq2[j] ) != lgth2 )
1617                 {
1618                         fprintf( stderr, "j = %d / %d\n", j, jcyc );
1619                         fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
1620                         exit( 1 );
1621                 }
1622         }
1623
1624         MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1625         MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1626
1627
1628 #if 1
1629
1630         if( sgap1 )
1631         {
1632                 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
1633                 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
1634                 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap2 );
1635                 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
1636         }
1637         else
1638         {
1639                 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1640                 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1641                 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1642                 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1643         }
1644
1645 #if 1
1646         for( i=0; i<lgth1; i++ )
1647         {
1648                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
1649                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
1650 //              fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
1651         }
1652         for( i=0; i<lgth2; i++ )
1653         {
1654                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
1655                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
1656 //              fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
1657         }
1658 #else
1659         for( i=0; i<lgth1; i++ )
1660         {
1661                 ogcp1[i] = 0.5 * fpenalty;
1662                 fgcp1[i] = 0.5 * fpenalty;
1663         }
1664         for( i=0; i<lgth2; i++ )
1665         {
1666                 ogcp2[i] = 0.5 * fpenalty;
1667                 fgcp2[i] = 0.5 * fpenalty;
1668         }
1669 #endif
1670
1671         gapinfo[0] = ogcp1;
1672         gapinfo[1] = fgcp1;
1673         gapinfo[2] = ogcp2;
1674         gapinfo[3] = fgcp2;
1675 #endif
1676
1677 #if 0
1678         fprintf( stdout, "in MSalignmm.c\n" );
1679         for( i=0; i<icyc; i++ )
1680         {
1681                 fprintf( stdout, ">%d of GROUP1\n", i );
1682                 fprintf( stdout, "%s\n", seq1[i] );
1683         }
1684         for( i=0; i<jcyc; i++ )
1685         {
1686                 fprintf( stdout, ">%d of GROUP2\n", i );
1687                 fprintf( stdout, "%s\n", seq2[i] );
1688         }
1689         fflush( stdout );
1690 #endif
1691
1692         wm = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, lgth1, lgth2, mseq1, mseq2, 0, gapinfo, chudanpt, chudanref, chudanres, headgp, tailgp );
1693 #ifdef enablemultithread
1694         if( chudanres && *chudanres )
1695         {
1696 //              fprintf( stderr, "\n\n## CHUUDAN!!! relay\n" );
1697                 *chudanres = 1;
1698                 freearrays( ogcp1, ogcp2, fgcp1, fgcp2, cpmx1, cpmx2, gapinfo, mseq1, mseq2 );
1699                 return( -1.0 );
1700         }
1701 #endif
1702
1703 #if DEBUG
1704                 fprintf( stderr, " seq1[0] = %s\n", seq1[0] );
1705                 fprintf( stderr, " seq2[0] = %s\n", seq2[0] );
1706                 fprintf( stderr, "mseq1[0] = %s\n", mseq1[0] );
1707                 fprintf( stderr, "mseq2[0] = %s\n", mseq2[0] );
1708 #endif
1709
1710 //      fprintf( stderr, "wm = %f\n", wm );
1711
1712
1713         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1714         for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
1715
1716         if( seqlen( seq1[0] ) != nglen1 )
1717         {
1718                 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
1719                 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
1720                 exit( 1 );
1721         }
1722         if( seqlen( seq2[0] ) != nglen2 )
1723         {
1724                 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
1725                 exit( 1 );
1726         }
1727
1728
1729         freearrays( ogcp1, ogcp2, fgcp1, fgcp2, cpmx1, cpmx2, gapinfo, mseq1, mseq2 );
1730
1731         lgth1 = strlen( seq1[0] );
1732         lgth2 = strlen( seq2[0] );
1733         for( i=0; i<icyc; i++ )
1734         {
1735                 if( strlen( seq1[i] ) != lgth1 )
1736                 {
1737                         fprintf( stderr, "i = %d / %d\n", i, icyc );
1738                         fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
1739                         exit( 1 );
1740                 }
1741         }
1742         for( j=0; j<jcyc; j++ )
1743         {
1744                 if( strlen( seq2[j] ) != lgth2 )
1745                 {
1746                         fprintf( stderr, "j = %d / %d\n", j, jcyc );
1747                         fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
1748                         exit( 1 );
1749                 }
1750         }
1751
1752         return( wm );
1753 }