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

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