binaries/src/mafft/core/fftFunctions.c

   1 #include "mltaln.h"
   2
   3 #define SEGMENTSIZE 150
   4 #define TMPTMPTMP 0
   5
   6 #define DEBUG 0
   7
   8 void keika( char *str, int current, int all )
   9 {
  10         if( current == 0 )
  11                 fprintf( stderr, "%s :         ", str );
  12
  13                 fprintf( stderr, "\b\b\b\b\b\b\b\b" );
  14                 fprintf( stderr, "%3d /%3d", current+1, all+1 );
  15
  16         if( current+1 == all )
  17                 fprintf( stderr, "\b\b\b\b\b\b\b\bdone.     \n" );
  18 }
  19
  20 double maxItch( double *soukan, int size )
  21 {
  22         int i;
  23         double value = 0.0;
  24         double cand;
  25         for( i=0; i<size; i++ )
  26                 if( ( cand = soukan[i] ) > value ) value = cand;
  27         return( value );
  28 }
  29
  30 void calcNaiseki( Fukusosuu *value, Fukusosuu *x, Fukusosuu *y )
  31 {
  32         value->R =  x->R * y->R + x->I * y->I;
  33         value->I = -x->R * y->I + x->I * y->R;
  34 }
  35
  36 Fukusosuu *AllocateFukusosuuVec( int l1 )
  37 {
  38         Fukusosuu *value;
  39         value = (Fukusosuu *)calloc( l1, sizeof( Fukusosuu ) );
  40         if( !value )
  41         {
  42                 fprintf( stderr, "Cannot allocate %d FukusosuuVec\n", l1 );
  43                 return( NULL );
  44         }
  45         return( value );
  46 }
  47
  48 Fukusosuu **AllocateFukusosuuMtx( int l1, int l2 )
  49 {
  50         Fukusosuu **value;
  51         int j;
  52 //      fprintf( stderr, "allocating %d x %d FukusosuuMtx\n", l1, l2 );
  53         value = (Fukusosuu **)calloc( l1+1, sizeof( Fukusosuu * ) );
  54         if( !value )
  55         {
  56                 fprintf( stderr, "Cannot allocate %d x %d FukusosuuVecMtx\n", l1, l2 );
  57                 exit( 1 );
  58         }
  59         for( j=0; j<l1; j++ )
  60         {
  61                 value[j] = AllocateFukusosuuVec( l2 );
  62                 if( !value[j] )
  63                 {
  64                         fprintf( stderr, "Cannot allocate %d x %d FukusosuuVecMtx\n", l1, l2 );
  65                         exit( 1 );
  66                 }
  67         }
  68         value[l1] = NULL;
  69         return( value );
  70 }
  71
  72 Fukusosuu ***AllocateFukusosuuCub( int l1, int l2, int l3 )
  73 {
  74         Fukusosuu ***value;
  75         int i;
  76         value = calloc( l1+1, sizeof( Fukusosuu ** ) );
  77         if( !value ) ErrorExit( "Cannot allocate Fukusosuu" );
  78         for( i=0; i<l1; i++ ) value[i] = AllocateFukusosuuMtx( l2, l3 );
  79         value[l1] = NULL;
  80         return( value );
  81 }
  82
  83 void FreeFukusosuuVec( Fukusosuu *vec )
  84 {
  85         free( (void *)vec );
  86 }
  87
  88 void FreeFukusosuuMtx( Fukusosuu **mtx )
  89 {
  90         int i;
  91
  92         for( i=0; mtx[i]; i++ )
  93                 free( (void *)mtx[i] );
  94         free( (void *)mtx );
  95 }
  96
  97 int getKouho( int *kouho, int nkouho, double *soukan, int nlen2 )
  98 {
  99         int i, j;
 100         int nlen4 = nlen2 / 2;
 101         double max;
 102         double tmp;
 103         int ikouho = 0; // by D.Mathog, iinoka?
 104         for( j=0; j<nkouho; j++ )
 105         {
 106                 max = -9999.9;
 107                 for( i=0; i<nlen2; i++ )
 108                 {
 109                         if( ( tmp = soukan[i] ) > max )
 110                         {
 111                                 ikouho = i;
 112                                 max = tmp;
 113                         }
 114                 }
 115 #if 0
 116                 if( max < 0.15 )
 117                 {
 118                         break;
 119                 }
 120 #endif
 121 #if 0
 122                 fprintf( stderr, "Kouho No.%d, pos=%d, score=%f, lag=%d\n", j, ikouho, soukan[ikouho], ikouho-nlen4 );
 123 #endif
 124                 soukan[ikouho] = -9999.9;
 125                 kouho[j] = ( ikouho - nlen4 );
 126         }
 127         return( j );
 128 }
 129
 130 void zurasu2( int lag, int    clus1, int    clus2,
 131                        char  **seq1, char  **seq2,
 132                                            char **aseq1, char **aseq2 )
 133 {
 134         int i;
 135 #if 0
 136         fprintf( stderr, "### lag = %d\n", lag );
 137 #endif
 138         if( lag > 0 )
 139         {
 140                 for( i=0; i<clus1; i++ ) aseq1[i] = seq1[i];
 141                 for( i=0; i<clus2; i++ ) aseq2[i] = seq2[i]+lag;
 142         }
 143         else
 144         {
 145                 for( i=0; i<clus1; i++ ) aseq1[i] = seq1[i]-lag;
 146                 for( i=0; i<clus2; i++ ) aseq2[i] = seq2[i];
 147         }
 148 }
 149
 150 void zurasu( int lag, int    clus1, int    clus2,
 151                       char  **seq1, char  **seq2,
 152                                           char **aseq1, char **aseq2 )
 153 {
 154         int i;
 155 #if DEBUG
 156         fprintf( stderr, "lag = %d\n", lag );
 157 #endif
 158         if( lag > 0 )
 159         {
 160                 for( i=0; i<clus1; i++ ) strcpy( aseq1[i], seq1[i] );
 161                 for( i=0; i<clus2; i++ ) strcpy( aseq2[i], seq2[i]+lag );
 162         }
 163         else
 164         {
 165                 for( i=0; i<clus1; i++ ) strcpy( aseq1[i], seq1[i]-lag );
 166                 for( i=0; i<clus2; i++ ) strcpy( aseq2[i], seq2[i] );
 167         }
 168 }
 169
 170
 171 int alignableReagion( int    clus1, int    clus2,
 172                                            char  **seq1, char  **seq2,
 173                                            double *eff1, double *eff2,
 174                                            Segment *seg )
 175 {
 176         int i, j, k;
 177         int status, starttmp = 0; // by D.Mathog, a gess
 178         double score;
 179         int value = 0;
 180         int len, maxlen;
 181         int length = 0; // by D.Mathog, a gess
 182         static double *stra = NULL;
 183         static int alloclen = 0;
 184         double totaleff;
 185         double cumscore;
 186         static double threshold;
 187         static double prf1[26], prf2[26];
 188         static int hat1[27], hat2[27];
 189         int pre1, pre2;
 190 #if 0
 191         char **seq1pt;
 192         char **seq2pt;
 193         double *eff1pt;
 194         double *eff2pt;
 195 #endif
 196
 197 #if 0
 198         fprintf( stderr, "### In alignableRegion, clus1=%d, clus2=%d \n", clus1, clus2 );
 199         fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
 200         fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
 201         fprintf( stderr, "eff1[0] = %f\n", eff1[0] );
 202         fprintf( stderr, "eff2[0] = %f\n", eff2[0] );
 203 #endif
 204
 205         len = MIN( strlen( seq1[0] ), strlen( seq2[0] ) );
 206         maxlen = MAX( strlen( seq1[0] ), strlen( seq2[0] ) ) + fftWinSize;
 207         if( alloclen < maxlen )
 208         {
 209                 if( alloclen )
 210                 {
 211                         FreeDoubleVec( stra );
 212                 }
 213                 else
 214                 {
 215                         threshold = (int)fftThreshold / 100.0 * 600.0 * fftWinSize;
 216                 }
 217                 stra = AllocateDoubleVec( maxlen );
 218                 alloclen = maxlen;
 219         }
 220
 221
 222         totaleff = 0.0;
 223         for( i=0; i<clus1; i++ ) for( j=0; j<clus2; j++ ) totaleff += eff1[i] * eff2[j];
 224         for( i=0; i<len; i++ )
 225         {
 226                 /* make prfs */
 227                 for( j=0; j<26; j++ )
 228                 {
 229                         prf1[j] = 0.0;
 230                         prf2[j] = 0.0;
 231                 }
 232 #if 0
 233                 seq1pt = seq1;
 234                 eff1pt = eff1;
 235                 j = clus1;
 236                 while( j-- ) prf1[amino_n[(*seq1pt++)[i]]] += *eff1pt++;
 237 #else
 238                 for( j=0; j<clus1; j++ ) prf1[amino_n[(int)seq1[j][i]]] += eff1[j];
 239 #endif
 240                 for( j=0; j<clus2; j++ ) prf2[amino_n[(int)seq2[j][i]]] += eff2[j];
 241
 242                 /* make hats */
 243                 pre1 = pre2 = 26;
 244                 for( j=25; j>=0; j-- )
 245                 {
 246                         if( prf1[j] )
 247                         {
 248                                 hat1[pre1] = j;
 249                                 pre1 = j;
 250                         }
 251                         if( prf2[j] )
 252                         {
 253                                 hat2[pre2] = j;
 254                                 pre2 = j;
 255                         }
 256                 }
 257                 hat1[pre1] = -1;
 258                 hat2[pre2] = -1;
 259
 260                 /* make site score */
 261                 stra[i] = 0.0;
 262                 for( k=hat1[26]; k!=-1; k=hat1[k] )
 263                         for( j=hat2[26]; j!=-1; j=hat2[j] )
 264 //                              stra[i] += n_dis[k][j] * prf1[k] * prf2[j];
 265                                 stra[i] += n_disFFT[k][j] * prf1[k] * prf2[j];
 266                 stra[i] /= totaleff;
 267         }
 268
 269         (seg+0)->skipForeward = 0;
 270         (seg+1)->skipBackward = 0;
 271         status = 0;
 272         cumscore = 0.0;
 273         score = 0.0;
 274         for( j=0; j<fftWinSize; j++ ) score += stra[j];
 275
 276         for( i=1; i<len-fftWinSize; i++ )
 277         {
 278                 score = score - stra[i-1] + stra[i+fftWinSize-1];
 279 #if TMPTMPTMP
 280                 fprintf( stderr, "%d %10.0f   ? %10.0f\n", i, score, threshold );
 281 #endif
 282
 283                 if( score > threshold )
 284                 {
 285 #if 0
 286                         seg->start = i;
 287                         seg->end = i;
 288                         seg->center = ( seg->start + seg->end + fftWinSize ) / 2 ;
 289                         seg->score = score;
 290                         status = 0;
 291                         value++;
 292 #else
 293                         if( !status )
 294                         {
 295                                 status = 1;
 296                                 starttmp = i;
 297                                 length = 0;
 298                                 cumscore = 0.0;
 299                         }
 300                         length++;
 301                         cumscore += score;
 302 #endif
 303                 }
 304                 if( score <= threshold || length > SEGMENTSIZE )
 305                 {
 306                         if( status )
 307                         {
 308                                 if( length > fftWinSize )
 309                                 {
 310                                         seg->start = starttmp;
 311                                         seg->end = i;
 312                                         seg->center = ( seg->start + seg->end + fftWinSize ) / 2 ;
 313                                         seg->score = cumscore;
 314 #if 0
 315                                         fprintf( stderr, "%d-%d length = %d, score = %f, value = %d\n", seg->start, seg->end, length, cumscore, value );
 316 #endif
 317                                         if( length > SEGMENTSIZE )
 318                                         {
 319                                                 (seg+0)->skipForeward = 1;
 320                                                 (seg+1)->skipBackward = 1;
 321                                         }
 322                                         else
 323                                         {
 324                                                 (seg+0)->skipForeward = 0;
 325                                                 (seg+1)->skipBackward = 0;
 326                                         }
 327                                         value++;
 328                                         seg++;
 329                                 }
 330                                 length = 0;
 331                                 cumscore = 0.0;
 332                                 status = 0;
 333                                 starttmp = i;
 334                                 if( value > MAXSEG - 3 ) ErrorExit( "TOO MANY SEGMENTS!");
 335                         }
 336                 }
 337         }
 338         if( status && length > fftWinSize )
 339         {
 340                 seg->end = i;
 341                 seg->start = starttmp;
 342                 seg->center = ( starttmp + i + fftWinSize ) / 2 ;
 343                 seg->score = cumscore;
 344 #if 0
 345 fprintf( stderr, "%d-%d length = %d\n", seg->start, seg->end, length );
 346 #endif
 347                 value++;
 348         }
 349 #if TMPTMPTMP
 350         exit( 0 );
 351 #endif
 352 //      fprintf( stderr, "returning %d\n", value );
 353         return( value );
 354 }
 355
 356 void blockAlign( int *cut1, int *cut2, double **ocrossscore, int *ncut )
 357 {
 358         int i, j, shift, cur1, cur2, count;
 359         static int result1[MAXSEG], result2[MAXSEG];
 360         static int ocut1[MAXSEG], ocut2[MAXSEG];
 361         double maximum;
 362         static double max[MAXSEG], maxi;
 363         static double point[MAXSEG], pointi;
 364         static double **crossscore = NULL;
 365         static int **track = NULL;
 366
 367         if( crossscore == NULL )
 368         {
 369                 crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
 370                 track = AllocateIntMtx( MAXSEG, MAXSEG );
 371         }
 372
 373 #if 0
 374         for( i=0; i<*ncut; i++ )
 375                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 376         for( i=0; i<*ncut; i++ )
 377         {
 378                 for( j=0; j<*ncut; j++ )
 379                         fprintf( stderr, "%#4.0f ", ocrossscore[i][j]/100 );
 380                 fprintf( stderr, "\n" );
 381         }
 382 #endif
 383
 384         for( i=0; i<*ncut; i++ ) for( j=0; j<*ncut; j++ )  /* mudadanaa */
 385                 crossscore[i][j] = ocrossscore[i][j];
 386         for( i=0; i<*ncut; i++ )
 387         {
 388                 ocut1[i] = cut1[i];
 389                 ocut2[i] = cut2[i];
 390         }
 391
 392         for( j=0; j<*ncut; j++ )
 393         {
 394                 max[j] = 0.0;
 395                 point[j] = 0;
 396         }
 397
 398         for( i=1; i<*ncut; i++ )
 399         {
 400                 maxi = 0.0;
 401                 pointi = 0;
 402                 for( j=1; j<*ncut; j++ )
 403                 {
 404                         maximum = crossscore[i-1][j-1];
 405                         track[i][j] = 0;
 406
 407                         if( maximum < max[j] + penalty )
 408                         {
 409                                 maximum = max[j] + penalty;
 410                                 track[i][j] = point[j] - i;
 411                         }
 412                         if( maximum < maxi + penalty )
 413                         {
 414                                 maximum = maxi + penalty;
 415                                 track[i][j] = j - pointi;
 416                         }
 417                         crossscore[i][j] += maximum;
 418
 419                         if( maxi < crossscore[i-1][j-1] )
 420                         {
 421                                 maxi = crossscore[i-1][j-1];
 422                                 pointi = j-1;
 423                         }
 424                         if( max[j] < crossscore[i-1][j-1] )
 425                         {
 426                                 max[j] = crossscore[i-1][j-1];
 427                                 point[j] = i-1;
 428                         }
 429                 }
 430         }
 431 #if 1
 432         for( i=0; i<*ncut; i++ )
 433         {
 434                 for( j=0; j<*ncut; j++ )
 435                         fprintf( stderr, "%3d ", track[i][j] );
 436                 fprintf( stderr, "\n" );
 437         }
 438 #endif
 439
 440
 441         result1[MAXSEG-1] = *ncut-1;
 442         result2[MAXSEG-1] = *ncut-1;
 443
 444         for( i=MAXSEG-1; i>=1; i-- )
 445         {
 446                 cur1 = result1[i];
 447                 cur2 = result2[i];
 448                 if( cur1 == 0 || cur2 == 0 ) break;
 449                 shift = track[cur1][cur2];
 450                 if( shift == 0 )
 451                 {
 452                         result1[i-1] = cur1 - 1;
 453                         result2[i-1] = cur2 - 1;
 454                         continue;
 455                 }
 456                 else if( shift > 0 )
 457                 {
 458                         result1[i-1] = cur1 - 1;
 459                         result2[i-1] = cur2 - shift;
 460                 }
 461                 else if( shift < 0 )
 462                 {
 463                         result1[i-1] = cur1 + shift;
 464                         result2[i-1] = cur2 - 1;
 465                 }
 466         }
 467
 468         count = 0;
 469         for( j=i; j<MAXSEG; j++ )
 470         {
 471                 if( ocrossscore[result1[j]][result2[j]] == 0.0 ) continue;
 472
 473                 if( result1[j] == result1[j-1] || result2[j] == result2[j-1] )
 474                         if( ocrossscore[result1[j]][result2[j]] > ocrossscore[result1[j-1]][result2[j-1]] )
 475                                 count--;
 476
 477                 cut1[count] = ocut1[result1[j]];
 478                 cut2[count] = ocut2[result2[j]];
 479                 count++;
 480         }
 481
 482         *ncut = count;
 483 #if 0
 484         for( i=0; i<*ncut; i++ )
 485                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 486 #endif
 487 }
 488
 489 static int permit( Segment *seg1, Segment *seg2 )
 490 {
 491         return( 0 );
 492         if( seg1->end >= seg2->start ) return( 0 );
 493         if( seg1->pair->end >= seg2->pair->start ) return( 0 );
 494         else return( 1 );
 495 }
 496
 497 void blockAlign2( int *cut1, int *cut2, Segment **seg1, Segment **seg2, double **ocrossscore, int *ncut )
 498 {
 499         int i, j, k, shift, cur1, cur2, count, klim;
 500         static int crossscoresize = 0;
 501         static int result1[MAXSEG], result2[MAXSEG];
 502         static int ocut1[MAXSEG], ocut2[MAXSEG];
 503         double maximum;
 504         static double **crossscore = NULL;
 505         static int **track = NULL;
 506         static double maxj, maxi;
 507         static int pointj, pointi;
 508
 509
 510     if( crossscoresize < *ncut+2 )
 511     {
 512         crossscoresize = *ncut+2;
 513                 if( fftkeika ) fprintf( stderr, "allocating crossscore and track, size = %d\n", crossscoresize );
 514                 if( track ) FreeIntMtx( track );
 515         if( crossscore ) FreeDoubleMtx( crossscore );
 516                 track = AllocateIntMtx( crossscoresize, crossscoresize );
 517         crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
 518     }
 519
 520 #if 0
 521         for( i=0; i<*ncut-2; i++ )
 522                 fprintf( stderr, "%d.start = %d, score = %f\n", i, seg1[i]->start, seg1[i]->score );
 523
 524         for( i=0; i<*ncut; i++ )
 525                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 526         for( i=0; i<*ncut; i++ )
 527         {
 528                 for( j=0; j<*ncut; j++ )
 529                         fprintf( stderr, "%#4.0f ", ocrossscore[i][j] );
 530                 fprintf( stderr, "\n" );
 531         }
 532 #endif
 533
 534         for( i=0; i<*ncut; i++ ) for( j=0; j<*ncut; j++ )  /* mudadanaa */
 535                 crossscore[i][j] = ocrossscore[i][j];
 536         for( i=0; i<*ncut; i++ )
 537         {
 538                 ocut1[i] = cut1[i];
 539                 ocut2[i] = cut2[i];
 540         }
 541
 542         for( i=1; i<*ncut; i++ )
 543         {
 544 #if 0
 545                 fprintf( stderr, "### i=%d/%d\n", i,*ncut );
 546 #endif
 547                 for( j=1; j<*ncut; j++ )
 548                 {
 549                         pointi = 0; maxi = 0.0;
 550                         klim = j-2;
 551                         for( k=0; k<klim; k++ )
 552                         {
 553 /*
 554                                 fprintf( stderr, "k=%d, i=%d\n", k, i );
 555 */
 556                                 if( k && k<*ncut-1 && j<*ncut-1 && !permit( seg1[k-1], seg1[j-1] ) ) continue;
 557                                 if( crossscore[i-1][k] > maxj )
 558                                 {
 559                                         pointi = k;
 560                                         maxi = crossscore[i-1][k];
 561                                 }
 562                         }
 563
 564                         pointj = 0; maxj = 0.0;
 565                         klim = i-2;
 566                         for( k=0; k<klim; k++ )
 567                         {
 568                                 if( k && k<*ncut-1 && i<*ncut-1 && !permit( seg2[k-1], seg2[i-1] ) ) continue;
 569                                 if( crossscore[k][j-1] > maxj )
 570                                 {
 571                                         pointj = k;
 572                                         maxj = crossscore[k][j-1];
 573                                 }
 574                         }
 575
 576                         maxi += penalty;
 577                         maxj += penalty;
 578
 579                         maximum = crossscore[i-1][j-1];
 580                         track[i][j] = 0;
 581
 582                         if( maximum < maxi )
 583                         {
 584                                 maximum = maxi ;
 585                                 track[i][j] = j - pointi;
 586                         }
 587
 588                         if( maximum < maxj )
 589                         {
 590                                 maximum = maxj ;
 591                                 track[i][j] = pointj - i;
 592                         }
 593
 594                         crossscore[i][j] += maximum;
 595                 }
 596         }
 597 #if 0
 598         for( i=0; i<*ncut; i++ )
 599         {
 600                 for( j=0; j<*ncut; j++ )
 601                         fprintf( stderr, "%3d ", track[i][j] );
 602                 fprintf( stderr, "\n" );
 603         }
 604 #endif
 605
 606
 607         result1[MAXSEG-1] = *ncut-1;
 608         result2[MAXSEG-1] = *ncut-1;
 609
 610         for( i=MAXSEG-1; i>=1; i-- )
 611         {
 612                 cur1 = result1[i];
 613                 cur2 = result2[i];
 614                 if( cur1 == 0 || cur2 == 0 ) break;
 615                 shift = track[cur1][cur2];
 616                 if( shift == 0 )
 617                 {
 618                         result1[i-1] = cur1 - 1;
 619                         result2[i-1] = cur2 - 1;
 620                         continue;
 621                 }
 622                 else if( shift > 0 )
 623                 {
 624                         result1[i-1] = cur1 - 1;
 625                         result2[i-1] = cur2 - shift;
 626                 }
 627                 else if( shift < 0 )
 628                 {
 629                         result1[i-1] = cur1 + shift;
 630                         result2[i-1] = cur2 - 1;
 631                 }
 632         }
 633
 634         count = 0;
 635         for( j=i; j<MAXSEG; j++ )
 636         {
 637                 if( ocrossscore[result1[j]][result2[j]] == 0.0 ) continue;
 638
 639                 if( result1[j] == result1[j-1] || result2[j] == result2[j-1] )
 640                         if( ocrossscore[result1[j]][result2[j]] > ocrossscore[result1[j-1]][result2[j-1]] )
 641                                 count--;
 642
 643                 cut1[count] = ocut1[result1[j]];
 644                 cut2[count] = ocut2[result2[j]];
 645
 646                 count++;
 647         }
 648
 649         *ncut = count;
 650 #if 0
 651         for( i=0; i<*ncut; i++ )
 652                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 653 #endif
 654 }
 655 void blockAlign3( int *cut1, int *cut2, Segment **seg1, Segment **seg2, double **ocrossscore, int *ncut )
 656 // memory complexity = O(n^3), time complexity = O(n^2)
 657 {
 658         int i, j, shift, cur1, cur2, count;
 659         static int crossscoresize = 0;
 660         static int jumpposi, *jumppos;
 661         static double jumpscorei, *jumpscore;
 662         static int result1[MAXSEG], result2[MAXSEG];
 663         static int ocut1[MAXSEG], ocut2[MAXSEG];
 664         double maximum;
 665         static double **crossscore = NULL;
 666         static int **track = NULL;
 667
 668     if( crossscoresize < *ncut+2 )
 669     {
 670         crossscoresize = *ncut+2;
 671                 if( fftkeika ) fprintf( stderr, "allocating crossscore and track, size = %d\n", crossscoresize );
 672                 if( track ) FreeIntMtx( track );
 673         if( crossscore ) FreeDoubleMtx( crossscore );
 674         if( jumppos ) FreeIntVec( jumppos );
 675         if( jumpscore ) FreeDoubleVec( jumpscore );
 676                 track = AllocateIntMtx( crossscoresize, crossscoresize );
 677         crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
 678         jumppos = AllocateIntVec( crossscoresize );
 679         jumpscore = AllocateDoubleVec( crossscoresize );
 680     }
 681
 682 #if 0
 683         for( i=0; i<*ncut-2; i++ )
 684                 fprintf( stderr, "%d.start = %d, score = %f\n", i, seg1[i]->start, seg1[i]->score );
 685
 686         for( i=0; i<*ncut; i++ )
 687                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 688         for( i=0; i<*ncut; i++ )
 689         {
 690                 for( j=0; j<*ncut; j++ )
 691                         fprintf( stderr, "%#4.0f ", ocrossscore[i][j] );
 692                 fprintf( stderr, "\n" );
 693         }
 694 #endif
 695
 696         for( i=0; i<*ncut; i++ ) for( j=0; j<*ncut; j++ )  /* mudadanaa */
 697                 crossscore[i][j] = ocrossscore[i][j];
 698         for( i=0; i<*ncut; i++ )
 699         {
 700                 ocut1[i] = cut1[i];
 701                 ocut2[i] = cut2[i];
 702         }
 703         for( j=0; j<*ncut; j++ )
 704         {
 705                 jumpscore[j] = -999.999;
 706                 jumppos[j] = -1;
 707         }
 708
 709         for( i=1; i<*ncut; i++ )
 710         {
 711
 712                 jumpscorei = -999.999;
 713                 jumpposi = -1;
 714
 715                 for( j=1; j<*ncut; j++ )
 716                 {
 717 #if 1
 718                         fprintf( stderr, "in blockalign3, ### i=%d, j=%d\n", i, j );
 719 #endif
 720
 721
 722 #if 0
 723                         for( k=0; k<j-2; k++ )
 724                         {
 725 /*
 726                                 fprintf( stderr, "k=%d, i=%d\n", k, i );
 727 */
 728                                 if( k && k<*ncut-1 && j<*ncut-1 && !permit( seg1[k-1], seg1[j-1] ) ) continue;
 729                                 if( crossscore[i-1][k] > maxj )
 730                                 {
 731                                         pointi = k;
 732                                         maxi = crossscore[i-1][k];
 733                                 }
 734                         }
 735
 736                         pointj = 0; maxj = 0.0;
 737                         for( k=0; k<i-2; k++ )
 738                         {
 739                                 if( k && k<*ncut-1 && i<*ncut-1 && !permit( seg2[k-1], seg2[i-1] ) ) continue;
 740                                 if( crossscore[k][j-1] > maxj )
 741                                 {
 742                                         pointj = k;
 743                                         maxj = crossscore[k][j-1];
 744                                 }
 745                         }
 746
 747
 748                         maxi += penalty;
 749                         maxj += penalty;
 750 #endif
 751                         maximum = crossscore[i-1][j-1];
 752                         track[i][j] = 0;
 753
 754                         if( maximum < jumpscorei && permit( seg1[jumpposi], seg1[i] ) )
 755                         {
 756                                 maximum = jumpscorei;
 757                                 track[i][j] = j - jumpposi;
 758                         }
 759
 760                         if( maximum < jumpscore[j] && permit( seg2[jumppos[j]], seg2[j] ) )
 761                         {
 762                                 maximum = jumpscore[j];
 763                                 track[i][j] = jumpscore[j] - i;
 764                         }
 765
 766                         crossscore[i][j] += maximum;
 767
 768                         if( jumpscorei < crossscore[i-1][j] )
 769                         {
 770                                 jumpscorei = crossscore[i-1][j];
 771                                 jumpposi = j;
 772                         }
 773
 774                         if( jumpscore[j] < crossscore[i][j-1] )
 775                         {
 776                                 jumpscore[j] = crossscore[i][j-1];
 777                                 jumppos[j] = i;
 778                         }
 779                 }
 780         }
 781 #if 0
 782         for( i=0; i<*ncut; i++ )
 783         {
 784                 for( j=0; j<*ncut; j++ )
 785                         fprintf( stderr, "%3d ", track[i][j] );
 786                 fprintf( stderr, "\n" );
 787         }
 788 #endif
 789
 790
 791         result1[MAXSEG-1] = *ncut-1;
 792         result2[MAXSEG-1] = *ncut-1;
 793
 794         for( i=MAXSEG-1; i>=1; i-- )
 795         {
 796                 cur1 = result1[i];
 797                 cur2 = result2[i];
 798                 if( cur1 == 0 || cur2 == 0 ) break;
 799                 shift = track[cur1][cur2];
 800                 if( shift == 0 )
 801                 {
 802                         result1[i-1] = cur1 - 1;
 803                         result2[i-1] = cur2 - 1;
 804                         continue;
 805                 }
 806                 else if( shift > 0 )
 807                 {
 808                         result1[i-1] = cur1 - 1;
 809                         result2[i-1] = cur2 - shift;
 810                 }
 811                 else if( shift < 0 )
 812                 {
 813                         result1[i-1] = cur1 + shift;
 814                         result2[i-1] = cur2 - 1;
 815                 }
 816         }
 817
 818         count = 0;
 819         for( j=i; j<MAXSEG; j++ )
 820         {
 821                 if( ocrossscore[result1[j]][result2[j]] == 0.0 ) continue;
 822
 823                 if( result1[j] == result1[j-1] || result2[j] == result2[j-1] )
 824                         if( ocrossscore[result1[j]][result2[j]] > ocrossscore[result1[j-1]][result2[j-1]] )
 825                                 count--;
 826
 827                 cut1[count] = ocut1[result1[j]];
 828                 cut2[count] = ocut2[result2[j]];
 829
 830                 count++;
 831         }
 832
 833         *ncut = count;
 834 #if 0
 835         for( i=0; i<*ncut; i++ )
 836                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 837 #endif
 838 }
 839