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

   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 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 *gapf1;
 599         static float *gapf2;
 600         static float *ogcp1g;
 601         static float *ogcp2g;
 602         static float *fgcp1g;
 603         static float *fgcp2g;
 604         static float *ogcp1;
 605         static float *ogcp2;
 606         static float *fgcp1;
 607         static float *fgcp2;
 608         static float **cpmx1;
 609         static float **cpmx2;
 610         static int **intwork;
 611         static float **floatwork;
 612         static int orlgth1 = 0, orlgth2 = 0;
 613         float fpenalty = (float)penalty;
 614         float tmppenal;
 615         float cumpenal;
 616         float *fgcp2pt;
 617         float *ogcp2pt;
 618         float fgcp1va;
 619         float ogcp1va;
 620
 621
 622
 623 #if 0
 624         fprintf( stderr, "####  eff in SA+++align\n" );
 625         fprintf( stderr, "####  seq1[0] = %s\n", seq1[0] );
 626         fprintf( stderr, "####  strlen( seq1[0] ) = %d\n", strlen( seq1[0] ) );
 627         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
 628 #endif
 629         if( orlgth1 == 0 )
 630         {
 631                 mseq1 = AllocateCharMtx( njob, 0 );
 632                 mseq2 = AllocateCharMtx( njob, 0 );
 633         }
 634
 635
 636         lgth1 = strlen( seq1[0] );
 637         lgth2 = strlen( seq2[0] );
 638 #if 0
 639         if( lgth1 == 0 || lgth2 == 0 )
 640         {
 641                 fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
 642         }
 643 #endif
 644
 645         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 646         {
 647                 int ll1, ll2;
 648
 649                 if( orlgth1 > 0 && orlgth2 > 0 )
 650                 {
 651                         FreeFloatVec( w1 );
 652                         FreeFloatVec( w2 );
 653                         FreeFloatVec( match );
 654                         FreeFloatVec( initverticalw );
 655                         FreeFloatVec( lastverticalw );
 656
 657                         FreeFloatVec( m );
 658                         FreeIntVec( mp );
 659
 660                         FreeCharMtx( mseq );
 661
 662                         FreeFloatVec( digf1 );
 663                         FreeFloatVec( digf2 );
 664                         FreeFloatVec( diaf1 );
 665                         FreeFloatVec( diaf2 );
 666                         FreeFloatVec( gapf1 );
 667                         FreeFloatVec( gapf2 );
 668                         FreeFloatVec( ogcp1 );
 669                         FreeFloatVec( ogcp2 );
 670                         FreeFloatVec( fgcp1 );
 671                         FreeFloatVec( fgcp2 );
 672                         FreeFloatVec( ogcp1g );
 673                         FreeFloatVec( ogcp2g );
 674                         FreeFloatVec( fgcp1g );
 675                         FreeFloatVec( fgcp2g );
 676
 677
 678                         FreeFloatMtx( cpmx1 );
 679                         FreeFloatMtx( cpmx2 );
 680
 681                         FreeFloatMtx( floatwork );
 682                         FreeIntMtx( intwork );
 683                 }
 684
 685                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 686                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 687
 688 #if DEBUG
 689                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 690 #endif
 691
 692                 w1 = AllocateFloatVec( ll2+2 );
 693                 w2 = AllocateFloatVec( ll2+2 );
 694                 match = AllocateFloatVec( ll2+2 );
 695
 696                 initverticalw = AllocateFloatVec( ll1+2 );
 697                 lastverticalw = AllocateFloatVec( ll1+2 );
 698
 699                 m = AllocateFloatVec( ll2+2 );
 700                 mp = AllocateIntVec( ll2+2 );
 701
 702                 mseq = AllocateCharMtx( njob, ll1+ll2 );
 703
 704                 digf1 = AllocateFloatVec( ll1+2 );
 705                 digf2 = AllocateFloatVec( ll2+2 );
 706                 diaf1 = AllocateFloatVec( ll1+2 );
 707                 diaf2 = AllocateFloatVec( ll2+2 );
 708                 gapf1 = AllocateFloatVec( ll1+2 );
 709                 gapf2 = AllocateFloatVec( ll2+2 );
 710                 ogcp1 = AllocateFloatVec( ll1+2 );
 711                 ogcp2 = AllocateFloatVec( ll2+2 );
 712                 fgcp1 = AllocateFloatVec( ll1+2 );
 713                 fgcp2 = AllocateFloatVec( ll2+2 );
 714                 ogcp1g = AllocateFloatVec( ll1+2 );
 715                 ogcp2g = AllocateFloatVec( ll2+2 );
 716                 fgcp1g = AllocateFloatVec( ll1+2 );
 717                 fgcp2g = AllocateFloatVec( ll2+2 );
 718
 719                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
 720                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
 721
 722 #if FASTMATCHCALC
 723                 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
 724                 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 );
 725 #else
 726                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
 727                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
 728 #endif
 729
 730 #if DEBUG
 731                 fprintf( stderr, "succeeded\n" );
 732 #endif
 733
 734                 orlgth1 = ll1 - 100;
 735                 orlgth2 = ll2 - 100;
 736         }
 737
 738
 739         for( i=0; i<icyc; i++ )
 740         {
 741                 mseq1[i] = mseq[i];
 742                 seq1[i][lgth1] = 0;
 743         }
 744         for( j=0; j<jcyc; j++ )
 745         {
 746                 mseq2[j] = mseq[icyc+j];
 747                 seq2[j][lgth2] = 0;
 748         }
 749
 750
 751         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
 752         {
 753                 int ll1, ll2;
 754
 755                 if( commonAlloc1 && commonAlloc2 )
 756                 {
 757                         FreeShortMtx( commonIP );
 758                 }
 759
 760                 ll1 = MAX( orlgth1, commonAlloc1 );
 761                 ll2 = MAX( orlgth2, commonAlloc2 );
 762
 763 #if DEBUG
 764                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
 765 #endif
 766
 767                 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
 768
 769 #if DEBUG
 770                 fprintf( stderr, "succeeded\n\n" );
 771 #endif
 772
 773                 commonAlloc1 = ll1;
 774                 commonAlloc2 = ll2;
 775         }
 776         ijp = commonIP;
 777
 778 #if 0
 779         {
 780                 float t = 0.0;
 781                 for( i=0; i<icyc; i++ )
 782                         t += eff1[i];
 783         fprintf( stderr, "## totaleff = %f\n", t );
 784         }
 785 #endif
 786
 787         cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
 788         cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
 789
 790         if( sgap1 )
 791         {
 792                 new_OpeningGapCount( ogcp1g, icyc, seq1, eff1, lgth1, sgap1 );
 793                 new_OpeningGapCount( ogcp2g, jcyc, seq2, eff2, lgth2, sgap2 );
 794                 new_FinalGapCount( fgcp1g, icyc, seq1, eff1, lgth1, egap1 );
 795                 new_FinalGapCount( fgcp2g, jcyc, seq2, eff2, lgth2, egap2 );
 796                 getdigapfreq_part( digf1, icyc, seq1, eff1, lgth1, sgap1, egap1 ); // sgap1 ha iranai ?
 797                 getdigapfreq_part( digf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 ); // sgap1 ha iranai ?
 798                 getdiaminofreq_part( diaf1, icyc, seq1, eff1, lgth1, sgap1, egap1 ); // sgap1 ha iranai ?
 799                 getdiaminofreq_part( diaf2, jcyc, seq2, eff2, lgth2, sgap1, egap2 ); // sgap1 ha iranai ?
 800                 getgapfreq( gapf1, icyc, seq1, eff1, lgth1 ); // atode
 801                 getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 ); // atode
 802         }
 803         else
 804         {
 805                 st_OpeningGapCount( ogcp1g, icyc, seq1, eff1, lgth1 );
 806                 st_OpeningGapCount( ogcp2g, jcyc, seq2, eff2, lgth2 );
 807                 st_FinalGapCount( fgcp1g, icyc, seq1, eff1, lgth1 );
 808                 st_FinalGapCount( fgcp2g, jcyc, seq2, eff2, lgth2 );
 809                 getdigapfreq_st( digf1, icyc, seq1, eff1, lgth1 );
 810                 getdigapfreq_st( digf2, jcyc, seq2, eff2, lgth2 );
 811                 getdiaminofreq_st( diaf1, icyc, seq1, eff1, lgth1 );
 812                 getdiaminofreq_st( diaf2, jcyc, seq2, eff2, lgth2 );
 813                 getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
 814                 getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
 815         }
 816
 817         for( i=0; i<lgth1; i++ )
 818         {
 819                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1g[i] ) * fpenalty;
 820                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1g[i] ) * fpenalty;
 821         }
 822         for( i=0; i<lgth2; i++ )
 823         {
 824                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2g[i] ) * fpenalty;
 825                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2g[i] ) * fpenalty;
 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[-1] * ( 1.0-fgcp2g[-1] ) * 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[-1] * ( 1.0-fgcp1g[-1] ) * 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-1]) + gapf2[j-1]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
 882                         tmppenal = ( (1.0-0.0)*(1.0-ogcp1g[0]+0.0) + 0.0*(1.0-0.0-0.0) ) * 0.5 * fpenalty; // mada
 883 //                      tmppenal = 0.5 * fpenalty;
 884 //                      tmppenal -= ( (1.0-0.0) * (1.0-diaf1[0]) + 0.0 ) * 0.5 * fpenalty; // 0.
 885 //                      tmppenal -= ( (1.0-gapf2[j-1]) * ogcp1g[i] + gapf2[j-1] ) * 0.5 * fpenalty;
 886 //                      fprintf( stderr, "0,0<-%d,%d, tmppenal 1 = %f\n", i, 0, tmppenal );
 887                         initverticalw[i] += tmppenal;
 888
 889 //                      tmppenal = ( (1.0-gapf2[j])*(1.0-fgcp1g[i-1]+ogcp1g[i]) + gapf2[j]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
 890                         tmppenal = ( (1.0-gapf2[0])*(1.0-fgcp1g[i-1]+ogcp1g[i]) + gapf2[0]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
 891 //                      tmppenal = 0.5 * fpenalty;
 892 //                      tmppenal -= ( (1.0-gapf2[0]) * (1.0-diaf1[i]) + gapf2[0] ) * 0.5 * fpenalty;
 893 //                      tmppenal -= ( (1.0-gapf2[j]) * fgcp1g[i-1] + gapf2[j] ) * 0.5 * fpenalty;
 894                         initverticalw[i] += tmppenal;
 895                         cumpenal += ( ogcp1g[i-0] + fgcp1g[i-1] ) * fpenalty * 0.5;
 896                         initverticalw[i] += cumpenal - fgcp1g[i-1] * fpenalty * 0.5;
 897 //                      fprintf( stderr, "0,0<-%d,%d, tmppenal 2 = %f, cumpenal=%f, fgcp1g[i-1]=%f, ogcp1g[i]=%f\n", i, 0, tmppenal, cumpenal, fgcp1g[i-1], ogcp1g[i] );
 898
 899                 }
 900                 cumpenal = 0.0;
 901                 for( j=1; j<lgth2+1; j++ )
 902                 {
 903 //                      currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
 904
 905 //                      tmppenal = ( (1.0-gapf1[i-1])*(1.0-ogcp2g[j]+fgcp2g[j-1]) + gapf1[i-1]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
 906                         tmppenal = ( (1.0-0.0)*(1.0-ogcp2g[0]+0.0) + 0.0*(1.0-0.0-0.0) ) * 0.5 * fpenalty; // mada
 907 //                      tmppenal = 0.5 * fpenalty;
 908 //                      tmppenal -= ( (1.0-0.0) * (1.0-diaf2[0]) + 0.0 ) * 0.5 * fpenalty; // 0.
 909 //                      tmppenal -= ( (1.0-gapf1[0]) * fgcp2g[j-1] + gapf1[0] ) * 0.5 * fpenalty;
 910 //                      fprintf( stderr, "0,0<-%d,%d, tmppenal 3 = %f\n", 0, j, tmppenal );
 911                         currentw[j] += tmppenal;
 912
 913 //                      tmppenal = ( (1.0-gapf1[i])*(1.0-fgcp2g[j-1]+ogcp2g[j]) + gapf1[i]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
 914                         tmppenal = ( (1.0-gapf1[0])*(1.0-fgcp2g[j-1]+ogcp2g[j]) + gapf1[0]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
 915 //                      tmppenal = 0.5 * fpenalty;
 916 //                      tmppenal -= ( (1.0-gapf1[0]) * (1.0-diaf2[j]) + gapf1[0] ) * 0.5 * fpenalty;
 917 //                      tmppenal -= ( (1.0-gapf1[0]) * ogcp2g[j] + gapf1[i-1] ) * 0.5 * fpenalty;
 918 //                      fprintf( stderr, "0,0<-%d,%d, tmppenal 4 = %f\n", 0, j, tmppenal );
 919                         currentw[j] += tmppenal;
 920                         cumpenal += ( ogcp2g[j-0] + fgcp2g[j-1] ) * fpenalty * 0.5;
 921                         currentw[j] += cumpenal - fgcp2g[j-1] * fpenalty * 0.5;
 922                 }
 923         }
 924 #if OUTGAP0TRY
 925         else
 926         {
 927                 for( j=1; j<lgth2+1; j++ )
 928                         currentw[j] -= offset * j / 2.0;
 929                 for( i=1; i<lgth1+1; i++ )
 930                         initverticalw[i] -= offset * i / 2.0;
 931         }
 932 #endif
 933
 934         m[0] = 0.0; // iranai
 935         for( j=1; j<lgth2+1; ++j )
 936         {
 937 //              m[j] = currentw[j-1] + ogcp1[1];
 938                 mp[j] = 0;
 939
 940                 tmppenal = ( (1.0-gapf2[j-1])*(1.0-ogcp1g[1]+fgcp1g[0]) + gapf2[j-1]*(1.0-digf1[1]-diaf1[1]) ) * 0.5 * fpenalty; // mada
 941 //              tmppenal = ( (1.0-gapf2[j-1])*(1.0-ogcp1g[i]+fgcp1g[i-1]) + gapf2[j-1]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
 942 //              tmppenal = 0.5 * fpenalty;
 943 //              tmppenal -= ( (1.0-0.0) * (1.0-0.0) + 0.0 ) * 0.5 * fpenalty;
 944 //              tmppenal -= ( (1.0-gapf2[-1]) * (1.0-diaf1[0]) + gapf2[-1] ) * 0.5 * fpenalty;
 945 //              if( tmppenal ) fprintf( stderr, "%c=%c, end j tmppenal=%f\n", seq1[0][0], seq2[0][j-1], tmppenal );
 946                 m[j] = currentw[j-1] + tmppenal;
 947 //              m[j] = currentw[j-1] + ogcp1[1];
 948         }
 949         if( lgth2 == 0 )
 950                 lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari
 951         else
 952                 lastverticalw[0] = currentw[lgth2-1];
 953
 954         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
 955
 956 #if XXXXXXX
 957 fprintf( stderr, "currentw = \n" );
 958 for( i=0; i<lgth1+1; i++ )
 959 {
 960         fprintf( stderr, "%5.2f ", currentw[i] );
 961 }
 962 fprintf( stderr, "\n" );
 963 fprintf( stderr, "initverticalw = \n" );
 964 for( i=0; i<lgth2+1; i++ )
 965 {
 966         fprintf( stderr, "%5.2f ", initverticalw[i] );
 967 }
 968 fprintf( stderr, "\n" );
 969 fprintf( stderr, "fcgp\n" );
 970 for( i=0; i<lgth1; i++ )
 971         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
 972 for( i=0; i<lgth2; i++ )
 973         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
 974 #endif
 975
 976         for( i=1; i<lasti; i++ )
 977         {
 978                 wtmp = previousw;
 979                 previousw = currentw;
 980                 currentw = wtmp;
 981
 982                 previousw[0] = initverticalw[i-1];
 983
 984                 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
 985 #if XXXXXXX
 986 fprintf( stderr, "\n" );
 987 fprintf( stderr, "i=%d\n", i );
 988 fprintf( stderr, "currentw = \n" );
 989 for( j=0; j<lgth2; j++ )
 990 {
 991         fprintf( stderr, "%5.2f ", currentw[j] );
 992 }
 993 fprintf( stderr, "\n" );
 994 #endif
 995                 if( localhom )
 996                 {
 997 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
 998 #if  0
 999                         imp_match_out_veadQ( currentw, i, lgth2 );
1000 #else
1001                         imp_match_out_veadQ( currentw, i, lgth2 );
1002 #endif
1003                 }
1004 #if XXXXXXX
1005 fprintf( stderr, "\n" );
1006 fprintf( stderr, "i=%d\n", i );
1007 fprintf( stderr, "currentw = \n" );
1008 for( j=0; j<lgth2; j++ )
1009 {
1010         fprintf( stderr, "%5.2f ", currentw[j] );
1011 }
1012 fprintf( stderr, "\n" );
1013 #endif
1014                 currentw[0] = initverticalw[i];
1015
1016
1017 //              mi = previousw[0] + ogcp2[1];  // machigai
1018                 tmppenal = ( (1.0-gapf1[i-1])*(1.0-ogcp2g[1]+fgcp2g[0]) + gapf1[i-1]*(1.0-digf2[1]-diaf2[1]) ) * 0.5 * fpenalty; // mada
1019 //              fprintf( stderr, "%c=%c, end i tmppenal=%f, ogcp2g[1]=%f\n", seq1[0][i-1], seq2[0][0], tmppenal, ogcp2g[1] );
1020                 mi = previousw[0] + tmppenal;
1021
1022                 mpi = 0;
1023                 ijppt = ijp[i] + 1;
1024                 mjpt = m + 1;
1025                 prept = previousw;
1026                 curpt = currentw + 1;
1027                 mpjpt = mp + 1;
1028                 fgcp2pt = fgcp2;
1029                 ogcp2pt = ogcp2 + 1;
1030                 fgcp1va = fgcp1[i-1];
1031                 ogcp1va = ogcp1[i];
1032                 lastj = lgth2+1;
1033                 for( j=1; j<lastj; j++ )
1034                 {
1035                         wm = *prept;
1036
1037                         g = ogcp1g[i] * ( 1.0-ogcp2g[j] ) * fpenalty * 0.5;
1038 //                      if( g ) fprintf( stderr, "match penal1=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1039                         wm += g;
1040
1041                         g = ogcp2g[j] * ( 1.0-ogcp1g[i] ) * fpenalty * 0.5;
1042 //                      if( g ) fprintf( stderr, "match penal2=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1043                         wm += g;
1044
1045                         g = fgcp1g[i-1] * ( 1.0-fgcp2g[j-1] ) * fpenalty * 0.5;
1046 //                      if( g ) fprintf( stderr, "match penal3=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1047                         wm += g;
1048
1049                         g = fgcp2g[j-1] * ( 1.0-fgcp1g[i-1] ) * fpenalty * 0.5;
1050 //                      if( g ) fprintf( stderr, "match penal4=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1051                         wm += g;
1052
1053 #if 0
1054                         g = fgcp2g[j-1] * ogcp1g[i] * fpenalty * 0.5;
1055                         if( g ) fprintf( stderr, "match m-penal5=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1056                         wm -= g;
1057
1058                         g = fgcp1g[i-1] * ogcp2g[j] * fpenalty * 0.5;
1059                         if( g ) fprintf( stderr, "match m-penal6=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1060                         wm -= g;
1061 #endif
1062
1063                         *ijppt = 0;
1064
1065 #if 0
1066                         fprintf( stderr, "%5.0f->", wm );
1067 #endif
1068 #if 0
1069                         fprintf( stderr, "%5.0f?", g );
1070 #endif
1071 //                      tmppenal = ( (1.0-gapf1[i])*(1.0-fgcp2g[j-1]) + gapf1[i]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
1072                         tmppenal = ( (1.0-gapf1[i])*(1.0-fgcp2g[j-1]+ogcp2g[j]) + gapf1[i]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
1073 //                      tmppenal = 0.5 * fpenalty;
1074 //                      tmppenal -= ( (1.0-gapf1[i])*(1.0-diaf2[j]) + gapf1[i] ) * 0.5 * fpenalty;
1075 //                      tmppenal -= ( (1.0-gapf1[i])*fgcp2g[j-1] + gapf1[i] ) * 0.5 * fpenalty;
1076 //                      tmppenal = *fgcp2pt-fpenalty*0.5*gapf1[i];
1077 //                      tmppenal = *fgcp2pt;
1078                         if( (g=mi+tmppenal) > wm )
1079                         {
1080 //                              fprintf( stderr, "jump i start=%f (j=%d, fgcp2g[j-1]=%f, digf2[j]=%f, diaf2[j]=%f), %c-%c\n", g-mi, j, fgcp2g[j-1], digf2[j], diaf2[j], seq1[0][i], seq2[0][j] );
1081                                 wm = g;
1082                                 *ijppt = -( j - mpi );
1083                         }
1084                         tmppenal = ( (1.0-gapf1[i-1])*(1.0-ogcp2g[j]+fgcp2g[j-1]) + gapf1[i-1]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
1085 //                      tmppenal = *ogcp2pt-fpenalty*0.5*(gapf2[j-1]+gapf1[i-1]);
1086 //                      tmppenal = 0.5 * fpenalty;
1087 //                      tmppenal -= ( (1.0-gapf1[i-1])*(1.0-diaf2[j]) + gapf1[i-1] ) * 0.5 * fpenalty;
1088 //                      tmppenal -= ( (1.0-gapf1[i-1])*ogcp2g[j] + gapf1[i-1] ) * 0.5 * fpenalty;
1089 //                      tmppenal = *prept+*ogcp2pt;
1090                         if( (g=*prept+tmppenal) >= mi )
1091                         {
1092 //                              fprintf( stderr, "jump i end=%f, %c-%c\n", g-*prept, seq1[0][i-1], seq2[0][j-1] );
1093                                 mi = g;
1094                                 mpi = j-1;
1095                         }
1096                         else if( j != 1 )
1097                         {
1098                                 mi += ( ogcp2g[j-0] + fgcp2g[j-1] ) * fpenalty * 0.5;
1099 //                              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-1] * fpenalty*0.5 );
1100                         }
1101 #if USE_PENALTY_EX
1102                                 mi += fpenalty_ex;
1103 #endif
1104
1105 #if 0
1106                         fprintf( stderr, "%5.0f?", g );
1107 #endif
1108                         tmppenal = ( (1.0-gapf2[j])*(1.0-fgcp1g[i-1]+ogcp1g[i]) + gapf2[j]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
1109 //                      tmppenal = ( (1.0-gapf2[j])*(1.0-fgcp1g[i-1]) + gapf2[j]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
1110 //                      tmppenal = 0.5 * fpenalty;
1111 //                      tmppenal -= ( (1.0-gapf2[j])*(1.0-diaf1[i]) + gapf2[j] ) * 0.5 * fpenalty;
1112 //                      tmppenal -= ( (1.0-gapf2[j])*fgcp1g[i-1] + gapf2[j] ) * 0.5 * fpenalty;
1113 //                      tmppenal = fgcp1va-fpenalty*0.5*gapf2[j];
1114 //                      tmppenal = fgcp1va;
1115                         if( (g=*mjpt+tmppenal) > wm )
1116                         {
1117 //                              fprintf( stderr, "jump j start=%f, %c-%c\n", g-*mjpt, seq1[0][i], seq2[0][j] );
1118                                 wm = g;
1119                                 *ijppt = +( i - *mpjpt );
1120                         }
1121                         tmppenal = ( (1.0-gapf2[j-1])*(1.0-ogcp1g[i]+fgcp1g[i-1]) + gapf2[j-1]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
1122 //                      tmppenal = ( (1.0-gapf2[j-1])*(1.0-ogcp1g[i]) + gapf2[j-1]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
1123 //                      tmppenal = ogcp1va-fpenalty*0.5*(gapf1[i-1]+gapf2[j-1]);
1124 //                      tmppenal = 0.5 * fpenalty;
1125 //                      tmppenal -= ( (1.0-gapf2[j-1]) * (1.0-diaf1[i]) + gapf2[j-1] ) * 0.5 * fpenalty;
1126 //                      tmppenal -= ( (1.0-gapf2[j-1]) * ogcp1g[i] + gapf2[j-1] ) * 0.5 * fpenalty;
1127 //                      tmppenal = 0.5 * fpenalty - ( (1.0-gapf2[j-1]) * (ogcp1g[i]) + gapf2[j-1] ) * ( 0.5 * fpenalty );
1128 //                      tmppenal = ogcp1va-fpenalty*0.5*gapf2[j-1];
1129 //                      tmppenal = ogcp1va;
1130                         if( (g=*prept+tmppenal) >= *mjpt )
1131                         {
1132 //                              fprintf( stderr, "jump j end=%f, %c-%c\n", g-*prept, seq1[0][i-1], seq2[0][j-1] );
1133                                 *mjpt = g;
1134                                 *mpjpt = i-1;
1135                         }
1136                         else if( i != 1 )
1137                         {
1138                                 m[j] += ( ogcp1g[i-0] + fgcp1g[i-1] ) * fpenalty * 0.5;
1139 //                              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-1] * fpenalty*0.5 );
1140                         }
1141 #if USE_PENALTY_EX
1142                         m[j] += fpenalty_ex;
1143 #endif
1144
1145 #if 0
1146                         fprintf( stderr, "%5.0f ", wm );
1147 #endif
1148                         *curpt++ += wm;
1149                         ijppt++;
1150                         mjpt++;
1151                         prept++;
1152                         mpjpt++;
1153                         fgcp2pt++;
1154                         ogcp2pt++;
1155                 }
1156                 lastverticalw[i] = currentw[lgth2-1];
1157         }
1158
1159 //      fprintf( stderr, "wm = %f\n", wm );
1160
1161 #if OUTGAP0TRY
1162         if( !outgap )
1163         {
1164                 for( j=1; j<lgth2+1; j++ )
1165                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1166                 for( i=1; i<lgth1+1; i++ )
1167                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1168         }
1169 #endif
1170
1171         /*
1172         fprintf( stderr, "\n" );
1173         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1174         fprintf( stderr, "#####\n" );
1175         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1176         fprintf( stderr, "====>" );
1177         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1178         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1179         */
1180         if( localhom )
1181         {
1182                 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1183         }
1184         else
1185                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1186
1187 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1188
1189         resultlen = strlen( mseq1[0] );
1190         if( alloclen < resultlen || resultlen > N )
1191         {
1192                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1193                 ErrorExit( "LENGTH OVER!\n" );
1194         }
1195
1196
1197         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1198         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1199         /*
1200         fprintf( stderr, "\n" );
1201         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1202         fprintf( stderr, "#####\n" );
1203         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1204         */
1205
1206         fprintf( stderr, "wm = %f\n", wm );
1207
1208         return( wm );
1209 }
1210
1211