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

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