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

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