binaries/src/mafft/core/MSalignmm.c.algLargeA

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