binaries/src/mafft/core/Lalign11.c.fast

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