binaries/src/mafft/core/partSalignmm.c

   1 #include "mltaln.h"
   2 #include "dp.h"
   3
   4 #define MACHIGAI 0
   5 #define OUTGAP0TRY 1
   6 #define DEBUG 0
   7 #define XXXXXXX    0
   8 #define USE_PENALTY_EX  0
   9 #define FASTMATCHCALC 1
  10
  11 #if 0
  12 static void st_OpeningGapCount( float *ogcp, int clus, char **seq, double *eff, int len )
  13 {
  14         int i, j, gc, gb;
  15         float feff;
  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         }
  31 }
  32
  33 static void st_FinalGapCount( float *fgcp, int clus, char **seq, double *eff, int len )
  34 {
  35         int i, j, gc, gb;
  36         float feff;
  37
  38         for( i=0; i<len; i++ ) fgcp[i] = 0.0;
  39         for( j=0; j<clus; j++ )
  40         {
  41                 feff = (float)eff[j];
  42                 gc = ( seq[j][0] == '-' );
  43                 for( i=1; i<len+1; i++ )
  44                 {
  45                         gb = gc;
  46                         gc = ( seq[j][i] == '-' );
  47                         {
  48                                 if( gb * !gc ) fgcp[i-1] += feff;
  49                         }
  50                 }
  51         }
  52 }
  53 #endif
  54
  55
  56
  57
  58
  59 static int impalloclen = 0;
  60 static float **impmtx = NULL;
  61 float part_imp_match_out_sc( int i1, int j1 )
  62 {
  63 //      fprintf( stderr, "impalloclen = %d\n", impalloclen );
  64 //      fprintf( stderr, "i1,j1=%d,%d -> impmtx=%f\n", i1, j1, impmtx[i1][j1] );
  65         return( impmtx[i1][j1] );
  66 #if 0
  67         if( i1 == l1 || j1 == l2 ) return( 0.0 );
  68         return( impmtx[i1+start1][j1+start2] );
  69 #endif
  70 }
  71 static void part_imp_match_out_vead_gapmap( float *imp, int i1, int lgth2, int start2, int *gapmap2 )
  72 {
  73 #if FASTMACHCALC
  74         float *pt = imp;
  75         int *gapmappt = gapmap2;
  76         while( lgth2-- )
  77                 *pt++ += impmtx[i1][start2+*gapmappt++];
  78 #else
  79         int j;
  80         for( j=0; j<lgth2; j++ )
  81         {
  82                 imp[j] += impmtx[i1][start2+gapmap2[j]];
  83         }
  84 #endif
  85 }
  86
  87 static void part_imp_match_out_vead_tate_gapmap( float *imp, int j1, int lgth1, int start1, int *gapmap1 )
  88 {
  89 #if FASTMACHCALC
  90         float *pt = imp;
  91         int *gapmappt = gapmap1;
  92         while( lgth1-- )
  93                 *pt++ = impmtx[start1+*gapmappt++][j1];
  94 #else
  95         int i;
  96         for( i=0; i<lgth1; i++ )
  97         {
  98                 imp[i] += impmtx[start1+gapmap1[i]][j1];
  99         }
 100 #endif
 101 }
 102
 103 void part_imp_match_init_strict( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, double *eff1_kozo, double *eff2_kozo, LocalHom ***localhom, int forscore )
 104 {
 105         int i, j, k1, k2, tmpint, start1, start2, end1, end2;
 106         double effij, effijx, effij_kozo;
 107         char *pt, *pt1, *pt2;
 108         LocalHom *tmpptr;
 109
 110         if( impalloclen <= lgth1 + 2 || impalloclen <= lgth2 + 2 )
 111         {
 112                 if( impmtx ) FreeFloatMtx( impmtx );
 113                 impalloclen = MAX( lgth1, lgth2 ) + 2;
 114                 impmtx = AllocateFloatMtx( impalloclen+100, impalloclen+100 );
 115         }
 116
 117
 118 #if 0
 119         fprintf( stderr, "eff1 in _init_strict = \n" );
 120         for( i=0; i<clus1; i++ )
 121                 fprintf( stderr, "eff1[] = %f\n", eff1[i] );
 122         for( i=0; i<clus2; i++ )
 123                 fprintf( stderr, "eff2[] = %f\n", eff2[i] );
 124 #endif
 125
 126         for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
 127                 impmtx[i][j] = 0.0;
 128         effijx = 1.0 * fastathreshold;
 129         for( i=0; i<clus1; i++ )
 130         {
 131                 for( j=0; j<clus2; j++ )
 132                 {
 133                         effij = eff1[i] * eff2[j] * effijx;
 134                         effij_kozo = eff1_kozo[i] * eff2_kozo[j] * effijx;
 135                         tmpptr = localhom[i][j];
 136                         while( tmpptr )
 137                         {
 138 //                              fprintf( stderr, "start1 = %d\n", tmpptr->start1 );
 139 //                              fprintf( stderr, "end1   = %d\n", tmpptr->end1   );
 140 //                              fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] );
 141 //                              fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] );
 142                                 pt = seq1[i];
 143                                 tmpint = -1;
 144                                 while( *pt != 0 )
 145                                 {
 146                                         if( *pt++ != '-' ) tmpint++;
 147                                         if( tmpint == tmpptr->start1 ) break;
 148                                 }
 149                                 start1 = (int)( pt - seq1[i] ) - 1;
 150
 151                                 if( tmpptr->start1 == tmpptr->end1 ) end1 = start1;
 152                                 else
 153                                 {
 154 #if MACHIGAI
 155                                         while( *pt != 0 )
 156                                         {
 157                                                 if( tmpint == tmpptr->end1 ) break;
 158                                                 if( *pt++ != '-' ) tmpint++;
 159                                         }
 160                                         end1 = (int)( pt - seq1[i] ) - 1;
 161 #else
 162                                         while( *pt != 0 )
 163                                         {
 164 //                                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
 165                                                 if( *pt++ != '-' ) tmpint++;
 166                                                 if( tmpint == tmpptr->end1 ) break;
 167                                         }
 168                                         end1 = (int)( pt - seq1[i] ) - 1;
 169 #endif
 170                                 }
 171
 172                                 pt = seq2[j];
 173                                 tmpint = -1;
 174                                 while( *pt != 0 )
 175                                 {
 176                                         if( *pt++ != '-' ) tmpint++;
 177                                         if( tmpint == tmpptr->start2 ) break;
 178                                 }
 179                                 start2 = (int)( pt - seq2[j] ) - 1;
 180                                 if( tmpptr->start2 == tmpptr->end2 ) end2 = start2;
 181                                 else
 182                                 {
 183 #if MACHIGAI
 184                                         while( *pt != 0 )
 185                                         {
 186                                                 if( tmpint == tmpptr->end2 ) break;
 187                                                 if( *pt++ != '-' ) tmpint++;
 188                                         }
 189                                         end2 = (int)( pt - seq2[j] ) - 1;
 190 #else
 191                                         while( *pt != 0 )
 192                                         {
 193                                                 if( *pt++ != '-' ) tmpint++;
 194                                                 if( tmpint == tmpptr->end2 ) break;
 195                                         }
 196                                         end2 = (int)( pt - seq2[j] ) - 1;
 197 #endif
 198                                 }
 199 //                              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] );
 200 //                              fprintf( stderr, "step 0\n" );
 201                                 if( end1 - start1 != end2 - start2 )
 202                                 {
 203 //                                      fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
 204                                 }
 205
 206                                 k1 = start1; k2 = start2;
 207                                 pt1 = seq1[i] + k1;
 208                                 pt2 = seq2[j] + k2;
 209                                 while( *pt1 && *pt2 )
 210                                 {
 211                                         if( *pt1 != '-' && *pt2 != '-' )
 212                                         {
 213 // ½Å¤ß¤òÆó½Å¤Ë¤«¤±¤Ê¤¤¤è¤¦¤ËÃí°Õ¤·¤Æ²¼¤µ¤¤¡£
 214 //                                              impmtx[k1][k2] += tmpptr->wimportance * fastathreshold;
 215 //                                              impmtx[k1][k2] += tmpptr->importance * effij;
 216 //                                              impmtx[k1][k2] += tmpptr->fimportance * effij;
 217                                                 if( tmpptr->korh == 'k' )
 218                                                         impmtx[k1][k2] += tmpptr->fimportance * effij_kozo;
 219                                                 else
 220                                                         impmtx[k1][k2] += tmpptr->fimportance * effij;
 221 //                                              fprintf( stderr, "k1=%d, k2=%d, impalloclen=%d\n", k1, k2, impalloclen );
 222 //                                              fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 223                                                 k1++; k2++;
 224                                                 pt1++; pt2++;
 225                                         }
 226                                         else if( *pt1 != '-' && *pt2 == '-' )
 227                                         {
 228 //                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 229                                                 k2++; pt2++;
 230                                         }
 231                                         else if( *pt1 == '-' && *pt2 != '-' )
 232                                         {
 233 //                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 234                                                 k1++; pt1++;
 235                                         }
 236                                         else if( *pt1 == '-' && *pt2 == '-' )
 237                                         {
 238 //                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 239                                                 k1++; pt1++;
 240                                                 k2++; pt2++;
 241                                         }
 242                                         if( k1 > end1 || k2 > end2 ) break;
 243                                 }
 244                                 tmpptr = tmpptr->next;
 245                         }
 246                 }
 247         }
 248 #if 0
 249         fprintf( stderr, "impmtx = \n" );
 250         for( k2=0; k2<lgth2; k2++ )
 251                 fprintf( stderr, "%6.3f ", (double)k2 );
 252         fprintf( stderr, "\n" );
 253         for( k1=0; k1<lgth1; k1++ )
 254         {
 255                 fprintf( stderr, "%d", k1 );
 256                 for( k2=0; k2<lgth2; k2++ )
 257                         fprintf( stderr, "%2.1f ", impmtx[k1][k2] );
 258                 fprintf( stderr, "\n" );
 259         }
 260         exit( 1 );
 261 #endif
 262 }
 263
 264
 265 void part_imp_rna( int nseq1, int nseq2, char **seq1, char **seq2, double *eff1, double *eff2, RNApair ***grouprna1, RNApair ***grouprna2, int *gapmap1, int *gapmap2, RNApair *additionalpair )
 266 {
 267         foldrna( nseq1, nseq2, seq1, seq2, eff1, eff2, grouprna1, grouprna2, impmtx, gapmap1, gapmap2, additionalpair );
 268 }
 269
 270
 271 void part_imp_match_init( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom )
 272 {
 273         int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2;
 274         static int impalloclen = 0;
 275         char *pt;
 276         static char *nocount1 = NULL;
 277         static char *nocount2 = NULL;
 278
 279         if( impalloclen < lgth1 || impalloclen < lgth2 )
 280         {
 281                 if( impmtx ) FreeFloatMtx( impmtx );
 282                 if( nocount1 ) free( nocount1 );
 283                 if( nocount2 ) free( nocount2 );
 284                 impalloclen = MAX( lgth1, lgth2 ) + 2;
 285                 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
 286                 nocount1 = AllocateCharVec( impalloclen );
 287                 nocount2 = AllocateCharVec( impalloclen );
 288                 impalloclen -= 2;
 289         }
 290
 291         for( i=0; i<lgth1; i++ )
 292         {
 293                 for( j=0; j<clus1; j++ )
 294                         if( seq1[j][i] == '-' ) break;
 295                 if( j != clus1 ) nocount1[i] = 1;
 296                 else                     nocount1[i] = 0;
 297         }
 298         for( i=0; i<lgth2; i++ )
 299         {
 300                 for( j=0; j<clus2; j++ )
 301                         if( seq2[j][i] == '-' ) break;
 302                 if( j != clus2 ) nocount2[i] = 1;
 303                 else                     nocount2[i] = 0;
 304         }
 305
 306 #if 0
 307 fprintf( stderr, "nocount2 =\n" );
 308 for( i = 0; i<impalloclen; i++ )
 309 {
 310         fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
 311 }
 312 #endif
 313
 314         for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
 315                 impmtx[i][j] = 0.0;
 316         for( i=0; i<clus1; i++ )
 317         {
 318                 fprintf( stderr, "i = %d, seq1 = %s\n", i, seq1[i] );
 319                 for( j=0; j<clus2; j++ )
 320                 {
 321                         fprintf( stderr, "start1 = %d\n", localhom[i][j]->start1 );
 322                         fprintf( stderr, "end1   = %d\n", localhom[i][j]->end1   );
 323                         fprintf( stderr, "j = %d, seq2 = %s\n", j, seq2[j] );
 324                         pt = seq1[i];
 325                         tmpint = -1;
 326                         while( *pt != 0 )
 327                         {
 328                                 if( *pt++ != '-' ) tmpint++;
 329                                 if( tmpint == localhom[i][j]->start1 ) break;
 330                         }
 331                         start1 = pt - seq1[i] - 1;
 332
 333                         while( *pt != 0 )
 334                         {
 335 //                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, localhom[i][j].end1, pt-seq1[i] );
 336                                 if( *pt++ != '-' ) tmpint++;
 337                                 if( tmpint == localhom[i][j]->end1 ) break;
 338                         }
 339                         end1 = pt - seq1[i] - 1;
 340
 341                         pt = seq2[j];
 342                         tmpint = -1;
 343                         while( *pt != 0 )
 344                         {
 345                                 if( *pt++ != '-' ) tmpint++;
 346                                 if( tmpint == localhom[i][j]->start2 ) break;
 347                         }
 348                         start2 = pt - seq2[j] - 1;
 349                         while( *pt != 0 )
 350                         {
 351                                 if( *pt++ != '-' ) tmpint++;
 352                                 if( tmpint == localhom[i][j]->end2 ) break;
 353                         }
 354                         end2 = pt - seq2[j] - 1;
 355 //                      fprintf( stderr, "start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
 356                         k1 = start1;
 357                         k2 = start2;
 358                         fprintf( stderr, "step 0\n" );
 359                         while( k1 <= end1 && k2 <= end2 )
 360                         {
 361 #if 0
 362                                 if( !nocount1[k1] && !nocount2[k2] )
 363                                         impmtx[k1][k2] += localhom[i][j].wimportance * eff1[i] * eff2[j];
 364                                 k1++; k2++;
 365 #else
 366                                 if( !nocount1[k1] && !nocount2[k2] )
 367                                         impmtx[k1][k2] += localhom[i][j]->wimportance * eff1[i] * eff2[j];
 368                                 k1++; k2++;
 369 #endif
 370                         }
 371
 372                         dif = ( end1 - start1 ) - ( end2 - start2 );
 373                         fprintf( stderr, "dif = %d\n", dif );
 374                         if( dif > 0 )
 375                         {
 376                                 do
 377                                 {
 378                                         fprintf( stderr, "dif = %d\n", dif );
 379                                         k1 = start1;
 380                                         k2 = start2 - dif;
 381                                         while( k1 <= end1 && k2 <= end2 )
 382                                         {
 383                                                 if( 0 <= k2 && start2 <= k2 && !nocount1[k1] && !nocount2[k2] )
 384                                                         impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
 385                                                 k1++; k2++;
 386                                         }
 387                                 }
 388                                 while( dif-- );
 389                         }
 390                         else
 391                         {
 392                                 do
 393                                 {
 394                                         k1 = start1 + dif;
 395                                         k2 = start2;
 396                                         while( k1 <= end1 )
 397                                         {
 398                                                 if( k1 >= 0 && k1 >= start1 && !nocount1[k1] && !nocount2[k2] )
 399                                                         impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
 400                                                 k1++; k2++;
 401                                         }
 402                                 }
 403                                 while( dif++ );
 404                         }
 405                 }
 406         }
 407 #if 0
 408         fprintf( stderr, "impmtx = \n" );
 409         for( k2=0; k2<lgth2; k2++ )
 410                 fprintf( stderr, "%6.3f ", (double)k2 );
 411         fprintf( stderr, "\n" );
 412         for( k1=0; k1<lgth1; k1++ )
 413         {
 414                 fprintf( stderr, "%d", k1 );
 415                 for( k2=0; k2<lgth2; k2++ )
 416                         fprintf( stderr, "%6.3f ", impmtx[k1][k2] );
 417                 fprintf( stderr, "\n" );
 418         }
 419         exit( 1 );
 420 #endif
 421 }
 422
 423 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
 424 {
 425 #if FASTMATCHCALC
 426         int j, l;
 427         float scarr[26];
 428         float **cpmxpd = floatwork;
 429         int **cpmxpdn = intwork;
 430         float *matchpt, *cpmxpdpt, **cpmxpdptpt;
 431         int *cpmxpdnpt, **cpmxpdnptpt;
 432         if( initialize )
 433         {
 434                 int count = 0;
 435                 for( j=0; j<lgth2; j++ )
 436                 {
 437                         count = 0;
 438                         for( l=0; l<26; l++ )
 439                         {
 440                                 if( cpmx2[l][j] )
 441                                 {
 442                                         cpmxpd[j][count] = cpmx2[l][j];
 443                                         cpmxpdn[j][count] = l;
 444                                         count++;
 445                                 }
 446                         }
 447                         cpmxpdn[j][count] = -1;
 448                 }
 449         }
 450
 451         {
 452                 for( l=0; l<26; l++ )
 453                 {
 454                         scarr[l] = 0.0;
 455                         for( j=0; j<26; j++ )
 456                                 scarr[l] += n_dis_consweight_multi[j][l] * cpmx1[j][i1];
 457 //                              scarr[l] += n_dis[j][l] * cpmx1[j][i1];
 458                 }
 459                 matchpt = match;
 460                 cpmxpdnptpt = cpmxpdn;
 461                 cpmxpdptpt = cpmxpd;
 462                 while( lgth2-- )
 463                 {
 464                         *matchpt = 0.0;
 465                         cpmxpdnpt = *cpmxpdnptpt++;
 466                         cpmxpdpt = *cpmxpdptpt++;
 467                         while( *cpmxpdnpt>-1 )
 468                                 *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;
 469                         matchpt++;
 470                 }
 471         }
 472 #else
 473         int j, k, l;
 474         float scarr[26];
 475         float **cpmxpd = floatwork;
 476         int **cpmxpdn = intwork;
 477         // simple
 478         if( initialize )
 479         {
 480                 int count = 0;
 481                 for( j=0; j<lgth2; j++ )
 482                 {
 483                         count = 0;
 484                         for( l=0; l<26; l++ )
 485                         {
 486                                 if( cpmx2[l][j] )
 487                                 {
 488                                         cpmxpd[count][j] = cpmx2[l][j];
 489                                         cpmxpdn[count][j] = l;
 490                                         count++;
 491                                 }
 492                         }
 493                         cpmxpdn[count][j] = -1;
 494                 }
 495         }
 496         for( l=0; l<26; l++ )
 497         {
 498                 scarr[l] = 0.0;
 499                 for( k=0; k<26; k++ )
 500                         scarr[l] += n_dis_consweight_multi[k][l] * cpmx1[k][i1];
 501 //                      scarr[l] += n_dis[k][l] * cpmx1[k][i1];
 502         }
 503         for( j=0; j<lgth2; j++ )
 504         {
 505                 match[j] = 0.0;
 506                 for( k=0; cpmxpdn[k][j]>-1; k++ )
 507                         match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
 508         }
 509 #endif
 510 }
 511
 512 static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
 513                                                 char **seq1, char **seq2,
 514                         char **mseq1, char **mseq2,
 515                         float **cpmx1, float **cpmx2,
 516                         int **ijp, int icyc, int jcyc,
 517                                                 int start1, int end1, int start2, int end2,
 518                                                 int *gapmap1, int *gapmap2 )
 519 {
 520         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
 521         char gap[] = "-";
 522         float wm;
 523         lgth1 = strlen( seq1[0] );
 524         lgth2 = strlen( seq2[0] );
 525
 526 #if 0
 527         for( i=0; i<lgth1; i++ )
 528         {
 529                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 530         }
 531 #endif
 532
 533         if( outgap == 1 )
 534                 ;
 535         else
 536         {
 537                 wm = lastverticalw[0];
 538                 for( i=0; i<lgth1; i++ )
 539                 {
 540                         if( lastverticalw[i] >= wm )
 541                         {
 542                                 wm = lastverticalw[i];
 543                                 iin = i; jin = lgth2-1;
 544                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 545                         }
 546                 }
 547                 for( j=0; j<lgth2; j++ )
 548                 {
 549                         if( lasthorizontalw[j] >= wm )
 550                         {
 551                                 wm = lasthorizontalw[j];
 552                                 iin = lgth1-1; jin = j;
 553                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 554                         }
 555                 }
 556         }
 557
 558     for( i=0; i<lgth1+1; i++ )
 559     {
 560         ijp[i][0] = i + 1;
 561     }
 562     for( j=0; j<lgth2+1; j++ )
 563     {
 564         ijp[0][j] = -( j + 1 );
 565     }
 566
 567         for( i=0; i<icyc; i++ )
 568         {
 569                 mseq1[i] += lgth1+lgth2;
 570                 *mseq1[i] = 0;
 571         }
 572         for( j=0; j<jcyc; j++ )
 573         {
 574                 mseq2[j] += lgth1+lgth2;
 575                 *mseq2[j] = 0;
 576         }
 577         iin = lgth1; jin = lgth2;
 578         *impwmpt = 0.0;
 579         for( k=0; k<=lgth1+lgth2; k++ )
 580         {
 581                 if( ijp[iin][jin] < 0 )
 582                 {
 583                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 584                 }
 585                 else if( ijp[iin][jin] > 0 )
 586                 {
 587                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 588                 }
 589                 else
 590                 {
 591                         ifi = iin-1; jfi = jin-1;
 592                 }
 593                 l = iin - ifi;
 594                 while( --l )
 595                 {
 596                         for( i=0; i<icyc; i++ )
 597                                 *--mseq1[i] = seq1[i][ifi+l];
 598                         for( j=0; j<jcyc; j++ )
 599                                 *--mseq2[j] = *gap;
 600                         k++;
 601                 }
 602                 l= jin - jfi;
 603                 while( --l )
 604                 {
 605                         for( i=0; i<icyc; i++ )
 606                                 *--mseq1[i] = *gap;
 607                         for( j=0; j<jcyc; j++ )
 608                                 *--mseq2[j] = seq2[j][jfi+l];
 609                         k++;
 610                 }
 611                 if( iin != lgth1 && jin != lgth2 ) // ??
 612                 {
 613                         *impwmpt += part_imp_match_out_sc( gapmap1[iin]+start1, gapmap2[jin]+start2 );
 614 //                      fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
 615                 }
 616                 if( iin <= 0 || jin <= 0 ) break;
 617                 for( i=0; i<icyc; i++ )
 618                         *--mseq1[i] = seq1[i][ifi];
 619                 for( j=0; j<jcyc; j++ )
 620                         *--mseq2[j] = seq2[j][jfi];
 621                 k++;
 622                 iin = ifi; jin = jfi;
 623         }
 624 }
 625 static float Atracking( float *lasthorizontalw, float *lastverticalw,
 626                                                 char **seq1, char **seq2,
 627                         char **mseq1, char **mseq2,
 628                         float **cpmx1, float **cpmx2,
 629                         int **ijp, int icyc, int jcyc )
 630 {
 631         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, lastk;
 632         char gap[] = "-";
 633         float wm = 0.0;
 634         lgth1 = strlen( seq1[0] );
 635         lgth2 = strlen( seq2[0] );
 636
 637 #if 0
 638         for( i=0; i<lgth1; i++ )
 639         {
 640                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 641         }
 642 #endif
 643
 644         if( outgap == 1 )
 645                 ;
 646         else
 647         {
 648                 wm = lastverticalw[0];
 649                 for( i=0; i<lgth1; i++ )
 650                 {
 651                         if( lastverticalw[i] >= wm )
 652                         {
 653                                 wm = lastverticalw[i];
 654                                 iin = i; jin = lgth2-1;
 655                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 656                         }
 657                 }
 658                 for( j=0; j<lgth2; j++ )
 659                 {
 660                         if( lasthorizontalw[j] >= wm )
 661                         {
 662                                 wm = lasthorizontalw[j];
 663                                 iin = lgth1-1; jin = j;
 664                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 665                         }
 666                 }
 667         }
 668
 669     for( i=0; i<lgth1+1; i++ )
 670     {
 671         ijp[i][0] = i + 1;
 672     }
 673     for( j=0; j<lgth2+1; j++ )
 674     {
 675         ijp[0][j] = -( j + 1 );
 676     }
 677
 678         for( i=0; i<icyc; i++ )
 679         {
 680                 mseq1[i] += lgth1+lgth2;
 681                 *mseq1[i] = 0;
 682         }
 683         for( j=0; j<jcyc; j++ )
 684         {
 685                 mseq2[j] += lgth1+lgth2;
 686                 *mseq2[j] = 0;
 687         }
 688         iin = lgth1; jin = lgth2;
 689         lastk = lgth1+lgth2;
 690         for( k=0; k<=lastk; k++ )
 691         {
 692                 if( ijp[iin][jin] < 0 )
 693                 {
 694                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 695                 }
 696                 else if( ijp[iin][jin] > 0 )
 697                 {
 698                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 699                 }
 700                 else
 701                 {
 702                         ifi = iin-1; jfi = jin-1;
 703                 }
 704                 l = iin - ifi;
 705                 while( --l )
 706                 {
 707                         for( i=0; i<icyc; i++ )
 708                                 *--mseq1[i] = seq1[i][ifi+l];
 709                         for( j=0; j<jcyc; j++ )
 710                                 *--mseq2[j] = *gap;
 711                         k++;
 712                 }
 713                 l= jin - jfi;
 714                 while( --l )
 715                 {
 716                         for( i=0; i<icyc; i++ )
 717                                 *--mseq1[i] = *gap;
 718                         for( j=0; j<jcyc; j++ )
 719                                 *--mseq2[j] = seq2[j][jfi+l];
 720                         k++;
 721                 }
 722                 if( iin <= 0 || jin <= 0 ) break;
 723                 for( i=0; i<icyc; i++ )
 724                         *--mseq1[i] = seq1[i][ifi];
 725                 for( j=0; j<jcyc; j++ )
 726                         *--mseq2[j] = seq2[j][jfi];
 727                 k++;
 728                 iin = ifi; jin = jfi;
 729         }
 730         return( 0.0 );
 731 }
 732
 733 float partA__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 )
 734 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 735 {
 736 //      int k;
 737         register int i, j;
 738         int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
 739         int lgth1, lgth2;
 740         int resultlen;
 741         float wm = 0.0;   /* int ?????? */
 742         float g;
 743         float *currentw, *previousw;
 744 #if 1
 745         float *wtmp;
 746         int *ijppt;
 747         float *mjpt, *prept, *curpt;
 748         int *mpjpt;
 749 #endif
 750         static float mi, *m;
 751         static int **ijp;
 752         static int mpi, *mp;
 753         static float *w1, *w2;
 754         static float *match;
 755         static float *initverticalw;    /* kufuu sureba iranai */
 756         static float *lastverticalw;    /* kufuu sureba iranai */
 757         static char **mseq1;
 758         static char **mseq2;
 759         static char **mseq;
 760         static float *ogcp1;
 761         static float *ogcp2;
 762         static float *fgcp1;
 763         static float *fgcp2;
 764         static float **cpmx1;
 765         static float **cpmx2;
 766         static int **intwork;
 767         static float **floatwork;
 768         static int orlgth1 = 0, orlgth2 = 0;
 769         float fpenalty = (float)penalty;
 770 #if USE_PENALTY_EX
 771         float fpenalty_ex = (float)penalty_ex;
 772 #endif
 773         float *fgcp2pt;
 774         float *ogcp2pt;
 775         float fgcp1va;
 776         float ogcp1va;
 777
 778
 779
 780 #if 0
 781         fprintf( stderr, "eff in SA+++align\n" );
 782         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
 783 #endif
 784         if( orlgth1 == 0 )
 785         {
 786                 mseq1 = AllocateCharMtx( njob, 0 );
 787                 mseq2 = AllocateCharMtx( njob, 0 );
 788         }
 789
 790
 791         lgth1 = strlen( seq1[0] );
 792         lgth2 = strlen( seq2[0] );
 793
 794         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 795         {
 796                 int ll1, ll2;
 797
 798                 if( orlgth1 > 0 && orlgth2 > 0 )
 799                 {
 800                         FreeFloatVec( w1 );
 801                         FreeFloatVec( w2 );
 802                         FreeFloatVec( match );
 803                         FreeFloatVec( initverticalw );
 804                         FreeFloatVec( lastverticalw );
 805
 806                         FreeFloatVec( m );
 807                         FreeIntVec( mp );
 808
 809                         FreeCharMtx( mseq );
 810
 811                         FreeFloatVec( ogcp1 );
 812                         FreeFloatVec( ogcp2 );
 813                         FreeFloatVec( fgcp1 );
 814                         FreeFloatVec( fgcp2 );
 815
 816
 817                         FreeFloatMtx( cpmx1 );
 818                         FreeFloatMtx( cpmx2 );
 819
 820                         FreeFloatMtx( floatwork );
 821                         FreeIntMtx( intwork );
 822                 }
 823
 824                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 825                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 826
 827 #if DEBUG
 828                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 829 #endif
 830
 831                 w1 = AllocateFloatVec( ll2+2 );
 832                 w2 = AllocateFloatVec( ll2+2 );
 833                 match = AllocateFloatVec( ll2+2 );
 834
 835                 initverticalw = AllocateFloatVec( ll1+2 );
 836                 lastverticalw = AllocateFloatVec( ll1+2 );
 837
 838                 m = AllocateFloatVec( ll2+2 );
 839                 mp = AllocateIntVec( ll2+2 );
 840
 841                 mseq = AllocateCharMtx( njob, ll1+ll2 );
 842
 843                 ogcp1 = AllocateFloatVec( ll1+2 );
 844                 ogcp2 = AllocateFloatVec( ll2+2 );
 845                 fgcp1 = AllocateFloatVec( ll1+2 );
 846                 fgcp2 = AllocateFloatVec( ll2+2 );
 847
 848                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
 849                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
 850
 851 #if FASTMATCHCALC
 852                 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
 853                 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
 854 #else
 855                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
 856                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
 857 #endif
 858
 859 #if DEBUG
 860                 fprintf( stderr, "succeeded\n" );
 861 #endif
 862
 863                 orlgth1 = ll1 - 100;
 864                 orlgth2 = ll2 - 100;
 865         }
 866
 867
 868         for( i=0; i<icyc; i++ ) mseq1[i] = mseq[i];
 869         for( j=0; j<jcyc; j++ ) mseq2[j] = mseq[icyc+j];
 870
 871
 872         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
 873         {
 874                 int ll1, ll2;
 875
 876                 if( commonAlloc1 && commonAlloc2 )
 877                 {
 878                         FreeIntMtx( commonIP );
 879                 }
 880
 881                 ll1 = MAX( orlgth1, commonAlloc1 );
 882                 ll2 = MAX( orlgth2, commonAlloc2 );
 883
 884 #if DEBUG
 885                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
 886 #endif
 887
 888                 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
 889
 890 #if DEBUG
 891                 fprintf( stderr, "succeeded\n\n" );
 892 #endif
 893
 894                 commonAlloc1 = ll1;
 895                 commonAlloc2 = ll2;
 896         }
 897         ijp = commonIP;
 898
 899         cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
 900         cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
 901
 902         if( sgap1 )
 903         {
 904                 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
 905                 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
 906                 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap1 );
 907                 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
 908         }
 909         else
 910         {
 911                 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
 912                 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
 913                 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
 914                 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
 915         }
 916
 917         for( i=0; i<lgth1; i++ )
 918         {
 919                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
 920                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
 921         }
 922         for( i=0; i<lgth2; i++ )
 923         {
 924                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
 925                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
 926         }
 927 #if 0
 928         for( i=0; i<lgth1; i++ )
 929                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
 930 #endif
 931
 932         currentw = w1;
 933         previousw = w2;
 934
 935
 936         match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
 937         if( localhom )
 938                 part_imp_match_out_vead_tate_gapmap( initverticalw, gapmap2[0]+start2, lgth1, start1, gapmap1 );
 939
 940
 941         match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
 942         if( localhom )
 943                 part_imp_match_out_vead_gapmap( currentw, gapmap1[0]+start1, lgth2, start2, gapmap2 );
 944 #if 0 // -> tbfast.c
 945         if( localhom )
 946                 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
 947
 948 #endif
 949
 950         if( outgap == 1 )
 951         {
 952                 for( i=1; i<lgth1+1; i++ )
 953                 {
 954                         initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
 955                 }
 956                 for( j=1; j<lgth2+1; j++ )
 957                 {
 958                         currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
 959                 }
 960         }
 961 #if OUTGAP0TRY
 962         else
 963         {
 964                 for( j=1; j<lgth2+1; j++ )
 965                         currentw[j] -= offset * j / 2.0;
 966                 for( i=1; i<lgth1+1; i++ )
 967                         initverticalw[i] -= offset * i / 2.0;
 968         }
 969 #endif
 970
 971         for( j=1; j<lgth2+1; ++j )
 972         {
 973                 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
 974         }
 975
 976         lastverticalw[0] = currentw[lgth2-1];
 977
 978         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
 979         lastj = lgth2+1;
 980
 981 #if XXXXXXX
 982 fprintf( stderr, "currentw = \n" );
 983 for( i=0; i<lgth1+1; i++ )
 984 {
 985         fprintf( stderr, "%5.2f ", currentw[i] );
 986 }
 987 fprintf( stderr, "\n" );
 988 fprintf( stderr, "initverticalw = \n" );
 989 for( i=0; i<lgth2+1; i++ )
 990 {
 991         fprintf( stderr, "%5.2f ", initverticalw[i] );
 992 }
 993 fprintf( stderr, "\n" );
 994 fprintf( stderr, "fcgp\n" );
 995 for( i=0; i<lgth1; i++ )
 996         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
 997 for( i=0; i<lgth2; i++ )
 998         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
 999 #endif
1000
1001         for( i=1; i<lasti; i++ )
1002         {
1003                 wtmp = previousw;
1004                 previousw = currentw;
1005                 currentw = wtmp;
1006
1007                 previousw[0] = initverticalw[i-1];
1008
1009                 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1010 #if XXXXXXX
1011 fprintf( stderr, "\n" );
1012 fprintf( stderr, "i=%d\n", i );
1013 fprintf( stderr, "currentw = \n" );
1014 for( j=0; j<lgth2; j++ )
1015 {
1016         fprintf( stderr, "%5.2f ", currentw[j] );
1017 }
1018 fprintf( stderr, "\n" );
1019 #endif
1020                 if( localhom )
1021                 {
1022 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1023 //                      imp_match_out_vead( currentw, i, lgth2 );
1024                         part_imp_match_out_vead_gapmap( currentw, gapmap1[i]+start1, lgth2, start2, gapmap2 );
1025                 }
1026 #if XXXXXXX
1027 fprintf( stderr, "\n" );
1028 fprintf( stderr, "i=%d\n", i );
1029 fprintf( stderr, "currentw = \n" );
1030 for( j=0; j<lgth2; j++ )
1031 {
1032         fprintf( stderr, "%5.2f ", currentw[j] );
1033 }
1034 fprintf( stderr, "\n" );
1035 #endif
1036                 currentw[0] = initverticalw[i];
1037
1038
1039                 mi = previousw[0] + ogcp2[1]; mpi = 0;
1040
1041                 ijppt = ijp[i] + 1;
1042                 mjpt = m + 1;
1043                 prept = previousw;
1044                 curpt = currentw + 1;
1045                 mpjpt = mp + 1;
1046                 fgcp2pt = fgcp2;
1047                 ogcp2pt = ogcp2+1;
1048                 fgcp1va = fgcp1[i-1];
1049                 ogcp1va = ogcp1[i];
1050                 for( j=1; j<lastj; j++ )
1051                 {
1052                         wm = *prept;
1053                         *ijppt = 0;
1054
1055 #if 0
1056                         fprintf( stderr, "%5.0f->", wm );
1057 #endif
1058                         g = mi + *fgcp2pt;
1059 #if 0
1060                         fprintf( stderr, "%5.0f?", g );
1061 #endif
1062                         if( g > wm )
1063                         {
1064                                 wm = g;
1065                                 *ijppt = -( j - mpi );
1066                         }
1067                         g = *prept + *ogcp2pt;
1068                         if( g >= mi )
1069                         {
1070                                 mi = g;
1071                                 mpi = j-1;
1072                         }
1073 #if USE_PENALTY_EX
1074                         mi += fpenalty_ex;
1075 #endif
1076
1077                         g = *mjpt + fgcp1va;
1078 #if 0
1079                         fprintf( stderr, "%5.0f?", g );
1080 #endif
1081                         if( g > wm )
1082                         {
1083                                 wm = g;
1084                                 *ijppt = +( i - *mpjpt );
1085                         }
1086                         g = *prept + ogcp1va;
1087                         if( g >= *mjpt )
1088                         {
1089                                 *mjpt = g;
1090                                 *mpjpt = i-1;
1091                         }
1092 #if USE_PENALTY_EX
1093                         m[j] += fpenalty_ex;
1094 #endif
1095
1096 #if 0
1097                         fprintf( stderr, "%5.0f ", wm );
1098 #endif
1099                         *curpt += wm;
1100                         ijppt++;
1101                         mjpt++;
1102                         prept++;
1103                         mpjpt++;
1104                         curpt++;
1105                         fgcp2pt++;
1106                         ogcp2pt++;
1107                 }
1108                 lastverticalw[i] = currentw[lgth2-1];
1109         }
1110
1111 #if OUTGAP0TRY
1112         if( !outgap )
1113         {
1114                 for( j=1; j<lgth2+1; j++ )
1115                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1116                 for( i=1; i<lgth1+1; i++ )
1117                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1118         }
1119 #endif
1120
1121         /*
1122         fprintf( stderr, "\n" );
1123         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1124         fprintf( stderr, "#####\n" );
1125         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1126         fprintf( stderr, "====>" );
1127         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1128         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1129         */
1130         if( localhom )
1131         {
1132                 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, start1, end1, start2, end2, gapmap1, gapmap2 );
1133         }
1134         else
1135                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1136
1137 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1138
1139         resultlen = strlen( mseq1[0] );
1140         if( alloclen < resultlen || resultlen > N )
1141         {
1142                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1143                 ErrorExit( "LENGTH OVER!\n" );
1144         }
1145
1146
1147         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1148         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1149         /*
1150         fprintf( stderr, "\n" );
1151         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1152         fprintf( stderr, "#####\n" );
1153         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1154         */
1155
1156
1157         return( wm );
1158 }
1159