binaries/src/mafft/core/partQalignmm.c.dame

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