binaries/src/mafft/core/Salignmm.c

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