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