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