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         char *gap;
 471         float wm;
 472         gap = newgapstr;
 473         lgth1 = strlen( seq1[0] );
 474         lgth2 = strlen( seq2[0] );
 475
 476 #if 0
 477         for( i=0; i<lgth1; i++ )
 478         {
 479                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 480         }
 481 #endif
 482
 483         if( outgap == 1 )
 484                 ;
 485         else
 486         {
 487                 wm = lastverticalw[0];
 488                 for( i=0; i<lgth1; i++ )
 489                 {
 490                         if( lastverticalw[i] >= wm )
 491                         {
 492                                 wm = lastverticalw[i];
 493                                 iin = i; jin = lgth2-1;
 494                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 495                         }
 496                 }
 497                 for( j=0; j<lgth2; j++ )
 498                 {
 499                         if( lasthorizontalw[j] >= wm )
 500                         {
 501                                 wm = lasthorizontalw[j];
 502                                 iin = lgth1-1; jin = j;
 503                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 504                         }
 505                 }
 506         }
 507
 508     for( i=0; i<lgth1+1; i++ )
 509     {
 510         ijp[i][0] = i + 1;
 511     }
 512     for( j=0; j<lgth2+1; j++ )
 513     {
 514         ijp[0][j] = -( j + 1 );
 515     }
 516
 517         for( i=0; i<icyc; i++ )
 518         {
 519                 mseq1[i] += lgth1+lgth2;
 520                 *mseq1[i] = 0;
 521         }
 522         for( j=0; j<jcyc; j++ )
 523         {
 524                 mseq2[j] += lgth1+lgth2;
 525                 *mseq2[j] = 0;
 526         }
 527         iin = lgth1; jin = lgth2;
 528         *impwmpt = 0.0;
 529         for( k=0; k<=lgth1+lgth2; k++ )
 530         {
 531                 if( ijp[iin][jin] < 0 )
 532                 {
 533                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 534                 }
 535                 else if( ijp[iin][jin] > 0 )
 536                 {
 537                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 538                 }
 539                 else
 540                 {
 541                         ifi = iin-1; jfi = jin-1;
 542                 }
 543                 l = iin - ifi;
 544                 while( --l )
 545                 {
 546                         for( i=0; i<icyc; i++ )
 547                                 *--mseq1[i] = seq1[i][ifi+l];
 548                         for( j=0; j<jcyc; j++ )
 549                                 *--mseq2[j] = *gap;
 550                         k++;
 551                 }
 552                 l= jin - jfi;
 553                 while( --l )
 554                 {
 555                         for( i=0; i<icyc; i++ )
 556                                 *--mseq1[i] = *gap;
 557                         for( j=0; j<jcyc; j++ )
 558                                 *--mseq2[j] = seq2[j][jfi+l];
 559                         k++;
 560                 }
 561                 if( iin != lgth1 && jin != lgth2 ) // ??
 562                 {
 563                         *impwmpt += part_imp_match_out_scQ( gapmap1[iin]+start1, gapmap2[jin]+start2 );
 564 //                      fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
 565                 }
 566                 if( iin <= 0 || jin <= 0 ) break;
 567                 for( i=0; i<icyc; i++ )
 568                         *--mseq1[i] = seq1[i][ifi];
 569                 for( j=0; j<jcyc; j++ )
 570                         *--mseq2[j] = seq2[j][jfi];
 571                 k++;
 572                 iin = ifi; jin = jfi;
 573         }
 574 }
 575 static float Atracking( float *lasthorizontalw, float *lastverticalw,
 576                                                 char **seq1, char **seq2,
 577                         char **mseq1, char **mseq2,
 578                         float **cpmx1, float **cpmx2,
 579                         int **ijp, int icyc, int jcyc )
 580 {
 581         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, lastk;
 582 //      char gap[] = "-";
 583         char *gap;
 584         float wm = 0.0;
 585         gap = newgapstr;
 586         lgth1 = strlen( seq1[0] );
 587         lgth2 = strlen( seq2[0] );
 588
 589 #if 0
 590         for( i=0; i<lgth1; i++ )
 591         {
 592                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 593         }
 594 #endif
 595
 596         if( outgap == 1 )
 597                 ;
 598         else
 599         {
 600                 wm = lastverticalw[0];
 601                 for( i=0; i<lgth1; i++ )
 602                 {
 603                         if( lastverticalw[i] >= wm )
 604                         {
 605                                 wm = lastverticalw[i];
 606                                 iin = i; jin = lgth2-1;
 607                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 608                         }
 609                 }
 610                 for( j=0; j<lgth2; j++ )
 611                 {
 612                         if( lasthorizontalw[j] >= wm )
 613                         {
 614                                 wm = lasthorizontalw[j];
 615                                 iin = lgth1-1; jin = j;
 616                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 617                         }
 618                 }
 619         }
 620
 621     for( i=0; i<lgth1+1; i++ )
 622     {
 623         ijp[i][0] = i + 1;
 624     }
 625     for( j=0; j<lgth2+1; j++ )
 626     {
 627         ijp[0][j] = -( j + 1 );
 628     }
 629
 630         for( i=0; i<icyc; i++ )
 631         {
 632                 mseq1[i] += lgth1+lgth2;
 633                 *mseq1[i] = 0;
 634         }
 635         for( j=0; j<jcyc; j++ )
 636         {
 637                 mseq2[j] += lgth1+lgth2;
 638                 *mseq2[j] = 0;
 639         }
 640         iin = lgth1; jin = lgth2;
 641         lastk = lgth1+lgth2;
 642         for( k=0; k<=lastk; k++ )
 643         {
 644                 if( ijp[iin][jin] < 0 )
 645                 {
 646                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 647                 }
 648                 else if( ijp[iin][jin] > 0 )
 649                 {
 650                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 651                 }
 652                 else
 653                 {
 654                         ifi = iin-1; jfi = jin-1;
 655                 }
 656                 l = iin - ifi;
 657                 while( --l )
 658                 {
 659                         for( i=0; i<icyc; i++ )
 660                                 *--mseq1[i] = seq1[i][ifi+l];
 661                         for( j=0; j<jcyc; j++ )
 662                                 *--mseq2[j] = *gap;
 663                         k++;
 664                 }
 665                 l= jin - jfi;
 666                 while( --l )
 667                 {
 668                         for( i=0; i<icyc; i++ )
 669                                 *--mseq1[i] = *gap;
 670                         for( j=0; j<jcyc; j++ )
 671                                 *--mseq2[j] = seq2[j][jfi+l];
 672                         k++;
 673                 }
 674                 if( iin <= 0 || jin <= 0 ) break;
 675                 for( i=0; i<icyc; i++ )
 676                         *--mseq1[i] = seq1[i][ifi];
 677                 for( j=0; j<jcyc; j++ )
 678                         *--mseq2[j] = seq2[j][jfi];
 679                 k++;
 680                 iin = ifi; jin = jfi;
 681         }
 682         return( 0.0 );
 683 }
 684
 685 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 )
 686 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 687 {
 688 //      int k;
 689         register int i, j;
 690         int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
 691         int lgth1, lgth2;
 692         int resultlen;
 693         float wm = 0.0;   /* int ?????? */
 694         float g;
 695         float *currentw, *previousw;
 696 #if 1
 697         float *wtmp;
 698         int *ijppt;
 699         float *mjpt, *prept, *curpt;
 700         int *mpjpt;
 701 #endif
 702         static float mi, *m;
 703         static int **ijp;
 704         static int mpi, *mp;
 705         static float *w1, *w2;
 706         static float *match;
 707         static float *initverticalw;    /* kufuu sureba iranai */
 708         static float *lastverticalw;    /* kufuu sureba iranai */
 709         static char **mseq1;
 710         static char **mseq2;
 711         static char **mseq;
 712         static float *digf1;
 713         static float *digf2;
 714         static float *diaf1;
 715         static float *diaf2;
 716         static float *gapz1;
 717         static float *gapz2;
 718         static float *gapf1;
 719         static float *gapf2;
 720         static float *ogcp1g;
 721         static float *ogcp2g;
 722         static float *fgcp1g;
 723         static float *fgcp2g;
 724         static float *og_h_dg_n1_p;
 725         static float *og_h_dg_n2_p;
 726         static float *fg_h_dg_n1_p;
 727         static float *fg_h_dg_n2_p;
 728         static float *og_t_fg_h_dg_n1_p;
 729         static float *og_t_fg_h_dg_n2_p;
 730         static float *fg_t_og_h_dg_n1_p;
 731         static float *fg_t_og_h_dg_n2_p;
 732         static float *gapz_n1;
 733         static float *gapz_n2;
 734         static float **cpmx1;
 735         static float **cpmx2;
 736         static int **intwork;
 737         static float **floatwork;
 738         static int orlgth1 = 0, orlgth2 = 0;
 739         float fpenalty = (float)penalty;
 740 #if USE_PENALTY_EX
 741         float fpenalty_ex = (float)penalty_ex;
 742 #endif
 743         float tmppenal;
 744         float *fg_t_og_h_dg_n2_p_pt;
 745         float *og_t_fg_h_dg_n2_p_pt;
 746         float *og_h_dg_n2_p_pt;
 747         float *fg_h_dg_n2_p_pt;
 748         float *gapz_n2_pt0;
 749         float *gapz_n2_pt1;
 750         float *fgcp2pt;
 751         float *ogcp2pt;
 752         float fg_t_og_h_dg_n1_p_va;
 753         float og_t_fg_h_dg_n1_p_va;
 754         float og_h_dg_n1_p_va;
 755         float fg_h_dg_n1_p_va;
 756         float gapz_n1_va0;
 757         float gapz_n1_va1;
 758         float fgcp1va;
 759         float ogcp1va;
 760
 761
 762
 763 #if 0
 764         fprintf( stderr, "eff in SA+++align\n" );
 765         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
 766 #endif
 767         if( orlgth1 == 0 )
 768         {
 769                 mseq1 = AllocateCharMtx( njob, 0 );
 770                 mseq2 = AllocateCharMtx( njob, 0 );
 771         }
 772
 773
 774         lgth1 = strlen( seq1[0] );
 775         lgth2 = strlen( seq2[0] );
 776
 777         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 778         {
 779                 int ll1, ll2;
 780
 781                 if( orlgth1 > 0 && orlgth2 > 0 )
 782                 {
 783                         FreeFloatVec( w1 );
 784                         FreeFloatVec( w2 );
 785                         FreeFloatVec( match );
 786                         FreeFloatVec( initverticalw );
 787                         FreeFloatVec( lastverticalw );
 788
 789                         FreeFloatVec( m );
 790                         FreeIntVec( mp );
 791
 792                         FreeCharMtx( mseq );
 793
 794                         FreeFloatVec( digf1 );
 795                         FreeFloatVec( digf2 );
 796                         FreeFloatVec( diaf1 );
 797                         FreeFloatVec( diaf2 );
 798                         FreeFloatVec( gapz1 );
 799                         FreeFloatVec( gapz2 );
 800                         FreeFloatVec( gapf1 );
 801                         FreeFloatVec( gapf2 );
 802                         FreeFloatVec( ogcp1g );
 803                         FreeFloatVec( ogcp2g );
 804                         FreeFloatVec( fgcp1g );
 805                         FreeFloatVec( fgcp2g );
 806                         FreeFloatVec( og_h_dg_n1_p );
 807                         FreeFloatVec( og_h_dg_n2_p );
 808                         FreeFloatVec( fg_h_dg_n1_p );
 809                         FreeFloatVec( fg_h_dg_n2_p );
 810                         FreeFloatVec( og_t_fg_h_dg_n1_p );
 811                         FreeFloatVec( og_t_fg_h_dg_n2_p );
 812                         FreeFloatVec( fg_t_og_h_dg_n1_p );
 813                         FreeFloatVec( fg_t_og_h_dg_n2_p );
 814                         FreeFloatVec( gapz_n1 );
 815                         FreeFloatVec( gapz_n2 );
 816
 817                         FreeFloatMtx( cpmx1 );
 818                         FreeFloatMtx( cpmx2 );
 819
 820                         FreeFloatMtx( floatwork );
 821                         FreeIntMtx( intwork );
 822                 }
 823
 824                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 825                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 826
 827 #if DEBUG
 828                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 829 #endif
 830
 831                 w1 = AllocateFloatVec( ll2+2 );
 832                 w2 = AllocateFloatVec( ll2+2 );
 833                 match = AllocateFloatVec( ll2+2 );
 834
 835                 initverticalw = AllocateFloatVec( ll1+2 );
 836                 lastverticalw = AllocateFloatVec( ll1+2 );
 837
 838                 m = AllocateFloatVec( ll2+2 );
 839                 mp = AllocateIntVec( ll2+2 );
 840
 841                 mseq = AllocateCharMtx( njob, ll1+ll2 );
 842
 843                 digf1 = AllocateFloatVec( ll1+2 );
 844                 digf2 = AllocateFloatVec( ll2+2 );
 845                 diaf1 = AllocateFloatVec( ll1+2 );
 846                 diaf2 = AllocateFloatVec( ll2+2 );
 847                 gapz1 = AllocateFloatVec( ll1+2 );
 848                 gapz2 = AllocateFloatVec( ll2+2 );
 849                 gapf1 = AllocateFloatVec( ll1+2 );
 850                 gapf2 = AllocateFloatVec( ll2+2 );
 851                 ogcp1g = AllocateFloatVec( ll1+2 );
 852                 ogcp2g = AllocateFloatVec( ll2+2 );
 853                 fgcp1g = AllocateFloatVec( ll1+2 );
 854                 fgcp2g = AllocateFloatVec( ll2+2 );
 855                 og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
 856                 og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
 857                 fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
 858                 fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
 859                 og_t_fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
 860                 og_t_fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
 861                 fg_t_og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
 862                 fg_t_og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
 863                 gapz_n1 = AllocateFloatVec( ll1+2 );
 864                 gapz_n2 = AllocateFloatVec( ll2+2 );
 865
 866                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
 867                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
 868
 869 #if FASTMATCHCALC
 870                 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
 871                 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
 872 #else
 873                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
 874                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
 875 #endif
 876
 877 #if DEBUG
 878                 fprintf( stderr, "succeeded\n" );
 879 #endif
 880
 881                 orlgth1 = ll1 - 100;
 882                 orlgth2 = ll2 - 100;
 883         }
 884
 885
 886         for( i=0; i<icyc; i++ ) mseq1[i] = mseq[i];
 887         for( j=0; j<jcyc; j++ ) mseq2[j] = mseq[icyc+j];
 888
 889
 890         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
 891         {
 892                 int ll1, ll2;
 893
 894                 if( commonAlloc1 && commonAlloc2 )
 895                 {
 896                         FreeIntMtx( commonIP );
 897                 }
 898
 899                 ll1 = MAX( orlgth1, commonAlloc1 );
 900                 ll2 = MAX( orlgth2, commonAlloc2 );
 901
 902 #if DEBUG
 903                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
 904 #endif
 905
 906                 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
 907
 908 #if DEBUG
 909                 fprintf( stderr, "succeeded\n\n" );
 910 #endif
 911
 912                 commonAlloc1 = ll1;
 913                 commonAlloc2 = ll2;
 914         }
 915         ijp = commonIP;
 916
 917         cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
 918         cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
 919
 920         if( sgap1 )
 921         {
 922                 new_OpeningGapCount_zure( ogcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
 923                 new_OpeningGapCount_zure( ogcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
 924                 new_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
 925                 new_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
 926                 getdigapfreq_part( digf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
 927                 getdigapfreq_part( digf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
 928                 getdiaminofreq_part( diaf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
 929                 getdiaminofreq_part( diaf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
 930                 getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
 931                 getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
 932                 getgapfreq_zure_part( gapz1, icyc, seq1, eff1, lgth1, sgap1 );
 933                 getgapfreq_zure_part( gapz2, jcyc, seq2, eff2, lgth2, sgap1 );
 934         }
 935         else
 936         {
 937                 st_OpeningGapCount( ogcp1g, icyc, seq1, eff1, lgth1 );
 938                 st_OpeningGapCount( ogcp2g, jcyc, seq2, eff2, lgth2 );
 939                 st_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1 );
 940                 st_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2 );
 941                 getdigapfreq_st( digf1, icyc, seq1, eff1, lgth1 );
 942                 getdigapfreq_st( digf2, jcyc, seq2, eff2, lgth2 );
 943                 getdiaminofreq_x( diaf1, icyc, seq1, eff1, lgth1 );
 944                 getdiaminofreq_x( diaf2, jcyc, seq2, eff2, lgth2 );
 945                 getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
 946                 getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
 947                 getgapfreq_zure( gapz1, icyc, seq1, eff1, lgth1 );
 948                 getgapfreq_zure( gapz2, jcyc, seq2, eff2, lgth2 );
 949         }
 950
 951 #if 1
 952         lastj = lgth2+2;
 953         for( i=0; i<lastj; i++ )
 954         {
 955                 og_h_dg_n2_p[i] = ( 1.0-ogcp2g[i]-digf2[i] ) * fpenalty * 0.5;
 956                 fg_h_dg_n2_p[i] = ( 1.0-fgcp2g[i]-digf2[i] ) * fpenalty * 0.5;
 957                 og_t_fg_h_dg_n2_p[i] = (1.0-ogcp2g[i]+fgcp2g[i]-digf2[i]) * 0.5 * fpenalty;
 958                 fg_t_og_h_dg_n2_p[i] = (1.0-fgcp2g[i]+ogcp2g[i]-digf2[i]) * 0.5 * fpenalty;
 959                 gapz_n2[i] = (1.0-gapz2[i]);
 960         }
 961         lastj = lgth1+2;
 962         for( i=0; i<lastj; i++ )
 963         {
 964                 og_h_dg_n1_p[i] = ( 1.0-ogcp1g[i]-digf1[i] ) * fpenalty * 0.5;
 965                 fg_h_dg_n1_p[i] = ( 1.0-fgcp1g[i]-digf1[i] ) * fpenalty * 0.5;
 966                 og_t_fg_h_dg_n1_p[i] = (1.0-ogcp1g[i]+fgcp1g[i]-digf1[i]) * 0.5 * fpenalty;
 967                 fg_t_og_h_dg_n1_p[i] = (1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) * 0.5 * fpenalty;
 968                 gapz_n1[i] = (1.0-gapz1[i]);
 969         }
 970 #endif
 971
 972         currentw = w1;
 973         previousw = w2;
 974
 975
 976         match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
 977         if( localhom )
 978                 part_imp_match_out_vead_tate_gapmapQ( initverticalw, gapmap2[0]+start2, lgth1, start1, gapmap1 );
 979
 980
 981         match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
 982         if( localhom )
 983                 part_imp_match_out_vead_gapmapQ( currentw, gapmap1[0]+start1, lgth2, start2, gapmap2 );
 984 #if 0 // -> tbfast.c
 985         if( localhom )
 986                 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
 987
 988 #endif
 989
 990         if( outgap == 1 )
 991         {
 992                 g = 0.0;
 993
 994                 g += ogcp1g[0] * og_h_dg_n2_p[0];
 995
 996                 g += ogcp2g[0] * og_h_dg_n1_p[0];
 997
 998                 g += fgcp1g[0] * fg_h_dg_n2_p[0];
 999
1000                 g += fgcp2g[0] * fg_h_dg_n1_p[0];
1001
1002                 initverticalw[0] += g;
1003                 currentw[0] += g;
1004
1005                 for( i=1; i<lgth1+1; i++ )
1006                 {
1007                         tmppenal = gapz_n2[0]*og_t_fg_h_dg_n1_p[0];
1008                         initverticalw[i] += tmppenal;
1009
1010                         tmppenal = gapz_n2[1]*fg_t_og_h_dg_n1_p[i];
1011                         initverticalw[i] += tmppenal;
1012
1013                 }
1014                 for( j=1; j<lgth2+1; j++ )
1015                 {
1016                         tmppenal = gapz_n1[0]*og_t_fg_h_dg_n2_p[0];
1017                         currentw[j] += tmppenal;
1018
1019                         tmppenal = gapz_n1[1]*fg_t_og_h_dg_n2_p[j];
1020                         currentw[j] += tmppenal;
1021                 }
1022         }
1023 #if OUTGAP0TRY
1024         else
1025         {
1026                 for( j=1; j<lgth2+1; j++ )
1027                                 currentw[j] -= offset * j / 2.0;
1028                 for( i=1; i<lgth1+1; i++ )
1029                         initverticalw[i] -= offset * i / 2.0;
1030         }
1031 #endif
1032
1033         m[0] = 0.0; // iranai
1034         for( j=1; j<lgth2+1; ++j )
1035         {
1036                 mp[j] = 0;
1037                 m[j] = currentw[j-1] + 10000 * fpenalty; //iinoka?
1038         }
1039         if( lgth2 == 0 )
1040                 lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari
1041         else
1042                 lastverticalw[0] = currentw[lgth2-1];
1043
1044         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1045         lastj = lgth2+1;
1046
1047 #if XXXXXXX
1048 fprintf( stderr, "currentw = \n" );
1049 for( i=0; i<lgth1+1; i++ )
1050 {
1051         fprintf( stderr, "%5.2f ", currentw[i] );
1052 }
1053 fprintf( stderr, "\n" );
1054 fprintf( stderr, "initverticalw = \n" );
1055 for( i=0; i<lgth2+1; i++ )
1056 {
1057         fprintf( stderr, "%5.2f ", initverticalw[i] );
1058 }
1059 fprintf( stderr, "\n" );
1060 fprintf( stderr, "fcgp\n" );
1061 for( i=0; i<lgth1; i++ )
1062         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1063 for( i=0; i<lgth2; i++ )
1064         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1065 #endif
1066
1067         for( i=1; i<lasti; i++ )
1068         {
1069                 wtmp = previousw;
1070                 previousw = currentw;
1071                 currentw = wtmp;
1072
1073                 previousw[0] = initverticalw[i-1];
1074
1075                 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1076 #if XXXXXXX
1077 fprintf( stderr, "\n" );
1078 fprintf( stderr, "i=%d\n", i );
1079 fprintf( stderr, "currentw = \n" );
1080 for( j=0; j<lgth2; j++ )
1081 {
1082         fprintf( stderr, "%5.2f ", currentw[j] );
1083 }
1084 fprintf( stderr, "\n" );
1085 #endif
1086                 if( localhom )
1087                 {
1088 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1089 //                      imp_match_out_vead( currentw, i, lgth2 );
1090                         part_imp_match_out_vead_gapmapQ( currentw, gapmap1[i]+start1, lgth2, start2, gapmap2 );
1091                 }
1092 #if XXXXXXX
1093 fprintf( stderr, "\n" );
1094 fprintf( stderr, "i=%d\n", i );
1095 fprintf( stderr, "currentw = \n" );
1096 for( j=0; j<lgth2; j++ )
1097 {
1098         fprintf( stderr, "%5.2f ", currentw[j] );
1099 }
1100 fprintf( stderr, "\n" );
1101 #endif
1102                 currentw[0] = initverticalw[i];
1103
1104                 mpi = 0;
1105                 mi = previousw[0] + 10000 * fpenalty;
1106
1107                 ijppt = ijp[i] + 1;
1108                 mjpt = m + 1;
1109                 prept = previousw;
1110                 curpt = currentw + 1;
1111                 mpjpt = mp + 1;
1112                 fg_t_og_h_dg_n2_p_pt = fg_t_og_h_dg_n2_p + 1;
1113                 og_t_fg_h_dg_n2_p_pt = og_t_fg_h_dg_n2_p + 1;
1114                 og_h_dg_n2_p_pt = og_h_dg_n2_p + 1;
1115                 fg_h_dg_n2_p_pt = fg_h_dg_n2_p + 1;
1116                 gapz_n2_pt0 = gapz_n2 + 1;
1117                 gapz_n2_pt1 = gapz_n2 + 2;
1118                 fgcp2pt = fgcp2g + 1;
1119                 ogcp2pt = ogcp2g + 1;
1120
1121                 fg_t_og_h_dg_n1_p_va = fg_t_og_h_dg_n1_p[i];
1122                 og_t_fg_h_dg_n1_p_va = og_t_fg_h_dg_n1_p[i];
1123                 og_h_dg_n1_p_va = og_h_dg_n1_p[i];
1124                 fg_h_dg_n1_p_va = fg_h_dg_n1_p[i];
1125                 gapz_n1_va0 = gapz_n1[i];
1126                 gapz_n1_va1 = gapz_n1[i+1];
1127                 fgcp1va = fgcp1g[i];
1128                 ogcp1va = ogcp1g[i];
1129
1130                 for( j=1; j<lastj; j++ )
1131                 {
1132                         wm = *prept;
1133
1134                         g = ogcp1va * *og_h_dg_n2_p_pt;
1135                         wm += g;
1136
1137                         g = *ogcp2pt * og_h_dg_n1_p_va;
1138                         wm += g;
1139
1140                         g = fgcp1va * *fg_h_dg_n2_p_pt;
1141                         wm += g;
1142
1143                         g = *fgcp2pt * fg_h_dg_n1_p_va;
1144                         wm += g;
1145
1146                         *ijppt = 0;
1147
1148
1149 #if 0
1150                         fprintf( stderr, "%5.0f->", wm );
1151 #endif
1152                         tmppenal = gapz_n1_va1 * *fg_t_og_h_dg_n2_p_pt;
1153 #if 0
1154                         fprintf( stderr, "%5.0f?", g );
1155 #endif
1156                         if( (g=mi+tmppenal) > wm )
1157                         {
1158                                 wm = g;
1159                                 *ijppt = -( j - mpi );
1160                         }
1161
1162                         tmppenal = gapz_n1_va0 * *og_t_fg_h_dg_n2_p_pt;
1163                         if( (g=*prept+tmppenal) >= mi )
1164                         {
1165                                 mi = g;
1166                                 mpi = j-1;
1167                         }
1168 #if USE_PENALTY_EX
1169                         mi += fpenalty_ex;
1170 #endif
1171
1172                         tmppenal = *gapz_n2_pt1 * fg_t_og_h_dg_n1_p_va;
1173 #if 0
1174                         fprintf( stderr, "%5.0f?", g );
1175 #endif
1176                         if( (g=*mjpt+tmppenal) > wm )
1177                         {
1178                                 wm = g;
1179                                 *ijppt = +( i - *mpjpt );
1180                         }
1181
1182                         tmppenal = *gapz_n2_pt0 * og_t_fg_h_dg_n1_p_va;
1183                         if( (g=*prept+tmppenal) >= *mjpt )
1184                         {
1185                                 *mjpt = g;
1186                                 *mpjpt = i-1;
1187                         }
1188 #if USE_PENALTY_EX
1189                         m[j] += fpenalty_ex;
1190 #endif
1191
1192 #if 0
1193                         fprintf( stderr, "%5.0f ", wm );
1194 #endif
1195                         *curpt++ += wm;
1196                         ijppt++;
1197                         mjpt++;
1198                         prept++;
1199                         mpjpt++;
1200                         fg_t_og_h_dg_n2_p_pt++;
1201                         og_t_fg_h_dg_n2_p_pt++;
1202                         og_h_dg_n2_p_pt++;
1203                         fg_h_dg_n2_p_pt++;
1204                         gapz_n2_pt0++;
1205                         gapz_n2_pt1++;
1206                         fgcp2pt++;
1207                         ogcp2pt++;
1208                 }
1209                 lastverticalw[i] = currentw[lgth2-1];
1210         }
1211
1212 #if OUTGAP0TRY
1213         if( !outgap )
1214         {
1215                 for( j=1; j<lgth2+1; j++ )
1216                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1217                 for( i=1; i<lgth1+1; i++ )
1218                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1219         }
1220 #endif
1221
1222         /*
1223         fprintf( stderr, "\n" );
1224         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1225         fprintf( stderr, "#####\n" );
1226         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1227         fprintf( stderr, "====>" );
1228         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1229         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1230         */
1231         if( localhom )
1232         {
1233                 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, start1, end1, start2, end2, gapmap1, gapmap2 );
1234         }
1235         else
1236                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1237
1238 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1239
1240         resultlen = strlen( mseq1[0] );
1241         if( alloclen < resultlen || resultlen > N )
1242         {
1243                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1244                 ErrorExit( "LENGTH OVER!\n" );
1245         }
1246
1247
1248         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1249         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1250         /*
1251         fprintf( stderr, "\n" );
1252         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1253         fprintf( stderr, "#####\n" );
1254         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1255         */
1256
1257 //      fprintf( stderr, "wm = %f\n", wm );
1258
1259
1260         return( wm );
1261 }
1262