binaries/src/mafft/core/Salignmm.c.060407

   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 static void OpeningGapCount( float *ogcp, int clus, char **seq, double *eff, int len )
  12 {
  13         int i, j, gc, gb;
  14         float feff;
  15
  16
  17         for( i=0; i<len; i++ ) ogcp[i] = 0.0;
  18         for( j=0; j<clus; j++ )
  19         {
  20                 feff = (float)eff[j];
  21                 gc = 0;
  22                 for( i=0; i<len; i++ )
  23                 {
  24                         gb = gc;
  25                         gc = ( seq[j][i] == '-' );
  26                         {
  27                                 if( !gb *  gc ) ogcp[i] += feff;
  28                         }
  29                 }
  30 //              totaleff+= eff[j];
  31         }
  32 //      for( i=0; i<len; i++ )
  33 //      {
  34 //              fprintf( stderr, "ogcp[%d] = %30.20f\n", i, ogcp[i] );
  35 //              ogcp[i] /= totaleff;
  36 //      }
  37 }
  38
  39 static void FinalGapCount( float *fgcp, int clus, char **seq, double *eff, int len )
  40 {
  41         int i, j, gc, gb;
  42         float feff;
  43
  44         for( i=0; i<len+1; i++ ) fgcp[i] = 0.0;
  45         for( j=0; j<clus; j++ )
  46         {
  47                 feff = (float)eff[j];
  48                 gc = ( seq[j][0] == '-' );
  49                 for( i=1; i<len+1; i++ )
  50                 {
  51                         gb = gc;
  52                         gc = ( seq[j][i] == '-' );
  53                         {
  54                                 if( gb * !gc ) fgcp[i-1] += feff;
  55                         }
  56                 }
  57 //              totaleff += eff[j];
  58 //              fprintf( stderr, "eff[%d] = %30.20f\n", j, eff[j] );
  59         }
  60 //      for( i=0; i<len; i++ )
  61 //      {
  62 //              fgcp[i] /= totaleff;
  63 //              fprintf( stderr, "fgcp[%d] = %30.20f\n", i, fgcp[i] );
  64 //      }
  65 }
  66 static float **impmtx = NULL;
  67 float imp_match_out_sc( int i1, int j1 )
  68 {
  69 //      fprintf( stderr, "imp+match = %f\n", impmtx[i1][j1] * fastathreshold );
  70 //      fprintf( stderr, "val = %f\n", impmtx[i1][j1] );
  71         return( impmtx[i1][j1] );
  72 }
  73 static void imp_match_out_vead_gapmap( float *imp, int i1, int lgth2, int *gapmap2 )
  74 {
  75 #if FASTMATCHCALC
  76         float *pt = impmtx[i1];
  77         int *gapmappt = gapmap2;
  78         while( lgth2-- )
  79                 *imp++ += pt[*gapmappt++];
  80 #else
  81         int j;
  82         float *pt = impmtx[i1];
  83         for( j=0; j<lgth2; j++ )
  84                 *imp++ += pt[gapmap2[j]];
  85 #endif
  86 }
  87
  88 static void imp_match_out_vead_tate_gapmap( float *imp, int j1, int lgth1, int *gapmap1 )
  89 {
  90 #if FASTMATCHCALC
  91         int *gapmappt = gapmap1;
  92         while( lgth1-- )
  93                 *imp++ += impmtx[*gapmappt++][j1];
  94 #else
  95         int i;
  96         for( i=0; i<lgth1; i++ )
  97                 *imp++ += impmtx[gapmap1[i]][j1];
  98 #endif
  99 }
 100
 101 static void imp_match_out_vead( float *imp, int i1, int lgth2 )
 102 {
 103 #if FASTMATCHCALC
 104         float *pt = impmtx[i1];
 105         while( lgth2-- )
 106                 *imp++ += *pt++;
 107 #else
 108         int j;
 109         float *pt = impmtx[i1];
 110         for( j=0; j<lgth2; j++ )
 111                 *imp++ += pt[j];
 112 #endif
 113 }
 114 static void imp_match_out_vead_tate( float *imp, int j1, int lgth1 )
 115 {
 116         int i;
 117         for( i=0; i<lgth1; i++ )
 118                 *imp++ += impmtx[i][j1];
 119 }
 120
 121 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 )
 122 {
 123         int i, j, k1, k2, tmpint, start1, start2, end1, end2;
 124         static int impalloclen = 0;
 125         float effij;
 126         double effijx;
 127         char *pt, *pt1, *pt2;
 128         static char *nocount1 = NULL;
 129         static char *nocount2 = NULL;
 130         LocalHom *tmpptr;
 131
 132         if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
 133         {
 134                 if( impmtx ) FreeFloatMtx( impmtx );
 135                 if( nocount1 ) free( nocount1 );
 136                 if( nocount2 ) free( nocount2 );
 137                 impalloclen = MAX( lgth1, lgth2 ) + 2;
 138                 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
 139                 nocount1 = AllocateCharVec( impalloclen );
 140                 nocount2 = AllocateCharVec( impalloclen );
 141         }
 142
 143         for( i=0; i<lgth1; i++ )
 144         {
 145                 for( j=0; j<clus1; j++ )
 146                         if( seq1[j][i] == '-' ) break;
 147                 if( j != clus1 ) nocount1[i] = 1;
 148                 else                     nocount1[i] = 0;
 149         }
 150         for( i=0; i<lgth2; i++ )
 151         {
 152                 for( j=0; j<clus2; j++ )
 153                         if( seq2[j][i] == '-' ) break;
 154                 if( j != clus2 ) nocount2[i] = 1;
 155                 else                     nocount2[i] = 0;
 156         }
 157
 158 #if 0
 159 fprintf( stderr, "nocount2 =\n" );
 160 for( i = 0; i<impalloclen; i++ )
 161 {
 162         fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
 163 }
 164 #endif
 165
 166
 167
 168 #if 0
 169         fprintf( stderr, "eff1 in _init_strict = \n" );
 170         for( i=0; i<clus1; i++ )
 171                 fprintf( stderr, "eff1[] = %f\n", eff1[i] );
 172         for( i=0; i<clus2; i++ )
 173                 fprintf( stderr, "eff2[] = %f\n", eff2[i] );
 174 #endif
 175
 176         for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
 177                 impmtx[i][j] = 0.0;
 178         effijx =  fastathreshold;
 179         for( i=0; i<clus1; i++ )
 180         {
 181                 for( j=0; j<clus2; j++ )
 182                 {
 183                         effij = (float)( eff1[i] * eff2[j] * effijx );
 184                         tmpptr = localhom[i][j];
 185                         while( tmpptr )
 186                         {
 187 //                              fprintf( stderr, "start1 = %d\n", tmpptr->start1 );
 188 //                              fprintf( stderr, "end1   = %d\n", tmpptr->end1   );
 189 //                              fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] );
 190 //                              fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] );
 191                                 pt = seq1[i];
 192                                 tmpint = -1;
 193                                 while( *pt != 0 )
 194                                 {
 195                                         if( *pt++ != '-' ) tmpint++;
 196                                         if( tmpint == tmpptr->start1 ) break;
 197                                 }
 198                                 start1 = pt - seq1[i] - 1;
 199
 200                                 if( tmpptr->start1 == tmpptr->end1 ) end1 = start1;
 201                                 else
 202                                 {
 203 #if MACHIGAI
 204                                         while( *pt != 0 )
 205                                         {
 206 //                                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
 207                                                 if( tmpint == tmpptr->end1 ) break;
 208                                                 if( *pt++ != '-' ) tmpint++;
 209                                         }
 210                                         end1 = pt - seq1[i] - 0;
 211 #else
 212                                         while( *pt != 0 )
 213                                         {
 214 //                                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
 215                                                 if( *pt++ != '-' ) tmpint++;
 216                                                 if( tmpint == tmpptr->end1 ) break;
 217                                         }
 218                                         end1 = pt - seq1[i] - 1;
 219 #endif
 220                                 }
 221
 222                                 pt = seq2[j];
 223                                 tmpint = -1;
 224                                 while( *pt != 0 )
 225                                 {
 226                                         if( *pt++ != '-' ) tmpint++;
 227                                         if( tmpint == tmpptr->start2 ) break;
 228                                 }
 229                                 start2 = pt - seq2[j] - 1;
 230                                 if( tmpptr->start2 == tmpptr->end2 ) end2 = start2;
 231                                 else
 232                                 {
 233 #if MACHIGAI
 234                                         while( *pt != 0 )
 235                                         {
 236                                                 if( tmpint == tmpptr->end2 ) break;
 237                                                 if( *pt++ != '-' ) tmpint++;
 238                                         }
 239                                         end2 = pt - seq2[j] - 0;
 240 #else
 241                                         while( *pt != 0 )
 242                                         {
 243                                                 if( *pt++ != '-' ) tmpint++;
 244                                                 if( tmpint == tmpptr->end2 ) break;
 245                                         }
 246                                         end2 = pt - seq2[j] - 1;
 247 #endif
 248                                 }
 249 //                              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] );
 250 //                              fprintf( stderr, "step 0\n" );
 251                                 if( end1 - start1 != end2 - start2 )
 252                                 {
 253 //                                      fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
 254                                 }
 255
 256 #if 1
 257                                 k1 = start1; k2 = start2;
 258                                 pt1 = seq1[i] + k1;
 259                                 pt2 = seq2[j] + k2;
 260                                 while( *pt1 && *pt2 )
 261                                 {
 262                                         if( *pt1 != '-' && *pt2 != '-' )
 263                                         {
 264 // ½Å¤ß¤òÆó½Å¤Ë¤«¤±¤Ê¤¤¤è¤¦¤ËÃí°Õ¤·¤Æ²¼¤µ¤¤¡£
 265 //                                              impmtx[k1][k2] += tmpptr->wimportance * fastathreshold;
 266 //                                              impmtx[k1][k2] += tmpptr->importance * effij;
 267                                                 impmtx[k1][k2] += tmpptr->fimportance * effij;
 268 //                                              fprintf( stderr, "#### impmtx[k1][k2] = %f, tmpptr->fimportance=%f, effij=%f\n", impmtx[k1][k2], tmpptr->fimportance, effij );
 269 //                                              fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 270 //                                              fprintf( stderr, "%d (%c) - %d (%c)  - %f\n", k1, *pt1, k2, *pt2, tmpptr->fimportance * effij );
 271                                                 k1++; k2++;
 272                                                 pt1++; pt2++;
 273                                         }
 274                                         else if( *pt1 != '-' && *pt2 == '-' )
 275                                         {
 276 //                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 277                                                 k2++; pt2++;
 278                                         }
 279                                         else if( *pt1 == '-' && *pt2 != '-' )
 280                                         {
 281 //                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 282                                                 k1++; pt1++;
 283                                         }
 284                                         else if( *pt1 == '-' && *pt2 == '-' )
 285                                         {
 286 //                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 287                                                 k1++; pt1++;
 288                                                 k2++; pt2++;
 289                                         }
 290                                         if( k1 > end1 || k2 > end2 ) break;
 291                                 }
 292 #else
 293                                 while( k1 <= end1 && k2 <= end2 )
 294                                 {
 295                                         fprintf( stderr, "k1,k2=%d,%d - ", k1, k2 );
 296                                         if( !nocount1[k1] && !nocount2[k2] )
 297                                         {
 298                                                 impmtx[k1][k2] += tmpptr->wimportance * eff1[i] * eff2[j]  * fastathreshold;
 299                                                 fprintf( stderr, "marked\n" );
 300                                         }
 301                                         else
 302                                                 fprintf( stderr, "no count\n" );
 303                                         k1++; k2++;
 304                                 }
 305 #endif
 306                                 tmpptr = tmpptr->next;
 307                         }
 308                 }
 309         }
 310 #if 0
 311         if( clus1 == 1 && clus2 == 6 )
 312         {
 313                 fprintf( stderr, "\n" );
 314                 fprintf( stderr, "seq1[0] =  %s\n", seq1[0] );
 315                 fprintf( stderr, "seq2[0] =  %s\n", seq2[0] );
 316                 fprintf( stderr, "impmtx = \n" );
 317                 for( k2=0; k2<lgth2; k2++ )
 318                         fprintf( stderr, "%6.3f ", (double)k2 );
 319                 fprintf( stderr, "\n" );
 320                 for( k1=0; k1<lgth1; k1++ )
 321                 {
 322                         fprintf( stderr, "%d ", k1 );
 323                         for( k2=0; k2<3; k2++ )
 324                                 fprintf( stderr, "%2.1f ", impmtx[k1][k2] );
 325                         fprintf( stderr, "\n" );
 326                 }
 327                 exit( 1 );
 328         }
 329 #endif
 330 }
 331
 332 #if 0
 333 void imp_match_init( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom )
 334 {
 335         int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2;
 336         static int impalloclen = 0;
 337         char *pt;
 338         int allgap;
 339         static char *nocount1 = NULL;
 340         static char *nocount2 = NULL;
 341
 342         if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
 343         {
 344                 if( impmtx ) FreeFloatMtx( impmtx );
 345                 if( nocount1 ) free( nocount1 );
 346                 if( nocount2 ) free( nocount2 );
 347                 impalloclen = MAX( lgth1, lgth2 ) + 2;
 348                 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
 349                 nocount1 = AllocateCharVec( impalloclen );
 350                 nocount2 = AllocateCharVec( impalloclen );
 351         }
 352
 353         for( i=0; i<lgth1; i++ )
 354         {
 355                 for( j=0; j<clus1; j++ )
 356                         if( seq1[j][i] == '-' ) break;
 357                 if( j != clus1 ) nocount1[i] = 1;
 358                 else                     nocount1[i] = 0;
 359         }
 360         for( i=0; i<lgth2; i++ )
 361         {
 362                 for( j=0; j<clus2; j++ )
 363                         if( seq2[j][i] == '-' ) break;
 364                 if( j != clus2 ) nocount2[i] = 1;
 365                 else                     nocount2[i] = 0;
 366         }
 367
 368 #if 0
 369 fprintf( stderr, "nocount2 =\n" );
 370 for( i = 0; i<impalloclen; i++ )
 371 {
 372         fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
 373 }
 374 #endif
 375
 376         for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
 377                 impmtx[i][j] = 0;
 378         for( i=0; i<clus1; i++ )
 379         {
 380                 fprintf( stderr, "i = %d, seq1 = %s\n", i, seq1[i] );
 381                 for( j=0; j<clus2; j++ )
 382                 {
 383                         fprintf( stderr, "start1 = %d\n", localhom[i][j]->start1 );
 384                         fprintf( stderr, "end1   = %d\n", localhom[i][j]->end1   );
 385                         fprintf( stderr, "j = %d, seq2 = %s\n", j, seq2[j] );
 386                         pt = seq1[i];
 387                         tmpint = -1;
 388                         while( *pt != 0 )
 389                         {
 390                                 if( *pt++ != '-' ) tmpint++;
 391                                 if( tmpint == localhom[i][j]->start1 ) break;
 392                         }
 393                         start1 = pt - seq1[i] - 1;
 394
 395                         while( *pt != 0 )
 396                         {
 397 //                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, localhom[i][j].end1, pt-seq1[i] );
 398                                 if( *pt++ != '-' ) tmpint++;
 399                                 if( tmpint == localhom[i][j]->end1 ) break;
 400                         }
 401                         end1 = pt - seq1[i] - 1;
 402
 403                         pt = seq2[j];
 404                         tmpint = -1;
 405                         while( *pt != 0 )
 406                         {
 407                                 if( *pt++ != '-' ) tmpint++;
 408                                 if( tmpint == localhom[i][j]->start2 ) break;
 409                         }
 410                         start2 = pt - seq2[j] - 1;
 411                         while( *pt != 0 )
 412                         {
 413                                 if( *pt++ != '-' ) tmpint++;
 414                                 if( tmpint == localhom[i][j]->end2 ) break;
 415                         }
 416                         end2 = pt - seq2[j] - 1;
 417 //                      fprintf( stderr, "start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
 418                         k1 = start1;
 419                         k2 = start2;
 420                         fprintf( stderr, "step 0\n" );
 421                         while( k1 <= end1 && k2 <= end2 )
 422                         {
 423 #if 0
 424                                 if( !nocount1[k1] && !nocount2[k2] )
 425                                         impmtx[k1][k2] += localhom[i][j].wimportance * eff1[i] * eff2[j];
 426                                 k1++; k2++;
 427 #else
 428                                 if( !nocount1[k1] && !nocount2[k2] )
 429                                         impmtx[k1][k2] += localhom[i][j]->wimportance * eff1[i] * eff2[j];
 430                                 k1++; k2++;
 431 #endif
 432                         }
 433
 434                         dif = ( end1 - start1 ) - ( end2 - start2 );
 435                         fprintf( stderr, "dif = %d\n", dif );
 436                         if( dif > 0 )
 437                         {
 438                                 do
 439                                 {
 440                                         fprintf( stderr, "dif = %d\n", dif );
 441                                         k1 = start1;
 442                                         k2 = start2 - dif;
 443                                         while( k1 <= end1 && k2 <= end2 )
 444                                         {
 445                                                 if( 0 <= k2 && start2 <= k2 && !nocount1[k1] && !nocount2[k2] )
 446                                                         impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
 447                                                 k1++; k2++;
 448                                         }
 449                                 }
 450                                 while( dif-- );
 451                         }
 452                         else
 453                         {
 454                                 do
 455                                 {
 456                                         k1 = start1 + dif;
 457                                         k2 = start2;
 458                                         while( k1 <= end1 )
 459                                         {
 460                                                 if( k1 >= 0 && k1 >= start1 && !nocount1[k1] && !nocount2[k2] )
 461                                                         impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
 462                                                 k1++; k2++;
 463                                         }
 464                                 }
 465                                 while( dif++ );
 466                         }
 467                 }
 468         }
 469 #if 0
 470         fprintf( stderr, "impmtx = \n" );
 471         for( k2=0; k2<lgth2; k2++ )
 472                 fprintf( stderr, "%6.3f ", (double)k2 );
 473         fprintf( stderr, "\n" );
 474         for( k1=0; k1<lgth1; k1++ )
 475         {
 476                 fprintf( stderr, "%d", k1 );
 477                 for( k2=0; k2<lgth2; k2++ )
 478                         fprintf( stderr, "%6.3f ", impmtx[k1][k2] );
 479                 fprintf( stderr, "\n" );
 480         }
 481 #endif
 482 }
 483 #endif
 484
 485 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
 486 {
 487 #if FASTMATCHCALC
 488         int j, l;
 489         float scarr[26];
 490         float **cpmxpd = floatwork;
 491         int **cpmxpdn = intwork;
 492         float *matchpt, *cpmxpdpt, **cpmxpdptpt;
 493         int *cpmxpdnpt, **cpmxpdnptpt;
 494         if( initialize )
 495         {
 496                 int count = 0;
 497                 for( j=0; j<lgth2; j++ )
 498                 {
 499                         count = 0;
 500                         for( l=0; l<26; l++ )
 501                         {
 502                                 if( cpmx2[l][j] )
 503                                 {
 504                                         cpmxpd[j][count] = cpmx2[l][j];
 505                                         cpmxpdn[j][count] = l;
 506                                         count++;
 507                                 }
 508                         }
 509                         cpmxpdn[j][count] = -1;
 510                 }
 511         }
 512
 513         {
 514                 for( l=0; l<26; l++ )
 515                 {
 516                         scarr[l] = 0.0;
 517                         for( j=0; j<26; j++ )
 518                                 scarr[l] += n_dis[j][l] * cpmx1[j][i1];
 519                 }
 520                 matchpt = match;
 521                 cpmxpdnptpt = cpmxpdn;
 522                 cpmxpdptpt = cpmxpd;
 523                 while( lgth2-- )
 524                 {
 525                         *matchpt = 0.0;
 526                         cpmxpdnpt = *cpmxpdnptpt++;
 527                         cpmxpdpt = *cpmxpdptpt++;
 528                         while( *cpmxpdnpt>-1 )
 529                                 *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;
 530                         matchpt++;
 531                 }
 532         }
 533 #else
 534         int j, k, l;
 535         float scarr[26];
 536         float **cpmxpd = floatwork;
 537         int **cpmxpdn = intwork;
 538 // simple
 539         if( initialize )
 540         {
 541                 int count = 0;
 542                 for( j=0; j<lgth2; j++ )
 543                 {
 544                         count = 0;
 545                         for( l=0; l<26; l++ )
 546                         {
 547                                 if( cpmx2[l][j] )
 548                                 {
 549                                         cpmxpd[count][j] = cpmx2[l][j];
 550                                         cpmxpdn[count][j] = l;
 551                                         count++;
 552                                 }
 553                         }
 554                         cpmxpdn[count][j] = -1;
 555                 }
 556         }
 557         for( l=0; l<26; l++ )
 558         {
 559                 scarr[l] = 0.0;
 560                 for( k=0; k<26; k++ )
 561                         scarr[l] += n_dis[k][l] * cpmx1[k][i1];
 562         }
 563         for( j=0; j<lgth2; j++ )
 564         {
 565                 match[j] = 0.0;
 566                 for( k=0; cpmxpdn[k][j]>-1; k++ )
 567                         match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
 568         }
 569 #endif
 570 }
 571
 572 static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
 573                                                 char **seq1, char **seq2,
 574                         char **mseq1, char **mseq2,
 575                         float **cpmx1, float **cpmx2,
 576                         short **ijp, int icyc, int jcyc )
 577 {
 578         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
 579         char gap[] = "-";
 580         float wm;
 581         lgth1 = strlen( seq1[0] );
 582         lgth2 = strlen( seq2[0] );
 583
 584 #if 0
 585         for( i=0; i<lgth1; i++ )
 586         {
 587                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 588         }
 589 #endif
 590
 591         if( outgap == 1 )
 592                 ;
 593         else
 594         {
 595                 wm = lastverticalw[0];
 596                 for( i=0; i<lgth1; i++ )
 597                 {
 598                         if( lastverticalw[i] >= wm )
 599                         {
 600                                 wm = lastverticalw[i];
 601                                 iin = i; jin = lgth2-1;
 602                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 603                         }
 604                 }
 605                 for( j=0; j<lgth2; j++ )
 606                 {
 607                         if( lasthorizontalw[j] >= wm )
 608                         {
 609                                 wm = lasthorizontalw[j];
 610                                 iin = lgth1-1; jin = j;
 611                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 612                         }
 613                 }
 614         }
 615
 616     for( i=0; i<lgth1+1; i++ )
 617     {
 618         ijp[i][0] = i + 1;
 619     }
 620     for( j=0; j<lgth2+1; j++ )
 621     {
 622         ijp[0][j] = -( j + 1 );
 623     }
 624
 625         for( i=0; i<icyc; i++ )
 626         {
 627                 mseq1[i] += lgth1+lgth2;
 628                 *mseq1[i] = 0;
 629         }
 630         for( j=0; j<jcyc; j++ )
 631         {
 632                 mseq2[j] += lgth1+lgth2;
 633                 *mseq2[j] = 0;
 634         }
 635         iin = lgth1; jin = lgth2;
 636         *impwmpt = 0.0;
 637         for( k=0; k<=lgth1+lgth2; k++ )
 638         {
 639                 if( ijp[iin][jin] < 0 )
 640                 {
 641                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 642                 }
 643                 else if( ijp[iin][jin] > 0 )
 644                 {
 645                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 646                 }
 647                 else
 648                 {
 649                         ifi = iin-1; jfi = jin-1;
 650                 }
 651                 l = iin - ifi;
 652                 while( --l )
 653                 {
 654                         for( i=0; i<icyc; i++ )
 655                                 *--mseq1[i] = seq1[i][ifi+l];
 656                         for( j=0; j<jcyc; j++ )
 657                                 *--mseq2[j] = *gap;
 658                         k++;
 659                 }
 660                 l= jin - jfi;
 661                 while( --l )
 662                 {
 663                         for( i=0; i<icyc; i++ )
 664                                 *--mseq1[i] = *gap;
 665                         for( j=0; j<jcyc; j++ )
 666                                 *--mseq2[j] = seq2[j][jfi+l];
 667                         k++;
 668                 }
 669                 if( iin == lgth1 || jin == lgth2 )
 670                         ;
 671                 else
 672                 {
 673                         *impwmpt += imp_match_out_sc( iin, jin );
 674
 675 //              fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
 676                 }
 677                 if( iin <= 0 || jin <= 0 ) break;
 678                 for( i=0; i<icyc; i++ )
 679                         *--mseq1[i] = seq1[i][ifi];
 680                 for( j=0; j<jcyc; j++ )
 681                         *--mseq2[j] = seq2[j][jfi];
 682                 k++;
 683                 iin = ifi; jin = jfi;
 684         }
 685 }
 686 static void Atracking_localhom_gapmap( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
 687                                                 char **seq1, char **seq2,
 688                         char **mseq1, char **mseq2,
 689                         float **cpmx1, float **cpmx2,
 690                         short **ijp, int icyc, int jcyc,
 691                                                 int *gapmap1, int *gapmap2 )
 692 {
 693         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
 694         char gap[] = "-";
 695         float wm;
 696         lgth1 = strlen( seq1[0] );
 697         lgth2 = strlen( seq2[0] );
 698
 699 #if 0
 700         for( i=0; i<lgth1; i++ )
 701         {
 702                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 703         }
 704 #endif
 705
 706         if( outgap == 1 )
 707                 ;
 708         else
 709         {
 710                 wm = lastverticalw[0];
 711                 for( i=0; i<lgth1; i++ )
 712                 {
 713                         if( lastverticalw[i] >= wm )
 714                         {
 715                                 wm = lastverticalw[i];
 716                                 iin = i; jin = lgth2-1;
 717                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 718                         }
 719                 }
 720                 for( j=0; j<lgth2; j++ )
 721                 {
 722                         if( lasthorizontalw[j] >= wm )
 723                         {
 724                                 wm = lasthorizontalw[j];
 725                                 iin = lgth1-1; jin = j;
 726                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 727                         }
 728                 }
 729         }
 730
 731     for( i=0; i<lgth1+1; i++ )
 732     {
 733         ijp[i][0] = i + 1;
 734     }
 735     for( j=0; j<lgth2+1; j++ )
 736     {
 737         ijp[0][j] = -( j + 1 );
 738     }
 739
 740         for( i=0; i<icyc; i++ )
 741         {
 742                 mseq1[i] += lgth1+lgth2;
 743                 *mseq1[i] = 0;
 744         }
 745         for( j=0; j<jcyc; j++ )
 746         {
 747                 mseq2[j] += lgth1+lgth2;
 748                 *mseq2[j] = 0;
 749         }
 750         iin = lgth1; jin = lgth2;
 751         *impwmpt = 0.0;
 752         for( k=0; k<=lgth1+lgth2; k++ )
 753         {
 754                 if( ijp[iin][jin] < 0 )
 755                 {
 756                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 757                 }
 758                 else if( ijp[iin][jin] > 0 )
 759                 {
 760                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 761                 }
 762                 else
 763                 {
 764                         ifi = iin-1; jfi = jin-1;
 765                 }
 766                 l = iin - ifi;
 767                 while( --l )
 768                 {
 769                         for( i=0; i<icyc; i++ )
 770                                 *--mseq1[i] = seq1[i][ifi+l];
 771                         for( j=0; j<jcyc; j++ )
 772                                 *--mseq2[j] = *gap;
 773                         k++;
 774                 }
 775                 l= jin - jfi;
 776                 while( --l )
 777                 {
 778                         for( i=0; i<icyc; i++ )
 779                                 *--mseq1[i] = *gap;
 780                         for( j=0; j<jcyc; j++ )
 781                                 *--mseq2[j] = seq2[j][jfi+l];
 782                         k++;
 783                 }
 784                 if( iin == lgth1 || jin == lgth2 )
 785                         ;
 786                 else
 787                 {
 788                         *impwmpt += imp_match_out_sc( gapmap1[iin], gapmap2[jin] );
 789
 790 //              fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
 791                 }
 792                 if( iin <= 0 || jin <= 0 ) break;
 793                 for( i=0; i<icyc; i++ )
 794                         *--mseq1[i] = seq1[i][ifi];
 795                 for( j=0; j<jcyc; j++ )
 796                         *--mseq2[j] = seq2[j][jfi];
 797                 k++;
 798                 iin = ifi; jin = jfi;
 799         }
 800 }
 801 static float Atracking( float *lasthorizontalw, float *lastverticalw,
 802                                                 char **seq1, char **seq2,
 803                         char **mseq1, char **mseq2,
 804                         float **cpmx1, float **cpmx2,
 805                         short **ijp, int icyc, int jcyc )
 806 {
 807         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
 808         char gap[] = "-";
 809         float wm;
 810         lgth1 = strlen( seq1[0] );
 811         lgth2 = strlen( seq2[0] );
 812
 813 #if 0
 814         for( i=0; i<lgth1; i++ )
 815         {
 816                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 817         }
 818 #endif
 819
 820         if( outgap == 1 )
 821                 ;
 822         else
 823         {
 824                 wm = lastverticalw[0];
 825                 for( i=0; i<lgth1; i++ )
 826                 {
 827                         if( lastverticalw[i] >= wm )
 828                         {
 829                                 wm = lastverticalw[i];
 830                                 iin = i; jin = lgth2-1;
 831                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 832                         }
 833                 }
 834                 for( j=0; j<lgth2; j++ )
 835                 {
 836                         if( lasthorizontalw[j] >= wm )
 837                         {
 838                                 wm = lasthorizontalw[j];
 839                                 iin = lgth1-1; jin = j;
 840                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 841                         }
 842                 }
 843         }
 844
 845     for( i=0; i<lgth1+1; i++ )
 846     {
 847         ijp[i][0] = i + 1;
 848     }
 849     for( j=0; j<lgth2+1; j++ )
 850     {
 851         ijp[0][j] = -( j + 1 );
 852     }
 853
 854         for( i=0; i<icyc; i++ )
 855         {
 856                 mseq1[i] += lgth1+lgth2;
 857                 *mseq1[i] = 0;
 858         }
 859         for( j=0; j<jcyc; j++ )
 860         {
 861                 mseq2[j] += lgth1+lgth2;
 862                 *mseq2[j] = 0;
 863         }
 864         iin = lgth1; jin = lgth2;
 865         for( k=0; k<=lgth1+lgth2; k++ )
 866         {
 867                 if( ijp[iin][jin] < 0 )
 868                 {
 869                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 870                 }
 871                 else if( ijp[iin][jin] > 0 )
 872                 {
 873                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 874                 }
 875                 else
 876                 {
 877                         ifi = iin-1; jfi = jin-1;
 878                 }
 879                 l = iin - ifi;
 880                 while( --l )
 881                 {
 882                         for( i=0; i<icyc; i++ )
 883                                 *--mseq1[i] = seq1[i][ifi+l];
 884                         for( j=0; j<jcyc; j++ )
 885                                 *--mseq2[j] = *gap;
 886                         k++;
 887                 }
 888                 l= jin - jfi;
 889                 while( --l )
 890                 {
 891                         for( i=0; i<icyc; i++ )
 892                                 *--mseq1[i] = *gap;
 893                         for( j=0; j<jcyc; j++ )
 894                                 *--mseq2[j] = seq2[j][jfi+l];
 895                         k++;
 896                 }
 897                 if( iin <= 0 || jin <= 0 ) break;
 898                 for( i=0; i<icyc; i++ )
 899                         *--mseq1[i] = seq1[i][ifi];
 900                 for( j=0; j<jcyc; j++ )
 901                         *--mseq2[j] = seq2[j][jfi];
 902                 k++;
 903                 iin = ifi; jin = jfi;
 904         }
 905         return( 0.0 );
 906 }
 907
 908 float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch )
 909 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 910 {
 911 //      int k;
 912         register int i, j;
 913         int lasti;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
 914         int lgth1, lgth2;
 915         int resultlen;
 916         float wm = 0.0;   /* int ?????? */
 917         float g;
 918         float *currentw, *previousw;
 919 //      float fpenalty = (float)penalty;
 920 #if USE_PENALTY_EX
 921         float fpenalty_ex = (float)penalty_ex;
 922 #endif
 923 #if 1
 924         float *wtmp;
 925         short *ijppt;
 926         float *mjpt, *prept, *curpt;
 927         int *mpjpt;
 928 #endif
 929         static float mi, *m;
 930         static short **ijp;
 931         static int mpi, *mp;
 932         static float *w1, *w2;
 933         static float *match;
 934         static float *initverticalw;    /* kufuu sureba iranai */
 935         static float *lastverticalw;    /* kufuu sureba iranai */
 936         static char **mseq1;
 937         static char **mseq2;
 938         static char **mseq;
 939         static float *ogcp1;
 940         static float *ogcp2;
 941         static float *fgcp1;
 942         static float *fgcp2;
 943         static float **cpmx1;
 944         static float **cpmx2;
 945         static int **intwork;
 946         static float **floatwork;
 947         static int orlgth1 = 0, orlgth2 = 0;
 948         float fpenalty = (float)penalty;
 949
 950
 951 #if 0
 952         fprintf( stderr, "eff in SA+++align\n" );
 953         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
 954 #endif
 955         if( orlgth1 == 0 )
 956         {
 957                 mseq1 = AllocateCharMtx( njob, 0 );
 958                 mseq2 = AllocateCharMtx( njob, 0 );
 959         }
 960
 961
 962         lgth1 = strlen( seq1[0] );
 963         lgth2 = strlen( seq2[0] );
 964
 965         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 966         {
 967                 int ll1, ll2;
 968
 969                 if( orlgth1 > 0 && orlgth2 > 0 )
 970                 {
 971                         FreeFloatVec( w1 );
 972                         FreeFloatVec( w2 );
 973                         FreeFloatVec( match );
 974                         FreeFloatVec( initverticalw );
 975                         FreeFloatVec( lastverticalw );
 976
 977                         FreeFloatVec( m );
 978                         FreeIntVec( mp );
 979
 980                         FreeCharMtx( mseq );
 981
 982                         FreeFloatVec( ogcp1 );
 983                         FreeFloatVec( ogcp2 );
 984                         FreeFloatVec( fgcp1 );
 985                         FreeFloatVec( fgcp2 );
 986
 987
 988                         FreeFloatMtx( cpmx1 );
 989                         FreeFloatMtx( cpmx2 );
 990
 991                         FreeFloatMtx( floatwork );
 992                         FreeIntMtx( intwork );
 993                 }
 994
 995                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 996                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 997
 998 #if DEBUG
 999                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1000 #endif
1001
1002                 w1 = AllocateFloatVec( ll2+2 );
1003                 w2 = AllocateFloatVec( ll2+2 );
1004                 match = AllocateFloatVec( ll2+2 );
1005
1006                 initverticalw = AllocateFloatVec( ll1+2 );
1007                 lastverticalw = AllocateFloatVec( ll1+2 );
1008
1009                 m = AllocateFloatVec( ll2+2 );
1010                 mp = AllocateIntVec( ll2+2 );
1011
1012                 mseq = AllocateCharMtx( njob, ll1+ll2 );
1013
1014                 ogcp1 = AllocateFloatVec( ll1+2 );
1015                 ogcp2 = AllocateFloatVec( ll2+2 );
1016                 fgcp1 = AllocateFloatVec( ll1+2 );
1017                 fgcp2 = AllocateFloatVec( ll2+2 );
1018
1019                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1020                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1021
1022 #if FASTMATCHCALC
1023                 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
1024                 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
1025 #else
1026                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
1027                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
1028 #endif
1029
1030 #if DEBUG
1031                 fprintf( stderr, "succeeded\n" );
1032 #endif
1033
1034                 orlgth1 = ll1 - 100;
1035                 orlgth2 = ll2 - 100;
1036         }
1037
1038
1039         for( i=0; i<icyc; i++ )
1040         {
1041                 mseq1[i] = mseq[i];
1042                 seq1[i][lgth1] = 0;
1043         }
1044         for( j=0; j<jcyc; j++ )
1045         {
1046                 mseq2[j] = mseq[icyc+j];
1047                 seq2[j][lgth2] = 0;
1048         }
1049
1050
1051         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1052         {
1053                 int ll1, ll2;
1054
1055                 if( commonAlloc1 && commonAlloc2 )
1056                 {
1057                         FreeShortMtx( commonIP );
1058                 }
1059
1060                 ll1 = MAX( orlgth1, commonAlloc1 );
1061                 ll2 = MAX( orlgth2, commonAlloc2 );
1062
1063 #if DEBUG
1064                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1065 #endif
1066
1067                 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
1068
1069 #if DEBUG
1070                 fprintf( stderr, "succeeded\n\n" );
1071 #endif
1072
1073                 commonAlloc1 = ll1;
1074                 commonAlloc2 = ll2;
1075         }
1076         ijp = commonIP;
1077
1078 #if 0
1079         {
1080                 float t = 0.0;
1081                 for( i=0; i<icyc; i++ )
1082                         t += eff1[i];
1083         fprintf( stderr, "## totaleff = %f\n", t );
1084         }
1085 #endif
1086         cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1087         cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1088
1089         OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1090         OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1091         FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1092         FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1093
1094         for( i=0; i<lgth1; i++ )
1095         {
1096                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
1097                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
1098         }
1099         for( i=0; i<lgth2; i++ )
1100         {
1101                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
1102                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
1103         }
1104 #if 0
1105         for( i=0; i<lgth1; i++ )
1106                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1107 #endif
1108
1109         currentw = w1;
1110         previousw = w2;
1111
1112         match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1113         if( localhom )
1114                 imp_match_out_vead_tate( initverticalw, 0, lgth1 ); // 060306
1115
1116         match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1117         if( localhom )
1118                 imp_match_out_vead( currentw, 0, lgth2 ); // 060306
1119 #if 0 // -> tbfast.c
1120         if( localhom )
1121                 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1122
1123 #endif
1124
1125         if( outgap == 1 )
1126         {
1127                 for( i=1; i<lgth1+1; i++ )
1128                 {
1129                         initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1130                 }
1131                 for( j=1; j<lgth2+1; j++ )
1132                 {
1133                         currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1134                 }
1135         }
1136 #if OUTGAP0TRY
1137         else
1138         {
1139                 for( j=1; j<lgth2+1; j++ )
1140                         currentw[j] -= offset * j / 2.0;
1141                 for( i=1; i<lgth1+1; i++ )
1142                         initverticalw[i] -= offset * i / 2.0;
1143         }
1144 #endif
1145
1146         for( j=1; j<lgth2+1; ++j )
1147         {
1148                 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1149         }
1150
1151         lastverticalw[0] = currentw[lgth2-1];
1152
1153         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1154
1155 #if XXXXXXX
1156 fprintf( stderr, "currentw = \n" );
1157 for( i=0; i<lgth1+1; i++ )
1158 {
1159         fprintf( stderr, "%5.2f ", currentw[i] );
1160 }
1161 fprintf( stderr, "\n" );
1162 fprintf( stderr, "initverticalw = \n" );
1163 for( i=0; i<lgth2+1; i++ )
1164 {
1165         fprintf( stderr, "%5.2f ", initverticalw[i] );
1166 }
1167 fprintf( stderr, "\n" );
1168 fprintf( stderr, "fcgp\n" );
1169 for( i=0; i<lgth1; i++ )
1170         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1171 for( i=0; i<lgth2; i++ )
1172         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1173 #endif
1174
1175         for( i=1; i<lasti; i++ )
1176         {
1177                 wtmp = previousw;
1178                 previousw = currentw;
1179                 currentw = wtmp;
1180
1181                 previousw[0] = initverticalw[i-1];
1182
1183                 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1184 #if XXXXXXX
1185 fprintf( stderr, "\n" );
1186 fprintf( stderr, "i=%d\n", i );
1187 fprintf( stderr, "currentw = \n" );
1188 for( j=0; j<lgth2; j++ )
1189 {
1190         fprintf( stderr, "%5.2f ", currentw[j] );
1191 }
1192 fprintf( stderr, "\n" );
1193 #endif
1194                 if( localhom )
1195                 {
1196 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1197 #if  0
1198                         imp_match_out_vead( currentw, i, lgth2 );
1199 #else
1200                         imp_match_out_vead( currentw, i, lgth2 );
1201 #endif
1202                 }
1203 #if XXXXXXX
1204 fprintf( stderr, "\n" );
1205 fprintf( stderr, "i=%d\n", i );
1206 fprintf( stderr, "currentw = \n" );
1207 for( j=0; j<lgth2; j++ )
1208 {
1209         fprintf( stderr, "%5.2f ", currentw[j] );
1210 }
1211 fprintf( stderr, "\n" );
1212 #endif
1213                 currentw[0] = initverticalw[i];
1214
1215 #if 1
1216
1217                 mi = previousw[0] + ogcp2[1]; mpi = 0;
1218
1219                 ijppt = ijp[i] + 1;
1220                 mjpt = m + 1;
1221                 prept = previousw;
1222                 curpt = currentw + 1;
1223                 mpjpt = mp + 1;
1224                 for( j=1; j<lgth2+1; j++ )
1225                 {
1226                         wm = *prept;
1227                         *ijppt = 0;
1228
1229 #if 0
1230                         fprintf( stderr, "%5.0f->", wm );
1231 #endif
1232                         g = mi + fgcp2[j-1];
1233 #if 0
1234                         fprintf( stderr, "%5.0f?", g );
1235 #endif
1236                         if( g > wm )
1237                         {
1238                                 wm = g;
1239                                 *ijppt = -( j - mpi );
1240                         }
1241                         g = *prept + ogcp2[j];
1242                         if( g >= mi )
1243                         {
1244                                 mi = g;
1245                                 mpi = j-1;
1246                         }
1247 #if USE_PENALTY_EX
1248                         mi += fpenalty_ex;
1249 #endif
1250
1251                         g = *mjpt + fgcp1[i-1];
1252 #if 0
1253                         fprintf( stderr, "%5.0f?", g );
1254 #endif
1255                         if( g > wm )
1256                         {
1257                                 wm = g;
1258                                 *ijppt = +( i - *mpjpt );
1259                         }
1260                         g = *prept + ogcp1[i];
1261                         if( g >= *mjpt )
1262                         {
1263                                 *mjpt = g;
1264                                 *mpjpt = i-1;
1265                         }
1266 #if USE_PENALTY_EX
1267                         m[j] += fpenalty_ex;
1268 #endif
1269
1270 #if 0
1271                         fprintf( stderr, "%5.0f ", wm );
1272 #endif
1273                         *curpt += wm;
1274                         ijppt++;
1275                         mjpt++;
1276                         prept++;
1277                         mpjpt++;
1278                         curpt++;
1279                 }
1280                 lastverticalw[i] = currentw[lgth2-1];
1281 #else
1282         floatncpy( previousw, currentw, lgth2+1 );
1283         previousw[0] = initverticalw[i-1];
1284
1285
1286         match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1287         currentw[0] = initverticalw[i];
1288
1289         mi = previousw[0] + ogcp2[1]; mpi = 0;
1290         for( j=1; j<lgth2+1; j++ )
1291         {
1292             wm = previousw[j-1];
1293             ijp[i][j] = 0;
1294
1295             g = fgcp2[j-1];
1296             if( mi + g > wm )
1297             {
1298                 wm = mi + g;
1299                 ijp[i][j] = -( j - mpi );
1300             }
1301             g = ogcp2[j];
1302             if( mi <= previousw[j-1] + g )
1303             {
1304                 mi = previousw[j-1] + g;
1305                 mpi = j-1;
1306             }
1307
1308             g = fgcp1[i-1];
1309             if( m[j] + g > wm )
1310             {
1311                 wm = m[j] + g;
1312                 ijp[i][j] = +( i - mp[j] );
1313             }
1314             g = ogcp1[i];
1315             if( m[j] <= previousw[j-1] + g )
1316             {
1317                 m[j] = previousw[j-1] + g ;
1318                 mp[j] = i-1;
1319             }
1320             currentw[j] += wm;
1321         }
1322         lastverticalw[i] = currentw[lgth2-1];
1323
1324 #endif
1325         }
1326
1327 //      fprintf( stderr, "wm = %f\n", wm );
1328
1329 #if OUTGAP0TRY
1330         if( !outgap )
1331         {
1332                 for( j=1; j<lgth2+1; j++ )
1333                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1334                 for( i=1; i<lgth1+1; i++ )
1335                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1336         }
1337 #endif
1338
1339         /*
1340         fprintf( stderr, "\n" );
1341         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1342         fprintf( stderr, "#####\n" );
1343         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1344         fprintf( stderr, "====>" );
1345         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1346         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1347         */
1348         if( localhom )
1349         {
1350                 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1351         }
1352         else
1353                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1354
1355 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1356
1357         resultlen = strlen( mseq1[0] );
1358         if( alloclen < resultlen || resultlen > N )
1359         {
1360                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1361                 ErrorExit( "LENGTH OVER!\n" );
1362         }
1363
1364
1365         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1366         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1367         /*
1368         fprintf( stderr, "\n" );
1369         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1370         fprintf( stderr, "#####\n" );
1371         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1372         */
1373
1374 //      fprintf( stderr, "wm = %f\n", wm );
1375
1376         return( wm );
1377 }
1378
1379 float A__align_gapmap( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, int *gapmap1, int *gapmap2 )
1380 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1381 {
1382 //      int k;
1383         register int i, j;
1384         int lasti;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
1385         int lgth1, lgth2;
1386         int resultlen;
1387         float wm = 0.0;   /* int ?????? */
1388         float g;
1389         float *currentw, *previousw;
1390 //      float fpenalty = (float)penalty;
1391 #if USE_PENALTY_EX
1392         float fpenalty_ex = (float)penalty_ex;
1393 #endif
1394 #if 1
1395         float *wtmp;
1396         short *ijppt;
1397         float *mjpt, *prept, *curpt;
1398         int *mpjpt;
1399 #endif
1400         static float mi, *m;
1401         static short **ijp;
1402         static int mpi, *mp;
1403         static float *w1, *w2;
1404         static float *match;
1405         static float *initverticalw;    /* kufuu sureba iranai */
1406         static float *lastverticalw;    /* kufuu sureba iranai */
1407         static char **mseq1;
1408         static char **mseq2;
1409         static char **mseq;
1410         static float *ogcp1;
1411         static float *ogcp2;
1412         static float *fgcp1;
1413         static float *fgcp2;
1414         static float **cpmx1;
1415         static float **cpmx2;
1416         static int **intwork;
1417         static float **floatwork;
1418         static int orlgth1 = 0, orlgth2 = 0;
1419         float *fgcp2pt;
1420         float *ogcp2pt;
1421         float fgcp1va;
1422         float ogcp1va;
1423
1424
1425 #if 0
1426         fprintf( stderr, "eff in SA+++align\n" );
1427         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1428 #endif
1429         if( orlgth1 == 0 )
1430         {
1431                 mseq1 = AllocateCharMtx( njob, 0 );
1432                 mseq2 = AllocateCharMtx( njob, 0 );
1433         }
1434
1435
1436         lgth1 = strlen( seq1[0] );
1437         lgth2 = strlen( seq2[0] );
1438
1439         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
1440         {
1441                 int ll1, ll2;
1442
1443                 if( orlgth1 > 0 && orlgth2 > 0 )
1444                 {
1445                         FreeFloatVec( w1 );
1446                         FreeFloatVec( w2 );
1447                         FreeFloatVec( match );
1448                         FreeFloatVec( initverticalw );
1449                         FreeFloatVec( lastverticalw );
1450
1451                         FreeFloatVec( m );
1452                         FreeIntVec( mp );
1453
1454                         FreeCharMtx( mseq );
1455
1456                         FreeFloatVec( ogcp1 );
1457                         FreeFloatVec( ogcp2 );
1458                         FreeFloatVec( fgcp1 );
1459                         FreeFloatVec( fgcp2 );
1460
1461
1462                         FreeFloatMtx( cpmx1 );
1463                         FreeFloatMtx( cpmx2 );
1464
1465                         FreeFloatMtx( floatwork );
1466                         FreeIntMtx( intwork );
1467                 }
1468
1469                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
1470                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
1471
1472 #if DEBUG
1473                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1474 #endif
1475
1476                 w1 = AllocateFloatVec( ll2+2 );
1477                 w2 = AllocateFloatVec( ll2+2 );
1478                 match = AllocateFloatVec( ll2+2 );
1479
1480                 initverticalw = AllocateFloatVec( ll1+2 );
1481                 lastverticalw = AllocateFloatVec( ll1+2 );
1482
1483                 m = AllocateFloatVec( ll2+2 );
1484                 mp = AllocateIntVec( ll2+2 );
1485
1486                 mseq = AllocateCharMtx( njob, ll1+ll2 );
1487
1488                 ogcp1 = AllocateFloatVec( ll1+2 );
1489                 ogcp2 = AllocateFloatVec( ll2+2 );
1490                 fgcp1 = AllocateFloatVec( ll1+2 );
1491                 fgcp2 = AllocateFloatVec( ll2+2 );
1492
1493                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1494                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1495
1496 #if FASTMATCHCALC
1497                 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
1498                 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
1499 #else
1500                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
1501                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
1502 #endif
1503
1504 #if DEBUG
1505                 fprintf( stderr, "succeeded\n" );
1506 #endif
1507
1508                 orlgth1 = ll1 - 100;
1509                 orlgth2 = ll2 - 100;
1510         }
1511
1512
1513         for( i=0; i<icyc; i++ )
1514         {
1515                 mseq1[i] = mseq[i];
1516                 seq1[i][lgth1] = 0;
1517         }
1518         for( j=0; j<jcyc; j++ )
1519         {
1520                 mseq2[j] = mseq[icyc+j];
1521                 seq2[j][lgth2] = 0;
1522         }
1523
1524
1525         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1526         {
1527                 int ll1, ll2;
1528
1529                 if( commonAlloc1 && commonAlloc2 )
1530                 {
1531                         FreeShortMtx( commonIP );
1532                 }
1533
1534                 ll1 = MAX( orlgth1, commonAlloc1 );
1535                 ll2 = MAX( orlgth2, commonAlloc2 );
1536
1537 #if DEBUG
1538                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1539 #endif
1540
1541                 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
1542
1543 #if DEBUG
1544                 fprintf( stderr, "succeeded\n\n" );
1545 #endif
1546
1547                 commonAlloc1 = ll1;
1548                 commonAlloc2 = ll2;
1549         }
1550         ijp = commonIP;
1551
1552         cpmx_calc_new( seq1, cpmx1, eff1, strlen( seq1[0] ), icyc );
1553         cpmx_calc_new( seq2, cpmx2, eff2, strlen( seq2[0] ), jcyc );
1554
1555         OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1  );
1556         OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1557         FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1558         FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1559
1560         for( i=0; i<lgth1; i++ )
1561         {
1562                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] );
1563                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] );
1564         }
1565         for( i=0; i<lgth2; i++ )
1566         {
1567                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] );
1568                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] );
1569         }
1570 #if 0
1571         for( i=0; i<lgth1; i++ )
1572                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1573 #endif
1574
1575         currentw = w1;
1576         previousw = w2;
1577
1578
1579         match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1580         if( localhom )
1581                 imp_match_out_vead_tate_gapmap( initverticalw, gapmap2[0], lgth1, gapmap1 ); // 060306
1582
1583
1584         match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1585         if( localhom )
1586                 imp_match_out_vead_gapmap( currentw, gapmap1[0], lgth2, gapmap2 ); // 060306
1587 #if 0 // -> tbfast.c
1588         if( localhom )
1589                 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1590
1591 #endif
1592
1593         if( outgap == 1 )
1594         {
1595                 for( i=1; i<lgth1+1; i++ )
1596                 {
1597                         initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1598                 }
1599                 for( j=1; j<lgth2+1; j++ )
1600                 {
1601                         currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1602                 }
1603         }
1604 #if OUTGAP0TRY
1605         else
1606         {
1607                 for( j=1; j<lgth2+1; j++ )
1608                         currentw[j] -= offset * j / 2.0;
1609                 for( i=1; i<lgth1+1; i++ )
1610                         initverticalw[i] -= offset * i / 2.0;
1611         }
1612 #endif
1613
1614         for( j=1; j<lgth2+1; ++j )
1615         {
1616                 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1617         }
1618
1619         lastverticalw[0] = currentw[lgth2-1];
1620
1621         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1622
1623 #if XXXXXXX
1624 fprintf( stderr, "currentw = \n" );
1625 for( i=0; i<lgth1+1; i++ )
1626 {
1627         fprintf( stderr, "%5.2f ", currentw[i] );
1628 }
1629 fprintf( stderr, "\n" );
1630 fprintf( stderr, "initverticalw = \n" );
1631 for( i=0; i<lgth2+1; i++ )
1632 {
1633         fprintf( stderr, "%5.2f ", initverticalw[i] );
1634 }
1635 fprintf( stderr, "\n" );
1636 fprintf( stderr, "fcgp\n" );
1637 for( i=0; i<lgth1; i++ )
1638         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1639 for( i=0; i<lgth2; i++ )
1640         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1641 #endif
1642
1643         for( i=1; i<lasti; i++ )
1644         {
1645                 wtmp = previousw;
1646                 previousw = currentw;
1647                 currentw = wtmp;
1648
1649                 previousw[0] = initverticalw[i-1];
1650
1651                 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1652 #if XXXXXXX
1653 fprintf( stderr, "\n" );
1654 fprintf( stderr, "i=%d\n", i );
1655 fprintf( stderr, "currentw = \n" );
1656 for( j=0; j<lgth2; j++ )
1657 {
1658         fprintf( stderr, "%5.2f ", currentw[j] );
1659 }
1660 fprintf( stderr, "\n" );
1661 #endif
1662                 if( localhom )
1663                 {
1664 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1665 #if  0
1666                         imp_match_out_vead( currentw, i, lgth2 );
1667 #else
1668                         imp_match_out_vead_gapmap( currentw, gapmap1[i], lgth2, gapmap2 );
1669 #endif
1670                 }
1671 #if XXXXXXX
1672 fprintf( stderr, "\n" );
1673 fprintf( stderr, "i=%d\n", i );
1674 fprintf( stderr, "currentw = \n" );
1675 for( j=0; j<lgth2; j++ )
1676 {
1677         fprintf( stderr, "%5.2f ", currentw[j] );
1678 }
1679 fprintf( stderr, "\n" );
1680 #endif
1681                 currentw[0] = initverticalw[i];
1682
1683
1684                 mi = previousw[0] + ogcp2[1]; mpi = 0;
1685
1686                 ijppt = ijp[i] + 1;
1687                 mjpt = m + 1;
1688                 prept = previousw;
1689                 curpt = currentw + 1;
1690                 mpjpt = mp + 1;
1691                 fgcp2pt = fgcp2;
1692                 ogcp2pt = ogcp2 + 1;
1693                 fgcp1va = fgcp1[i-1];
1694                 ogcp1va = ogcp1[i];
1695                 for( j=1; j<lgth2+1; j++ )
1696                 {
1697                         wm = *prept;
1698                         *ijppt = 0;
1699
1700 #if 0
1701                         fprintf( stderr, "%5.0f->", wm );
1702 #endif
1703                         g = mi + *fgcp2pt;
1704 #if 0
1705                         fprintf( stderr, "%5.0f?", g );
1706 #endif
1707                         if( g > wm )
1708                         {
1709                                 wm = g;
1710                                 *ijppt = -( j - mpi );
1711                         }
1712                         g = *prept + *ogcp2pt;
1713                         if( g >= mi )
1714                         {
1715                                 mi = g;
1716                                 mpi = j-1;
1717                         }
1718 #if USE_PENALTY_EX
1719                         mi += fpenalty_ex;
1720 #endif
1721
1722                         g = *mjpt + fgcp1va;
1723 #if 0
1724                         fprintf( stderr, "%5.0f?", g );
1725 #endif
1726                         if( g > wm )
1727                         {
1728                                 wm = g;
1729                                 *ijppt = +( i - *mpjpt );
1730                         }
1731                         g = *prept + ogcp1va;
1732                         if( g >= *mjpt )
1733                         {
1734                                 *mjpt = g;
1735                                 *mpjpt = i-1;
1736                         }
1737 #if USE_PENALTY_EX
1738                         m[j] += fpenalty_ex;
1739 #endif
1740
1741 #if 0
1742                         fprintf( stderr, "%5.0f ", wm );
1743 #endif
1744                         *curpt += wm;
1745                         ijppt++;
1746                         mjpt++;
1747                         prept++;
1748                         mpjpt++;
1749                         curpt++;
1750                         fgcp2pt++;
1751                         ogcp2pt++;
1752                 }
1753                 lastverticalw[i] = currentw[lgth2-1];
1754         }
1755
1756 #if OUTGAP0TRY
1757         if( !outgap )
1758         {
1759                 for( j=1; j<lgth2+1; j++ )
1760                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1761                 for( i=1; i<lgth1+1; i++ )
1762                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1763         }
1764 #endif
1765
1766         /*
1767         fprintf( stderr, "\n" );
1768         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1769         fprintf( stderr, "#####\n" );
1770         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1771         fprintf( stderr, "====>" );
1772         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1773         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1774         */
1775         if( localhom )
1776         {
1777                 Atracking_localhom_gapmap( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, gapmap1, gapmap2 );
1778         }
1779         else
1780                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1781
1782 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1783
1784         resultlen = strlen( mseq1[0] );
1785         if( alloclen < resultlen || resultlen > N )
1786         {
1787                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1788                 ErrorExit( "LENGTH OVER!\n" );
1789         }
1790
1791
1792         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1793         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1794         /*
1795         fprintf( stderr, "\n" );
1796         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1797         fprintf( stderr, "#####\n" );
1798         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1799         */
1800
1801 //      fprintf( stderr, "wm = %f\n", wm );
1802
1803         return( wm );
1804 }
1805
1806 float translate_and_Calign( char **mseq1, char **mseq2, double *effarr1, double *effarr2, int clus1, int clus2, int alloclen )
1807 {
1808     int i;
1809     float wm;
1810     char **result;
1811     char *seq, **aseq;
1812     double *effarr;
1813     int nseq;
1814         int resultlen;
1815
1816     if     ( clus1 == 1 && clus2 != 1 )
1817     {
1818         seq = mseq1[0]; aseq = mseq2; effarr = effarr2; nseq = clus2+1;
1819 #if 0
1820                 printf( "effarr in transl... = \n" );
1821                 for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
1822 #endif
1823     }
1824     else if( clus1 != 1 && clus2 == 1 )
1825     {
1826         seq = mseq2[0]; aseq = mseq1; effarr = effarr1; nseq = clus1+1;
1827     }
1828     else ErrorExit( "ERROR in translate_and_Calign" );
1829
1830     result = Calignm1( &wm, aseq, seq, effarr, nseq-2, 0 );
1831
1832         resultlen = strlen( result[0] );
1833         if( alloclen < resultlen || resultlen > N )
1834         {
1835                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1836                 ErrorExit( "LENGTH OVER!\n" );
1837         }
1838     for( i=0; i<nseq-1; i++ ) strcpy( aseq[i], result[i] );
1839     strcpy( seq, result[nseq-1] );
1840
1841     return( 0.0 );
1842 }