Mercurial > repos > nick > duplex

comparison mafft/core/fftFunctions.c @ 18:e4d75f9efb90 draft

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