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