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

   1 #include "mltaln.h"
   2 #include "dp.h"
   3
   4 #define DEBUG 1
   5 #define XXXXXXX    0
   6 #define USE_PENALTY_EX  0
   7
   8 #define DPTANNI 100
   9
  10 static void OpeningGapCount( double *ogcp, int clus, char **seq, double *eff )
  11 {
  12         int i, j, gc, gb;
  13         int len = strlen( seq[0] );
  14         double totaleff = 0.0;
  15
  16         for( i=0; i<len; i++ ) ogcp[i] = 0.0;
  17         for( j=0; j<clus; j++ )
  18         {
  19                 gc = 0;
  20                 for( i=0; i<len; i++ )
  21                 {
  22                         gb = gc;
  23                         gc = ( seq[j][i] == '-' );
  24                         {
  25                                 if( !gb *  gc ) ogcp[i] += eff[j];
  26                         }
  27                 }
  28                 totaleff+= eff[j];
  29         }
  30         for( i=0; i<len; i++ )
  31                 ogcp[i] /= totaleff;
  32 }
  33
  34 static void FinalGapCount( double *fgcp, int clus, char **seq, double *eff )
  35 {
  36         int i, j, gc, gb;
  37         int len = strlen( seq[0] );
  38         double totaleff = 0.0;
  39
  40         for( i=0; i<len; i++ ) fgcp[i] = 0.0;
  41         for( j=0; j<clus; j++ )
  42         {
  43                 gc = ( seq[j][0] == '-' );
  44                 for( i=1; i<len+1; i++ )
  45                 {
  46                         gb = gc;
  47                         gc = ( seq[j][i] == '-' );
  48                         {
  49                                 if( gb * !gc ) fgcp[i-1] += eff[j];
  50                         }
  51                 }
  52                 totaleff += eff[j];
  53         }
  54         for( i=0; i<len; i++ )
  55                 fgcp[i] /= totaleff;
  56 }
  57
  58 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
  59 {
  60         int j, k, l;
  61         float scarr[26];
  62         float **cpmxpd = floatwork;
  63         int **cpmxpdn = intwork;
  64         int count = 0;
  65
  66         if( initialize )
  67         {
  68                 for( j=0; j<lgth2; j++ )
  69                 {
  70                         count = 0;
  71                         for( l=0; l<26; l++ )
  72                         {
  73                                 if( cpmx2[j][l] )
  74                                 {
  75                                         cpmxpd[j][count] = cpmx2[j][l];
  76                                         cpmxpdn[j][count] = l;
  77                                         count++;
  78                                 }
  79                         }
  80                         cpmxpdn[j][count] = -1;
  81                 }
  82         }
  83
  84         for( l=0; l<26; l++ )
  85         {
  86                 scarr[l] = 0.0;
  87                 for( k=0; k<26; k++ )
  88                 {
  89                         scarr[l] += n_dis[k][l] * cpmx1[i1][k];
  90                 }
  91         }
  92 #if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
  93         {
  94                 float *fpt, **fptpt, *fpt2;
  95                 int *ipt, **iptpt;
  96                 fpt2 = match;
  97                 iptpt = cpmxpdn;
  98                 fptpt = cpmxpd;
  99                 while( lgth2-- )
 100                 {
 101                         *fpt2 = 0.0;
 102                         ipt=*iptpt,fpt=*fptpt;
 103                         while( *ipt > -1 )
 104                                 *fpt2 += scarr[*ipt++] * *fpt++;
 105                         fpt2++,iptpt++,fptpt++;
 106                 }
 107         }
 108 #else
 109         for( j=0; j<lgth2; j++ )
 110         {
 111                 match[j] = 0.0;
 112                 for( k=0; cpmxpdn[j][k]>-1; k++ )
 113                         match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
 114         }
 115 #endif
 116 }
 117
 118 static float Atracking( float *lasthorizontalw, float *lastverticalw,
 119                                                 char **seq1, char **seq2,
 120                         char **mseq1, char **mseq2,
 121                         float **cpmx1, float **cpmx2,
 122                         short **ijp, int icyc, int jcyc,
 123                                                 int ist, int ien, int jst, int jen )
 124 {
 125         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, klim;
 126         char gap[] = "-";
 127         float wm;
 128         lgth1 = ien-ist+1;
 129         lgth2 = jen-jst+1;
 130
 131 #if 1
 132         for( i=0; i<lgth1; i++ )
 133         {
 134                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 135         }
 136 #endif
 137
 138     for( i=0; i<lgth1+1; i++ )
 139     {
 140         ijp[i][0] = i + 1;
 141     }
 142     for( j=0; j<lgth2+1; j++ )
 143     {
 144         ijp[0][j] = -( j + 1 );
 145     }
 146
 147         for( i=0; i<icyc; i++ )
 148         {
 149                 mseq1[i] += lgth1+lgth2;
 150                 *mseq1[i] = 0;
 151         }
 152         for( j=0; j<jcyc; j++ )
 153         {
 154                 mseq2[j] += lgth1+lgth2;
 155                 *mseq2[j] = 0;
 156         }
 157         iin = lgth1; jin = lgth2;
 158         klim = lgth1+lgth2;
 159         for( k=0; k<=klim; k++ )
 160         {
 161                 if( ijp[iin][jin] < 0 )
 162                 {
 163                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 164                 }
 165                 else if( ijp[iin][jin] > 0 )
 166                 {
 167                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 168                 }
 169                 else
 170                 {
 171                         ifi = iin-1; jfi = jin-1;
 172                 }
 173                 l = iin - ifi;
 174                 while( --l )
 175                 {
 176                         for( i=0; i<icyc; i++ )
 177                                 *--mseq1[i] = seq1[i][ist+ifi+l];
 178                         for( j=0; j<jcyc; j++ )
 179                                 *--mseq2[j] = *gap;
 180                         k++;
 181                 }
 182                 l= jin - jfi;
 183                 while( --l )
 184                 {
 185                         for( i=0; i<icyc; i++ )
 186                                 *--mseq1[i] = *gap;
 187                         for( j=0; j<jcyc; j++ )
 188                                 *--mseq2[j] = seq2[j][jst+jfi+l];
 189                         k++;
 190                 }
 191                 if( iin <= 0 || jin <= 0 ) break;
 192                 for( i=0; i<icyc; i++ )
 193                         *--mseq1[i] = seq1[i][ist+ifi];
 194                 for( j=0; j<jcyc; j++ )
 195                         *--mseq2[j] = seq2[j][jst+jfi];
 196                 k++;
 197                 iin = ifi; jin = jfi;
 198
 199                 fprintf( stderr, "in Atracking, mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
 200         }
 201         fprintf( stderr, "in Atracking (owari), mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
 202         return( 0.0 );
 203 }
 204
 205 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  )
 206 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 207 {
 208 //      int k;
 209         register int i, j;
 210         int ll1, ll2;
 211         int lasti, lastj;
 212         int resultlen;
 213         float wm;   /* int ?????? */
 214         float g;
 215         float *currentw, *previousw;
 216         float fpenalty = (float)penalty;
 217         float fpenalty_ex = (float)penalty_ex;
 218 #if 1
 219         float *wtmp;
 220         short *ijppt;
 221         float *mjpt, *prept, *curpt;
 222         int *mpjpt;
 223 #endif
 224         float mi, *m;
 225         short **ijp;
 226         int mpi, *mp;
 227         float *w1, *w2;
 228         float *match;
 229         float *initverticalw;    /* kufuu sureba iranai */
 230         float *lastverticalw;    /* kufuu sureba iranai */
 231         int **intwork;
 232         float **floatwork;
 233         short **shortmtx;
 234         float **WMMTX;
 235         float dumfl;
 236         int lgth1, lgth2;
 237         static char **aseq1 = NULL;
 238         static char **aseq2 = NULL;
 239         static char **aseq1bk, **aseq2bk;
 240
 241         if( !aseq1 )
 242         {
 243                 aseq1 = AllocateCharMtx( icyc, 0 );
 244                 aseq2 = AllocateCharMtx( jcyc, 0 );
 245         }
 246
 247         lgth1 = ien-ist+1;
 248         lgth2 = jen-jst+1;
 249
 250         fprintf( stderr, "seq1[0]+ist = %s\n", seq1[0]+ist );
 251         fprintf( stderr, "seq2[0]+jst = %s\n", seq2[0]+jst );
 252
 253         fprintf( stderr, "ist,ien = %d,%d, lgth1=%d\n", ist, ien, lgth1 );
 254         fprintf( stderr, "jst,jen = %d,%d, lgth2=%d\n", jst, jen, lgth2 );
 255
 256
 257         ll1 = ( (int)(1.3*lgth1) ) + 100;
 258         ll2 = ( (int)(1.3*lgth2) ) + 100;
 259
 260         aseq1bk = AllocateCharMtx( icyc, lgth1+lgth2+100 );
 261         aseq2bk = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
 262         for( i=0; i<icyc; i++ ) aseq1[i] = aseq1bk[i];
 263         for( i=0; i<jcyc; i++ ) aseq2[i] = aseq2bk[i];
 264
 265         w1 = AllocateFloatVec( ll2+2 );
 266         w2 = AllocateFloatVec( ll2+2 );
 267         match = AllocateFloatVec( ll2+2 );
 268
 269         initverticalw = AllocateFloatVec( ll1+2 );
 270         lastverticalw = AllocateFloatVec( ll1+2 );
 271
 272         m = AllocateFloatVec( ll2+2 );
 273         mp = AllocateIntVec( ll2+2 );
 274
 275         floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
 276         intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
 277
 278
 279         shortmtx = AllocateShortMtx( ll1, ll2 );
 280         WMMTX = AllocateFloatMtx( ll1, ll2 );
 281
 282         ijp = shortmtx;
 283
 284         currentw = w1;
 285         previousw = w2;
 286
 287         fprintf( stderr, "calling match_calc\n" );
 288         match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
 289
 290         match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
 291
 292         WMMTX[0][0] = initverticalw[0];
 293         for( i=1; i<lgth1+1; i++ )
 294         {
 295                 initverticalw[i] += fpenalty;
 296                 WMMTX[i][0] = initverticalw[i];
 297         }
 298         for( j=1; j<lgth2+1; j++ )
 299         {
 300                 currentw[j] += fpenalty;
 301                 WMMTX[0][j] = currentw[j];
 302         }
 303
 304         for( j=1; j<lgth2+1; ++j )
 305         {
 306                 m[j] = currentw[j-1];
 307         }
 308
 309         lastverticalw[0] = currentw[lgth2-1];
 310
 311         fprintf( stderr, "entering to loop\n" );
 312
 313         lasti = lgth1+1;
 314         for( i=1; i<lasti; i++ )
 315         {
 316
 317                 fprintf( stderr, "i=%d\n", i );
 318                 wtmp = previousw;
 319                 previousw = currentw;
 320                 currentw = wtmp;
 321
 322                 previousw[0] = initverticalw[i-1];
 323
 324                 fprintf( stderr, "calling match_calc, ist=%d, lgth2=%d\n", ist, lgth2 );
 325                 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
 326                 fprintf( stderr, "done\n" );
 327                 currentw[0] = initverticalw[i];
 328
 329                 mi = previousw[0]; mpi = 0;
 330
 331                 ijppt = ijp[i] + 1;
 332                 mjpt = m + 1;
 333                 prept = previousw;
 334                 curpt = currentw + 1;
 335                 mpjpt = mp + 1;
 336                 lastj = lgth2+1;
 337                 for( j=1; j<lastj; j++ )
 338                 {
 339                         fprintf( stderr, "j=%d\n", j );
 340                         wm = *prept;
 341                         *ijppt = 0;
 342
 343 #if 0
 344                         fprintf( stderr, "%5.0f->", wm );
 345 #endif
 346                         g = mi + fpenalty;
 347 #if 0
 348                         fprintf( stderr, "%5.0f?", g );
 349 #endif
 350                         if( g > wm )
 351                         {
 352                                 wm = g;
 353                                 *ijppt = -( j - mpi );
 354                         }
 355                         g = *prept;
 356                         if( g >= mi )
 357                         {
 358                                 mi = g;
 359                                 mpi = j-1;
 360                         }
 361 #if USE_PENALTY_EX
 362                         mi += fpenalty_ex;
 363 #endif
 364
 365                         g = *mjpt + fpenalty;
 366 #if 0
 367                         fprintf( stderr, "%5.0f?", g );
 368 #endif
 369                         if( g > wm )
 370                         {
 371                                 wm = g;
 372                                 *ijppt = +( i - *mpjpt );
 373                         }
 374                         g = *prept;
 375                         if( g >= *mjpt )
 376                         {
 377                                 *mjpt = g;
 378                                 *mpjpt = i-1;
 379                         }
 380 #if USE_PENALTY_EX
 381                         m[j] += fpenalty_ex;
 382 #endif
 383
 384 #if 0
 385                         fprintf( stderr, "%5.0f ", wm );
 386 #endif
 387                         *curpt += wm;
 388
 389                         WMMTX[i][j] = *curpt;
 390
 391                         ijppt++;
 392                         mjpt++;
 393                         prept++;
 394                         mpjpt++;
 395                         curpt++;
 396                         fprintf( stderr, "j=%d end\n", j );
 397                 }
 398                 lastverticalw[i] = currentw[lgth2-1];
 399                 fprintf( stderr, "i=%d end\n", i );
 400         }
 401
 402 #if 1
 403     for( i=0; i<lgth1; i++ )
 404     {
 405         for( j=0; j<lgth2; j++ )
 406         {
 407             fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
 408         }
 409         fprintf( stderr, "\n" );
 410     }
 411         fprintf( stderr, "\n" );
 412 #endif
 413
 414         Atracking( currentw, lastverticalw, seq1, seq2, aseq1, aseq2, cpmx1+ist, cpmx2+jst, ijp, icyc, jcyc, ist, ien, jst, jen );
 415
 416         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
 417         for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
 418
 419         fprintf( stderr, "in _tanni, aseq1 = %s\n", aseq1[0] );
 420         fprintf( stderr, "in _tanni, mseq1 = %s\n", mseq1[0] );
 421
 422         FreeFloatVec( w1 );
 423         FreeFloatVec( w2 );
 424         FreeFloatVec( match );
 425         FreeFloatVec( initverticalw );
 426         FreeFloatVec( lastverticalw );
 427
 428         FreeFloatVec( m );
 429         FreeIntVec( mp );
 430
 431
 432         FreeFloatMtx( floatwork );
 433         FreeIntMtx( intwork );
 434
 435         FreeShortMtx( shortmtx );
 436         FreeFloatMtx( WMMTX );
 437
 438         FreeCharMtx( aseq1bk );
 439         FreeCharMtx( aseq2bk );
 440
 441         return( wm );
 442 }
 443
 444 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 )
 445 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 446 {
 447 //      int k;
 448         register int i, j;
 449         char **aseq1, **aseq2;
 450         int ll1, ll2;
 451         int lasti, lastj, imid, jmid;
 452         int resultlen;
 453         float wm;   /* int ?????? */
 454         float g;
 455         float *currentw, *previousw;
 456         float fpenalty = (float)penalty;
 457         float fpenalty_ex = (float)penalty_ex;
 458 #if 1
 459         float *wtmp;
 460         short *ijppt;
 461         float *mjpt, *prept, *curpt;
 462         int *mpjpt;
 463 #endif
 464         float mi, *m;
 465         short **ijp;
 466         int mpi, *mp;
 467         float *w1, *w2;
 468         float *match;
 469         float *initverticalw;    /* kufuu sureba iranai */
 470         float *lastverticalw;    /* kufuu sureba iranai */
 471         int **intwork;
 472         float **floatwork;
 473         short **shortmtx;
 474         float **WMMTX;
 475         float dumfl;
 476         int lgth1, lgth2;
 477         float maxwm;
 478
 479         depth++;
 480
 481         lgth1 = ien-ist+1;
 482         lgth2 = jen-jst+1;
 483
 484         fprintf( stderr, "==== MSalign (%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, ist, ien, jst, jen );
 485         if( lgth2 <= 0 )
 486         {
 487                 fprintf( stderr, "==== Jimei\n" );
 488                 for( i=0; i<icyc; i++ )
 489                 {
 490                         strncpy( mseq1[i], seq1[i]+ist, lgth1 );
 491                         mseq1[i][lgth1] = 0;
 492                 }
 493                 for( i=0; i<jcyc; i++ )
 494                 {
 495                         mseq2[i][0] = 0;
 496                         for( j=0; j<lgth1; j++ )
 497                                 strcat( mseq2[i], "-" );
 498                 }
 499
 500                 fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
 501                 fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
 502
 503                 return( (float)offset * lgth1 );
 504         }
 505
 506
 507         aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
 508         aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
 509
 510         if( lgth1 < DPTANNI )
 511         {
 512                 fprintf( stderr, "==== Going to _tanni\n" );
 513                 wm = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2 );
 514
 515
 516                 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
 517                 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
 518
 519                 fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
 520                 fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
 521
 522                 fprintf( stderr, "freeing aseq\n" );
 523                 FreeCharMtx( aseq1 );
 524                 FreeCharMtx( aseq2 );
 525
 526                 return( wm );
 527         }
 528         fprintf( stderr, "Trying to divide the mtx\n" );
 529
 530
 531         ll1 = ( (int)(1.3*lgth1) ) + 100;
 532         ll2 = ( (int)(1.3*lgth2) ) + 100;
 533
 534         fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
 535
 536         w1 = AllocateFloatVec( ll2+2 );
 537         w2 = AllocateFloatVec( ll2+2 );
 538         match = AllocateFloatVec( ll2+2 );
 539
 540         initverticalw = AllocateFloatVec( ll1+2 );
 541         lastverticalw = AllocateFloatVec( ll1+2 );
 542
 543         m = AllocateFloatVec( ll2+2 );
 544         mp = AllocateIntVec( ll2+2 );
 545
 546         floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
 547         intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
 548
 549 #if DEBUG
 550         fprintf( stderr, "succeeded\n" );
 551 #endif
 552
 553
 554 #if DEBUG
 555         fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
 556 #endif
 557
 558         shortmtx = AllocateShortMtx( ll1, ll2 );
 559         WMMTX = AllocateFloatMtx( ll1, ll2 );
 560
 561 #if DEBUG
 562         fprintf( stderr, "succeeded\n\n" );
 563 #endif
 564
 565         ijp = shortmtx;
 566
 567         currentw = w1;
 568         previousw = w2;
 569
 570         match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
 571
 572         match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
 573
 574         WMMTX[0][0] = initverticalw[0];
 575         for( i=1; i<lgth1+1; i++ )
 576         {
 577                 initverticalw[i] += fpenalty;
 578                 WMMTX[i][0] = initverticalw[i];
 579         }
 580         for( j=1; j<lgth2+1; j++ )
 581         {
 582                 currentw[j] += fpenalty;
 583                 WMMTX[0][j] = currentw[j];
 584         }
 585
 586         for( j=1; j<lgth2+1; ++j )
 587         {
 588                 m[j] = currentw[j-1];
 589         }
 590
 591         lastverticalw[0] = currentw[lgth2-1];
 592
 593         imid = lgth1 / 2;
 594
 595         lasti = lgth1+1;
 596 //      for( i=1; i<lasti; i++ )
 597         for( i=1; i<=imid; i++ )
 598         {
 599                 wtmp = previousw;
 600                 previousw = currentw;
 601                 currentw = wtmp;
 602
 603                 previousw[0] = initverticalw[i-1];
 604
 605                 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
 606                 currentw[0] = initverticalw[i];
 607
 608                 mi = previousw[0]; mpi = 0;
 609
 610                 ijppt = ijp[i] + 1;
 611                 mjpt = m + 1;
 612                 prept = previousw;
 613                 curpt = currentw + 1;
 614                 mpjpt = mp + 1;
 615                 lastj = lgth2+1;
 616                 for( j=1; j<lastj; j++ )
 617                 {
 618                         wm = *prept;
 619                         *ijppt = 0;
 620
 621 #if 0
 622                         fprintf( stderr, "%5.0f->", wm );
 623 #endif
 624                         g = mi + fpenalty;
 625 #if 0
 626                         fprintf( stderr, "%5.0f?", g );
 627 #endif
 628                         if( g > wm )
 629                         {
 630                                 wm = g;
 631                                 *ijppt = -( j - mpi );
 632                         }
 633                         g = *prept;
 634                         if( g >= mi )
 635                         {
 636                                 mi = g;
 637                                 mpi = j-1;
 638                         }
 639 #if USE_PENALTY_EX
 640                         mi += fpenalty_ex;
 641 #endif
 642
 643                         g = *mjpt + fpenalty;
 644 #if 0
 645                         fprintf( stderr, "%5.0f?", g );
 646 #endif
 647                         if( g > wm )
 648                         {
 649                                 wm = g;
 650                                 *ijppt = +( i - *mpjpt );
 651                         }
 652                         g = *prept;
 653                         if( g >= *mjpt )
 654                         {
 655                                 *mjpt = g;
 656                                 *mpjpt = i-1;
 657                         }
 658 #if USE_PENALTY_EX
 659                         m[j] += fpenalty_ex;
 660 #endif
 661
 662 #if 0
 663                         fprintf( stderr, "%5.0f ", wm );
 664 #endif
 665                         *curpt += wm;
 666
 667                         WMMTX[i][j] = *curpt;
 668
 669                         ijppt++;
 670                         mjpt++;
 671                         prept++;
 672                         mpjpt++;
 673                         curpt++;
 674                 }
 675                 lastverticalw[i] = currentw[lgth2-1];
 676         }
 677
 678 #if 0
 679     for( i=0; i<lgth1; i++ )
 680     {
 681         for( j=0; j<lgth2; j++ )
 682         {
 683             fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
 684         }
 685         fprintf( stderr, "\n" );
 686     }
 687         fprintf( stderr, "\n" );
 688 #endif
 689
 690 // gyakudp
 691
 692         match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
 693         match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
 694
 695         for( i=0; i<lgth1-1; i++ )
 696         {
 697                 initverticalw[i] += fpenalty;
 698                 WMMTX[i][lgth2-1] += fpenalty;
 699         }
 700         for( j=0; j<lgth2-1; j++ )
 701         {
 702                 currentw[j] += fpenalty;
 703                 WMMTX[lgth1-1][j] += fpenalty;
 704         }
 705
 706         for( j=lgth2-1; j>0; --j )
 707         {
 708                 m[j-1] = currentw[j];
 709         }
 710
 711 //      for( j=0; j<lgth2; j++ ) m[j] = 0.0;
 712         // m[lgth2-1] ha irunoka?
 713
 714         for( i=lgth1-2; i>=imid; i-- )
 715 //      for( i=lgth1-2; i>-1; i-- )
 716         {
 717                 wtmp = previousw;
 718                 previousw = currentw;
 719                 currentw = wtmp;
 720                 previousw[lgth2-1] = initverticalw[i+1];
 721 //              match_calc( currentw, seq1, seq2, i, lgth2 );
 722                 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
 723
 724                 currentw[lgth2-1] = initverticalw[i];
 725
 726                 mi = previousw[lgth2-1];
 727
 728                 mjpt = m + lgth2 - 2;
 729                 prept = previousw + lgth2 - 1;
 730                 curpt = currentw + lgth2 - 2;
 731                 mpjpt = mp + lgth2 - 2;
 732
 733                 for( j=lgth2-2; j>-1; j-- )
 734                 {
 735                         wm = *prept;
 736                         g = mi + fpenalty;
 737                         if( g > wm ) wm = g;
 738
 739                         g = *prept;
 740                         if( g >= mi ) mi = g;
 741
 742 #if USE_PENALTY_EX
 743                         mi += fpenalty_ex;
 744 #endif
 745
 746                         g = *mjpt + fpenalty;
 747                         if( g > wm ) wm = g;
 748
 749                         g = *prept;
 750                         if( g >= *mjpt ) *mjpt = g;
 751
 752 #if USE_PENALTY_EX
 753                         m[j] += fpenalty_ex;
 754 #endif
 755
 756                         WMMTX[i][j] += wm;
 757                         *curpt += wm;
 758
 759                         mjpt--;
 760                         prept--;
 761                         mpjpt--;
 762                         curpt--;
 763                 }
 764         }
 765
 766 #if 0
 767     for( i=0; i<lgth1; i++ )
 768     {
 769         for( j=0; j<lgth2; j++ )
 770         {
 771             fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
 772         }
 773         fprintf( stderr, "\n" );
 774     }
 775 #endif
 776
 777 //      Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
 778
 779 #if 0 // irukamo
 780         resultlen = strlen( mseq1[0] );
 781         if( alloclen < resultlen || resultlen > N )
 782         {
 783                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
 784                 ErrorExit( "LENGTH OVER!\n" );
 785         }
 786 #endif
 787
 788
 789         maxwm = -999999999.9;
 790         jmid = -100;
 791         for( j=0; j<lgth2; j++ )
 792         {
 793                 wm = WMMTX[imid][j];
 794                 if( wm > maxwm )
 795                 {
 796                         jmid = j;
 797                         maxwm = wm;
 798                 }
 799         }
 800         // gap no tochu de kirrareru kanousei ari.
 801
 802         fprintf( stderr, "wm = %f\n", wm );
 803         fprintf( stderr, "imid = %d\n", imid );
 804         fprintf( stderr, "jmid = %d\n", jmid );
 805
 806         FreeFloatVec( w1 );
 807         FreeFloatVec( w2 );
 808         FreeFloatVec( match );
 809         FreeFloatVec( initverticalw );
 810         FreeFloatVec( lastverticalw );
 811
 812         FreeFloatVec( m );
 813         FreeIntVec( mp );
 814
 815         FreeFloatMtx( floatwork );
 816         FreeIntMtx( intwork );
 817
 818         FreeShortMtx( shortmtx );
 819         FreeFloatMtx( WMMTX );
 820
 821
 822         fprintf( stderr, "==== calling myself (first)\n" );
 823         MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ist+imid, jst, jst+jmid, alloclen, aseq1, aseq2, depth );
 824         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
 825         for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
 826
 827         fprintf( stderr, "====(f) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
 828         fprintf( stderr, "====(f) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
 829
 830
 831         fprintf( stderr, "==== calling myself (second)\n" );
 832         MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist+imid+1, ien, jst+jmid+1, jen, alloclen, aseq1, aseq2, depth );
 833
 834         for( i=0; i<icyc; i++ ) strcat( mseq1[i], aseq1[i] );
 835         for( i=0; i<jcyc; i++ ) strcat( mseq2[i], aseq2[i] );
 836
 837         fprintf( stderr, "====(s) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
 838         fprintf( stderr, "====(s) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
 839
 840         fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
 841         fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
 842
 843         FreeCharMtx( aseq1 );
 844         FreeCharMtx( aseq2 );
 845
 846         return( wm );
 847 }
 848
 849
 850
 851 float MSalignmm( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen )
 852 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 853 {
 854 //      int k;
 855         int i, j;
 856         int ll1, ll2;
 857         int lgth1, lgth2;
 858         float wm;   /* int ?????? */
 859         static char **mseq1 = NULL;
 860         static char **mseq2 = NULL;
 861         char **mseq;
 862         double *ogcp1;
 863         double *ogcp2;
 864         double *fgcp1;
 865         double *fgcp2;
 866         float **cpmx1;
 867         float **cpmx2;
 868
 869 #if 1
 870         fprintf( stderr, "eff in SA+++align\n" );
 871         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
 872 #endif
 873         if( mseq1 == NULL )
 874         {
 875                 mseq1 = AllocateCharMtx( njob, 0 );
 876                 mseq2 = AllocateCharMtx( njob, 0 );
 877         }
 878
 879
 880
 881
 882         lgth1 = strlen( seq1[0] );
 883         lgth2 = strlen( seq2[0] );
 884
 885         ll1 = ( (int)(1.3*lgth1) ) + 100;
 886         ll2 = ( (int)(1.3*lgth2) ) + 100;
 887
 888         mseq = AllocateCharMtx( njob, ll1+ll2 );
 889
 890         ogcp1 = AllocateDoubleVec( ll1+2 );
 891         ogcp2 = AllocateDoubleVec( ll2+2 );
 892         fgcp1 = AllocateDoubleVec( ll1+2 );
 893         fgcp2 = AllocateDoubleVec( ll2+2 );
 894
 895         cpmx1 = AllocateFloatMtx( ll1+2, 26 );
 896         cpmx2 = AllocateFloatMtx( ll2+2, 26 );
 897
 898         for( i=0; i<icyc; i++ ) mseq1[i] = mseq[i];
 899         for( j=0; j<jcyc; j++ ) mseq2[j] = mseq[icyc+j];
 900
 901
 902         MScpmx_calc_new( seq1, cpmx1, eff1, strlen( seq1[0] ), icyc );
 903         MScpmx_calc_new( seq2, cpmx2, eff2, strlen( seq2[0] ), jcyc );
 904
 905 #if 0
 906         OpeningGapCount( ogcp1, icyc, seq1, eff1 );
 907         OpeningGapCount( ogcp2, jcyc, seq2, eff2 );
 908         FinalGapCount( fgcp1, icyc, seq1, eff1 );
 909         FinalGapCount( fgcp2, jcyc, seq2, eff2 );
 910
 911         for( i=0; i<lgth1; i++ )
 912         {
 913                 ogcp1[i] = 1.0 - ogcp1[i];
 914                 fgcp1[i] = 1.0 - fgcp1[i];
 915         }
 916         for( i=0; i<lgth2; i++ )
 917         {
 918                 ogcp2[i] = 1.0 - ogcp2[i];
 919                 fgcp2[i] = 1.0 - fgcp2[i];
 920         }
 921 #endif
 922
 923         wm = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0 );
 924
 925         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
 926         for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
 927
 928
 929         FreeDoubleVec( ogcp1 );
 930         FreeDoubleVec( ogcp2 );
 931         FreeDoubleVec( fgcp1 );
 932         FreeFloatMtx( cpmx1 );
 933         FreeFloatMtx( cpmx2 );
 934         FreeDoubleVec( fgcp2 );
 935
 936         FreeCharMtx( mseq );
 937
 938         return( wm );
 939 }