binaries/src/mafft/core/Qalignmm.c.er

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