binaries/src/mafft/core/partSalignmm.c

   1 #include "mltaln.h"
   2 #include "dp.h"
   3
   4 #define MACHIGAI 0
   5 #define OUTGAP0TRY 1
   6 #define DEBUG 0
   7 #define XXXXXXX    0
   8 #define USE_PENALTY_EX  0
   9 #define FASTMATCHCALC 1
  10
  11 #if 0
  12 static void st_OpeningGapCount( double *ogcp, int clus, char **seq, double *eff, int len )
  13 {
  14         int i, j, gc, gb;
  15         double feff;
  16
  17         for( i=0; i<len; i++ ) ogcp[i] = 0.0;
  18         for( j=0; j<clus; j++ )
  19         {
  20                 feff = (double)eff[j];
  21                 gc = 0;
  22                 for( i=0; i<len; i++ )
  23                 {
  24                         gb = gc;
  25                         gc = ( seq[j][i] == '-' );
  26                         {
  27                                 if( !gb *  gc ) ogcp[i] += feff;
  28                         }
  29                 }
  30         }
  31 }
  32
  33 static void st_FinalGapCount( double *fgcp, int clus, char **seq, double *eff, int len )
  34 {
  35         int i, j, gc, gb;
  36         double feff;
  37
  38         for( i=0; i<len; i++ ) fgcp[i] = 0.0;
  39         for( j=0; j<clus; j++ )
  40         {
  41                 feff = (double)eff[j];
  42                 gc = ( seq[j][0] == '-' );
  43                 for( i=1; i<len+1; i++ )
  44                 {
  45                         gb = gc;
  46                         gc = ( seq[j][i] == '-' );
  47                         {
  48                                 if( gb * !gc ) fgcp[i-1] += feff;
  49                         }
  50                 }
  51         }
  52 }
  53 #endif
  54
  55
  56
  57
  58
  59 static TLS int impalloclen = 0;
  60 static TLS double **impmtx = NULL;
  61 double part_imp_match_out_sc( int i1, int j1 )
  62 {
  63 //      fprintf( stderr, "impalloclen = %d\n", impalloclen );
  64 //      fprintf( stderr, "i1,j1=%d,%d -> impmtx=%f\n", i1, j1, impmtx[i1][j1] );
  65         return( impmtx[i1][j1] );
  66 #if 0
  67         if( i1 == l1 || j1 == l2 ) return( 0.0 );
  68         return( impmtx[i1+start1][j1+start2] );
  69 #endif
  70 }
  71 static void part_imp_match_out_vead_gapmap( double *imp, int i1, int lgth2, int start2, int *gapmap2 )
  72 {
  73 #if FASTMACHCALC
  74         double *pt = imp;
  75         int *gapmappt = gapmap2;
  76         while( lgth2-- )
  77                 *pt++ += impmtx[i1][start2+*gapmappt++];
  78 #else
  79         int j;
  80         for( j=0; j<lgth2; j++ )
  81         {
  82                 imp[j] += impmtx[i1][start2+gapmap2[j]];
  83         }
  84 #endif
  85 }
  86
  87 static void part_imp_match_out_vead_tate_gapmap( double *imp, int j1, int lgth1, int start1, int *gapmap1 )
  88 {
  89 #if FASTMACHCALC
  90         double *pt = imp;
  91         int *gapmappt = gapmap1;
  92         while( lgth1-- )
  93                 *pt++ = impmtx[start1+*gapmappt++][j1];
  94 #else
  95         int i;
  96         for( i=0; i<lgth1; i++ )
  97         {
  98                 imp[i] += impmtx[start1+gapmap1[i]][j1];
  99         }
 100 #endif
 101 }
 102
 103 #if 1
 104 void part_imp_match_init_strict( double *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, double *eff1_kozo, double *eff2_kozo, LocalHom ***localhom, char *swaplist, int forscore, int *orinum1, int *orinum2 )
 105 {
 106 //      int i, j, k1, k2, tmpint, start1, start2, end1, end2;
 107 //      double effij;
 108 //      double effij_kozo;
 109 //      double effijx;
 110 //      char *pt, *pt1, *pt2;
 111 //      static TLS char *nocount1 = NULL;
 112 //      static TLS char *nocount2 = NULL;
 113 //      LocalHom *tmpptr;
 114
 115         if( seq1 == NULL )
 116         {
 117                 if( impmtx ) FreeFloatMtx( impmtx );
 118                 impmtx = NULL;
 119 //              if( nocount1 ) free( nocount1 );
 120 //              nocount1 = NULL;
 121 //              if( nocount2 ) free( nocount2 );
 122 //              nocount2 = NULL;
 123
 124                 return;
 125         }
 126
 127         if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
 128         {
 129                 if( impmtx ) FreeFloatMtx( impmtx );
 130 //              if( nocount1 ) free( nocount1 );
 131 //              if( nocount2 ) free( nocount2 );
 132                 impalloclen = MAX( lgth1, lgth2 ) + 2;
 133                 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
 134 //              nocount1 = AllocateCharVec( impalloclen );
 135 //              nocount2 = AllocateCharVec( impalloclen );
 136         }
 137
 138         fillimp( impmtx, imp, clus1, clus2, lgth1, lgth2, seq1, seq2, eff1, eff2, eff1_kozo, eff2_kozo, localhom, swaplist, forscore, orinum1, orinum2 );
 139 }
 140 #else
 141 #endif
 142
 143
 144 void part_imp_rna( int nseq1, int nseq2, char **seq1, char **seq2, double *eff1, double *eff2, RNApair ***grouprna1, RNApair ***grouprna2, int *gapmap1, int *gapmap2, RNApair *additionalpair )
 145 {
 146         foldrna( nseq1, nseq2, seq1, seq2, eff1, eff2, grouprna1, grouprna2, impmtx, gapmap1, gapmap2, additionalpair );
 147 }
 148
 149
 150
 151
 152 static void match_calc( double *match, double **cpmx1, double **cpmx2, int i1, int lgth2, double **doublework, int **intwork, int initialize )
 153 {
 154 #if FASTMATCHCALC
 155         int j, l;
 156 //      double scarr[26];
 157         double **cpmxpd = doublework;
 158         int **cpmxpdn = intwork;
 159         double *matchpt, *cpmxpdpt, **cpmxpdptpt;
 160         int *cpmxpdnpt, **cpmxpdnptpt;
 161         double *scarr;
 162         scarr = calloc( nalphabets, sizeof( double ) );
 163         if( initialize )
 164         {
 165                 int count = 0;
 166                 for( j=0; j<lgth2; j++ )
 167                 {
 168                         count = 0;
 169                         for( l=0; l<nalphabets; l++ )
 170                         {
 171                                 if( cpmx2[l][j] )
 172                                 {
 173                                         cpmxpd[j][count] = cpmx2[l][j];
 174                                         cpmxpdn[j][count] = l;
 175                                         count++;
 176                                 }
 177                         }
 178                         cpmxpdn[j][count] = -1;
 179                 }
 180         }
 181
 182         {
 183                 for( l=0; l<nalphabets; l++ )
 184                 {
 185                         scarr[l] = 0.0;
 186                         for( j=0; j<nalphabets; j++ )
 187                                 scarr[l] += n_dis_consweight_multi[j][l] * cpmx1[j][i1];
 188 //                              scarr[l] += n_dis[j][l] * cpmx1[j][i1];
 189                 }
 190                 matchpt = match;
 191                 cpmxpdnptpt = cpmxpdn;
 192                 cpmxpdptpt = cpmxpd;
 193                 while( lgth2-- )
 194                 {
 195                         *matchpt = 0.0;
 196                         cpmxpdnpt = *cpmxpdnptpt++;
 197                         cpmxpdpt = *cpmxpdptpt++;
 198                         while( *cpmxpdnpt>-1 )
 199                                 *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;
 200                         matchpt++;
 201                 }
 202         }
 203         free( scarr );
 204 #else
 205         int j, k, l;
 206 //      double scarr[26];
 207         double **cpmxpd = doublework;
 208         int **cpmxpdn = intwork;
 209         double *scarr;
 210         scarr = calloc( nalphabets, sizeof( double ) );
 211         // simple
 212         if( initialize )
 213         {
 214                 int count = 0;
 215                 for( j=0; j<lgth2; j++ )
 216                 {
 217                         count = 0;
 218                         for( l=0; l<nalphabets; l++ )
 219                         {
 220                                 if( cpmx2[l][j] )
 221                                 {
 222                                         cpmxpd[count][j] = cpmx2[l][j];
 223                                         cpmxpdn[count][j] = l;
 224                                         count++;
 225                                 }
 226                         }
 227                         cpmxpdn[count][j] = -1;
 228                 }
 229         }
 230         for( l=0; l<nalphabets; l++ )
 231         {
 232                 scarr[l] = 0.0;
 233                 for( k=0; k<nalphabets; k++ )
 234                         scarr[l] += n_dis_consweight_multi[k][l] * cpmx1[k][i1];
 235 //                      scarr[l] += n_dis[k][l] * cpmx1[k][i1];
 236         }
 237         for( j=0; j<lgth2; j++ )
 238         {
 239                 match[j] = 0.0;
 240                 for( k=0; cpmxpdn[k][j]>-1; k++ )
 241                         match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
 242         }
 243 #endif
 244 }
 245
 246
 247 static void fillzero( double *s, int l )
 248 {
 249         while( l-- ) *s++ = 0.0;
 250 }
 251
 252
 253 static void match_calc_del( int **which, double ***matrices, double *match, int n1, char **seq1, double *eff1, int n2, char **seq2, double *eff2, int i1, int lgth2, int mid, int nmask, int *mask1, int *mask2 )
 254 {
 255 // osoi!
 256         int i, j, k, m;
 257         int c1, c2;
 258 //      fprintf( stderr, "\nmatch_calc_dynamicmtx... %d", i1 );
 259 //      fprintf( stderr, "\nseq1[0]=%s\n", seq1[0] );
 260 //      fprintf( stderr, "\nseq2[0]=%s\n", seq2[0] );
 261 //      for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
 262 //      {
 263 //              if( flip ) reporterr( "in match_calc_slow, which[%d][%d] = %d\n", j, i, which[j][i] );
 264 //              else       reporterr( "in match_calc_slow, which[%d][%d] = %d\n", i, j, which[i][j] );
 265 //      }
 266         for( k=0; k<lgth2; k++ )
 267         {
 268                 for( m=0; m<nmask; m++ )
 269                 {
 270                         i = mask1[m];
 271                         j = mask2[m];
 272 //                      reporterr( "Deleting %d-%d (c=%d)\n", i, j, mid );
 273 //                      if( k==0 ) fprintf( stderr, "pairoffset[%d][%d] = %f\n", i, j, po );
 274                         c1 = amino_n[(int)seq1[i][i1]];
 275                         c2 = amino_n[(int)seq2[j][k]];
 276 //                      reporterr( "k=%d, c1=%d, c2=%d, seq1[i][i1]=%c, seq2[%d][%d]=%c\n", k, c1, c2, seq1[i][i1], j, k, seq2[j][k] );
 277                         if( seq1[i][i1] == '-' || seq2[j][k] == '-' ) continue;
 278                         if( c1 < 0 || c2 < 0 ) continue;
 279 //                      fprintf( stderr, "c1=%d, c2=%d\n", c1, c2 );
 280 //                      fprintf( stderr, "match[k] = %f -> ", match[k], mid );
 281                         match[k] -= matrices[mid][c1][c2] * eff1[i] * eff2[j];
 282 //                      fprintf( stderr, "match[k] = %f (mid=%d)\n", match[k], mid );
 283                 }
 284         }
 285 //      fprintf( stderr, "done\n" );
 286         return;
 287 }
 288
 289 static void match_calc_add( double **scoreingmtx, double *match, double **cpmx1, double **cpmx2, int i1, int lgth2, double **doublework, int **intwork, int initialize )
 290 {
 291 #if FASTMATCHCALC
 292 //      fprintf( stderr, "\nmatch_calc... %d", i1 );
 293         int j, l;
 294 //      double scarr[26];
 295         double **cpmxpd = doublework;
 296         int **cpmxpdn = intwork;
 297         double *matchpt, *cpmxpdpt, **cpmxpdptpt;
 298         int *cpmxpdnpt, **cpmxpdnptpt;
 299         double *scarr;
 300         scarr = calloc( nalphabets, sizeof( double ) );
 301         if( initialize )
 302         {
 303                 int count = 0;
 304                 for( j=0; j<lgth2; j++ )
 305                 {
 306                         count = 0;
 307                         for( l=0; l<nalphabets; l++ )
 308                         {
 309                                 if( cpmx2[l][j] )
 310                                 {
 311                                         cpmxpd[j][count] = cpmx2[l][j];
 312                                         cpmxpdn[j][count] = l;
 313                                         count++;
 314                                 }
 315                         }
 316                         cpmxpdn[j][count] = -1;
 317                 }
 318         }
 319
 320         {
 321                 for( l=0; l<nalphabets; l++ )
 322                 {
 323                         scarr[l] = 0.0;
 324                         for( j=0; j<nalphabets; j++ )
 325 //                              scarr[l] += n_dis[j][l] * cpmx1[j][i1];
 326 //                              scarr[l] += n_dis_consweight_multi[j][l] * cpmx1[j][i1];
 327                                 scarr[l] += scoreingmtx[j][l] * cpmx1[j][i1];
 328                 }
 329                 matchpt = match;
 330                 cpmxpdnptpt = cpmxpdn;
 331                 cpmxpdptpt = cpmxpd;
 332                 while( lgth2-- )
 333                 {
 334 //                      *matchpt = 0.0;
 335                         cpmxpdnpt = *cpmxpdnptpt++;
 336                         cpmxpdpt = *cpmxpdptpt++;
 337                         while( *cpmxpdnpt>-1 )
 338                                 *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;
 339                         matchpt++;
 340                 }
 341         }
 342         free( scarr );
 343 //      fprintf( stderr, "done\n" );
 344 #else
 345         int j, k, l;
 346 //      double scarr[26];
 347         double **cpmxpd = doublework;
 348         int **cpmxpdn = intwork;
 349         double *scarr;
 350         scarr = calloc( nalphabets, sizeof( double ) );
 351 // simple
 352         if( initialize )
 353         {
 354                 int count = 0;
 355                 for( j=0; j<lgth2; j++ )
 356                 {
 357                         count = 0;
 358                         for( l=0; l<nalphabets; l++ )
 359                         {
 360                                 if( cpmx2[l][j] )
 361                                 {
 362                                         cpmxpd[count][j] = cpmx2[l][j];
 363                                         cpmxpdn[count][j] = l;
 364                                         count++;
 365                                 }
 366                         }
 367                         cpmxpdn[count][j] = -1;
 368                 }
 369         }
 370         for( l=0; l<nalphabets; l++ )
 371         {
 372                 scarr[l] = 0.0;
 373                 for( k=0; k<nalphabets; k++ )
 374 //                      scarr[l] += n_dis[k][l] * cpmx1[k][i1];
 375 //                      scarr[l] += n_dis_consweight_multi[k][l] * cpmx1[k][i1];
 376                         scarr[l] += scoreingmtx[k][l] * cpmx1[k][i1];
 377         }
 378         for( j=0; j<lgth2; j++ )
 379         {
 380                 match[j] = 0.0;
 381                 for( k=0; cpmxpdn[k][j]>-1; k++ )
 382                         match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
 383         }
 384         free( scarr );
 385 #endif
 386 }
 387
 388 static void Atracking_localhom( double *impwmpt, double *lasthorizontalw, double *lastverticalw,
 389                                                 char **seq1, char **seq2,
 390                         char **mseq1, char **mseq2,
 391                         int **ijp, int icyc, int jcyc,
 392                                                 int start1, int end1, int start2, int end2,
 393                                                 int *gapmap1, int *gapmap2,
 394                                                 int *warpis, int *warpjs, int warpbase )
 395 {
 396         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
 397 //      char gap[] = "-";
 398         char *gap;
 399         double wm;
 400         gap = newgapstr;
 401         lgth1 = strlen( seq1[0] );
 402         lgth2 = strlen( seq2[0] );
 403
 404 #if 0
 405         for( i=0; i<lgth1; i++ )
 406         {
 407                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 408         }
 409 #endif
 410
 411         if( outgap == 1 )
 412                 ;
 413         else
 414         {
 415                 wm = lastverticalw[0];
 416                 for( i=0; i<lgth1; i++ )
 417                 {
 418                         if( lastverticalw[i] >= wm )
 419                         {
 420                                 wm = lastverticalw[i];
 421                                 iin = i; jin = lgth2-1;
 422                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 423                         }
 424                 }
 425                 for( j=0; j<lgth2; j++ )
 426                 {
 427                         if( lasthorizontalw[j] >= wm )
 428                         {
 429                                 wm = lasthorizontalw[j];
 430                                 iin = lgth1-1; jin = j;
 431                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 432                         }
 433                 }
 434         }
 435
 436     for( i=0; i<lgth1+1; i++ )
 437     {
 438         ijp[i][0] = i + 1;
 439     }
 440     for( j=0; j<lgth2+1; j++ )
 441     {
 442         ijp[0][j] = -( j + 1 );
 443     }
 444
 445         for( i=0; i<icyc; i++ )
 446         {
 447                 mseq1[i] += lgth1+lgth2;
 448                 *mseq1[i] = 0;
 449         }
 450         for( j=0; j<jcyc; j++ )
 451         {
 452                 mseq2[j] += lgth1+lgth2;
 453                 *mseq2[j] = 0;
 454         }
 455         iin = lgth1; jin = lgth2;
 456         *impwmpt = 0.0;
 457         limk = lgth1+lgth2+1;
 458         for( k=0; k<limk; k++ )
 459         {
 460
 461                 if( ijp[iin][jin] >= warpbase )
 462                 {
 463                         ifi = warpis[ijp[iin][jin]-warpbase];
 464                         jfi = warpjs[ijp[iin][jin]-warpbase];
 465                 }
 466                 else if( ijp[iin][jin] < 0 )
 467                 {
 468                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 469                 }
 470                 else if( ijp[iin][jin] > 0 )
 471                 {
 472                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 473                 }
 474                 else
 475                 {
 476                         ifi = iin-1; jfi = jin-1;
 477                 }
 478                 if( ifi == -warpbase && jfi == -warpbase )
 479                 {
 480                         l = iin;
 481                         while( --l >= 0 )
 482                         {
 483                                 for( i=0; i<icyc; i++ )
 484                                         *--mseq1[i] = seq1[i][l];
 485                                 for( j=0; j<jcyc; j++ )
 486                                         *--mseq2[j] = *gap;
 487                                 k++;
 488                         }
 489                         l= jin;
 490                         while( --l >= 0 )
 491                         {
 492                                 for( i=0; i<icyc; i++ )
 493                                         *--mseq1[i] = *gap;
 494                                 for( j=0; j<jcyc; j++ )
 495                                         *--mseq2[j] = seq2[j][l];
 496                                 k++;
 497                         }
 498                         break;
 499                 }
 500                 else
 501                 {
 502                         l = iin - ifi;
 503                         while( --l )
 504                         {
 505                                 for( i=0; i<icyc; i++ )
 506                                         *--mseq1[i] = seq1[i][ifi+l];
 507                                 for( j=0; j<jcyc; j++ )
 508                                         *--mseq2[j] = *gap;
 509                                 k++;
 510                         }
 511                         l= jin - jfi;
 512                         while( --l )
 513                         {
 514                                 for( i=0; i<icyc; i++ )
 515                                         *--mseq1[i] = *gap;
 516                                 for( j=0; j<jcyc; j++ )
 517                                         *--mseq2[j] = seq2[j][jfi+l];
 518                                 k++;
 519                         }
 520                 }
 521                 if( iin != lgth1 && jin != lgth2 ) // ??
 522                 {
 523                         *impwmpt += (double)part_imp_match_out_sc( gapmap1[iin]+start1, gapmap2[jin]+start2 );
 524 //                      fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
 525                 }
 526                 if( iin <= 0 || jin <= 0 ) break;
 527                 for( i=0; i<icyc; i++ )
 528                         *--mseq1[i] = seq1[i][ifi];
 529                 for( j=0; j<jcyc; j++ )
 530                         *--mseq2[j] = seq2[j][jfi];
 531                 k++;
 532                 iin = ifi; jin = jfi;
 533         }
 534 }
 535
 536 static double Atracking( double *lasthorizontalw, double *lastverticalw,
 537                                                 char **seq1, char **seq2,
 538                         char **mseq1, char **mseq2,
 539                         int **ijp, int icyc, int jcyc,
 540                                                 int *warpis, int *warpjs, int warpbase )
 541 {
 542         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, lastk, limk;
 543 //      char gap[] = "-";
 544         char *gap;
 545         gap = newgapstr;
 546         double wm = 0.0;
 547         lgth1 = strlen( seq1[0] );
 548         lgth2 = strlen( seq2[0] );
 549
 550 #if 0
 551         for( i=0; i<lgth1; i++ )
 552         {
 553                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 554         }
 555 #endif
 556
 557         if( outgap == 1 )
 558                 ;
 559         else
 560         {
 561                 wm = lastverticalw[0];
 562                 for( i=0; i<lgth1; i++ )
 563                 {
 564                         if( lastverticalw[i] >= wm )
 565                         {
 566                                 wm = lastverticalw[i];
 567                                 iin = i; jin = lgth2-1;
 568                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 569                         }
 570                 }
 571                 for( j=0; j<lgth2; j++ )
 572                 {
 573                         if( lasthorizontalw[j] >= wm )
 574                         {
 575                                 wm = lasthorizontalw[j];
 576                                 iin = lgth1-1; jin = j;
 577                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 578                         }
 579                 }
 580         }
 581
 582     for( i=0; i<lgth1+1; i++ )
 583     {
 584         ijp[i][0] = i + 1;
 585     }
 586     for( j=0; j<lgth2+1; j++ )
 587     {
 588         ijp[0][j] = -( j + 1 );
 589     }
 590
 591         for( i=0; i<icyc; i++ )
 592         {
 593                 mseq1[i] += lgth1+lgth2;
 594                 *mseq1[i] = 0;
 595         }
 596         for( j=0; j<jcyc; j++ )
 597         {
 598                 mseq2[j] += lgth1+lgth2;
 599                 *mseq2[j] = 0;
 600         }
 601         iin = lgth1; jin = lgth2;
 602         lastk = lgth1+lgth2;
 603         limk = lgth1+lgth2+1;
 604         for( k=0; k<limk; k++ )
 605         {
 606                 if( ijp[iin][jin] >= warpbase )
 607                 {
 608                         ifi = warpis[ijp[iin][jin]-warpbase];
 609                         jfi = warpjs[ijp[iin][jin]-warpbase];
 610                 }
 611                 else if( ijp[iin][jin] < 0 )
 612                 {
 613                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 614                 }
 615                 else if( ijp[iin][jin] > 0 )
 616                 {
 617                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 618                 }
 619                 else
 620                 {
 621                         ifi = iin-1; jfi = jin-1;
 622                 }
 623                 if( ifi == -warpbase && jfi == -warpbase )
 624                 {
 625                         l = iin;
 626                         while( --l >= 0 )
 627                         {
 628                                 for( i=0; i<icyc; i++ )
 629                                         *--mseq1[i] = seq1[i][l];
 630                                 for( j=0; j<jcyc; j++ )
 631                                         *--mseq2[j] = *gap;
 632                                 k++;
 633                         }
 634                         l= jin;
 635                         while( --l >= 0 )
 636                         {
 637                                 for( i=0; i<icyc; i++ )
 638                                         *--mseq1[i] = *gap;
 639                                 for( j=0; j<jcyc; j++ )
 640                                         *--mseq2[j] = seq2[j][l];
 641                                 k++;
 642                         }
 643                         break;
 644                 }
 645                 else
 646                 {
 647                         l = iin - ifi;
 648                         while( --l )
 649                         {
 650                                 for( i=0; i<icyc; i++ )
 651                                         *--mseq1[i] = seq1[i][ifi+l];
 652                                 for( j=0; j<jcyc; j++ )
 653                                         *--mseq2[j] = *gap;
 654                                 k++;
 655                         }
 656                         l= jin - jfi;
 657                         while( --l )
 658                         {
 659                                 for( i=0; i<icyc; i++ )
 660                                         *--mseq1[i] = *gap;
 661                                 for( j=0; j<jcyc; j++ )
 662                                         *--mseq2[j] = seq2[j][jfi+l];
 663                                 k++;
 664                         }
 665                 }
 666
 667                 if( iin <= 0 || jin <= 0 ) break;
 668                 for( i=0; i<icyc; i++ )
 669                         *--mseq1[i] = seq1[i][ifi];
 670                 for( j=0; j<jcyc; j++ )
 671                         *--mseq2[j] = seq2[j][jfi];
 672                 k++;
 673                 iin = ifi; jin = jfi;
 674         }
 675         return( 0.0 );
 676 }
 677
 678 double partA__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, double *impmatch, int start1, int end1, int start2, int end2, int *gapmap1, int *gapmap2, char *sgap1, char *sgap2, char *egap1, char *egap2, int *chudanpt, int chudanref, int *chudanres )
 679 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 680 {
 681 //      int k;
 682         register int i, j;
 683         int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
 684         int lgth1, lgth2;
 685         int resultlen;
 686         double wm = 0.0;   /* int ?????? */
 687         double g;
 688         double *currentw, *previousw;
 689 #if 1
 690         double *wtmp;
 691         int *ijppt;
 692         double *mjpt, *prept, *curpt;
 693         int *mpjpt;
 694 #endif
 695         static TLS double mi, *m;
 696         static TLS int **ijp;
 697         static TLS int mpi, *mp;
 698         static TLS double *w1, *w2;
 699         static TLS double *match;
 700         static TLS double *initverticalw;    /* kufuu sureba iranai */
 701         static TLS double *lastverticalw;    /* kufuu sureba iranai */
 702         static TLS char **mseq1;
 703         static TLS char **mseq2;
 704         static TLS char **mseq;
 705         static TLS double *ogcp1;
 706         static TLS double *ogcp2;
 707         static TLS double *fgcp1;
 708         static TLS double *fgcp2;
 709         static TLS double **cpmx1;
 710         static TLS double **cpmx2;
 711         static TLS double *gapfreq1;
 712         static TLS double *gapfreq2;
 713         static TLS int **intwork;
 714         static TLS double **doublework;
 715         static TLS int orlgth1 = 0, orlgth2 = 0;
 716         double fpenalty = (double)penalty;
 717         double fpenalty_shift = (double)penalty_shift;
 718 #if USE_PENALTY_EX
 719         double fpenalty_ex = (double)penalty_ex;
 720 #endif
 721         double *fgcp2pt;
 722         double *ogcp2pt;
 723         double fgcp1va;
 724         double ogcp1va;
 725         double *gf2pt;
 726         double *gf2ptpre;
 727         double gf1va;
 728         double gf1vapre;
 729         double headgapfreq1;
 730         double headgapfreq2;
 731
 732         int *warpis = NULL;
 733         int *warpjs = NULL;
 734         int *warpi = NULL;
 735         int *warpj = NULL;
 736         int *prevwarpi = NULL;
 737         int *prevwarpj = NULL;
 738         double *wmrecords = NULL;
 739         double *prevwmrecords = NULL;
 740         int warpn = 0;
 741         int warpbase;
 742         double curm = 0.0;
 743         double *wmrecordspt, *wmrecords1pt, *prevwmrecordspt;
 744         int *warpipt, *warpjpt;
 745
 746
 747         if( seq1 == NULL )
 748         {
 749                 if( orlgth1 )
 750                 {
 751 //                      fprintf( stderr, "## Freeing local arrays in A__align\n" );
 752                         orlgth1 = 0;
 753                         orlgth2 = 0;
 754
 755                         part_imp_match_init_strict( NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, NULL );
 756
 757                         free( mseq1 );
 758                         free( mseq2 );
 759                         FreeFloatVec( w1 );
 760                         FreeFloatVec( w2 );
 761                         FreeFloatVec( match );
 762                         FreeFloatVec( initverticalw );
 763                         FreeFloatVec( lastverticalw );
 764
 765                         FreeFloatVec( m );
 766                         FreeIntVec( mp );
 767
 768                         FreeCharMtx( mseq );
 769
 770                         FreeFloatVec( ogcp1 );
 771                         FreeFloatVec( ogcp2 );
 772                         FreeFloatVec( fgcp1 );
 773                         FreeFloatVec( fgcp2 );
 774
 775
 776                         FreeFloatMtx( cpmx1 );
 777                         FreeFloatMtx( cpmx2 );
 778
 779                         FreeFloatVec( gapfreq1 );
 780                         FreeFloatVec( gapfreq2 );
 781
 782                         FreeFloatMtx( doublework );
 783                         FreeIntMtx( intwork );
 784
 785                 }
 786                 else
 787                 {
 788 //                      fprintf( stderr, "## Not allocated\n" );
 789                 }
 790                 return( 0.0 );
 791         }
 792 //      fprintf( stderr, "IN partA__align\n" );
 793
 794
 795         lgth1 = strlen( seq1[0] );
 796         lgth2 = strlen( seq2[0] );
 797 #if 1
 798 //      if( lgth1 == 0 ) fprintf( stderr, "WARNING: lgth1=0 in partA__align\n" );
 799 //      if( lgth2 == 0 ) fprintf( stderr, "WARNING: lgth2=0 in partA__align\n" );
 800
 801         if( lgth1 == 0 && lgth2 == 0 )
 802                 return( 0.0 );
 803
 804         if( lgth1 == 0 )
 805         {
 806                 for( i=0; i<icyc; i++ )
 807                 {
 808                         j = lgth2;
 809                         seq1[i][j] = 0;
 810                         while( j ) seq1[i][--j] = *newgapstr;
 811 //                      fprintf( stderr, "seq1[i] = %s\n", seq1[i] );
 812                 }
 813                 return( 0.0 );
 814         }
 815
 816         if( lgth2 == 0 )
 817         {
 818                 for( i=0; i<jcyc; i++ )
 819                 {
 820                         j = lgth1;
 821                         seq2[i][j] = 0;
 822                         while( j ) seq2[i][--j] = *newgapstr;
 823 //                      fprintf( stderr, "seq2[i] = %s\n", seq2[i] );
 824                 }
 825                 return( 0.0 );
 826         }
 827 #endif
 828
 829         warpbase = lgth1 + lgth2;
 830         warpis = NULL;
 831         warpjs = NULL;
 832         warpn = 0;
 833
 834
 835
 836
 837         if( trywarp )
 838         {
 839                 if( outgap == 0 )
 840                 {
 841                         fprintf( stderr, "At present, outgap must be 1.\n" );
 842                         exit( 1 );
 843                 }
 844                 wmrecords = AllocateFloatVec( lgth2+1 );
 845                 warpi = AllocateIntVec( lgth2+1 );
 846                 warpj = AllocateIntVec( lgth2+1 );
 847                 prevwmrecords = AllocateFloatVec( lgth2+1 );
 848                 prevwarpi = AllocateIntVec( lgth2+1 );
 849                 prevwarpj = AllocateIntVec( lgth2+1 );
 850                 for( i=0; i<lgth2+1; i++ ) wmrecords[i] = 0.0;
 851                 for( i=0; i<lgth2+1; i++ ) prevwmrecords[i] = 0.0;
 852                 for( i=0; i<lgth2+1; i++ ) prevwarpi[i] = -warpbase;
 853                 for( i=0; i<lgth2+1; i++ ) prevwarpj[i] = -warpbase;
 854                 for( i=0; i<lgth2+1; i++ ) warpi[i] = -warpbase;
 855                 for( i=0; i<lgth2+1; i++ ) warpj[i] = -warpbase;
 856         }
 857
 858 #if 0
 859         fprintf( stderr, "eff in SA+++align\n" );
 860         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
 861 #endif
 862         if( orlgth1 == 0 )
 863         {
 864                 mseq1 = AllocateCharMtx( njob, 0 );
 865                 mseq2 = AllocateCharMtx( njob, 0 );
 866         }
 867
 868
 869
 870
 871         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 872         {
 873                 int ll1, ll2;
 874
 875                 if( orlgth1 > 0 && orlgth2 > 0 )
 876                 {
 877                         FreeFloatVec( w1 );
 878                         FreeFloatVec( w2 );
 879                         FreeFloatVec( match );
 880                         FreeFloatVec( initverticalw );
 881                         FreeFloatVec( lastverticalw );
 882
 883                         FreeFloatVec( m );
 884                         FreeIntVec( mp );
 885
 886                         FreeCharMtx( mseq );
 887
 888                         FreeFloatVec( ogcp1 );
 889                         FreeFloatVec( ogcp2 );
 890                         FreeFloatVec( fgcp1 );
 891                         FreeFloatVec( fgcp2 );
 892
 893
 894                         FreeFloatMtx( cpmx1 );
 895                         FreeFloatMtx( cpmx2 );
 896
 897                         FreeFloatVec( gapfreq1 );
 898                         FreeFloatVec( gapfreq2 );
 899
 900                         FreeFloatMtx( doublework );
 901                         FreeIntMtx( intwork );
 902                 }
 903
 904                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 905                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 906
 907 #if DEBUG
 908                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 909 #endif
 910
 911                 w1 = AllocateFloatVec( ll2+2 );
 912                 w2 = AllocateFloatVec( ll2+2 );
 913                 match = AllocateFloatVec( ll2+2 );
 914
 915                 initverticalw = AllocateFloatVec( ll1+2 );
 916                 lastverticalw = AllocateFloatVec( ll1+2 );
 917
 918                 m = AllocateFloatVec( ll2+2 );
 919                 mp = AllocateIntVec( ll2+2 );
 920
 921                 mseq = AllocateCharMtx( njob, ll1+ll2 );
 922
 923                 ogcp1 = AllocateFloatVec( ll1+2 );
 924                 ogcp2 = AllocateFloatVec( ll2+2 );
 925                 fgcp1 = AllocateFloatVec( ll1+2 );
 926                 fgcp2 = AllocateFloatVec( ll2+2 );
 927
 928                 cpmx1 = AllocateFloatMtx( nalphabets, ll1+2 );
 929                 cpmx2 = AllocateFloatMtx( nalphabets, ll2+2 );
 930
 931                 gapfreq1 = AllocateFloatVec( ll1+2 );
 932                 gapfreq2 = AllocateFloatVec( ll2+2 );
 933
 934 #if FASTMATCHCALC
 935                 doublework = AllocateFloatMtx( MAX( ll1, ll2 )+2, nalphabets );
 936                 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, nalphabets );
 937 #else
 938                 doublework = AllocateFloatMtx( nalphabets, MAX( ll1, ll2 )+2 );
 939                 intwork = AllocateIntMtx( nalphabets, MAX( ll1, ll2 )+2 );
 940 #endif
 941
 942 #if DEBUG
 943                 fprintf( stderr, "succeeded\n" );
 944 #endif
 945
 946                 orlgth1 = ll1 - 100;
 947                 orlgth2 = ll2 - 100;
 948         }
 949
 950
 951         for( i=0; i<icyc; i++ ) mseq1[i] = mseq[i];
 952         for( j=0; j<jcyc; j++ ) mseq2[j] = mseq[icyc+j];
 953
 954
 955         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
 956         {
 957                 int ll1, ll2;
 958
 959                 if( commonAlloc1 && commonAlloc2 )
 960                 {
 961                         FreeIntMtx( commonIP );
 962                 }
 963
 964                 ll1 = MAX( orlgth1, commonAlloc1 );
 965                 ll2 = MAX( orlgth2, commonAlloc2 );
 966
 967 #if DEBUG
 968                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
 969 #endif
 970
 971                 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
 972
 973 #if DEBUG
 974                 fprintf( stderr, "succeeded\n\n" );
 975 #endif
 976
 977                 commonAlloc1 = ll1;
 978                 commonAlloc2 = ll2;
 979         }
 980         ijp = commonIP;
 981
 982         cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
 983         cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
 984
 985         if( sgap1 )
 986         {
 987                 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
 988                 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
 989                 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap1 );
 990                 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
 991                 outgapcount( &headgapfreq1, icyc, sgap1, eff1 );
 992                 outgapcount( &headgapfreq2, jcyc, sgap2, eff2 );
 993                 outgapcount( gapfreq1+lgth1, icyc, egap1, eff1 );
 994                 outgapcount( gapfreq2+lgth2, jcyc, egap2, eff2 );
 995         }
 996         else
 997         {
 998                 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
 999                 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1000                 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1001                 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1002                 headgapfreq1 = 0.0;
1003                 headgapfreq2 = 0.0;
1004                 gapfreq1[lgth1] = 0.0;
1005                 gapfreq2[lgth2] = 0.0;
1006         }
1007
1008         if( legacygapcost == 0 )
1009         {
1010                 gapcountf( gapfreq1, seq1, icyc, eff1, lgth1 );
1011                 gapcountf( gapfreq2, seq2, jcyc, eff2, lgth2 );
1012                 for( i=0; i<lgth1+1; i++ ) gapfreq1[i] = 1.0 - gapfreq1[i];
1013                 for( i=0; i<lgth2+1; i++ ) gapfreq2[i] = 1.0 - gapfreq2[i];
1014                 headgapfreq1 = 1.0 - headgapfreq1;
1015                 headgapfreq2 = 1.0 - headgapfreq2;
1016         }
1017         else
1018         {
1019                 for( i=0; i<lgth1+1; i++ ) gapfreq1[i] = 1.0;
1020                 for( i=0; i<lgth2+1; i++ ) gapfreq2[i] = 1.0;
1021                 headgapfreq1 = 1.0;
1022                 headgapfreq2 = 1.0;
1023         }
1024
1025         for( i=0; i<lgth1; i++ )
1026         {
1027                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty * ( gapfreq1[i] );
1028                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty * ( gapfreq1[i] );
1029         }
1030         for( i=0; i<lgth2; i++ )
1031         {
1032                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty * ( gapfreq2[i] );
1033                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty * ( gapfreq2[i] );
1034         }
1035 #if 0
1036         for( i=0; i<lgth1; i++ )
1037                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1038 #endif
1039
1040         currentw = w1;
1041         previousw = w2;
1042
1043
1044         match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, doublework, intwork, 1 );
1045         if( localhom )
1046                 part_imp_match_out_vead_tate_gapmap( initverticalw, gapmap2[0]+start2, lgth1, start1, gapmap1 );
1047
1048
1049         match_calc( currentw, cpmx1, cpmx2, 0, lgth2, doublework, intwork, 1 );
1050         if( localhom )
1051                 part_imp_match_out_vead_gapmap( currentw, gapmap1[0]+start1, lgth2, start2, gapmap2 );
1052 #if 0 // -> tbfast.c
1053         if( localhom )
1054                 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1055
1056 #endif
1057
1058         if( outgap == 1 )
1059         {
1060                 for( i=1; i<lgth1+1; i++ )
1061                 {
1062 //                      initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1063                         initverticalw[i] += ( ogcp1[0] * headgapfreq2 + fgcp1[i-1] * gapfreq2[0] ) ;
1064                 }
1065                 for( j=1; j<lgth2+1; j++ )
1066                 {
1067 //                      currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1068                         currentw[j] += ( ogcp2[0] * headgapfreq1 + fgcp2[j-1] * gapfreq1[0] ) ;
1069                 }
1070         }
1071 #if OUTGAP0TRY
1072         else
1073         {
1074                 for( j=1; j<lgth2+1; j++ )
1075                         currentw[j] -= offset * j / 2.0;
1076                 for( i=1; i<lgth1+1; i++ )
1077                         initverticalw[i] -= offset * i / 2.0;
1078         }
1079 #endif
1080
1081         for( j=1; j<lgth2+1; ++j )
1082         {
1083 //              m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1084                 m[j] = currentw[j-1] + ogcp1[1] * gapfreq2[j-1]; mp[j] = 0;;
1085         }
1086
1087         lastverticalw[0] = currentw[lgth2-1];
1088
1089         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1090         lastj = lgth2+1;
1091
1092 #if XXXXXXX
1093 fprintf( stderr, "currentw = \n" );
1094 for( i=0; i<lgth1+1; i++ )
1095 {
1096         fprintf( stderr, "%5.2f ", currentw[i] );
1097 }
1098 fprintf( stderr, "\n" );
1099 fprintf( stderr, "initverticalw = \n" );
1100 for( i=0; i<lgth2+1; i++ )
1101 {
1102         fprintf( stderr, "%5.2f ", initverticalw[i] );
1103 }
1104 fprintf( stderr, "\n" );
1105 fprintf( stderr, "fcgp\n" );
1106 for( i=0; i<lgth1; i++ )
1107         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1108 for( i=0; i<lgth2; i++ )
1109         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1110 #endif
1111
1112         for( i=1; i<lasti; i++ )
1113         {
1114
1115 #ifdef enablemultithread
1116 //              fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
1117                 if( chudanpt && *chudanpt != chudanref )
1118                 {
1119 //                      fprintf( stderr, "\n\n## CHUUDAN!!! i\n" );
1120                         *chudanres = 1;
1121                         return( -1.0 );
1122                 }
1123 #endif
1124
1125                 wtmp = previousw;
1126                 previousw = currentw;
1127                 currentw = wtmp;
1128
1129                 previousw[0] = initverticalw[i-1];
1130
1131                 match_calc( currentw, cpmx1, cpmx2, i, lgth2, doublework, intwork, 0 );
1132 #if 0
1133 fprintf( stderr, "\n" );
1134 fprintf( stderr, "i=%d\n", i );
1135 fprintf( stderr, "currentw before imp = \n" );
1136 for( j=0; j<lgth2; j++ )
1137 {
1138         fprintf( stderr, "%5.2f ", currentw[j] );
1139 }
1140 fprintf( stderr, "\n" );
1141 #endif
1142                 if( localhom )
1143                 {
1144 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1145 //                      imp_match_out_vead( currentw, i, lgth2 );
1146                         part_imp_match_out_vead_gapmap( currentw, gapmap1[i]+start1, lgth2, start2, gapmap2 );
1147                 }
1148 #if 0
1149 fprintf( stderr, "specificity = 0\n" );
1150 fprintf( stderr, "i=%d\n", i );
1151 fprintf( stderr, "currentw = \n" );
1152 for( j=0; j<lgth2; j++ )
1153 {
1154         fprintf( stderr, "%5.2f ", currentw[j] );
1155 }
1156 fprintf( stderr, "\n" );
1157 #endif
1158                 currentw[0] = initverticalw[i];
1159
1160
1161 //              mi = previousw[0] + ogcp2[1]; mpi = 0;
1162                 mi = previousw[0] + ogcp2[1] * gapfreq1[i-1]; mpi=0;
1163
1164                 ijppt = ijp[i] + 1;
1165                 mjpt = m + 1;
1166                 prept = previousw;
1167                 curpt = currentw + 1;
1168                 mpjpt = mp + 1;
1169                 fgcp2pt = fgcp2;
1170                 ogcp2pt = ogcp2+1;
1171                 fgcp1va = fgcp1[i-1];
1172                 ogcp1va = ogcp1[i];
1173                 gf1va = gapfreq1[i];
1174                 gf1vapre = gapfreq1[i-1];
1175                 gf2pt = gapfreq2+1;
1176                 gf2ptpre = gapfreq2;
1177
1178                 if( trywarp )
1179                 {
1180                         prevwmrecordspt = prevwmrecords;
1181                         wmrecordspt = wmrecords+1;
1182                         wmrecords1pt = wmrecords;
1183                         warpipt = warpi + 1;
1184                         warpjpt = warpj + 1;
1185                 }
1186
1187                 for( j=1; j<lastj; j++ )
1188                 {
1189 #ifdef xxxenablemultithread
1190 //                      fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
1191                         if( chudanpt && *chudanpt != chudanref )
1192                         {
1193 //                              fprintf( stderr, "\n\n## CHUUDAN!!! j\n" );
1194                                 *chudanres = 1;
1195                                 return( -1.0 );
1196                         }
1197 #endif
1198                         wm = *prept;
1199                         *ijppt = 0;
1200
1201 #if 0
1202                         fprintf( stderr, "%5.0f->", wm );
1203 #endif
1204 //                      g = mi + *fgcp2pt * gapfreq1[i];
1205                         if( (g = mi + *fgcp2pt * gf1va) > wm )
1206                         {
1207                                 wm = g;
1208                                 *ijppt = -( j - mpi );
1209                         }
1210 //                      g = *prept + *ogcp2pt * gapfreq1[i-1];
1211                         if( (g = *prept + *ogcp2pt * gf1vapre) >= mi )
1212                         {
1213                                 mi = g;
1214                                 mpi = j-1;
1215                         }
1216 #if USE_PENALTY_EX
1217                         mi += fpenalty_ex;
1218 #endif
1219
1220 //                      g = *mjpt + fgcp1va * gapfreq2[j];
1221                         if( (g = *mjpt + fgcp1va * *gf2pt) > wm )
1222                         {
1223                                 wm = g;
1224                                 *ijppt = +( i - *mpjpt );
1225                         }
1226 //                      g = *prept + ogcp1va * gapfreq2[j-1];
1227                         if( (g = *prept + ogcp1va * *gf2ptpre) >= *mjpt )
1228                         {
1229                                 *mjpt = g;
1230                                 *mpjpt = i-1;
1231                         }
1232 #if USE_PENALTY_EX
1233                         m[j] += fpenalty_ex;
1234 #endif
1235                         if( trywarp )
1236                         {
1237 #if USE_PENALTY_EX
1238                                 if( ( g=*prevwmrecordspt++ + fpenalty_shift + fpenalty_ex * ( i - prevwarpi[j-1] + j - prevwarpj[j-1] ) ) > wm ) // naka ha osokute kamawanai
1239 #else
1240                                 if( ( g=*prevwmrecordspt++ + fpenalty_shift ) > wm ) // naka ha osokute kamawanai
1241 #endif
1242                                 {
1243 //                                      fprintf( stderr, "WARP in partA__align\n" );
1244                                         if( warpn && prevwarpi[j-1] == warpis[warpn-1] && prevwarpj[j-1] == warpjs[warpn-1] )
1245                                         {
1246                                                 *ijppt = warpbase + warpn - 1;
1247                                         }
1248                                         else
1249                                         {
1250                                                 *ijppt = warpbase + warpn;
1251                                                 warpis = realloc( warpis, sizeof(int) * ( warpn+1 ) );
1252                                                 warpjs = realloc( warpjs, sizeof(int) * ( warpn+1 ) );
1253                                                 warpis[warpn] = prevwarpi[j-1];
1254                                                 warpjs[warpn] = prevwarpj[j-1];
1255                                                 warpn++;
1256                                         }
1257                                         wm = g;
1258                                 }
1259                                 curm = *curpt + wm;
1260                                 if( *wmrecords1pt > *wmrecordspt )
1261                                 {
1262                                         *wmrecordspt = *wmrecords1pt;
1263                                         *warpipt  = *(warpipt-1);
1264                                         *warpjpt  = *(warpjpt-1);
1265                                 }
1266                                 if( curm > *wmrecordspt )
1267                                 {
1268                                         *wmrecordspt = curm;
1269                                         *warpipt = i;
1270                                         *warpjpt = j;
1271                                 }
1272                                 wmrecordspt++;
1273                                 wmrecords1pt++;
1274                                 warpipt++;
1275                                 warpjpt++;
1276                         }
1277
1278 #if 0
1279                         fprintf( stderr, "%5.0f ", wm );
1280 #endif
1281                         *curpt += wm;
1282                         ijppt++;
1283                         mjpt++;
1284                         prept++;
1285                         mpjpt++;
1286                         curpt++;
1287                         fgcp2pt++;
1288                         ogcp2pt++;
1289                         gf2ptpre++;
1290                         gf2pt++;
1291
1292                 }
1293                 lastverticalw[i] = currentw[lgth2-1];
1294
1295                 if( trywarp )
1296                 {
1297                         fltncpy( prevwmrecords, wmrecords, lastj );
1298                         intncpy( prevwarpi, warpi, lastj );
1299                         intncpy( prevwarpj, warpj, lastj );
1300                 }
1301         }
1302         if( trywarp )
1303         {
1304 //              fprintf( stderr, "wm = %f\n", wm );
1305 //              fprintf( stderr, "warpn = %d\n", warpn );
1306                 free( wmrecords );
1307                 free( prevwmrecords );
1308                 free( warpi );
1309                 free( warpj );
1310                 free( prevwarpi );
1311                 free( prevwarpj );
1312         }
1313
1314 #if OUTGAP0TRY
1315         if( !outgap )
1316         {
1317                 for( j=1; j<lgth2+1; j++ )
1318                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1319                 for( i=1; i<lgth1+1; i++ )
1320                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1321         }
1322 #endif
1323
1324         /*
1325         fprintf( stderr, "\n" );
1326         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1327         fprintf( stderr, "#####\n" );
1328         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1329         fprintf( stderr, "====>" );
1330         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1331         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1332         */
1333         if( localhom )
1334         {
1335                 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, start1, end1, start2, end2, gapmap1, gapmap2, warpis, warpjs, warpbase );
1336         }
1337         else
1338                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, warpis, warpjs, warpbase );
1339
1340         if( warpis ) free( warpis );
1341         if( warpjs ) free( warpjs );
1342
1343 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1344
1345         resultlen = strlen( mseq1[0] );
1346         if( alloclen < resultlen || resultlen > N )
1347         {
1348                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1349                 ErrorExit( "LENGTH OVER!\n" );
1350         }
1351
1352
1353         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1354         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1355         /*
1356         fprintf( stderr, "\n" );
1357         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1358         fprintf( stderr, "#####\n" );
1359         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1360         */
1361
1362
1363         return( wm );
1364 }
1365
1366 double partA__align_variousdist( int **which, double ***matrices, double **n_dynamicmtx, char **seq1, char **seq2, double *eff1, double *eff2, double **eff1s, double **eff2s, int icyc, int jcyc, int alloclen, LocalHom ***localhom, double *impmatch, int start1, int end1, int start2, int end2, int *gapmap1, int *gapmap2, char *sgap1, char *sgap2, char *egap1, char *egap2, int *chudanpt, int chudanref, int *chudanres )
1367 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1368 {
1369 //      int k;
1370         register int i, j, c;
1371         int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
1372         int lgth1, lgth2;
1373         int resultlen;
1374         double wm = 0.0;   /* int ?????? */
1375         double g;
1376         double *currentw, *previousw;
1377 #if 1
1378         double *wtmp;
1379         int *ijppt;
1380         double *mjpt, *prept, *curpt;
1381         int *mpjpt;
1382 #endif
1383         static TLS double mi, *m;
1384         static TLS int **ijp;
1385         static TLS int mpi, *mp;
1386         static TLS double *w1, *w2;
1387         static TLS double *match;
1388         static TLS double *initverticalw;    /* kufuu sureba iranai */
1389         static TLS double *lastverticalw;    /* kufuu sureba iranai */
1390         static TLS char **mseq1;
1391         static TLS char **mseq2;
1392         static TLS char **mseq;
1393         static TLS double *ogcp1;
1394         static TLS double *ogcp2;
1395         static TLS double *fgcp1;
1396         static TLS double *fgcp2;
1397         static TLS double ***cpmx1s;
1398         static TLS double ***cpmx2s;
1399         static TLS double *gapfreq1;
1400         static TLS double *gapfreq2;
1401         static TLS int ***intwork;
1402         static TLS double ***doublework;
1403         static TLS int orlgth1 = 0, orlgth2 = 0;
1404         double fpenalty = (double)penalty;
1405         double fpenalty_shift = (double)penalty_shift;
1406 #if USE_PENALTY_EX
1407         double fpenalty_ex = (double)penalty_ex;
1408 #endif
1409         double *fgcp2pt;
1410         double *ogcp2pt;
1411         double fgcp1va;
1412         double ogcp1va;
1413         double *gf2pt;
1414         double *gf2ptpre;
1415         double gf1va;
1416         double gf1vapre;
1417         double headgapfreq1;
1418         double headgapfreq2;
1419
1420         int *warpis = NULL;
1421         int *warpjs = NULL;
1422         int *warpi = NULL;
1423         int *warpj = NULL;
1424         int *prevwarpi = NULL;
1425         int *prevwarpj = NULL;
1426         double *wmrecords = NULL;
1427         double *prevwmrecords = NULL;
1428         int warpn = 0;
1429         int warpbase;
1430         double curm = 0.0;
1431         double *wmrecordspt, *wmrecords1pt, *prevwmrecordspt;
1432         int *warpipt, *warpjpt;
1433         int *nmask, **masklist1, **masklist2;
1434
1435
1436         if( seq1 == NULL )
1437         {
1438                 if( orlgth1 )
1439                 {
1440 //                      fprintf( stderr, "## Freeing local arrays in A__align\n" );
1441                         orlgth1 = 0;
1442                         orlgth2 = 0;
1443
1444                         part_imp_match_init_strict( NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, NULL );
1445
1446                         free( mseq1 );
1447                         free( mseq2 );
1448                         FreeFloatVec( w1 );
1449                         FreeFloatVec( w2 );
1450                         FreeFloatVec( match );
1451                         FreeFloatVec( initverticalw );
1452                         FreeFloatVec( lastverticalw );
1453
1454                         FreeFloatVec( m );
1455                         FreeIntVec( mp );
1456
1457                         FreeCharMtx( mseq );
1458
1459                         FreeFloatVec( ogcp1 );
1460                         FreeFloatVec( ogcp2 );
1461                         FreeFloatVec( fgcp1 );
1462                         FreeFloatVec( fgcp2 );
1463
1464
1465                         FreeFloatCub( cpmx1s );
1466                         FreeFloatCub( cpmx2s );
1467
1468                         FreeFloatVec( gapfreq1 );
1469                         FreeFloatVec( gapfreq2 );
1470
1471                         FreeFloatCub( doublework );
1472                         FreeIntCub( intwork );
1473
1474                 }
1475                 else
1476                 {
1477 //                      fprintf( stderr, "## Not allocated\n" );
1478                 }
1479                 return( 0.0 );
1480         }
1481
1482         masklist1 = AllocateIntMtx( maxdistclass, 0 );
1483         masklist2 = AllocateIntMtx( maxdistclass, 0 );
1484         nmask = calloc( maxdistclass, sizeof( int ) );
1485
1486         for( c=0; c<maxdistclass; c++ )
1487         {
1488                 for( i=0; i<icyc; i++ ) for( j=0; j<jcyc; j++ )
1489                 {
1490                         if( eff1s[c][i] * eff2s[c][j] != 0.0 )
1491                         {
1492
1493                                 if( c != which[i][j] )
1494                                 {
1495                                         masklist1[c] = realloc( masklist1[c], sizeof( int ) * (nmask[c]+1) );
1496                                         masklist2[c] = realloc( masklist2[c], sizeof( int ) * (nmask[c]+1) );
1497
1498                                         masklist1[c][nmask[c]] = i;
1499                                         masklist2[c][nmask[c]] = j;
1500                                         nmask[c]++;
1501                                 }
1502                         }
1503                 }
1504         }
1505         for( c=0; c<maxdistclass; c++ ) if( nmask[c] ) break;
1506         if( c<maxdistclass ) reporterr( "Found a complex grouping. This step may be a bit slow.\n" );
1507
1508         lgth1 = strlen( seq1[0] );
1509         lgth2 = strlen( seq2[0] );
1510 #if 1
1511 //      if( lgth1 == 0 ) fprintf( stderr, "WARNING: lgth1=0 in partA__align\n" );
1512 //      if( lgth2 == 0 ) fprintf( stderr, "WARNING: lgth2=0 in partA__align\n" );
1513
1514         if( lgth1 == 0 && lgth2 == 0 )
1515                 return( 0.0 );
1516
1517         if( lgth1 == 0 )
1518         {
1519                 for( i=0; i<icyc; i++ )
1520                 {
1521                         j = lgth2;
1522                         seq1[i][j] = 0;
1523                         while( j ) seq1[i][--j] = *newgapstr;
1524 //                      fprintf( stderr, "seq1[i] = %s\n", seq1[i] );
1525                 }
1526                 return( 0.0 );
1527         }
1528
1529         if( lgth2 == 0 )
1530         {
1531                 for( i=0; i<jcyc; i++ )
1532                 {
1533                         j = lgth1;
1534                         seq2[i][j] = 0;
1535                         while( j ) seq2[i][--j] = *newgapstr;
1536 //                      fprintf( stderr, "seq2[i] = %s\n", seq2[i] );
1537                 }
1538                 return( 0.0 );
1539         }
1540 #endif
1541
1542         warpbase = lgth1 + lgth2;
1543         warpis = NULL;
1544         warpjs = NULL;
1545         warpn = 0;
1546
1547
1548         if( trywarp )
1549         {
1550 //              fprintf( stderr, "IN partA__align_variousdist\n" );
1551                 if( outgap == 0 )
1552                 {
1553                         fprintf( stderr, "At present, outgap must be 1 to allow shift.\n" );
1554                         exit( 1 );
1555                 }
1556                 wmrecords = AllocateFloatVec( lgth2+1 );
1557                 warpi = AllocateIntVec( lgth2+1 );
1558                 warpj = AllocateIntVec( lgth2+1 );
1559                 prevwmrecords = AllocateFloatVec( lgth2+1 );
1560                 prevwarpi = AllocateIntVec( lgth2+1 );
1561                 prevwarpj = AllocateIntVec( lgth2+1 );
1562                 for( i=0; i<lgth2+1; i++ ) wmrecords[i] = 0.0;
1563                 for( i=0; i<lgth2+1; i++ ) prevwmrecords[i] = 0.0;
1564                 for( i=0; i<lgth2+1; i++ ) prevwarpi[i] = -warpbase;
1565                 for( i=0; i<lgth2+1; i++ ) prevwarpj[i] = -warpbase;
1566                 for( i=0; i<lgth2+1; i++ ) warpi[i] = -warpbase;
1567                 for( i=0; i<lgth2+1; i++ ) warpj[i] = -warpbase;
1568         }
1569
1570 #if 0
1571         fprintf( stderr, "eff in SA+++align\n" );
1572         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1573 #endif
1574         if( orlgth1 == 0 )
1575         {
1576                 mseq1 = AllocateCharMtx( njob, 0 );
1577                 mseq2 = AllocateCharMtx( njob, 0 );
1578         }
1579
1580
1581
1582
1583         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
1584         {
1585                 int ll1, ll2;
1586
1587                 if( orlgth1 > 0 && orlgth2 > 0 )
1588                 {
1589                         FreeFloatVec( w1 );
1590                         FreeFloatVec( w2 );
1591                         FreeFloatVec( match );
1592                         FreeFloatVec( initverticalw );
1593                         FreeFloatVec( lastverticalw );
1594
1595                         FreeFloatVec( m );
1596                         FreeIntVec( mp );
1597
1598                         FreeCharMtx( mseq );
1599
1600                         FreeFloatVec( ogcp1 );
1601                         FreeFloatVec( ogcp2 );
1602                         FreeFloatVec( fgcp1 );
1603                         FreeFloatVec( fgcp2 );
1604
1605
1606                         FreeFloatCub( cpmx1s );
1607                         FreeFloatCub( cpmx2s );
1608
1609                         FreeFloatVec( gapfreq1 );
1610                         FreeFloatVec( gapfreq2 );
1611
1612                         FreeFloatCub( doublework );
1613                         FreeIntCub( intwork );
1614                 }
1615
1616                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
1617                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
1618
1619 #if DEBUG
1620                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1621 #endif
1622
1623                 w1 = AllocateFloatVec( ll2+2 );
1624                 w2 = AllocateFloatVec( ll2+2 );
1625                 match = AllocateFloatVec( ll2+2 );
1626
1627                 initverticalw = AllocateFloatVec( ll1+2 );
1628                 lastverticalw = AllocateFloatVec( ll1+2 );
1629
1630                 m = AllocateFloatVec( ll2+2 );
1631                 mp = AllocateIntVec( ll2+2 );
1632
1633                 mseq = AllocateCharMtx( njob, ll1+ll2 );
1634
1635                 ogcp1 = AllocateFloatVec( ll1+2 );
1636                 ogcp2 = AllocateFloatVec( ll2+2 );
1637                 fgcp1 = AllocateFloatVec( ll1+2 );
1638                 fgcp2 = AllocateFloatVec( ll2+2 );
1639
1640                 cpmx1s = AllocateFloatCub( maxdistclass, nalphabets, ll1+2 );
1641                 cpmx2s = AllocateFloatCub( maxdistclass, nalphabets, ll2+2 );
1642
1643                 gapfreq1 = AllocateFloatVec( ll1+2 );
1644                 gapfreq2 = AllocateFloatVec( ll2+2 );
1645
1646                 doublework = AllocateFloatCub( maxdistclass, MAX( ll1, ll2 )+2, nalphabets );
1647                 intwork = AllocateIntCub( maxdistclass, MAX( ll1, ll2 )+2, nalphabets );
1648
1649 #if DEBUG
1650                 fprintf( stderr, "succeeded\n" );
1651 #endif
1652
1653                 orlgth1 = ll1 - 100;
1654                 orlgth2 = ll2 - 100;
1655         }
1656
1657
1658         for( i=0; i<icyc; i++ ) mseq1[i] = mseq[i];
1659         for( j=0; j<jcyc; j++ ) mseq2[j] = mseq[icyc+j];
1660
1661
1662         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1663         {
1664                 int ll1, ll2;
1665
1666                 if( commonAlloc1 && commonAlloc2 )
1667                 {
1668                         FreeIntMtx( commonIP );
1669                 }
1670
1671                 ll1 = MAX( orlgth1, commonAlloc1 );
1672                 ll2 = MAX( orlgth2, commonAlloc2 );
1673
1674 #if DEBUG
1675                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1676 #endif
1677
1678                 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
1679
1680 #if DEBUG
1681                 fprintf( stderr, "succeeded\n\n" );
1682 #endif
1683
1684                 commonAlloc1 = ll1;
1685                 commonAlloc2 = ll2;
1686         }
1687         ijp = commonIP;
1688
1689 //      cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1690 //      cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1691         for( c=0; c<maxdistclass; c++ )
1692         {
1693                 cpmx_calc_new( seq1, cpmx1s[c], eff1s[c], lgth1, icyc );
1694                 cpmx_calc_new( seq2, cpmx2s[c], eff2s[c], lgth2, jcyc );
1695         }
1696
1697         if( sgap1 )
1698         {
1699                 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
1700                 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
1701                 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap1 );
1702                 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
1703                 outgapcount( &headgapfreq1, icyc, sgap1, eff1 );
1704                 outgapcount( &headgapfreq2, jcyc, sgap2, eff2 );
1705                 outgapcount( gapfreq1+lgth1, icyc, egap1, eff1 );
1706                 outgapcount( gapfreq2+lgth2, jcyc, egap2, eff2 );
1707         }
1708         else
1709         {
1710                 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1711                 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1712                 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1713                 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1714                 headgapfreq1 = 0.0;
1715                 headgapfreq2 = 0.0;
1716                 gapfreq1[lgth1] = 0.0;
1717                 gapfreq2[lgth2] = 0.0;
1718         }
1719
1720         if( legacygapcost == 0 )
1721         {
1722                 gapcountf( gapfreq1, seq1, icyc, eff1, lgth1 );
1723                 gapcountf( gapfreq2, seq2, jcyc, eff2, lgth2 );
1724                 for( i=0; i<lgth1+1; i++ ) gapfreq1[i] = 1.0 - gapfreq1[i];
1725                 for( i=0; i<lgth2+1; i++ ) gapfreq2[i] = 1.0 - gapfreq2[i];
1726                 headgapfreq1 = 1.0 - headgapfreq1;
1727                 headgapfreq2 = 1.0 - headgapfreq2;
1728         }
1729         else
1730         {
1731                 for( i=0; i<lgth1+1; i++ ) gapfreq1[i] = 1.0;
1732                 for( i=0; i<lgth2+1; i++ ) gapfreq2[i] = 1.0;
1733                 headgapfreq1 = 1.0;
1734                 headgapfreq2 = 1.0;
1735         }
1736
1737         for( i=0; i<lgth1; i++ )
1738         {
1739                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty * ( gapfreq1[i] );
1740                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty * ( gapfreq1[i] );
1741         }
1742         for( i=0; i<lgth2; i++ )
1743         {
1744                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty * ( gapfreq2[i] );
1745                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty * ( gapfreq2[i] );
1746         }
1747 #if 0
1748         for( i=0; i<lgth1; i++ )
1749                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1750 #endif
1751
1752         currentw = w1;
1753         previousw = w2;
1754
1755
1756 //      match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, doublework, intwork, 1 );
1757         fillzero( initverticalw, lgth1 );
1758         for( c=0; c<maxdistclass; c++ )
1759         {
1760                 match_calc_add( matrices[c], initverticalw, cpmx2s[c], cpmx1s[c], 0, lgth1, doublework[c], intwork[c], 1 );
1761                 if( nmask[c] ) match_calc_del( which, matrices, initverticalw, jcyc, seq2, eff2, icyc, seq1, eff1, 0, lgth1, c, nmask[c], masklist2[c], masklist1[c] );
1762         }
1763
1764         if( localhom )
1765                 part_imp_match_out_vead_tate_gapmap( initverticalw, gapmap2[0]+start2, lgth1, start1, gapmap1 );
1766
1767
1768 //      match_calc( currentw, cpmx1, cpmx2, 0, lgth2, doublework, intwork, 1 );
1769         fillzero( currentw, lgth2 );
1770         for( c=0; c<maxdistclass; c++ )
1771         {
1772                 match_calc_add( matrices[c], currentw, cpmx1s[c], cpmx2s[c], 0, lgth2, doublework[c], intwork[c], 1 );
1773                 if( nmask[c] ) match_calc_del( which, matrices, currentw, icyc, seq1, eff1, jcyc, seq2, eff2, 0, lgth2, c, nmask[c], masklist1[c], masklist2[c] );
1774         }
1775         if( localhom )
1776                 part_imp_match_out_vead_gapmap( currentw, gapmap1[0]+start1, lgth2, start2, gapmap2 );
1777 #if 0 // -> tbfast.c
1778         if( localhom )
1779                 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1780
1781 #endif
1782
1783         if( outgap == 1 )
1784         {
1785                 for( i=1; i<lgth1+1; i++ )
1786                 {
1787 //                      initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1788                         initverticalw[i] += ( ogcp1[0] * headgapfreq2 + fgcp1[i-1] * gapfreq2[0] ) ;
1789                 }
1790                 for( j=1; j<lgth2+1; j++ )
1791                 {
1792 //                      currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1793                         currentw[j] += ( ogcp2[0] * headgapfreq1 + fgcp2[j-1] * gapfreq1[0] ) ;
1794                 }
1795         }
1796 #if OUTGAP0TRY
1797         else
1798         {
1799                 for( j=1; j<lgth2+1; j++ )
1800                         currentw[j] -= offset * j / 2.0;
1801                 for( i=1; i<lgth1+1; i++ )
1802                         initverticalw[i] -= offset * i / 2.0;
1803         }
1804 #endif
1805
1806         for( j=1; j<lgth2+1; ++j )
1807         {
1808 //              m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1809                 m[j] = currentw[j-1] + ogcp1[1] * gapfreq2[j-1]; mp[j] = 0;;
1810         }
1811
1812         lastverticalw[0] = currentw[lgth2-1];
1813
1814         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1815         lastj = lgth2+1;
1816
1817 #if XXXXXXX
1818 fprintf( stderr, "currentw = \n" );
1819 for( i=0; i<lgth1+1; i++ )
1820 {
1821         fprintf( stderr, "%5.2f ", currentw[i] );
1822 }
1823 fprintf( stderr, "\n" );
1824 fprintf( stderr, "initverticalw = \n" );
1825 for( i=0; i<lgth2+1; i++ )
1826 {
1827         fprintf( stderr, "%5.2f ", initverticalw[i] );
1828 }
1829 fprintf( stderr, "\n" );
1830 fprintf( stderr, "fcgp\n" );
1831 for( i=0; i<lgth1; i++ )
1832         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1833 for( i=0; i<lgth2; i++ )
1834         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1835 #endif
1836
1837         for( i=1; i<lasti; i++ )
1838         {
1839
1840 #ifdef enablemultithread
1841 //              fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
1842                 if( chudanpt && *chudanpt != chudanref )
1843                 {
1844 //                      fprintf( stderr, "\n\n## CHUUDAN!!! i\n" );
1845                         *chudanres = 1;
1846                         if( masklist1 ) freeintmtx( masklist1, maxdistclass ); masklist1 = NULL;
1847                         if( masklist2 ) freeintmtx( masklist2, maxdistclass ); masklist2 = NULL;
1848                         if( nmask ) free( nmask ); nmask = NULL;
1849                         return( -1.0 );
1850                 }
1851 #endif
1852
1853                 wtmp = previousw;
1854                 previousw = currentw;
1855                 currentw = wtmp;
1856
1857                 previousw[0] = initverticalw[i-1];
1858
1859 //              match_calc( currentw, cpmx1, cpmx2, i, lgth2, doublework, intwork, 0 );
1860                 fillzero( currentw, lgth2 );
1861                 for( c=0; c<maxdistclass; c++ )
1862                 {
1863                         match_calc_add( matrices[c], currentw, cpmx1s[c], cpmx2s[c], i, lgth2, doublework[c], intwork[c], 0 );
1864                         if( nmask[c] ) match_calc_del( which, matrices, currentw, icyc, seq1, eff1, jcyc, seq2, eff2, i, lgth2, c, nmask[c], masklist1[c], masklist2[c] );
1865                 }
1866 #if 0
1867 fprintf( stderr, "\n" );
1868 fprintf( stderr, "i=%d\n", i );
1869 fprintf( stderr, "currentw before imp = \n" );
1870 for( j=0; j<lgth2; j++ )
1871 {
1872         fprintf( stderr, "%5.2f ", currentw[j] );
1873 }
1874 fprintf( stderr, "\n" );
1875 #endif
1876                 if( localhom )
1877                 {
1878 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1879 //                      imp_match_out_vead( currentw, i, lgth2 );
1880                         part_imp_match_out_vead_gapmap( currentw, gapmap1[i]+start1, lgth2, start2, gapmap2 );
1881                 }
1882 #if 0
1883 fprintf( stderr, "specificity = %f\n", specificityconsideration );
1884 fprintf( stderr, "i=%d\n", i );
1885 fprintf( stderr, "currentw = \n" );
1886 for( j=0; j<lgth2; j++ )
1887 {
1888         fprintf( stderr, "%5.2f ", currentw[j] );
1889 }
1890 fprintf( stderr, "\n" );
1891 #endif
1892                 currentw[0] = initverticalw[i];
1893
1894
1895 //              mi = previousw[0] + ogcp2[1]; mpi = 0;
1896                 mi = previousw[0] + ogcp2[1] * gapfreq1[i-1]; mpi=0;
1897
1898                 ijppt = ijp[i] + 1;
1899                 mjpt = m + 1;
1900                 prept = previousw;
1901                 curpt = currentw + 1;
1902                 mpjpt = mp + 1;
1903                 fgcp2pt = fgcp2;
1904                 ogcp2pt = ogcp2+1;
1905                 fgcp1va = fgcp1[i-1];
1906                 ogcp1va = ogcp1[i];
1907                 gf1va = gapfreq1[i];
1908                 gf1vapre = gapfreq1[i-1];
1909                 gf2pt = gapfreq2+1;
1910                 gf2ptpre = gapfreq2;
1911
1912                 if( trywarp )
1913                 {
1914                         prevwmrecordspt = prevwmrecords;
1915                         wmrecordspt = wmrecords+1;
1916                         wmrecords1pt = wmrecords;
1917                         warpipt = warpi + 1;
1918                         warpjpt = warpj + 1;
1919                 }
1920
1921                 for( j=1; j<lastj; j++ )
1922                 {
1923 #ifdef xxxenablemultithread
1924 //                      fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
1925                         if( chudanpt && *chudanpt != chudanref )
1926                         {
1927 //                              fprintf( stderr, "\n\n## CHUUDAN!!! j\n" );
1928                                 *chudanres = 1;
1929                                 if( masklist1 ) freeintmtx( masklist1, maxdistclass ); masklist1 = NULL;
1930                                 if( masklist2 ) freeintmtx( masklist2, maxdistclass ); masklist2 = NULL;
1931                                 if( nmask ) free( nmask ); nmask = NULL;
1932                                 return( -1.0 );
1933                         }
1934 #endif
1935                         wm = *prept;
1936                         *ijppt = 0;
1937
1938 #if 0
1939                         fprintf( stderr, "%5.0f->", wm );
1940 #endif
1941 //                      g = mi + *fgcp2pt * gapfreq1[i];
1942                         if( (g = mi + *fgcp2pt * gf1va) > wm )
1943                         {
1944                                 wm = g;
1945                                 *ijppt = -( j - mpi );
1946                         }
1947 //                      g = *prept + *ogcp2pt * gapfreq1[i-1];
1948                         if( (g = *prept + *ogcp2pt * gf1vapre) >= mi )
1949                         {
1950                                 mi = g;
1951                                 mpi = j-1;
1952                         }
1953 #if USE_PENALTY_EX
1954                         mi += fpenalty_ex;
1955 #endif
1956
1957 //                      g = *mjpt + fgcp1va * gapfreq2[j];
1958                         if( (g = *mjpt + fgcp1va * *gf2pt) > wm )
1959                         {
1960                                 wm = g;
1961                                 *ijppt = +( i - *mpjpt );
1962                         }
1963 //                      g = *prept + ogcp1va * gapfreq2[j-1];
1964                         if( (g = *prept + ogcp1va * *gf2ptpre) >= *mjpt )
1965                         {
1966                                 *mjpt = g;
1967                                 *mpjpt = i-1;
1968                         }
1969 #if USE_PENALTY_EX
1970                         m[j] += fpenalty_ex;
1971 #endif
1972                         if( trywarp )
1973                         {
1974 #if USE_PENALTY_EX
1975                                 if( ( g=*prevwmrecordspt++ + fpenalty_shift + fpenalty_ex * ( i - prevwarpi[j-1] + j - prevwarpj[j-1] ) ) > wm ) // naka ha osokute kamawanai
1976 #else
1977                                 if( ( g=*prevwmrecordspt++ + fpenalty_shift ) > wm ) // naka ha osokute kamawanai
1978 #endif
1979                                 {
1980                                         if( warpn && prevwarpi[j-1] == warpis[warpn-1] && prevwarpj[j-1] == warpjs[warpn-1] )
1981                                         {
1982                                                 *ijppt = warpbase + warpn - 1;
1983                                         }
1984                                         else
1985                                         {
1986                                                 *ijppt = warpbase + warpn;
1987                                                 warpis = realloc( warpis, sizeof(int) * ( warpn+1 ) );
1988                                                 warpjs = realloc( warpjs, sizeof(int) * ( warpn+1 ) );
1989                                                 warpis[warpn] = prevwarpi[j-1];
1990                                                 warpjs[warpn] = prevwarpj[j-1];
1991                                                 warpn++;
1992                                         }
1993                                         wm = g;
1994                                 }
1995
1996                                 curm = *curpt + wm;
1997                                 if( *wmrecords1pt > *wmrecordspt )
1998                                 {
1999                                         *wmrecordspt = *wmrecords1pt;
2000                                         *warpipt  = *(warpipt-1);
2001                                         *warpjpt  = *(warpjpt-1);
2002                                 }
2003                                 if( curm > *wmrecordspt )
2004                                 {
2005                                         *wmrecordspt = curm;
2006                                         *warpipt = i;
2007                                         *warpjpt = j;
2008                                 }
2009                                 wmrecordspt++;
2010                                 wmrecords1pt++;
2011                                 warpipt++;
2012                                 warpjpt++;
2013                         }
2014
2015 #if 0
2016                         fprintf( stderr, "%5.0f ", wm );
2017 #endif
2018                         *curpt += wm;
2019                         ijppt++;
2020                         mjpt++;
2021                         prept++;
2022                         mpjpt++;
2023                         curpt++;
2024                         fgcp2pt++;
2025                         ogcp2pt++;
2026                         gf2ptpre++;
2027                         gf2pt++;
2028
2029                 }
2030                 lastverticalw[i] = currentw[lgth2-1];
2031
2032                 if( trywarp )
2033                 {
2034                         fltncpy( prevwmrecords, wmrecords, lastj );
2035                         intncpy( prevwarpi, warpi, lastj );
2036                         intncpy( prevwarpj, warpj, lastj );
2037                 }
2038         }
2039         if( trywarp )
2040         {
2041 //              fprintf( stderr, "wm = %f\n", wm );
2042 //              fprintf( stderr, "warpn = %d\n", warpn );
2043                 free( wmrecords );
2044                 free( prevwmrecords );
2045                 free( warpi );
2046                 free( warpj );
2047                 free( prevwarpi );
2048                 free( prevwarpj );
2049         }
2050
2051 #if OUTGAP0TRY
2052         if( !outgap )
2053         {
2054                 for( j=1; j<lgth2+1; j++ )
2055                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
2056                 for( i=1; i<lgth1+1; i++ )
2057                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
2058         }
2059 #endif
2060
2061         /*
2062         fprintf( stderr, "\n" );
2063         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
2064         fprintf( stderr, "#####\n" );
2065         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
2066         fprintf( stderr, "====>" );
2067         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
2068         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
2069         */
2070         if( localhom )
2071         {
2072                 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, start1, end1, start2, end2, gapmap1, gapmap2, warpis, warpjs, warpbase );
2073         }
2074         else
2075                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, warpis, warpjs, warpbase );
2076
2077         if( warpis ) free( warpis );
2078         if( warpjs ) free( warpjs );
2079
2080 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
2081
2082         resultlen = strlen( mseq1[0] );
2083         if( alloclen < resultlen || resultlen > N )
2084         {
2085                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
2086                 ErrorExit( "LENGTH OVER!\n" );
2087         }
2088
2089
2090         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2091         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
2092         /*
2093         fprintf( stderr, "\n" );
2094         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
2095         fprintf( stderr, "#####\n" );
2096         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
2097         */
2098
2099         if( masklist1 ) freeintmtx( masklist1, maxdistclass ); masklist1 = NULL;
2100         if( masklist2 ) freeintmtx( masklist2, maxdistclass ); masklist2 = NULL;
2101         if( nmask ) free( nmask ); nmask = NULL;
2102
2103         return( wm );
2104 }