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