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

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