website/archive/binaries/mac/src/mafft/core/MSalignmm.c.bkbk

   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 10
  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 #if STOREWM
 518         static char ttt1[50000];
 519         static char ttt2[50000];
 520 #endif
 521
 522
 523         ogcp1 = gapinfo[0] + ist;
 524         fgcp1 = gapinfo[1] + ist;
 525         ogcp2 = gapinfo[2] + jst;
 526         fgcp2 = gapinfo[3] + jst;
 527
 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 )
 546         {
 547                 fprintf( stderr, "==== jimei\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( (float)offset * lgth1 );
 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         for( j=1; j<lgth2+1; ++j )
 670         {
 671                 m[j] = currentw[j-1] + fpenalty * ogcp1[1];
 672 //              m[j] = currentw[j-1];
 673                 mp[j] = 0;
 674         }
 675
 676         lastverticalw[0] = currentw[lgth2-1];
 677
 678         imid = lgth1 * 0.5;
 679
 680         jumpi = 0; // atode kawaru.
 681         lasti = lgth1+1;
 682 #if STOREWM
 683         for( i=1; i<lasti; i++ )
 684 #else
 685         for( i=1; i<=imid; i++ )
 686 #endif
 687         {
 688                 wtmp = previousw;
 689                 previousw = currentw;
 690                 currentw = wtmp;
 691
 692                 previousw[0] = initverticalw[i-1];
 693
 694                 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
 695                 currentw[0] = initverticalw[i];
 696
 697
 698                 mi = previousw[0] + fpenalty * ogcp2[1];
 699 //              mi = previousw[0];
 700                 mpi = 0;
 701
 702 //              ijppt = ijp[i] + 1;
 703                 mjpt = m + 1;
 704                 prept = previousw;
 705                 curpt = currentw + 1;
 706                 mpjpt = mp + 1;
 707
 708
 709                 lastj = lgth2+1;
 710                 for( j=1; j<lastj; j++ )
 711                 {
 712
 713                         wm = *prept;
 714
 715 #if 0
 716                         fprintf( stderr, "%5.0f->", wm );
 717 #endif
 718                         g = mi + fpenalty * fgcp2[j-1];
 719 //                      g = mi + fpenalty;
 720 #if 0
 721                         fprintf( stderr, "%5.0f?", g );
 722 #endif
 723                         if( g > wm )
 724                         {
 725                                 wm = g;
 726 //                              *ijppt = -( j - mpi );
 727                         }
 728                         g = *prept + fpenalty * ogcp2[j];
 729 //                      g = *prept;
 730                         if( g >= mi )
 731                         {
 732                                 mi = g;
 733                                 mpi = j-1;
 734                         }
 735 #if USE_PENALTY_EX
 736                         mi += fpenalty_ex;
 737 #endif
 738
 739                         g = *mjpt + fpenalty  * fgcp1[i-1];
 740 //                      g = *mjpt + fpenalty;
 741 #if 0
 742                         fprintf( stderr, "%5.0f?", g );
 743 #endif
 744                         if( g > wm )
 745                         {
 746                                 wm = g;
 747 //                              *ijppt = +( i - *mpjpt );
 748                         }
 749
 750
 751                         g = *prept + fpenalty * ogcp1[i];
 752 //                      g = *prept;
 753                         if( g >= *mjpt )
 754                         {
 755                                 *mjpt = g;
 756                                 *mpjpt = i-1;
 757                         }
 758 #if USE_PENALTY_EX
 759                         m[j] += fpenalty_ex;
 760 #endif
 761
 762 #if 0
 763                         fprintf( stderr, "%5.0f ", wm );
 764 #endif
 765                         *curpt += wm;
 766
 767 #if STOREWM
 768                         WMMTX[i][j] = *curpt;
 769                         WMMTX2[i][j] = *mjpt;
 770 #endif
 771
 772                         if( i == imid ) //muda
 773                         {
 774                                 jumpbackj[j] = *mpjpt; // muda atode matomeru
 775                                 jumpbacki[j] = mpi; // muda atode matomeru
 776 //                              fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
 777 //                              fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
 778                                 midw[j] = *curpt;
 779                                 midm[j] = *mjpt;
 780                                 midn[j] = mi;
 781                         }
 782
 783                         mjpt++;
 784                         prept++;
 785                         mpjpt++;
 786                         curpt++;
 787
 788                 }
 789                 lastverticalw[i] = currentw[lgth2-1];
 790 #if 0  // ue
 791                 if( i == imid )
 792                 {
 793                         for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
 794                         for( j=0; j<lgth2; j++ ) midm[j] = m[j];
 795                 }
 796 #endif
 797         }
 798 //      for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
 799 //      for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
 800
 801 #if 0
 802     for( i=0; i<lgth1; i++ )
 803     {
 804         for( j=0; j<lgth2; j++ )
 805         {
 806             fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
 807         }
 808         fprintf( stderr, "\n" );
 809     }
 810         fprintf( stderr, "\n" );
 811         fprintf( stderr, "WMMTX2 = \n" );
 812     for( i=0; i<lgth1; i++ )
 813     {
 814         for( j=0; j<lgth2; j++ )
 815         {
 816             fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
 817         }
 818         fprintf( stderr, "\n" );
 819     }
 820         fprintf( stderr, "\n" );
 821 #endif
 822
 823 // gyakudp
 824
 825         match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
 826         match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
 827
 828         for( i=0; i<lgth1-1; i++ )
 829         {
 830                 initverticalw[i] += fpenalty * ( fgcp1[lgth1-1] + ogcp1[i+1] );
 831 //              initverticalw[i] += fpenalty;
 832         }
 833         for( j=0; j<lgth2-1; j++ )
 834         {
 835                 currentw[j] += fpenalty * ( fgcp2[lgth2-1] + ogcp2[j+1] );
 836 //              currentw[j] += fpenalty;
 837         }
 838
 839 #if STOREWM
 840         for( i=0; i<lgth1-1; i++ )
 841         {
 842                 WMMTX[i][lgth2-1] += fpenalty * ( fgcp1[lgth1-1] + ogcp1[i+1] );
 843                 fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
 844         }
 845         for( j=0; j<lgth2-1; j++ )
 846         {
 847                 WMMTX[lgth1-1][j] += fpenalty * ( fgcp2[lgth2-1] + ogcp2[j+1] );
 848                 fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
 849         }
 850 #endif
 851
 852         for( j=lgth2-1; j>0; --j )
 853         {
 854                 m[j-1] = currentw[j] + fpenalty * fgcp2[lgth2-2];
 855 //              m[j-1] = currentw[j];
 856                 mp[j] = lgth1-1;
 857         }
 858
 859 //      for( j=0; j<lgth2; j++ ) m[j] = 0.0;
 860         // m[lgth2-1] ha irunoka?
 861
 862 //      for( i=lgth1-2; i>=imid; i-- )
 863         for( i=lgth1-2; i>-1; i-- )
 864         {
 865                 wtmp = previousw;
 866                 previousw = currentw;
 867                 currentw = wtmp;
 868                 previousw[lgth2-1] = initverticalw[i+1];
 869 //              match_calc( currentw, seq1, seq2, i, lgth2 );
 870                 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
 871
 872                 currentw[lgth2-1] = initverticalw[i];
 873
 874                 mi = previousw[lgth2-1] + fpenalty * fgcp2[lgth2-2];
 875 //              mi = previousw[lgth2-1];
 876                 mpi = lgth2 - 1;
 877
 878                 mjpt = m + lgth2 - 2;
 879                 prept = previousw + lgth2 - 1;
 880                 curpt = currentw + lgth2 - 2;
 881                 mpjpt = mp + lgth2 - 2;
 882
 883
 884                 for( j=lgth2-2; j>-1; j-- )
 885                 {
 886                         wm = *prept;
 887                         ijpi = i+1;
 888                         ijpj = j+1;
 889
 890                         g = mi + fpenalty * ogcp2[j+1];
 891 //                      g = mi + fpenalty;
 892                         if( g > wm )
 893                         {
 894                                 wm = g;
 895                                 ijpj = mpi;
 896                         }
 897
 898                         g = *prept + fpenalty * fgcp2[j];
 899 //                      g = *prept;
 900                         if( g >= mi )
 901                         {
 902 //                              fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
 903                                 mi = g;
 904                                 mpi = j + 1;
 905                         }
 906
 907 #if USE_PENALTY_EX
 908                         mi += fpenalty_ex;
 909 #endif
 910
 911 //                      fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
 912                         g = *mjpt + fpenalty * ogcp1[i+1];
 913 //                      g = *mjpt + fpenalty;
 914                         if( g > wm )
 915                         {
 916                                 wm = g;
 917                                 ijpi = *mpjpt;
 918                         }
 919
 920 //                      if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
 921                         g = *prept + fpenalty * fgcp1[i];
 922 //                      g = *prept;
 923                         if( g >= *mjpt )
 924                         {
 925                                 *mjpt = g;
 926                                 *mpjpt = i + 1;
 927                         }
 928
 929 #if USE_PENALTY_EX
 930                         m[j] += fpenalty_ex;
 931 #endif
 932
 933                         if( i == jumpi || i == imid - 1 )
 934                         {
 935                                 jumpforwi[j] = ijpi; //muda
 936                                 jumpforwj[j] = ijpj; //muda
 937 //                              fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
 938 //                              fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
 939                         }
 940                         if( i == imid ) // muda
 941                         {
 942                                 midw[j] += wm;
 943 //                              midm[j+1] += *mjpt + fpenalty; //??????
 944                                 midm[j+1] += *mjpt; //??????
 945                         }
 946                         if( i == imid - 1 )
 947                         {
 948 //                              midn[j] += mi + fpenalty;  //????
 949                                 midn[j] += mi;  //????
 950                         }
 951 #if STOREWM
 952                         WMMTX[i][j] += wm;
 953 //                      WMMTX2[i][j+1] += *mjpt + fpenalty;
 954                         WMMTX2[i][j+1] += *mjpt;
 955 #endif
 956                         *curpt += wm;
 957
 958                         mjpt--;
 959                         prept--;
 960                         mpjpt--;
 961                         curpt--;
 962                 }
 963
 964                 if( i == imid - 1 )
 965                 {
 966                         maxwm = midw[1];
 967 //                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
 968                         for( j=2; j<lgth2-1; j++ )
 969                         {
 970                                 wm = midw[j];
 971                                 if( wm > maxwm )
 972                                 {
 973                                         jmid = j;
 974                                         maxwm = wm;
 975                                         jumpi = imid-1;
 976                                         jumpj = jmid-1;
 977                                 }
 978 //                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
 979                         }
 980                         for( j=1; j<lgth2; j++ )
 981                         {
 982                                 wm = midm[j];
 983                                 if( wm > maxwm )
 984                                 {
 985                                         jmid = j;
 986                                         maxwm = wm;
 987                                         jumpi = jumpbackj[jmid];
 988                                         jumpj = jmid-1;
 989                                 }
 990 //                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
 991                         }
 992                         for( j=2; j<lgth2; j++ )
 993                         {
 994                                 wm = midn[j-1];
 995                                 if( wm > maxwm )
 996                                 {
 997                                         jmid = j-1;
 998                                         maxwm = wm;
 999                                         jumpi = imid-1;
1000                                         jumpj = jumpbacki[jmid];
1001                                 }
1002 //                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1003                         }
1004
1005                         if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1006
1007 #if 0
1008                         fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
1009                         wm = midw[jmid];
1010                         jumpi = imid-1;
1011                         jumpj = jmid-1;
1012                         if( midn[jmid-1] > wm )
1013                         {
1014                                 jumpi = imid-1;
1015                                 jumpj = jumpbacki[jmid];
1016                                 wm = midn[jmid-1];
1017                                 fprintf( stderr, "rejump (n)\n" );
1018                         }
1019                         if( midm[jmid] > wm )
1020                         {
1021                                 jumpi = jumpbackj[jmid];
1022                                 jumpj = jmid-1;
1023                                 wm = midm[jmid];
1024                                 fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
1025                         }
1026 #endif
1027
1028                         fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
1029                         fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
1030 #if STOREWM
1031                         fprintf( stderr, "imid = %d\n", imid );
1032                         fprintf( stderr, "midn = \n" );
1033                         for( j=0; j<lgth2; j++ )
1034                         {
1035                                 fprintf( stderr, "% 7.1f ", midn[j] );
1036                         }
1037                         fprintf( stderr, "\n" );
1038                         fprintf( stderr, "midw = \n" );
1039                         for( j=0; j<lgth2; j++ )
1040                         {
1041                                 fprintf( stderr, "% 7.1f ", midw[j] );
1042                         }
1043                         fprintf( stderr, "\n" );
1044                         fprintf( stderr, "midm = \n" );
1045                         for( j=0; j<lgth2; j++ )
1046                         {
1047                                 fprintf( stderr, "% 7.1f ", midm[j] );
1048                         }
1049                         fprintf( stderr, "\n" );
1050 #endif
1051 //                      fprintf( stderr, "maxwm = %f\n", maxwm );
1052                 }
1053
1054                 if( i == jumpi ) //saki?
1055                 {
1056 //                      fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
1057 //                      fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
1058                         imid = jumpforwi[jumpj];
1059                         jmid = jumpforwj[jumpj];
1060                         fprintf( stderr, "imid -> %d\n", imid );
1061                         fprintf( stderr, "jmid -> %d\n", jmid );
1062 #if STOREWM
1063 #else
1064                         break;
1065 #endif
1066                 }
1067         }
1068 //      fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1069 //      fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1070
1071 //      for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1072 //      for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1073 //      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1074 //      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1075
1076
1077 #if STOREWM
1078         fprintf( stderr, "WMMTX = \n" );
1079     for( i=0; i<lgth1; i++ )
1080     {
1081         fprintf( stderr, "%d ", i );
1082         for( j=0; j<lgth2; j++ )
1083         {
1084             fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
1085         }
1086         fprintf( stderr, "\n" );
1087     }
1088         fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
1089     for( i=0; i<lgth1; i++ )
1090     {
1091         fprintf( stderr, "%d ", i );
1092         for( j=0; j<lgth2; j++ )
1093         {
1094             fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
1095         }
1096         fprintf( stderr, "\n" );
1097     }
1098
1099         fprintf( stderr, "jumpbacki = \n" );
1100         for( j=0; j<lgth2; j++ )
1101         {
1102                 fprintf( stderr, "% 7d ", jumpbacki[j] );
1103         }
1104         fprintf( stderr, "\n" );
1105         fprintf( stderr, "jumpbackj = \n" );
1106         for( j=0; j<lgth2; j++ )
1107         {
1108                 fprintf( stderr, "% 7d ", jumpbackj[j] );
1109         }
1110         fprintf( stderr, "\n" );
1111         fprintf( stderr, "jumpforwi = \n" );
1112         for( j=0; j<lgth2; j++ )
1113         {
1114                 fprintf( stderr, "% 7d ", jumpforwi[j] );
1115         }
1116         fprintf( stderr, "\n" );
1117         fprintf( stderr, "jumpforwj = \n" );
1118         for( j=0; j<lgth2; j++ )
1119         {
1120                 fprintf( stderr, "% 7d ", jumpforwj[j] );
1121         }
1122         fprintf( stderr, "\n" );
1123
1124
1125 #endif
1126
1127 //      Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1128
1129 #if 0 // irukamo
1130         resultlen = strlen( mseq1[0] );
1131         if( alloclen < resultlen || resultlen > N )
1132         {
1133                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1134                 ErrorExit( "LENGTH OVER!\n" );
1135         }
1136 #endif
1137
1138
1139
1140 #if 0
1141         fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
1142         fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
1143
1144         fprintf( stderr, "imid = %d\n", imid );
1145         fprintf( stderr, "jmid = %d\n", jmid );
1146 #endif
1147
1148
1149         FreeFloatVec( w1 );
1150         FreeFloatVec( w2 );
1151         FreeFloatVec( initverticalw );
1152         FreeFloatVec( lastverticalw );
1153         FreeFloatVec( midw );
1154         FreeFloatVec( midm );
1155         FreeFloatVec( midn );
1156
1157         FreeIntVec( jumpbacki );
1158         FreeIntVec( jumpbackj );
1159         FreeIntVec( jumpforwi );
1160         FreeIntVec( jumpforwj );
1161         FreeIntVec( jumpdummi );
1162         FreeIntVec( jumpdummj );
1163
1164         FreeFloatVec( m );
1165         FreeIntVec( mp );
1166
1167         FreeFloatMtx( floatwork );
1168         FreeIntMtx( intwork );
1169
1170 #if STOREWM
1171         FreeFloatMtx( WMMTX );
1172         FreeFloatMtx( WMMTX2 );
1173 #endif
1174
1175 //      fprintf( stderr, "==== calling myself (first)\n" );
1176
1177         fprintf( stderr, "jumpi = %d\n", jumpi );
1178         value = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ist+jumpi, jst, jst+jumpj, alloclen, aseq1, aseq2, depth, gapinfo );
1179         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
1180         for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
1181
1182 //      fprintf( stderr, "====(f) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
1183 //      fprintf( stderr, "====(f) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
1184
1185         len = strlen( mseq1[0] );
1186 //      fprintf( stderr, "len = %d\n", len );
1187         l = jmid - jumpj - 1;
1188 //      fprintf( stderr, "l=%d\n", l );
1189         if( l > 0 )
1190         {
1191                 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
1192                 for( i=0; i<icyc; i++ )
1193                 {
1194                         strcat( mseq1[i], gaps );
1195                         mseq1[i][len+l] = 0;
1196                 }
1197                 for( j=0; j<jcyc; j++ )
1198                 {
1199                         strncat( mseq2[j], seq2[j]+jst+jumpj+1, l );
1200                         mseq2[j][len+l] = 0;
1201                 }
1202                 fprintf( stderr, "penalizing (2) .. %f(%d), %f(%d)\n", ogcp2[jumpj+1], jumpj+1, fgcp2[jmid-1], jmid-1 );
1203                 value += fpenalty * ( ogcp2[jumpj+1] + fgcp2[jmid-1] );
1204 //              value += fpenalty;
1205         }
1206         len = strlen( mseq1[0] );
1207         l = imid - jumpi - 1;
1208 //      fprintf( stderr, "l=%d\n", l );
1209         if( l > 0 )
1210         {
1211                 for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
1212                 for( i=0; i<icyc; i++ )
1213                 {
1214                         strncat( mseq1[i], seq1[i]+ist+jumpi+1, l );
1215                         mseq1[i][len+l] = 0;
1216                 }
1217                 for( j=0; j<jcyc; j++ )
1218                 {
1219                         strcat( mseq2[j], gaps );
1220                         mseq2[j][len+l] = 0;
1221                 }
1222
1223 //              for( i=0; i<lgth1; i++ ) fprintf( stderr, "ogcp1[%d] = %f\n", i, ogcp1[i] );
1224 //              for( i=0; i<lgth1; i++ ) fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
1225
1226
1227                 fprintf( stderr, "penalizing (1) .. ogcp1[%d] = %f, fgcp1[%d] = %f\n", jumpi+1, ogcp1[jumpi+1], imid-1, fgcp1[imid-1] );
1228                 value += fpenalty * ( ogcp1[jumpi+1] + fgcp1[imid-1] );
1229 //              value += fpenalty;
1230         }
1231 #if 0
1232         for( i=0; i<icyc; i++ ) fprintf( stderr, "after gapfill mseq1[%d]=%s\n", i, mseq1[i] );
1233         for( i=0; i<jcyc; i++ ) fprintf( stderr, "after gapfill mseq2[%d]=%s\n", i, mseq2[i] );
1234 #endif
1235
1236 //      fprintf( stderr, "==== calling myself (second)\n" );
1237
1238         value += MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist+imid, ien, jst+jmid, jen, alloclen, aseq1, aseq2, depth, gapinfo );
1239
1240
1241
1242         if( value - maxwm > 1 || maxwm - value > 1 )
1243         {
1244                 fprintf( stderr, "WARNING value  = %f, but maxwm = %f penalty = %f\n", value, maxwm, fpenalty );
1245 #if 1
1246                 for( i=0; i<icyc; i++ )
1247                 {
1248                         fprintf( stderr, ">1-%d\n%s\n", i, mseq1[i] );
1249                         fprintf( stderr, "%s\n", aseq1[i] );
1250                 }
1251                 for( i=0; i<jcyc; i++ )
1252                 {
1253                         fprintf( stderr, ">2-%d\n%s\n", i, mseq2[i] );
1254                         fprintf( stderr, "%s\n", aseq2[i] );
1255                 }
1256 #endif
1257
1258 //              exit( 1 );
1259         }
1260         else
1261         {
1262                 fprintf( stderr, "value = %.0f, maxwm = %.0f -> ok\n", value, maxwm );
1263         }
1264
1265         for( i=0; i<icyc; i++ ) strcat( mseq1[i], aseq1[i] );
1266         for( i=0; i<jcyc; i++ ) strcat( mseq2[i], aseq2[i] );
1267
1268 //      fprintf( stderr, "====(s) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
1269 //      fprintf( stderr, "====(s) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
1270
1271         free( gaps );
1272         FreeCharMtx( aseq1 );
1273         FreeCharMtx( aseq2 );
1274
1275         return( value );
1276 }
1277
1278
1279
1280 float MSalignmm( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen )
1281 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1282 {
1283 //      int k;
1284         int i, j;
1285         int ll1, ll2;
1286         int lgth1, lgth2;
1287         float wm;   /* int ?????? */
1288         char **mseq1;
1289         char **mseq2;
1290         char **mseq;
1291         float *ogcp1;
1292         float *ogcp2;
1293         float *fgcp1;
1294         float *fgcp2;
1295         float **cpmx1;
1296         float **cpmx2;
1297         float **gapinfo;
1298
1299 #if 0
1300         fprintf( stderr, "eff in SA+++align\n" );
1301         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1302 #endif
1303
1304         lgth1 = strlen( seq1[0] );
1305         lgth2 = strlen( seq2[0] );
1306
1307         ll1 = ( (int)(1.3*lgth1) ) + 100;
1308         ll2 = ( (int)(1.3*lgth2) ) + 100;
1309
1310         mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
1311         mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
1312
1313         gapinfo = AllocateFloatMtx( 4, 0 );
1314         ogcp1 = AllocateFloatVec( ll1+2 );
1315         ogcp2 = AllocateFloatVec( ll2+2 );
1316         fgcp1 = AllocateFloatVec( ll1+2 );
1317         fgcp2 = AllocateFloatVec( ll2+2 );
1318
1319
1320         cpmx1 = AllocateFloatMtx( ll1+2, 27 );
1321         cpmx2 = AllocateFloatMtx( ll2+2, 27 );
1322
1323         for( i=0; i<icyc; i++ )
1324         {
1325                 seq1[i][lgth1] = 0;
1326         }
1327         for( j=0; j<jcyc; j++ )
1328         {
1329                 seq2[j][lgth2] = 0;
1330         }
1331
1332         MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1333         MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1334
1335
1336 #if 1
1337
1338         OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1339         OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1340         FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1341         FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1342
1343 #if 1
1344         for( i=0; i<lgth1; i++ )
1345         {
1346                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] );
1347                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] );
1348 //              fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
1349         }
1350         for( i=0; i<lgth2; i++ )
1351         {
1352                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] );
1353                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] );
1354 //              fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
1355         }
1356 #else
1357         for( i=0; i<lgth1; i++ )
1358         {
1359                 ogcp1[i] = 0.5;
1360                 fgcp1[i] = 0.5;
1361         }
1362         for( i=0; i<lgth2; i++ )
1363         {
1364                 ogcp2[i] = 0.5;
1365                 fgcp2[i] = 0.5;
1366         }
1367 #endif
1368
1369         gapinfo[0] = ogcp1;
1370         gapinfo[1] = fgcp1;
1371         gapinfo[2] = ogcp2;
1372         gapinfo[3] = fgcp2;
1373 #endif
1374
1375 #if 0
1376         fprintf( stdout, "in MSalignmm.c\n" );
1377         for( i=0; i<icyc; i++ )
1378         {
1379                 fprintf( stdout, ">%d of GROUP1\n", i );
1380                 fprintf( stdout, "%s\n", seq1[i] );
1381         }
1382         for( i=0; i<jcyc; i++ )
1383         {
1384                 fprintf( stdout, ">%d of GROUP2\n", i );
1385                 fprintf( stdout, "%s\n", seq2[i] );
1386         }
1387         fflush( stdout );
1388 #endif
1389
1390         wm = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo );
1391
1392 //      fprintf( stderr, "wm = %f\n", wm );
1393
1394         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1395         for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
1396
1397
1398         FreeFloatVec( ogcp1 );
1399         FreeFloatVec( ogcp2 );
1400         FreeFloatVec( fgcp1 );
1401         FreeFloatVec( fgcp2 );
1402         FreeFloatMtx( cpmx1 );
1403         FreeFloatMtx( cpmx2 );
1404         free( (void *)gapinfo );
1405
1406         FreeCharMtx( mseq1 );
1407         FreeCharMtx( mseq2 );
1408
1409         return( wm );
1410 }