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

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