binaries/src/mafft/core/Ralignmm.c

   1 #include "mltaln.h"
   2 #include "dp.h"
   3
   4 #define MACHIGAI 0
   5 #define OUTGAP0TRY 1
   6 #define DEBUG 0
   7 #define XXXXXXX    0
   8 #define USE_PENALTY_EX  0
   9 #define FASTMATCHCALC 1
  10
  11
  12 static float **impmtx = NULL;
  13 #if 1 // tditeration to naiveRscore_imp de tsukawareru.
  14 float imp_match_out_scR( 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_veadR( 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_tateR( 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_strictR( 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                         int **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_scR( 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                         int **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 R__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         int *ijppt;
 581         float *mjpt, *prept, *curpt;
 582         int *mpjpt;
 583 #endif
 584         static float mi, *m;
 585         static int **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         float kyokaipenal;
 623
 624
 625
 626 #if 0
 627         fprintf( stderr, "####  seq1[0] = %s\n", seq1[0] );
 628         fprintf( stderr, "####  seq2[0] = %s\n", seq2[0] );
 629 #endif
 630         if( orlgth1 == 0 )
 631         {
 632                 mseq1 = AllocateCharMtx( njob, 0 );
 633                 mseq2 = AllocateCharMtx( njob, 0 );
 634         }
 635
 636
 637         lgth1 = strlen( seq1[0] );
 638         lgth2 = strlen( seq2[0] );
 639 #if 0
 640         if( lgth1 == 0 || lgth2 == 0 )
 641         {
 642                 fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
 643         }
 644 #endif
 645
 646         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 647         {
 648                 int ll1, ll2;
 649
 650                 if( orlgth1 > 0 && orlgth2 > 0 )
 651                 {
 652                         FreeFloatVec( w1 );
 653                         FreeFloatVec( w2 );
 654                         FreeFloatVec( match );
 655                         FreeFloatVec( initverticalw );
 656                         FreeFloatVec( lastverticalw );
 657
 658                         FreeFloatVec( m );
 659                         FreeIntVec( mp );
 660
 661                         FreeCharMtx( mseq );
 662
 663                         FreeFloatVec( digf1 );
 664                         FreeFloatVec( digf2 );
 665                         FreeFloatVec( diaf1 );
 666                         FreeFloatVec( diaf2 );
 667                         FreeFloatVec( gapz1 );
 668                         FreeFloatVec( gapz2 );
 669                         FreeFloatVec( gapf1 );
 670                         FreeFloatVec( gapf2 );
 671                         FreeFloatVec( ogcp1 );
 672                         FreeFloatVec( ogcp2 );
 673                         FreeFloatVec( fgcp1 );
 674                         FreeFloatVec( fgcp2 );
 675                         FreeFloatVec( ogcp1g );
 676                         FreeFloatVec( ogcp2g );
 677                         FreeFloatVec( fgcp1g );
 678                         FreeFloatVec( fgcp2g );
 679
 680
 681                         FreeFloatMtx( cpmx1 );
 682                         FreeFloatMtx( cpmx2 );
 683
 684                         FreeFloatMtx( floatwork );
 685                         FreeIntMtx( intwork );
 686                 }
 687
 688                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 689                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 690
 691 #if DEBUG
 692                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 693 #endif
 694
 695                 w1 = AllocateFloatVec( ll2+2 );
 696                 w2 = AllocateFloatVec( ll2+2 );
 697                 match = AllocateFloatVec( ll2+2 );
 698
 699                 initverticalw = AllocateFloatVec( ll1+2 );
 700                 lastverticalw = AllocateFloatVec( ll1+2 );
 701
 702                 m = AllocateFloatVec( ll2+2 );
 703                 mp = AllocateIntVec( ll2+2 );
 704
 705                 mseq = AllocateCharMtx( njob, ll1+ll2 );
 706
 707                 digf1 = AllocateFloatVec( ll1+2 );
 708                 digf2 = AllocateFloatVec( ll2+2 );
 709                 diaf1 = AllocateFloatVec( ll1+2 );
 710                 diaf2 = AllocateFloatVec( ll2+2 );
 711                 gapz1 = AllocateFloatVec( ll1+2 );
 712                 gapz2 = AllocateFloatVec( ll2+2 );
 713                 gapf1 = AllocateFloatVec( ll1+2 );
 714                 gapf2 = AllocateFloatVec( ll2+2 );
 715                 ogcp1 = AllocateFloatVec( ll1+2 );
 716                 ogcp2 = AllocateFloatVec( ll2+2 );
 717                 fgcp1 = AllocateFloatVec( ll1+2 );
 718                 fgcp2 = AllocateFloatVec( ll2+2 );
 719                 ogcp1g = AllocateFloatVec( ll1+2 );
 720                 ogcp2g = AllocateFloatVec( ll2+2 );
 721                 fgcp1g = AllocateFloatVec( ll1+2 );
 722                 fgcp2g = AllocateFloatVec( ll2+2 );
 723
 724                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
 725                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
 726
 727 #if FASTMATCHCALC
 728                 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
 729                 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 );
 730 #else
 731                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
 732                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
 733 #endif
 734
 735 #if DEBUG
 736                 fprintf( stderr, "succeeded\n" );
 737 #endif
 738
 739                 orlgth1 = ll1 - 100;
 740                 orlgth2 = ll2 - 100;
 741         }
 742
 743
 744         for( i=0; i<icyc; i++ )
 745         {
 746                 mseq1[i] = mseq[i];
 747                 seq1[i][lgth1] = 0;
 748         }
 749         for( j=0; j<jcyc; j++ )
 750         {
 751                 mseq2[j] = mseq[icyc+j];
 752                 seq2[j][lgth2] = 0;
 753         }
 754
 755
 756         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
 757         {
 758                 int ll1, ll2;
 759
 760                 if( commonAlloc1 && commonAlloc2 )
 761                 {
 762                         FreeIntMtx( commonIP );
 763                 }
 764
 765                 ll1 = MAX( orlgth1, commonAlloc1 );
 766                 ll2 = MAX( orlgth2, commonAlloc2 );
 767
 768 #if DEBUG
 769                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
 770 #endif
 771
 772                 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
 773
 774 #if DEBUG
 775                 fprintf( stderr, "succeeded\n\n" );
 776 #endif
 777
 778                 commonAlloc1 = ll1;
 779                 commonAlloc2 = ll2;
 780         }
 781         ijp = commonIP;
 782
 783 #if 0
 784         {
 785                 float t = 0.0;
 786                 for( i=0; i<icyc; i++ )
 787                         t += eff1[i];
 788         fprintf( stderr, "## totaleff = %f\n", t );
 789         }
 790 #endif
 791
 792         cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
 793         cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
 794
 795         if( sgap1 )
 796         {
 797                 new_OpeningGapCount_zure( ogcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
 798                 new_OpeningGapCount_zure( ogcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
 799                 new_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
 800                 new_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
 801                 getdigapfreq_part( digf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
 802                 getdigapfreq_part( digf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
 803                 getdiaminofreq_part( diaf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
 804                 getdiaminofreq_part( diaf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
 805                 getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
 806                 getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
 807                 getgapfreq_zure_part( gapz1, icyc, seq1, eff1, lgth1, sgap1 );
 808                 getgapfreq_zure_part( gapz2, jcyc, seq2, eff2, lgth2, sgap1 );
 809         }
 810         else
 811         {
 812                 st_OpeningGapCount( ogcp1g, icyc, seq1, eff1, lgth1 );
 813                 st_OpeningGapCount( ogcp2g, jcyc, seq2, eff2, lgth2 );
 814                 st_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1 );
 815                 st_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2 );
 816                 getdigapfreq_st( digf1, icyc, seq1, eff1, lgth1 );
 817                 getdigapfreq_st( digf2, jcyc, seq2, eff2, lgth2 );
 818                 getdiaminofreq_x( diaf1, icyc, seq1, eff1, lgth1 );
 819                 getdiaminofreq_x( diaf2, jcyc, seq2, eff2, lgth2 );
 820                 getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
 821                 getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
 822                 getgapfreq_zure( gapz1, icyc, seq1, eff1, lgth1 );
 823                 getgapfreq_zure( gapz2, jcyc, seq2, eff2, lgth2 );
 824         }
 825
 826 #if 0
 827         for( i=0; i<lgth1; i++ )
 828                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
 829 #endif
 830
 831         currentw = w1;
 832         previousw = w2;
 833
 834         match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
 835         if( localhom )
 836                 imp_match_out_vead_tateR( initverticalw, 0, lgth1 ); // 060306
 837
 838         match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
 839         if( localhom )
 840                 imp_match_out_veadR( currentw, 0, lgth2 ); // 060306
 841
 842
 843 #if 0 // -> tbfast.c
 844         if( localhom )
 845                 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
 846
 847 #endif
 848
 849
 850
 851         kyokaipenal = 0.0;
 852         if( outgap == 1 )
 853         {
 854                 g = 0.0;
 855
 856                 g += ogcp1g[0] * ( 1.0-ogcp2g[0] ) * fpenalty * 0.5;
 857 //              if( g ) fprintf( stderr, "init-match penal1=%f, %c-%c\n", g, seq1[0][0], seq2[0][0] );
 858
 859                 g += ogcp2g[0] * ( 1.0-ogcp1g[0] ) * fpenalty * 0.5;
 860 //              if( g ) fprintf( stderr, "init-match penal2=%f, %c-%c\n", g, seq1[0][0], seq2[0][0] );
 861
 862                 g += fgcp1g[0] * ( 1.0-fgcp2g[0] ) * fpenalty * 0.5;
 863 //              if( g ) fprintf( stderr, "match penal1=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
 864
 865                 g += fgcp2g[0] * ( 1.0-fgcp1g[0] ) * fpenalty * 0.5;
 866 //              if( g ) fprintf( stderr, "match penal2=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
 867
 868                 kyokaipenal = g;
 869                 initverticalw[0] += g;
 870                 currentw[0] += g;
 871
 872                 cumpenal = 0.0;
 873                 for( i=1; i<lgth1+1; i++ )
 874                 {
 875 //                      initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
 876
 877 //      Q               tmppenal = ( (1.0-gapz2[0])*(1.0-ogcp1g[0]+fgcp1g[0]-digf1[0]) ) * 0.5 * fpenalty; // mada
 878 //                      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
 879                         tmppenal = ( (1.0-gapz2[0])*(1.0-ogcp1g[0]+0.0) + gapz2[0]*(1.0-digf1[0] - diaf1[0]) ) * 0.5 * fpenalty; // mada
 880 //                      tmppenal = ( (1.0-0.0)*(1.0-ogcp1g[0]+0.0) + 0.0*(1.0-0.0-0.0) ) * 0.5 * fpenalty; // mada
 881 //                      tmppenal = 0.5 * fpenalty;
 882 //                      tmppenal -= ( (1.0-0.0) * (1.0-diaf1[0]) + 0.0 ) * 0.5 * fpenalty; // 0.
 883 //                      tmppenal -= ( (1.0-gapf2[j-1]) * ogcp1g[i] + gapf2[j-1] ) * 0.5 * fpenalty;
 884 //                      fprintf( stderr, "0,0<-%d,%d, tmppenal 1 = %f\n", i, 0, tmppenal );
 885                         initverticalw[i] += tmppenal;
 886
 887 //      Q               tmppenal = ( (1.0-gapz2[1])*(1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada
 888 //                      tmppenal = ( (1.0-gapf2[j])*(1.0-fgcp1g[i]+ogcp1g[i]) + gapf2[j]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
 889                         tmppenal = ( (1.0-gapz2[1])*(1.0-fgcp1g[i]+ogcp1g[i]) + gapz2[1]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
 890 //                      tmppenal = ( (1.0-gapf2[0])*(1.0-fgcp1g[i]+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] + gapf2[j] ) * 0.5 * fpenalty;
 894                         initverticalw[i] += tmppenal;
 895                         cumpenal += ( ogcp1g[i-0] + fgcp1g[i] ) * fpenalty * 0.5;
 896                         initverticalw[i] += cumpenal - fgcp1g[i] * fpenalty * 0.5; // honto ha iru
 897 //                      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] );
 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-gapz1[0])*(1.0-ogcp2g[0]+fgcp2g[0]-digf2[0]) ) * 0.5 * fpenalty; // mada
 906 //                      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
 907                         tmppenal = ( (1.0-gapz1[0])*(1.0-ogcp2g[0]+0.0) + gapz1[0]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
 908 //                      tmppenal = ( (1.0-0.0)*(1.0-ogcp2g[0]+0.0) + 0.0*(1.0-0.0-0.0) ) * 0.5 * fpenalty; // mada
 909 //                      tmppenal = 0.5 * fpenalty;
 910 //                      tmppenal -= ( (1.0-0.0) * (1.0-diaf2[0]) + 0.0 ) * 0.5 * fpenalty; // 0.
 911 //                      tmppenal -= ( (1.0-gapf1[0]) * fgcp2g[j] + gapf1[0] ) * 0.5 * fpenalty;
 912 //                      fprintf( stderr, "0,0<-%d,%d, tmppenal 3 = %f\n", 0, j, tmppenal );
 913                         currentw[j] += tmppenal;
 914
 915 //                      tmppenal = ( (1.0-gapz1[1])*(1.0-fgcp2g[j]+ogcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada
 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                 mp[j] = 0;
 942                 m[j] = currentw[j-1] + 10000 * fpenalty; //iinoka?
 943         }
 944         if( lgth2 == 0 )
 945                 lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari
 946         else
 947                 lastverticalw[0] = currentw[lgth2-1];
 948
 949         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
 950
 951 #if XXXXXXX
 952 fprintf( stderr, "currentw = \n" );
 953 for( i=0; i<lgth1+1; i++ )
 954 {
 955         fprintf( stderr, "%5.2f ", currentw[i] );
 956 }
 957 fprintf( stderr, "\n" );
 958 fprintf( stderr, "initverticalw = \n" );
 959 for( i=0; i<lgth2+1; i++ )
 960 {
 961         fprintf( stderr, "%5.2f ", initverticalw[i] );
 962 }
 963 fprintf( stderr, "\n" );
 964 fprintf( stderr, "fcgp\n" );
 965 for( i=0; i<lgth1; i++ )
 966         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
 967 for( i=0; i<lgth2; i++ )
 968         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
 969 #endif
 970
 971         for( i=1; i<lasti; i++ )
 972         {
 973                 wtmp = previousw;
 974                 previousw = currentw;
 975                 currentw = wtmp;
 976
 977                 previousw[0] = initverticalw[i-1];
 978
 979                 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
 980 #if XXXXXXX
 981 fprintf( stderr, "\n" );
 982 fprintf( stderr, "i=%d\n", i );
 983 fprintf( stderr, "currentw = \n" );
 984 for( j=0; j<lgth2; j++ )
 985 {
 986         fprintf( stderr, "%5.2f ", currentw[j] );
 987 }
 988 fprintf( stderr, "\n" );
 989 #endif
 990                 if( localhom )
 991                 {
 992 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
 993 #if  0
 994                         imp_match_out_veadR( currentw, i, lgth2 );
 995 #else
 996                         imp_match_out_veadR( currentw, i, lgth2 );
 997 #endif
 998                 }
 999 #if XXXXXXX
1000 fprintf( stderr, "\n" );
1001 fprintf( stderr, "i=%d\n", i );
1002 fprintf( stderr, "currentw = \n" );
1003 for( j=0; j<lgth2; j++ )
1004 {
1005         fprintf( stderr, "%5.2f ", currentw[j] );
1006 }
1007 fprintf( stderr, "\n" );
1008 #endif
1009                 currentw[0] = initverticalw[i];
1010
1011                 mpi = 0;
1012                 mi = previousw[0] + 10000 * fpenalty;
1013
1014                 ijppt = ijp[i] + 1;
1015                 mjpt = m + 1;
1016                 prept = previousw;
1017                 curpt = currentw + 1;
1018                 mpjpt = mp + 1;
1019                 fgcp2pt = fgcp2;
1020                 ogcp2pt = ogcp2 + 1;
1021                 fgcp1va = fgcp1[i-1];
1022                 ogcp1va = ogcp1[i];
1023                 lastj = lgth2+1;
1024                 for( j=1; j<lastj; j++ )
1025                 {
1026                         wm = *prept;
1027
1028                         g = ogcp1g[i] * ( 1.0-ogcp2g[j] ) * fpenalty * 0.5;
1029 //                      if( g && i==j ) fprintf( stderr, "match penal1=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1030                         wm += g;
1031
1032                         g = ogcp2g[j] * ( 1.0-ogcp1g[i] ) * fpenalty * 0.5;
1033 //                      if( g && i==j ) fprintf( stderr, "match penal2=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1034                         wm += g;
1035
1036                         g = fgcp1g[i] * ( 1.0-fgcp2g[j] ) * fpenalty * 0.5;
1037 //                      if( g && i==j ) fprintf( stderr, "match penal3=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1038                         wm += g;
1039
1040                         g = fgcp2g[j] * ( 1.0-fgcp1g[i] ) * fpenalty * 0.5;
1041 //                      if( g && i==j ) fprintf( stderr, "match penal4=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1042                         wm += g;
1043
1044 #if 0
1045                         g = fgcp2g[j] * ogcp1g[i] * fpenalty * 0.5;
1046                         if( g ) fprintf( stderr, "match m-penal5=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1047                         wm -= g;
1048
1049                         g = fgcp1g[i] * ogcp2g[j] * fpenalty * 0.5;
1050                         if( g ) fprintf( stderr, "match m-penal6=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1051                         wm -= g;
1052 #endif
1053
1054                         *ijppt = 0;
1055
1056 #if 0
1057                         fprintf( stderr, "%5.0f->", wm );
1058 #endif
1059 #if 0
1060                         fprintf( stderr, "%5.0f?", g );
1061 #endif
1062 //                      tmppenal = ( (1.0-gapz1[i+1])*(1.0-fgcp2g[j]+ogcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada
1063 //                      tmppenal = ( (1.0-gapf1[i])*(1.0-fgcp2g[j]) + gapf1[i]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
1064                         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
1065 //                      tmppenal = ( (1.0-gapf1[i])*(1.0-fgcp2g[j]+ogcp2g[j]) + gapf1[i]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
1066 //                      tmppenal = 0.5 * fpenalty;
1067 //                      tmppenal -= ( (1.0-gapf1[i])*(1.0-diaf2[j]) + gapf1[i] ) * 0.5 * fpenalty;
1068 //                      tmppenal -= ( (1.0-gapf1[i])*fgcp2g[j] + gapf1[i] ) * 0.5 * fpenalty;
1069 //                      tmppenal = *fgcp2pt-fpenalty*0.5*gapf1[i];
1070 //                      tmppenal = *fgcp2pt;
1071                         if( (g=mi+tmppenal) > wm )
1072                         {
1073 //                              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] );
1074                                 wm = g;
1075                                 *ijppt = -( j - mpi );
1076                         }
1077 //                      tmppenal = ( (1.0-gapz1[i])*(1.0-ogcp2g[j]+fgcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada
1078                         tmppenal = ( (1.0-gapz1[i])*(1.0-ogcp2g[j]+fgcp2g[j]) + gapz1[i]*(1.0-digf2[j]-diaf2[j]) ) * 0.5 * fpenalty; // mada
1079 //                      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
1080 //                      tmppenal = *ogcp2pt-fpenalty*0.5*(gapf2[j-1]+gapf1[i-1]);
1081 //                      tmppenal = 0.5 * fpenalty;
1082 //                      tmppenal -= ( (1.0-gapf1[i-1])*(1.0-diaf2[j]) + gapf1[i-1] ) * 0.5 * fpenalty;
1083 //                      tmppenal -= ( (1.0-gapf1[i-1])*ogcp2g[j] + gapf1[i-1] ) * 0.5 * fpenalty;
1084 //                      tmppenal = *prept+*ogcp2pt;
1085                         if( (g=*prept+tmppenal) >= mi )
1086                         {
1087 //                              fprintf( stderr, "jump i end=%f, %c-%c\n", g-*prept, seq1[0][i-1], seq2[0][j-1] );
1088                                 mi = g;
1089                                 mpi = j-1;
1090                         }
1091                         else if( j != 1 )
1092                         {
1093                                 mi += ( ogcp2g[j-0] + fgcp2g[j] ) * fpenalty * 0.5; // CHUUI honto ha iru
1094 //                              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 );
1095                         }
1096 #if USE_PENALTY_EX
1097                                 mi += fpenalty_ex;
1098 #endif
1099
1100 #if 0
1101                         fprintf( stderr, "%5.0f?", g );
1102 #endif
1103 //                      tmppenal = ( (1.0-gapz2[j+1])*(1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada
1104                         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
1105 //                      tmppenal = ( (1.0-gapf2[j])*(1.0-fgcp1g[i]+ogcp1g[i]) + gapf2[j]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
1106 //                      tmppenal = ( (1.0-gapf2[j])*(1.0-fgcp1g[i]) + gapf2[j]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
1107 //                      tmppenal = 0.5 * fpenalty;
1108 //                      tmppenal -= ( (1.0-gapf2[j])*(1.0-diaf1[i]) + gapf2[j] ) * 0.5 * fpenalty;
1109 //                      tmppenal -= ( (1.0-gapf2[j])*fgcp1g[i] + gapf2[j] ) * 0.5 * fpenalty;
1110 //                      tmppenal = fgcp1va-fpenalty*0.5*gapf2[j];
1111 //                      tmppenal = fgcp1va;
1112                         if( (g=*mjpt+tmppenal) > wm )
1113                         {
1114 //                              if( seq1[0][i] == 'Y' && seq2[0][j] == 'B' )
1115 //                                      fprintf( stderr, "jump j start=%f, %c-%c\n", g-*mjpt, seq1[0][i], seq2[0][j] );
1116                                 wm = g;
1117                                 *ijppt = +( i - *mpjpt );
1118                         }
1119 //                      tmppenal = ( (1.0-gapz2[j])*(1.0-ogcp1g[i]+fgcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada
1120                         tmppenal = ( (1.0-gapz2[j])*(1.0-ogcp1g[i]+fgcp1g[i]) + gapz2[j]*(1.0-digf1[i]-diaf1[i]) ) * 0.5 * fpenalty; // mada
1121 //                      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
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 //                              if( seq1[0][i] == 'T' && seq2[0][j] == 'T' )
1133 //                                      fprintf( stderr, "jump j end=%f, %c-%c\n", g-*prept, seq1[0][i-1], seq2[0][j-1] );
1134                                 *mjpt = g;
1135                                 *mpjpt = i-1;
1136                         }
1137                         else if( i != 1 )
1138                         {
1139                                 m[j] += ( ogcp1g[i-0] + fgcp1g[i] ) * fpenalty * 0.5; // CHUUI honto ha iru
1140 //                              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 );
1141                         }
1142 #if USE_PENALTY_EX
1143                         m[j] += fpenalty_ex;
1144 #endif
1145
1146 #if 0
1147                         fprintf( stderr, "%5.0f ", wm );
1148 #endif
1149                         *curpt++ += wm;
1150                         ijppt++;
1151                         mjpt++;
1152                         prept++;
1153                         mpjpt++;
1154                         fgcp2pt++;
1155                         ogcp2pt++;
1156                 }
1157                 lastverticalw[i] = currentw[lgth2-1];
1158         }
1159
1160 //      fprintf( stderr, "wm = %f\n", wm );
1161
1162 #if OUTGAP0TRY
1163         if( !outgap )
1164         {
1165                 for( j=1; j<lgth2+1; j++ )
1166                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1167                 for( i=1; i<lgth1+1; i++ )
1168                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1169         }
1170 #endif
1171
1172         /*
1173         fprintf( stderr, "\n" );
1174         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1175         fprintf( stderr, "#####\n" );
1176         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1177         fprintf( stderr, "====>" );
1178         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1179         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1180         */
1181         if( localhom )
1182         {
1183                 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1184         }
1185         else
1186                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1187
1188 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1189
1190         resultlen = strlen( mseq1[0] );
1191         if( alloclen < resultlen || resultlen > N )
1192         {
1193                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1194                 ErrorExit( "LENGTH OVER!\n" );
1195         }
1196
1197
1198         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1199         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1200         /*
1201         fprintf( stderr, "\n" );
1202         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1203         fprintf( stderr, "#####\n" );
1204         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1205         */
1206
1207         fprintf( stderr, "wm = %f\n", wm );
1208
1209
1210         return( wm );
1211 }
1212
1213