source: branches/lib/GDE/MAFFT/mafft-7.055-with-extensions/core/fftFunctions.c

Last change on this file was 10371, checked in by aboeckma, 12 years ago

updated mafft version. Added extensions (no svn ignore, yet)
File size: 16.3 KB
Line 
1#include "mltaln.h"
2
3#define SEGMENTSIZE 150
4#define TMPTMPTMP 0
5
6#define DEBUG 0
7
8void 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
20double 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
30void 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
36Fukusosuu *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       
48Fukusosuu **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
72Fukusosuu ***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
83void FreeFukusosuuVec( Fukusosuu *vec )
84{
85        free( (void *)vec );
86}
87
88void 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
97int 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
130void 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
150void 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
171int 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                if( stra ) FreeDoubleVec( stra ); stra = NULL;
210                if( prf1 ) FreeDoubleVec( prf1 ); prf1 = NULL;
211                if( prf2 ) FreeDoubleVec( prf2 ); prf2 = NULL;
212                if( hat1 ) FreeIntVec( hat1 ); hat1 = NULL;
213                if( hat2 ) 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
365fprintf( 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
377static 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
385void 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                        if( result1 ) free( result1 ); result1 = NULL;
404                        if( result2 ) free( result2 ); result2 = NULL;
405                        if( ocut1 ) free( ocut1 ); ocut1 = NULL;
406                        if( ocut2 ) free( ocut2 ); ocut2 = NULL;
407                        if( track ) FreeIntMtx( track ); track = NULL;
408                if( crossscore ) FreeDoubleMtx( crossscore ); crossscore = NULL;
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
567void 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
Note: See TracBrowser for help on using the repository browser.