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 TLS double *stra = NULL;
 183         static TLS int alloclen = 0;
 184         double totaleff;
 185         double cumscore;
 186         static TLS double threshold;
 187         static TLS double *prf1 = NULL;
 188         static TLS double *prf2 = NULL;
 189         static TLS int *hat1 = NULL;
 190         static TLS int *hat2 = NULL;
 191         int pre1, pre2;
 192 #if 0
 193         char **seq1pt;
 194         char **seq2pt;
 195         double *eff1pt;
 196         double *eff2pt;
 197 #endif
 198
 199 #if 0
 200         fprintf( stderr, "### In alignableRegion, clus1=%d, clus2=%d \n", clus1, clus2 );
 201         fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
 202         fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
 203         fprintf( stderr, "eff1[0] = %f\n", eff1[0] );
 204         fprintf( stderr, "eff2[0] = %f\n", eff2[0] );
 205 #endif
 206
 207         if( clus1 == 0 )
 208         {
 209                 FreeDoubleVec( stra ); stra = NULL;
 210                 FreeDoubleVec( prf1 ); prf1 = NULL;
 211                 FreeDoubleVec( prf2 ); prf2 = NULL;
 212                 FreeIntVec( hat1 ); hat1 = NULL;
 213                 FreeIntVec( hat2 ); hat2 = NULL;
 214                 return( 0 );
 215         }
 216
 217         if( prf1 == NULL )
 218         {
 219                 prf1 = AllocateDoubleVec( 26 );
 220                 prf2 = AllocateDoubleVec( 26 );
 221                 hat1 = AllocateIntVec( 27 );
 222                 hat2 = AllocateIntVec( 27 );
 223         }
 224
 225         len = MIN( strlen( seq1[0] ), strlen( seq2[0] ) );
 226         maxlen = MAX( strlen( seq1[0] ), strlen( seq2[0] ) ) + fftWinSize;
 227         if( alloclen < maxlen )
 228         {
 229                 if( alloclen )
 230                 {
 231                         FreeDoubleVec( stra );
 232                 }
 233                 else
 234                 {
 235                         threshold = (int)fftThreshold / 100.0 * 600.0 * fftWinSize;
 236                 }
 237                 stra = AllocateDoubleVec( maxlen );
 238                 alloclen = maxlen;
 239         }
 240
 241
 242         totaleff = 0.0;
 243         for( i=0; i<clus1; i++ ) for( j=0; j<clus2; j++ ) totaleff += eff1[i] * eff2[j];
 244         for( i=0; i<len; i++ )
 245         {
 246                 /* make prfs */
 247                 for( j=0; j<26; j++ )
 248                 {
 249                         prf1[j] = 0.0;
 250                         prf2[j] = 0.0;
 251                 }
 252 #if 0
 253                 seq1pt = seq1;
 254                 eff1pt = eff1;
 255                 j = clus1;
 256                 while( j-- ) prf1[amino_n[(*seq1pt++)[i]]] += *eff1pt++;
 257 #else
 258                 for( j=0; j<clus1; j++ ) prf1[amino_n[(int)seq1[j][i]]] += eff1[j];
 259 #endif
 260                 for( j=0; j<clus2; j++ ) prf2[amino_n[(int)seq2[j][i]]] += eff2[j];
 261
 262                 /* make hats */
 263                 pre1 = pre2 = 26;
 264                 for( j=25; j>=0; j-- )
 265                 {
 266                         if( prf1[j] )
 267                         {
 268                                 hat1[pre1] = j;
 269                                 pre1 = j;
 270                         }
 271                         if( prf2[j] )
 272                         {
 273                                 hat2[pre2] = j;
 274                                 pre2 = j;
 275                         }
 276                 }
 277                 hat1[pre1] = -1;
 278                 hat2[pre2] = -1;
 279
 280                 /* make site score */
 281                 stra[i] = 0.0;
 282                 for( k=hat1[26]; k!=-1; k=hat1[k] )
 283                         for( j=hat2[26]; j!=-1; j=hat2[j] )
 284 //                              stra[i] += n_dis[k][j] * prf1[k] * prf2[j];
 285                                 stra[i] += n_disFFT[k][j] * prf1[k] * prf2[j];
 286                 stra[i] /= totaleff;
 287         }
 288
 289         (seg+0)->skipForeward = 0;
 290         (seg+1)->skipBackward = 0;
 291         status = 0;
 292         cumscore = 0.0;
 293         score = 0.0;
 294         for( j=0; j<fftWinSize; j++ ) score += stra[j];
 295
 296         for( i=1; i<len-fftWinSize; i++ )
 297         {
 298                 score = score - stra[i-1] + stra[i+fftWinSize-1];
 299 #if TMPTMPTMP
 300                 fprintf( stderr, "%d %10.0f   ? %10.0f\n", i, score, threshold );
 301 #endif
 302
 303                 if( score > threshold )
 304                 {
 305 #if 0
 306                         seg->start = i;
 307                         seg->end = i;
 308                         seg->center = ( seg->start + seg->end + fftWinSize ) / 2 ;
 309                         seg->score = score;
 310                         status = 0;
 311                         value++;
 312 #else
 313                         if( !status )
 314                         {
 315                                 status = 1;
 316                                 starttmp = i;
 317                                 length = 0;
 318                                 cumscore = 0.0;
 319                         }
 320                         length++;
 321                         cumscore += score;
 322 #endif
 323                 }
 324                 if( score <= threshold || length > SEGMENTSIZE )
 325                 {
 326                         if( status )
 327                         {
 328                                 if( length > fftWinSize )
 329                                 {
 330                                         seg->start = starttmp;
 331                                         seg->end = i;
 332                                         seg->center = ( seg->start + seg->end + fftWinSize ) / 2 ;
 333                                         seg->score = cumscore;
 334 #if 0
 335                                         fprintf( stderr, "%d-%d length = %d, score = %f, value = %d\n", seg->start, seg->end, length, cumscore, value );
 336 #endif
 337                                         if( length > SEGMENTSIZE )
 338                                         {
 339                                                 (seg+0)->skipForeward = 1;
 340                                                 (seg+1)->skipBackward = 1;
 341                                         }
 342                                         else
 343                                         {
 344                                                 (seg+0)->skipForeward = 0;
 345                                                 (seg+1)->skipBackward = 0;
 346                                         }
 347                                         value++;
 348                                         seg++;
 349                                 }
 350                                 length = 0;
 351                                 cumscore = 0.0;
 352                                 status = 0;
 353                                 starttmp = i;
 354                                 if( value > MAXSEG - 3 ) ErrorExit( "TOO MANY SEGMENTS!");
 355                         }
 356                 }
 357         }
 358         if( status && length > fftWinSize )
 359         {
 360                 seg->end = i;
 361                 seg->start = starttmp;
 362                 seg->center = ( starttmp + i + fftWinSize ) / 2 ;
 363                 seg->score = cumscore;
 364 #if 0
 365 fprintf( stderr, "%d-%d length = %d\n", seg->start, seg->end, length );
 366 #endif
 367                 value++;
 368         }
 369 #if TMPTMPTMP
 370         exit( 0 );
 371 #endif
 372 //      fprintf( stderr, "returning %d\n", value );
 373         return( value );
 374 }
 375
 376
 377 static int permit( Segment *seg1, Segment *seg2 )
 378 {
 379         return( 0 );
 380         if( seg1->end >= seg2->start ) return( 0 );
 381         if( seg1->pair->end >= seg2->pair->start ) return( 0 );
 382         else return( 1 );
 383 }
 384
 385 void blockAlign2( int *cut1, int *cut2, Segment **seg1, Segment **seg2, double **ocrossscore, int *ncut )
 386 {
 387         int i, j, k, shift, cur1, cur2, count, klim;
 388         static TLS int crossscoresize = 0;
 389         static TLS int *result1 = NULL;
 390         static TLS int *result2 = NULL;
 391         static TLS int *ocut1 = NULL;
 392         static TLS int *ocut2 = NULL;
 393         double maximum;
 394         static TLS double **crossscore = NULL;
 395         static TLS int **track = NULL;
 396         static TLS double maxj, maxi;
 397         static TLS int pointj, pointi;
 398
 399         if( cut1 == NULL)
 400         {
 401                 if( result1 )
 402                 {
 403                         free( result1 );
 404                         free( result2 );
 405                         free( ocut1 );
 406                         free( ocut2 );
 407                         FreeIntMtx( track );
 408                 FreeDoubleMtx( crossscore );
 409                 }
 410                 return;
 411         }
 412
 413         if( result1 == NULL )
 414         {
 415                 result1 = AllocateIntVec( MAXSEG );
 416                 result2 = AllocateIntVec( MAXSEG );
 417                 ocut1 = AllocateIntVec( MAXSEG );
 418                 ocut2 = AllocateIntVec( MAXSEG );
 419         }
 420
 421     if( crossscoresize < *ncut+2 )
 422     {
 423         crossscoresize = *ncut+2;
 424                 if( fftkeika ) fprintf( stderr, "allocating crossscore and track, size = %d\n", crossscoresize );
 425                 if( track ) FreeIntMtx( track );
 426         if( crossscore ) FreeDoubleMtx( crossscore );
 427                 track = AllocateIntMtx( crossscoresize, crossscoresize );
 428         crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
 429     }
 430
 431 #if 0
 432         for( i=0; i<*ncut-2; i++ )
 433                 fprintf( stderr, "%d.start = %d, score = %f\n", i, seg1[i]->start, seg1[i]->score );
 434
 435         for( i=0; i<*ncut; i++ )
 436                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 437         for( i=0; i<*ncut; i++ )
 438         {
 439                 for( j=0; j<*ncut; j++ )
 440                         fprintf( stderr, "%#4.0f ", ocrossscore[i][j] );
 441                 fprintf( stderr, "\n" );
 442         }
 443 #endif
 444
 445         for( i=0; i<*ncut; i++ ) for( j=0; j<*ncut; j++ )  /* mudadanaa */
 446                 crossscore[i][j] = ocrossscore[i][j];
 447         for( i=0; i<*ncut; i++ )
 448         {
 449                 ocut1[i] = cut1[i];
 450                 ocut2[i] = cut2[i];
 451         }
 452
 453         for( i=1; i<*ncut; i++ )
 454         {
 455 #if 0
 456                 fprintf( stderr, "### i=%d/%d\n", i,*ncut );
 457 #endif
 458                 for( j=1; j<*ncut; j++ )
 459                 {
 460                         pointi = 0; maxi = 0.0;
 461                         klim = j-2;
 462                         for( k=0; k<klim; k++ )
 463                         {
 464 /*
 465                                 fprintf( stderr, "k=%d, i=%d\n", k, i );
 466 */
 467                                 if( k && k<*ncut-1 && j<*ncut-1 && !permit( seg1[k-1], seg1[j-1] ) ) continue;
 468                                 if( crossscore[i-1][k] > maxj )
 469                                 {
 470                                         pointi = k;
 471                                         maxi = crossscore[i-1][k];
 472                                 }
 473                         }
 474
 475                         pointj = 0; maxj = 0.0;
 476                         klim = i-2;
 477                         for( k=0; k<klim; k++ )
 478                         {
 479                                 if( k && k<*ncut-1 && i<*ncut-1 && !permit( seg2[k-1], seg2[i-1] ) ) continue;
 480                                 if( crossscore[k][j-1] > maxj )
 481                                 {
 482                                         pointj = k;
 483                                         maxj = crossscore[k][j-1];
 484                                 }
 485                         }
 486
 487                         maxi += penalty;
 488                         maxj += penalty;
 489
 490                         maximum = crossscore[i-1][j-1];
 491                         track[i][j] = 0;
 492
 493                         if( maximum < maxi )
 494                         {
 495                                 maximum = maxi ;
 496                                 track[i][j] = j - pointi;
 497                         }
 498
 499                         if( maximum < maxj )
 500                         {
 501                                 maximum = maxj ;
 502                                 track[i][j] = pointj - i;
 503                         }
 504
 505                         crossscore[i][j] += maximum;
 506                 }
 507         }
 508 #if 0
 509         for( i=0; i<*ncut; i++ )
 510         {
 511                 for( j=0; j<*ncut; j++ )
 512                         fprintf( stderr, "%3d ", track[i][j] );
 513                 fprintf( stderr, "\n" );
 514         }
 515 #endif
 516
 517
 518         result1[MAXSEG-1] = *ncut-1;
 519         result2[MAXSEG-1] = *ncut-1;
 520
 521         for( i=MAXSEG-1; i>=1; i-- )
 522         {
 523                 cur1 = result1[i];
 524                 cur2 = result2[i];
 525                 if( cur1 == 0 || cur2 == 0 ) break;
 526                 shift = track[cur1][cur2];
 527                 if( shift == 0 )
 528                 {
 529                         result1[i-1] = cur1 - 1;
 530                         result2[i-1] = cur2 - 1;
 531                         continue;
 532                 }
 533                 else if( shift > 0 )
 534                 {
 535                         result1[i-1] = cur1 - 1;
 536                         result2[i-1] = cur2 - shift;
 537                 }
 538                 else if( shift < 0 )
 539                 {
 540                         result1[i-1] = cur1 + shift;
 541                         result2[i-1] = cur2 - 1;
 542                 }
 543         }
 544
 545         count = 0;
 546         for( j=i; j<MAXSEG; j++ )
 547         {
 548                 if( ocrossscore[result1[j]][result2[j]] == 0.0 ) continue;
 549
 550                 if( result1[j] == result1[j-1] || result2[j] == result2[j-1] )
 551                         if( ocrossscore[result1[j]][result2[j]] > ocrossscore[result1[j-1]][result2[j-1]] )
 552                                 count--;
 553
 554                 cut1[count] = ocut1[result1[j]];
 555                 cut2[count] = ocut2[result2[j]];
 556
 557                 count++;
 558         }
 559
 560         *ncut = count;
 561 #if 0
 562         for( i=0; i<*ncut; i++ )
 563                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 564 #endif
 565 }
 566
 567 void blockAlign3( int *cut1, int *cut2, Segment **seg1, Segment **seg2, double **ocrossscore, int *ncut )
 568 // memory complexity = O(n^3), time complexity = O(n^2)
 569 {
 570         int i, j, shift, cur1, cur2, count;
 571         static TLS int crossscoresize = 0;
 572         static TLS int jumpposi, *jumppos;
 573         static TLS double jumpscorei, *jumpscore;
 574         static TLS int *result1 = NULL;
 575         static TLS int *result2 = NULL;
 576         static TLS int *ocut1 = NULL;
 577         static TLS int *ocut2 = NULL;
 578         double maximum;
 579         static TLS double **crossscore = NULL;
 580         static TLS int **track = NULL;
 581
 582         if( result1 == NULL )
 583         {
 584                 result1 = AllocateIntVec( MAXSEG );
 585                 result2 = AllocateIntVec( MAXSEG );
 586                 ocut1 = AllocateIntVec( MAXSEG );
 587                 ocut2 = AllocateIntVec( MAXSEG );
 588         }
 589     if( crossscoresize < *ncut+2 )
 590     {
 591         crossscoresize = *ncut+2;
 592                 if( fftkeika ) fprintf( stderr, "allocating crossscore and track, size = %d\n", crossscoresize );
 593                 if( track ) FreeIntMtx( track );
 594         if( crossscore ) FreeDoubleMtx( crossscore );
 595         if( jumppos ) FreeIntVec( jumppos );
 596         if( jumpscore ) FreeDoubleVec( jumpscore );
 597                 track = AllocateIntMtx( crossscoresize, crossscoresize );
 598         crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
 599         jumppos = AllocateIntVec( crossscoresize );
 600         jumpscore = AllocateDoubleVec( crossscoresize );
 601     }
 602
 603 #if 0
 604         for( i=0; i<*ncut-2; i++ )
 605                 fprintf( stderr, "%d.start = %d, score = %f\n", i, seg1[i]->start, seg1[i]->score );
 606
 607         for( i=0; i<*ncut; i++ )
 608                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 609         for( i=0; i<*ncut; i++ )
 610         {
 611                 for( j=0; j<*ncut; j++ )
 612                         fprintf( stderr, "%#4.0f ", ocrossscore[i][j] );
 613                 fprintf( stderr, "\n" );
 614         }
 615 #endif
 616
 617         for( i=0; i<*ncut; i++ ) for( j=0; j<*ncut; j++ )  /* mudadanaa */
 618                 crossscore[i][j] = ocrossscore[i][j];
 619         for( i=0; i<*ncut; i++ )
 620         {
 621                 ocut1[i] = cut1[i];
 622                 ocut2[i] = cut2[i];
 623         }
 624         for( j=0; j<*ncut; j++ )
 625         {
 626                 jumpscore[j] = -999.999;
 627                 jumppos[j] = -1;
 628         }
 629
 630         for( i=1; i<*ncut; i++ )
 631         {
 632
 633                 jumpscorei = -999.999;
 634                 jumpposi = -1;
 635
 636                 for( j=1; j<*ncut; j++ )
 637                 {
 638 #if 1
 639                         fprintf( stderr, "in blockalign3, ### i=%d, j=%d\n", i, j );
 640 #endif
 641
 642
 643 #if 0
 644                         for( k=0; k<j-2; k++ )
 645                         {
 646 /*
 647                                 fprintf( stderr, "k=%d, i=%d\n", k, i );
 648 */
 649                                 if( k && k<*ncut-1 && j<*ncut-1 && !permit( seg1[k-1], seg1[j-1] ) ) continue;
 650                                 if( crossscore[i-1][k] > maxj )
 651                                 {
 652                                         pointi = k;
 653                                         maxi = crossscore[i-1][k];
 654                                 }
 655                         }
 656
 657                         pointj = 0; maxj = 0.0;
 658                         for( k=0; k<i-2; k++ )
 659                         {
 660                                 if( k && k<*ncut-1 && i<*ncut-1 && !permit( seg2[k-1], seg2[i-1] ) ) continue;
 661                                 if( crossscore[k][j-1] > maxj )
 662                                 {
 663                                         pointj = k;
 664                                         maxj = crossscore[k][j-1];
 665                                 }
 666                         }
 667
 668
 669                         maxi += penalty;
 670                         maxj += penalty;
 671 #endif
 672                         maximum = crossscore[i-1][j-1];
 673                         track[i][j] = 0;
 674
 675                         if( maximum < jumpscorei && permit( seg1[jumpposi], seg1[i] ) )
 676                         {
 677                                 maximum = jumpscorei;
 678                                 track[i][j] = j - jumpposi;
 679                         }
 680
 681                         if( maximum < jumpscore[j] && permit( seg2[jumppos[j]], seg2[j] ) )
 682                         {
 683                                 maximum = jumpscore[j];
 684                                 track[i][j] = jumpscore[j] - i;
 685                         }
 686
 687                         crossscore[i][j] += maximum;
 688
 689                         if( jumpscorei < crossscore[i-1][j] )
 690                         {
 691                                 jumpscorei = crossscore[i-1][j];
 692                                 jumpposi = j;
 693                         }
 694
 695                         if( jumpscore[j] < crossscore[i][j-1] )
 696                         {
 697                                 jumpscore[j] = crossscore[i][j-1];
 698                                 jumppos[j] = i;
 699                         }
 700                 }
 701         }
 702 #if 0
 703         for( i=0; i<*ncut; i++ )
 704         {
 705                 for( j=0; j<*ncut; j++ )
 706                         fprintf( stderr, "%3d ", track[i][j] );
 707                 fprintf( stderr, "\n" );
 708         }
 709 #endif
 710
 711
 712         result1[MAXSEG-1] = *ncut-1;
 713         result2[MAXSEG-1] = *ncut-1;
 714
 715         for( i=MAXSEG-1; i>=1; i-- )
 716         {
 717                 cur1 = result1[i];
 718                 cur2 = result2[i];
 719                 if( cur1 == 0 || cur2 == 0 ) break;
 720                 shift = track[cur1][cur2];
 721                 if( shift == 0 )
 722                 {
 723                         result1[i-1] = cur1 - 1;
 724                         result2[i-1] = cur2 - 1;
 725                         continue;
 726                 }
 727                 else if( shift > 0 )
 728                 {
 729                         result1[i-1] = cur1 - 1;
 730                         result2[i-1] = cur2 - shift;
 731                 }
 732                 else if( shift < 0 )
 733                 {
 734                         result1[i-1] = cur1 + shift;
 735                         result2[i-1] = cur2 - 1;
 736                 }
 737         }
 738
 739         count = 0;
 740         for( j=i; j<MAXSEG; j++ )
 741         {
 742                 if( ocrossscore[result1[j]][result2[j]] == 0.0 ) continue;
 743
 744                 if( result1[j] == result1[j-1] || result2[j] == result2[j-1] )
 745                         if( ocrossscore[result1[j]][result2[j]] > ocrossscore[result1[j-1]][result2[j-1]] )
 746                                 count--;
 747
 748                 cut1[count] = ocut1[result1[j]];
 749                 cut2[count] = ocut2[result2[j]];
 750
 751                 count++;
 752         }
 753
 754         *ncut = count;
 755 #if 0
 756         for( i=0; i<*ncut; i++ )
 757                 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
 758 #endif
 759 }
 760