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