source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/core/Falign.c

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

updated mafft version. Added extensions (no svn ignore, yet)

File size: 60.4 KB
Line 
1#include "mltaln.h"
2
3#if 0
4static FILE *fftfp;
5#endif
6static TLS int n20or4or2;
7
8#define KEIKA 0
9#define RND   0
10#define DEBUG 0
11
12#if RND // by D.Mathog
13static void generateRndSeq( char *seq, int len )
14{
15        while( len-- )
16#if 1
17                *seq++ = (int)( rnd() * n20or4or2 );
18#else
19                *seq++ = (int)1;
20#endif
21}
22#endif
23
24static void vec_init( Fukusosuu *result, int nlen )
25{
26        while( nlen-- )
27        {
28                result->R = result->I = 0.0;
29                result++;
30        }
31}
32
33#if 0 // by D.Mathog
34static void vec_init2( Fukusosuu **result, char *seq, double eff, int st, int ed )
35{
36        int i;
37        for( i=st; i<ed; i++ )
38                result[(int)*seq++][i].R += eff;
39}
40#endif
41
42static void seq_vec_2( Fukusosuu *result, double *score, double incr, char *seq )
43{
44        static TLS int n;
45        for( ; *seq; result++ )
46        {
47                n = amino_n[(int)*seq++];
48                if( n < 20 && n >= 0 ) result->R += incr * score[n];
49#if 0
50                fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n,  score[n], incr * score[n], result->R );
51#endif
52        }
53}
54
55static void seq_vec_3( Fukusosuu **result, double incr, char *seq )
56{
57        int i;
58        int n;
59        for( i=0; *seq; i++ )
60        {
61                n = amino_n[(int)*seq++];
62                if( n < n20or4or2 && n >= 0 ) result[n][i].R += incr;
63        }
64}
65
66static void seq_vec_5( Fukusosuu *result, double *score1, double *score2, double incr, char *seq )
67{
68        int n;
69        for( ; *seq; result++ )
70        {
71                n = amino_n[(int)*seq++];
72                if( n > 20 ) continue;
73                result->R += incr * score1[n];
74                result->I += incr * score2[n];
75#if 0
76                fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n,  score[n], incr * score[n], result->R );
77#endif
78        }
79}
80
81
82static void seq_vec_4( Fukusosuu *result, double incr, char *seq )
83{
84        char s;
85        for( ; *seq; result++ )
86        {
87                s = *seq++;
88                if( s == 'a' )
89                        result->R += incr;
90                else if( s == 't' )
91                        result->R -= incr;
92                else if( s == 'g' )
93                        result->I += incr;
94                else if( s == 'c' )
95                        result->I -= incr;
96        }
97}
98
99#if 0 // by D.Mathog
100static void seq_vec( Fukusosuu *result, char query, double incr, char *seq )
101{
102#if 0
103        int bk = nlen;
104#endif
105        while( *seq )
106        {
107                if( *seq++ == query ) result->R += incr;
108                result++;
109#if 0
110fprintf( stderr, "i = %d result->R = %f\n", bk-nlen, (result-1)->R );
111#endif
112        }
113}
114
115static int checkRepeat( int num, int *cutpos )
116{
117        int tmp, buf;
118
119        buf = *cutpos;
120        while( num-- )
121        {
122                if( ( tmp = *cutpos++ ) < buf ) return( 1 );
123                buf = tmp;
124        }
125        return( 0 );
126}
127
128static int segcmp( void *ptr1, void *ptr2 )
129{
130        int diff;
131        Segment **seg1 = (Segment **)ptr1;
132        Segment **seg2 = (Segment **)ptr2;
133#if 0
134        return( (*seg1)->center - (*seg2)->center );
135#else
136        diff = (*seg1)->center - (*seg2)->center;
137        if( diff ) return( diff );
138
139        diff = (*seg1)->start - (*seg2)->start;
140        if( diff ) return( diff );
141
142        diff = (*seg1)->end - (*seg2)->end;
143        if( diff ) return( diff );
144
145        fprintf( stderr, "USE STABLE SORT !!\n" );
146        exit( 1 );
147        return( 0 );
148#endif
149}
150#endif
151
152
153static void mymergesort( int first, int last, Segment **seg )
154{
155        int middle;
156        static TLS int i, j, k, p;
157        static TLS int allo = 0;
158        static TLS Segment **work = NULL;
159
160        if( seg == NULL )
161        {
162                if( work ) free( work ); 
163                work = NULL;
164                return;
165        }
166
167        if( last > allo )
168        {
169                allo = last;
170                if( work ) free( work );
171                work = (Segment **)calloc( allo / 2 + 1, sizeof( Segment *) );
172        }
173
174        if( first < last )
175        {
176                middle = ( first + last ) / 2;
177                mymergesort( first, middle, seg );
178                mymergesort( middle+1, last, seg );
179                p = 0;
180                for( i=first; i<=middle; i++ ) work[p++] = seg[i];
181                i = middle + 1; j = 0; k = first;
182                while( i <= last && j < p )
183                {
184                        if( work[j]->center <= seg[i]->center ) 
185                                seg[k++] = work[j++];
186                        else
187                                seg[k++] = seg[i++];
188                }
189                while( j < p ) seg[k++] = work[j++];
190        }
191}
192
193
194double Fgetlag( char  **seq1, char  **seq2, 
195                            double *eff1, double *eff2, 
196                            int    clus1, int    clus2,
197                            int alloclen )
198{
199        int i, j, k, l, m;
200        int nlen, nlen2, nlen4;
201        static TLS int crossscoresize = 0;
202        static TLS char **tmpseq1 = NULL;
203        static TLS char **tmpseq2 = NULL;
204        static TLS char **tmpptr1 = NULL;
205        static TLS char **tmpptr2 = NULL;
206        static TLS char **tmpres1 = NULL;
207        static TLS char **tmpres2 = NULL;
208        static TLS char **result1 = NULL;
209        static TLS char **result2 = NULL;
210#if RND
211        static TLS char **rndseq1 = NULL;
212        static TLS char **rndseq2 = NULL;
213#endif
214        static TLS Fukusosuu **seqVector1 = NULL;
215        static TLS Fukusosuu **seqVector2 = NULL;
216        static TLS Fukusosuu **naiseki = NULL;   
217        static TLS Fukusosuu *naisekiNoWa = NULL; 
218        static TLS double *soukan = NULL;
219        static TLS double **crossscore = NULL;
220        int nlentmp;
221        static TLS int *kouho = NULL;
222        static TLS Segment *segment = NULL;
223        static TLS Segment *segment1 = NULL;
224        static TLS Segment *segment2 = NULL;
225        static TLS Segment **sortedseg1 = NULL;
226        static TLS Segment **sortedseg2 = NULL;
227        static TLS int *cut1 = NULL;
228        static TLS int *cut2 = NULL;
229        static TLS int localalloclen = 0;
230        int lag;
231        int tmpint;
232        int count, count0;
233        int len1, len2;
234        int totallen;
235        float dumfl = 0.0;
236        int headgp, tailgp;
237
238        len1 = strlen( seq1[0] );
239        len2 = strlen( seq2[0] );
240        nlentmp = MAX( len1, len2 );
241
242        nlen = 1;
243        while( nlentmp >= nlen ) nlen <<= 1;
244#if 0
245        fprintf( stderr, "###   nlen    = %d\n", nlen );
246#endif
247
248        nlen2 = nlen/2; nlen4 = nlen2 / 2;
249
250#if DEBUG
251        fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
252        fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
253#endif
254
255        if( !localalloclen )
256        {
257                kouho = AllocateIntVec( NKOUHO );
258                cut1 = AllocateIntVec( MAXSEG );
259                cut2 = AllocateIntVec( MAXSEG );
260                tmpptr1 = AllocateCharMtx( njob, 0 );
261                tmpptr2 = AllocateCharMtx( njob, 0 );
262                result1 = AllocateCharMtx( njob, alloclen );
263                result2 = AllocateCharMtx( njob, alloclen );
264                tmpres1 = AllocateCharMtx( njob, alloclen );
265                tmpres2 = AllocateCharMtx( njob, alloclen );
266//              crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
267                segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
268                segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
269                segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
270                sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
271                sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
272                if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
273                        ErrorExit( "Allocation error\n" );
274
275                if     ( scoremtx == -1 ) n20or4or2 = 4;
276                else if( fftscore == 1  ) n20or4or2 = 2;
277                else                      n20or4or2 = 20;
278        }
279        if( localalloclen < nlen )
280        {
281                if( localalloclen )
282                {
283#if 1
284                        FreeFukusosuuMtx ( seqVector1 );
285                        FreeFukusosuuMtx ( seqVector2 );
286                        FreeFukusosuuVec( naisekiNoWa );
287                        FreeFukusosuuMtx( naiseki );
288                        FreeDoubleVec( soukan );
289                        FreeCharMtx( tmpseq1 );
290                        FreeCharMtx( tmpseq2 );
291#endif
292#if RND
293                        FreeCharMtx( rndseq1 );
294                        FreeCharMtx( rndseq2 );
295#endif
296                }
297
298
299                tmpseq1 = AllocateCharMtx( njob, nlen );
300                tmpseq2 = AllocateCharMtx( njob, nlen );
301                naisekiNoWa = AllocateFukusosuuVec( nlen );
302                naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
303                seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
304                seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
305                soukan = AllocateDoubleVec( nlen+1 );
306
307#if RND
308                rndseq1 = AllocateCharMtx( njob, nlen );
309                rndseq2 = AllocateCharMtx( njob, nlen );
310                for( i=0; i<njob; i++ )
311                {
312                        generateRndSeq( rndseq1[i], nlen );
313                        generateRndSeq( rndseq2[i], nlen );
314                }
315#endif
316                localalloclen = nlen;
317        }
318       
319        for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
320        for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
321
322#if 0
323fftfp = fopen( "input_of_Falign", "w" );
324fprintf( fftfp, "nlen = %d\n", nlen );
325fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
326for( i=0; i<clus1; i++ )
327        fprintf( fftfp, "%s\n", seq1[i] );
328fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
329for( i=0; i<clus2; i++ )
330        fprintf( fftfp, "%s\n", seq2[i] );
331fclose( fftfp );
332system( "less input_of_Falign < /dev/tty > /dev/tty" );
333#endif
334
335        if( fftkeika ) fprintf( stderr,  " FFT ... " );
336
337        for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
338        if( fftscore && scoremtx != -1 )
339        {
340                for( i=0; i<clus1; i++ )
341                {
342                        seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
343                        seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
344                }
345        }
346        else
347        {
348#if 0
349                for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
350                        seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
351#else
352                for( i=0; i<clus1; i++ )
353                        seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
354#endif
355        }
356#if RND
357        for( i=0; i<clus1; i++ )
358        {
359                vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
360        }
361#endif
362#if 0
363fftfp = fopen( "seqVec", "w" );
364fprintf( fftfp, "before transform\n" );
365for( k=0; k<n20or4or2; k++ )
366{
367   fprintf( fftfp, "nlen=%d\n", nlen );
368   fprintf( fftfp, "%c\n", amino[k] );
369   for( l=0; l<nlen; l++ )
370   fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
371}
372fclose( fftfp );
373system( "less seqVec < /dev/tty > /dev/tty" );
374#endif
375
376        for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
377        if( fftscore && scoremtx != -1 )
378        {
379                for( i=0; i<clus2; i++ )
380                {
381                        seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
382                        seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
383                }
384        }
385        else
386        {
387#if 0
388                for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
389                        seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
390#else
391                for( i=0; i<clus2; i++ )
392                        seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
393#endif
394        }
395#if RND
396        for( i=0; i<clus2; i++ )
397        {
398                vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
399        }
400#endif
401
402#if 0
403fftfp = fopen( "seqVec2", "w" );
404fprintf( fftfp, "before fft\n" );
405for( k=0; k<n20or4or2; k++ )
406{
407   fprintf( fftfp, "%c\n", amino[k] );
408   for( l=0; l<nlen; l++ )
409   fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
410}
411fclose( fftfp );
412system( "less seqVec2 < /dev/tty > /dev/tty" );
413#endif
414
415        for( j=0; j<n20or4or2; j++ )
416        {
417                fft( nlen, seqVector2[j], 0 );
418                fft( nlen, seqVector1[j], 0 );
419        }
420#if 0
421fftfp = fopen( "seqVec2", "w" );
422fprintf( fftfp, "#after fft\n" );
423for( k=0; k<n20or4or2; k++ )
424{
425   fprintf( fftfp, "#%c\n", amino[k] );
426   for( l=0; l<nlen; l++ )
427   fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
428}
429fclose( fftfp );
430system( "less seqVec2 < /dev/tty > /dev/tty" );
431#endif
432
433        for( k=0; k<n20or4or2; k++ ) 
434        {
435                for( l=0; l<nlen; l++ ) 
436                        calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
437        }
438        for( l=0; l<nlen; l++ ) 
439        {
440                naisekiNoWa[l].R = 0.0;
441                naisekiNoWa[l].I = 0.0;
442                for( k=0; k<n20or4or2; k++ ) 
443                {
444                        naisekiNoWa[l].R += naiseki[k][l].R;
445                        naisekiNoWa[l].I += naiseki[k][l].I;
446                }
447        }
448
449#if 0
450fftfp = fopen( "naisekiNoWa", "w" );
451fprintf( fftfp, "#Before fft\n" );
452for( l=0; l<nlen; l++ )
453        fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
454fclose( fftfp );
455system( "less naisekiNoWa < /dev/tty > /dev/tty " );
456#endif
457
458        fft( -nlen, naisekiNoWa, 0 );
459
460        for( m=0; m<=nlen2; m++ ) 
461                soukan[m] = naisekiNoWa[nlen2-m].R;
462        for( m=nlen2+1; m<nlen; m++ ) 
463                soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
464
465#if 0
466fftfp = fopen( "naisekiNoWa", "w" );
467fprintf( fftfp, "#After fft\n" );
468for( l=0; l<nlen; l++ )
469        fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R );
470fclose( fftfp );
471fftfp = fopen( "list.plot", "w"  );
472fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
473fclose( fftfp );
474system( "/usr/bin/gnuplot list.plot &" );
475#endif
476#if 0
477fprintf( stderr, "frt write start\n" );
478fftfp = fopen( "frt", "w" );
479for( l=0; l<nlen; l++ )
480        fprintf( fftfp, "%d  %f\n", l-nlen2, soukan[l] );
481fclose( fftfp );
482system( "less frt < /dev/tty > /dev/tty" );
483#if 0
484fftfp = fopen( "list.plot", "w"  );
485fprintf( fftfp, "plot 'frt'\n pause +1" );
486fclose( fftfp );
487system( "/usr/bin/gnuplot list.plot" );
488#endif
489#endif
490
491
492        getKouho( kouho, NKOUHO, soukan, nlen );
493
494#if 0
495        for( i=0; i<NKOUHO; i++ )
496        {
497                fprintf( stdout, "kouho[%d] = %d\n", i, kouho[i] );
498        }
499#endif
500
501#if KEIKA
502        fprintf( stderr, "Searching anchors ... " );
503#endif
504        count = 0;
505
506
507
508#define CAND 0
509#if CAND
510        fftfp = fopen( "cand", "w" );
511        fclose( fftfp );
512#endif
513
514        for( k=0; k<NKOUHO; k++ ) 
515        {
516
517                lag = kouho[k];
518                zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
519#if CAND
520                fftfp = fopen( "cand", "a" );
521                fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
522                fprintf( fftfp, "%s\n", tmpptr1[0] );
523                fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
524                fprintf( fftfp, "%s\n", tmpptr2[0] );
525                fprintf( fftfp, ">\n", k+1, lag );
526                fclose( fftfp );
527#endif
528                tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
529               
530                if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
531
532
533                if( tmpint == 0 ) break; // 060430 iinoka ?
534                while( tmpint-- > 0 )
535                {
536                        if( lag > 0 )
537                        {
538                                segment1[count].start  = segment[count].start ;
539                                segment1[count].end    = segment[count].end   ;
540                                segment1[count].center = segment[count].center;
541                                segment1[count].score  = segment[count].score;
542
543                                segment2[count].start  = segment[count].start  + lag;
544                                segment2[count].end    = segment[count].end    + lag;
545                                segment2[count].center = segment[count].center + lag;
546                                segment2[count].score  = segment[count].score       ;
547                        }
548                        else
549                        {
550                                segment1[count].start  = segment[count].start  - lag;
551                                segment1[count].end    = segment[count].end    - lag;
552                                segment1[count].center = segment[count].center - lag;
553                                segment1[count].score  = segment[count].score       ;
554
555                                segment2[count].start  = segment[count].start ;
556                                segment2[count].end    = segment[count].end   ;
557                                segment2[count].center = segment[count].center;
558                                segment2[count].score  = segment[count].score ;
559                        }
560#if 0
561                        fprintf( stderr, "Goukaku=%dko\n", tmpint );
562                        fprintf( stderr, "in 1 %d\n", segment1[count].center );
563                        fprintf( stderr, "in 2 %d\n", segment2[count].center );
564#endif
565                        segment1[count].pair = &segment2[count];
566                        segment2[count].pair = &segment1[count];
567                        count++;
568#if 0
569                        fprintf( stderr, "count=%d\n", count );
570#endif
571                }
572        }
573
574#if 1
575        fprintf( stderr, "done. (%d anchors)\r", count );
576#endif
577        if( !count && fftNoAnchStop )
578                ErrorExit( "Cannot detect anchor!" );
579#if 0
580        fprintf( stdout, "RESULT before sort:\n" );
581        for( l=0; l<count+1; l++ )
582        {
583                fprintf( stdout, "cut[%d]=%d, ", l, segment1[l].center );
584                fprintf( stdout, "%d score = %f\n", segment2[l].center, segment1[l].score );
585        }
586        exit( 1 );
587#endif
588
589#if KEIKA
590        fprintf( stderr, "Aligning anchors ... " );
591#endif
592        for( i=0; i<count; i++ )
593        {
594                sortedseg1[i] = &segment1[i];
595                sortedseg2[i] = &segment2[i];
596        }
597
598        {
599                mymergesort( 0, count-1, sortedseg1 ); 
600                mymergesort( 0, count-1, sortedseg2 ); 
601                for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
602                for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
603
604                if( crossscoresize < count+2 )
605                {
606                        crossscoresize = count+2;
607                        fprintf( stderr, "####################################################################################################################################allocating crossscore, size = %d\n", crossscoresize );
608                        if( crossscore ) FreeDoubleMtx( crossscore );
609                        crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
610                }
611
612                for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
613                        crossscore[i][j] = 0.0;
614                for( i=0; i<count; i++ )
615                {
616                        crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
617                        cut1[i+1] = sortedseg1[i]->center;
618                        cut2[i+1] = sortedseg2[i]->center;
619                }
620
621#if DEBUG
622                fprintf( stderr, "AFTER SORT\n" );
623                for( i=0; i<count; i++ ) fprintf( stderr, "%d, %d\n", segment1[i].start, segment2[i].start );
624#endif
625
626                crossscore[0][0] = 10000000.0;
627                cut1[0] = 0; 
628                cut2[0] = 0;
629                crossscore[count+1][count+1] = 10000000.0;
630                cut1[count+1] = len1;
631                cut2[count+1] = len2;
632                count += 2;
633                count0 = count;
634
635                blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
636        }
637        if( fftkeika )
638        {
639                if( count0 > count )
640                {
641                        fprintf( stderr, "REPEAT!? \n" ); 
642                        if( fftRepeatStop ) exit( 1 );
643                }
644#if KEIKA
645                else 
646                        fprintf( stderr, "done\n" );
647                        fprintf( stderr, "done. (%d anchors)\n", count );
648#endif
649        }
650
651#if 0
652        fftfp = fopen( "fft", "a" );
653        fprintf( fftfp, "RESULT after sort:\n" );
654        for( l=0; l<count; l++ )
655        {
656                fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
657                fprintf( fftfp, "%d\n", segment2[l].center );
658        }
659        fclose( fftfp );
660#endif
661
662#if 0
663        fftfp = fopen( "fft", "a" );
664        fprintf( fftfp, "RESULT after sort:\n" );
665        for( l=0; l<count; l++ )
666        {
667                fprintf( fftfp, "cut : %d %d\n", cut1[l], cut2[l] );
668        }
669        fclose( fftfp );
670#endif
671
672#if KEIKA
673        fprintf( trap_g, "Devided to %d segments\n", count-1 );
674        fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
675#endif
676
677        totallen = 0;
678        for( j=0; j<clus1; j++ ) result1[j][0] = 0;
679        for( j=0; j<clus2; j++ ) result2[j][0] = 0;
680        for( i=0; i<count-1; i++ )
681        {
682                if( i == 0 ) headgp = outgap; else headgp = 1;
683                if( i == count-2 ) tailgp = outgap; else tailgp = 1;
684               
685#if DEBUG
686                fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
687#else
688#if KEIKA
689                fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
690#endif
691#endif
692                for( j=0; j<clus1; j++ )
693                {
694                        strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
695                        tmpres1[j][cut1[i+1]-cut1[i]] = 0;
696                }
697                for( j=0; j<clus2; j++ )
698                {
699                        strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
700                        tmpres2[j][cut2[i+1]-cut2[i]] = 0;
701                }
702                switch( alg )
703                {
704                        case( 'a' ):
705                                Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
706                                break;
707                        case( 'M' ):
708                                        MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, headgp, tailgp );
709                                break;
710                        case( 'A' ):
711                                if( clus1 == 1 && clus2 == 1 )
712                                        G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
713                                else
714                                        A__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, headgp, tailgp );
715                                break;
716                        case( 'H' ):
717                                if( clus1 == 1 && clus2 == 1 )
718                                        G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
719                                else
720                                        H__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
721                                break;
722                        case( 'Q' ):
723                                if( clus1 == 1 && clus2 == 1 )
724                                        G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
725                                else
726                                        Q__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
727                                break;
728                        default:
729                                fprintf( stderr, "alg = %c\n", alg );
730                                ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
731                                break;
732                }
733
734                nlen = strlen( tmpres1[0] );
735                if( totallen + nlen > alloclen ) ErrorExit( "LENGTH OVER in Falign\n " );
736                for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
737                for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
738                totallen += nlen;
739#if 0
740                fprintf( stderr, "%4d\r", totallen );
741                fprintf( stderr, "\n\n" );
742                for( j=0; j<clus1; j++ )
743                {
744                        fprintf( stderr, "%s\n", tmpres1[j] );
745                }
746                fprintf( stderr, "-------\n" );
747                for( j=0; j<clus2; j++ )
748                {
749                        fprintf( stderr, "%s\n", tmpres2[j] );
750                }
751#endif
752        }
753#if KEIKA
754        fprintf( stderr, "DP ... done   \n" );
755#endif
756
757        for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
758        for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
759#if 0
760        for( j=0; j<clus1; j++ )
761        {
762                fprintf( stderr, "%s\n", result1[j] );
763        }
764        fprintf( stderr, "- - - - - - - - - - -\n" );
765        for( j=0; j<clus2; j++ )
766        {
767                fprintf( stderr, "%s\n", result2[j] );
768        }
769#endif
770        return( 0.0 );
771}
772
773
774
775float Falign( char  **seq1, char  **seq2, 
776                          double *eff1, double *eff2, 
777                          int    clus1, int    clus2,
778                          int alloclen, int *fftlog,
779                          int *chudanpt, int chudanref, int *chudanres )
780{
781        int i, j, k, l, m, maxk;
782        int nlen, nlen2, nlen4;
783        static TLS int prevalloclen = 0;
784        static TLS int crossscoresize = 0;
785        static TLS char **tmpseq1 = NULL;
786        static TLS char **tmpseq2 = NULL;
787        static TLS char **tmpptr1 = NULL;
788        static TLS char **tmpptr2 = NULL;
789        static TLS char **tmpres1 = NULL;
790        static TLS char **tmpres2 = NULL;
791        static TLS char **result1 = NULL;
792        static TLS char **result2 = NULL;
793#if RND
794        static TLS char **rndseq1 = NULL;
795        static TLS char **rndseq2 = NULL;
796#endif
797        static TLS Fukusosuu **seqVector1 = NULL;
798        static TLS Fukusosuu **seqVector2 = NULL;
799        static TLS Fukusosuu **naiseki = NULL;   
800        static TLS Fukusosuu *naisekiNoWa = NULL; 
801        static TLS double *soukan = NULL;
802        static TLS double **crossscore = NULL;
803        int nlentmp;
804        static TLS int *kouho = NULL;
805        static TLS Segment *segment = NULL;
806        static TLS Segment *segment1 = NULL;
807        static TLS Segment *segment2 = NULL;
808        static TLS Segment **sortedseg1 = NULL;
809        static TLS Segment **sortedseg2 = NULL;
810        static TLS int *cut1 = NULL;
811        static TLS int *cut2 = NULL;
812        static TLS char *sgap1, *egap1, *sgap2, *egap2;
813        static TLS int localalloclen = 0;
814        int lag;
815        int tmpint;
816        int count, count0;
817        int len1, len2;
818        int totallen;
819        float totalscore;
820        float dumfl = 0.0;
821        int headgp, tailgp;
822
823
824        if( seq1 == NULL )
825        {
826                if( result1 ) 
827                {
828//                      fprintf( stderr, "Freeing localarrays in Falign\n" );
829                        localalloclen = 0;
830                        mymergesort( 0, 0, NULL );
831                        alignableReagion( 0, 0, NULL, NULL, NULL, NULL, NULL );
832                        fft( 0, NULL, 1 );
833                        A__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
834                        G__align11( NULL, NULL, 0, 0, 0 );
835                        blockAlign2( NULL, NULL, NULL, NULL, NULL, NULL );
836                        if( crossscore ) FreeDoubleMtx( crossscore );
837                        FreeCharMtx( result1 ); result1 = NULL;
838                        FreeCharMtx( result2 );
839                        FreeCharMtx( tmpres1 );
840                        FreeCharMtx( tmpres2 );
841                        FreeCharMtx( tmpseq1 );
842                        FreeCharMtx( tmpseq2 );
843                        free( sgap1 );
844                        free( egap1 );
845                        free( sgap2 );
846                        free( egap2 );
847                        free( kouho );
848                        free( cut1 );
849                        free( cut2 );
850                        free( tmpptr1 );
851                        free( tmpptr2 );
852                        free( segment );
853                        free( segment1 );
854                        free( segment2 );
855                        free( sortedseg1 );
856                        free( sortedseg2 );
857                        if( !kobetsubunkatsu )
858                        {
859                                FreeFukusosuuMtx ( seqVector1 );
860                                FreeFukusosuuMtx ( seqVector2 );
861                                FreeFukusosuuVec( naisekiNoWa );
862                                FreeFukusosuuMtx( naiseki );
863                                FreeDoubleVec( soukan );
864                        }
865                }
866                else
867                {
868//                      fprintf( stderr, "Did not allocate localarrays in Falign\n" );
869                }
870
871                return( 0.0 );
872        }
873
874        len1 = strlen( seq1[0] );
875        len2 = strlen( seq2[0] );
876        nlentmp = MAX( len1, len2 );
877
878        nlen = 1;
879        while( nlentmp >= nlen ) nlen <<= 1;
880#if 0
881        fprintf( stderr, "###   nlen    = %d\n", nlen );
882#endif
883
884        nlen2 = nlen/2; nlen4 = nlen2 / 2;
885
886#if DEBUG
887        fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
888        fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
889#endif
890
891        if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
892        {
893                if( prevalloclen )
894                {
895                        FreeCharMtx( result1 );
896                        FreeCharMtx( result2 );
897                        FreeCharMtx( tmpres1 );
898                        FreeCharMtx( tmpres2 );
899                }
900//              fprintf( stderr, "\n\n\nreallocating ...\n" );
901                result1 = AllocateCharMtx( njob, alloclen );
902                result2 = AllocateCharMtx( njob, alloclen );
903                tmpres1 = AllocateCharMtx( njob, alloclen );
904                tmpres2 = AllocateCharMtx( njob, alloclen );
905                prevalloclen = alloclen;
906        }
907        if( !localalloclen )
908        {
909                sgap1 = AllocateCharVec( njob );
910                egap1 = AllocateCharVec( njob );
911                sgap2 = AllocateCharVec( njob );
912                egap2 = AllocateCharVec( njob );
913                kouho = AllocateIntVec( NKOUHO );
914                cut1 = AllocateIntVec( MAXSEG );
915                cut2 = AllocateIntVec( MAXSEG );
916                tmpptr1 = AllocateCharMtx( njob, 0 );
917                tmpptr2 = AllocateCharMtx( njob, 0 );
918//              crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
919                segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
920                segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
921                segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
922                sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
923                sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
924                if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
925                        ErrorExit( "Allocation error\n" );
926
927                if     ( scoremtx == -1 ) n20or4or2 = 1;
928                else if( fftscore )       n20or4or2 = 1;
929                else                      n20or4or2 = 20;
930        }
931        if( localalloclen < nlen )
932        {
933                if( localalloclen )
934                {
935#if 1
936                        if( !kobetsubunkatsu )
937                        {
938                                FreeFukusosuuMtx ( seqVector1 );
939                                FreeFukusosuuMtx ( seqVector2 );
940                                FreeFukusosuuVec( naisekiNoWa );
941                                FreeFukusosuuMtx( naiseki );
942                                FreeDoubleVec( soukan );
943                        }
944                        FreeCharMtx( tmpseq1 );
945                        FreeCharMtx( tmpseq2 );
946#endif
947#if RND
948                        FreeCharMtx( rndseq1 );
949                        FreeCharMtx( rndseq2 );
950#endif
951                }
952
953                tmpseq1 = AllocateCharMtx( njob, nlen );
954                tmpseq2 = AllocateCharMtx( njob, nlen );
955                if( !kobetsubunkatsu )
956                {
957                        naisekiNoWa = AllocateFukusosuuVec( nlen );
958                        naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
959                        seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
960                        seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
961                        soukan = AllocateDoubleVec( nlen+1 );
962                }
963#if RND
964                rndseq1 = AllocateCharMtx( njob, nlen );
965                rndseq2 = AllocateCharMtx( njob, nlen );
966                for( i=0; i<njob; i++ )
967                {
968                        generateRndSeq( rndseq1[i], nlen );
969                        generateRndSeq( rndseq2[i], nlen );
970                }
971#endif
972                localalloclen = nlen;
973        }
974       
975        for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
976        for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
977
978#if 0
979fftfp = fopen( "input_of_Falign", "w" );
980fprintf( fftfp, "nlen = %d\n", nlen );
981fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
982for( i=0; i<clus1; i++ )
983        fprintf( fftfp, "%s\n", seq1[i] );
984fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
985for( i=0; i<clus2; i++ )
986        fprintf( fftfp, "%s\n", seq2[i] );
987fclose( fftfp );
988system( "less input_of_Falign < /dev/tty > /dev/tty" );
989#endif
990        if( !kobetsubunkatsu )
991        {
992                if( fftkeika ) fprintf( stderr,  " FFT ... " );
993
994                for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
995                if( fftscore && scoremtx != -1 )
996                {
997                        for( i=0; i<clus1; i++ )
998                        {
999#if 1
1000                                seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
1001#else
1002                                seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
1003                                seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
1004#endif
1005                        }
1006                }
1007                else
1008                {
1009#if 0
1010                        for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
1011                                seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
1012#else
1013                        for( i=0; i<clus1; i++ )
1014                                seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
1015#endif
1016                }
1017#if RND
1018                for( i=0; i<clus1; i++ )
1019                {
1020                        vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
1021                }
1022#endif
1023#if 0
1024fftfp = fopen( "seqVec", "w" );
1025fprintf( fftfp, "before transform\n" );
1026for( k=0; k<n20or4or2; k++ )
1027{
1028   fprintf( fftfp, "nlen=%d\n", nlen );
1029   fprintf( fftfp, "%c\n", amino[k] );
1030   for( l=0; l<nlen; l++ )
1031   fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
1032}
1033fclose( fftfp );
1034system( "less seqVec < /dev/tty > /dev/tty" );
1035#endif
1036
1037                for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
1038                if( fftscore && scoremtx != -1 )
1039                {
1040                        for( i=0; i<clus2; i++ )
1041                        {
1042#if 1
1043                                seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
1044#else
1045                                seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
1046                                seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
1047#endif
1048                        }
1049                }
1050                else
1051                {
1052#if 0
1053                        for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
1054                                seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
1055#else
1056                        for( i=0; i<clus2; i++ )
1057                                seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
1058#endif
1059                }
1060#if RND
1061                for( i=0; i<clus2; i++ )
1062                {
1063                        vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
1064                }
1065#endif
1066
1067#if 0
1068fftfp = fopen( "seqVec2", "w" );
1069fprintf( fftfp, "before fft\n" );
1070for( k=0; k<n20or4or2; k++ )
1071{
1072   fprintf( fftfp, "%c\n", amino[k] );
1073   for( l=0; l<nlen; l++ )
1074   fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1075}
1076fclose( fftfp );
1077system( "less seqVec2 < /dev/tty > /dev/tty" );
1078#endif
1079
1080                for( j=0; j<n20or4or2; j++ )
1081                {
1082                        fft( nlen, seqVector2[j], 0 );
1083                        fft( nlen, seqVector1[j], 0 );
1084                }
1085#if 0
1086fftfp = fopen( "seqVec2", "w" );
1087fprintf( fftfp, "#after fft\n" );
1088for( k=0; k<n20or4or2; k++ )
1089{
1090   fprintf( fftfp, "#%c\n", amino[k] );
1091   for( l=0; l<nlen; l++ )
1092           fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1093}
1094fclose( fftfp );
1095system( "less seqVec2 < /dev/tty > /dev/tty" );
1096#endif
1097
1098                for( k=0; k<n20or4or2; k++ ) 
1099                {
1100                        for( l=0; l<nlen; l++ ) 
1101                                calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
1102                }
1103                for( l=0; l<nlen; l++ ) 
1104                {
1105                        naisekiNoWa[l].R = 0.0;
1106                        naisekiNoWa[l].I = 0.0;
1107                        for( k=0; k<n20or4or2; k++ ) 
1108                        {
1109                                naisekiNoWa[l].R += naiseki[k][l].R;
1110                                naisekiNoWa[l].I += naiseki[k][l].I;
1111                        }
1112                }
1113       
1114#if 0
1115        fftfp = fopen( "naisekiNoWa", "w" );
1116        fprintf( fftfp, "#Before fft\n" );
1117        for( l=0; l<nlen; l++ )
1118                fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
1119        fclose( fftfp );
1120        system( "less naisekiNoWa < /dev/tty > /dev/tty " );
1121#endif
1122
1123                fft( -nlen, naisekiNoWa, 0 );
1124       
1125                for( m=0; m<=nlen2; m++ ) 
1126                        soukan[m] = naisekiNoWa[nlen2-m].R;
1127                for( m=nlen2+1; m<nlen; m++ ) 
1128                        soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
1129
1130#if 0
1131        fftfp = fopen( "naisekiNoWa", "w" );
1132        fprintf( fftfp, "#After fft\n" );
1133        for( l=0; l<nlen; l++ )
1134                fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R );
1135        fclose( fftfp );
1136        fftfp = fopen( "list.plot", "w"  );
1137        fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
1138        fclose( fftfp );
1139        system( "/usr/bin/gnuplot list.plot &" );
1140#endif
1141#if 0
1142        fprintf( stderr, "soukan\n" );
1143        for( l=0; l<nlen; l++ )
1144                fprintf( stderr, "%d  %f\n", l-nlen2, soukan[l] );
1145#if 0
1146        fftfp = fopen( "list.plot", "w"  );
1147        fprintf( fftfp, "plot 'frt'\n pause +1" );
1148        fclose( fftfp );
1149        system( "/usr/bin/gnuplot list.plot" );
1150#endif
1151#endif
1152
1153
1154                getKouho( kouho, NKOUHO, soukan, nlen );
1155
1156#if 0
1157                for( i=0; i<NKOUHO; i++ )
1158                {
1159                        fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
1160                }
1161#endif
1162        }
1163
1164#if KEIKA
1165        fprintf( stderr, "Searching anchors ... " );
1166#endif
1167        count = 0;
1168
1169
1170
1171#define CAND 0
1172#if CAND
1173        fftfp = fopen( "cand", "w" );
1174        fclose( fftfp );
1175#endif
1176        if( kobetsubunkatsu )
1177        {
1178                maxk = 1;
1179                kouho[0] = 0;
1180        }
1181        else
1182        {
1183                maxk = NKOUHO;
1184        }
1185
1186        for( k=0; k<maxk; k++ ) 
1187        {
1188                lag = kouho[k];
1189                if( lag <= -len1 || len2 <= lag ) continue;
1190                zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
1191#if CAND
1192                fftfp = fopen( "cand", "a" );
1193                fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1194                fprintf( fftfp, "%s\n", tmpptr1[0] );
1195                fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1196                fprintf( fftfp, "%s\n", tmpptr2[0] );
1197                fprintf( fftfp, ">\n", k+1, lag );
1198                fclose( fftfp );
1199#endif
1200
1201//              fprintf( stderr, "lag = %d\n", lag );
1202                tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
1203
1204//              if( lag == -50 ) exit( 1 );
1205               
1206                if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
1207
1208
1209                if( tmpint == 0 ) break; // 060430 iinoka ?
1210                while( tmpint-- > 0 )
1211                {
1212#if 0
1213                        if( segment[count].end - segment[count].start < fftWinSize )
1214                        {
1215                                count++;
1216                                continue;
1217                        }
1218#endif
1219                        if( lag > 0 )
1220                        {
1221                                segment1[count].start  = segment[count].start ;
1222                                segment1[count].end    = segment[count].end   ;
1223                                segment1[count].center = segment[count].center;
1224                                segment1[count].score  = segment[count].score;
1225
1226                                segment2[count].start  = segment[count].start  + lag;
1227                                segment2[count].end    = segment[count].end    + lag;
1228                                segment2[count].center = segment[count].center + lag;
1229                                segment2[count].score  = segment[count].score       ;
1230                        }
1231                        else
1232                        {
1233                                segment1[count].start  = segment[count].start  - lag;
1234                                segment1[count].end    = segment[count].end    - lag;
1235                                segment1[count].center = segment[count].center - lag;
1236                                segment1[count].score  = segment[count].score       ;
1237
1238                                segment2[count].start  = segment[count].start ;
1239                                segment2[count].end    = segment[count].end   ;
1240                                segment2[count].center = segment[count].center;
1241                                segment2[count].score  = segment[count].score ;
1242                        }
1243#if 0
1244                        fprintf( stderr, "in 1 %d\n", segment1[count].center );
1245                        fprintf( stderr, "in 2 %d\n", segment2[count].center );
1246#endif
1247                        segment1[count].pair = &segment2[count];
1248                        segment2[count].pair = &segment1[count];
1249                        count++;
1250                }
1251        }
1252#if 0
1253        if( !kobetsubunkatsu && fftkeika )
1254                fprintf( stderr, "%d anchors found\r", count );
1255#endif
1256        if( !count && fftNoAnchStop )
1257                ErrorExit( "Cannot detect anchor!" );
1258#if 0
1259        fprintf( stderr, "RESULT before sort:\n" );
1260        for( l=0; l<count+1; l++ )
1261        {
1262                fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
1263                fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
1264        }
1265#endif
1266
1267#if KEIKA
1268        fprintf( stderr, "done. (%d anchors)\n", count );
1269        fprintf( stderr, "Aligning anchors ... " );
1270#endif
1271        for( i=0; i<count; i++ )
1272        {
1273                sortedseg1[i] = &segment1[i];
1274                sortedseg2[i] = &segment2[i];
1275        }
1276#if 0
1277        tmpsort( count, sortedseg1 );
1278        tmpsort( count, sortedseg2 );
1279        qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
1280        qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
1281#else
1282        mymergesort( 0, count-1, sortedseg1 ); 
1283        mymergesort( 0, count-1, sortedseg2 ); 
1284#endif
1285        for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
1286        for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
1287
1288
1289        if( kobetsubunkatsu )
1290        {
1291                for( i=0; i<count; i++ )
1292            {
1293                        cut1[i+1] = sortedseg1[i]->center;
1294                        cut2[i+1] = sortedseg2[i]->center;
1295                }
1296                cut1[0] = 0;
1297                cut2[0] = 0;
1298                cut1[count+1] = len1;
1299                cut2[count+1] = len2;
1300                count += 2;
1301        }
1302        else
1303        {
1304                if( crossscoresize < count+2 )
1305                {
1306                        crossscoresize = count+2;
1307#if 1
1308                        if( fftkeika ) fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
1309#endif
1310                        if( crossscore ) FreeDoubleMtx( crossscore );
1311                        crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
1312                }
1313                for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
1314                        crossscore[i][j] = 0.0;
1315                for( i=0; i<count; i++ )
1316                {
1317                        crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
1318                        cut1[i+1] = sortedseg1[i]->center;
1319                        cut2[i+1] = sortedseg2[i]->center;
1320                }
1321
1322#if 0
1323                fprintf( stderr, "AFTER SORT\n" );
1324                for( i=0; i<count+1; i++ ) fprintf( stderr, "%d, %d\n", cut1[i], cut2[i] );
1325                fprintf( stderr, "crossscore = \n" );
1326                for( i=0; i<count+1; i++ )
1327                {
1328                        for( j=0; j<count+1; j++ )
1329                                fprintf( stderr, "%.0f ", crossscore[i][j] );
1330                        fprintf( stderr, "\n" );
1331                }
1332#endif
1333
1334                crossscore[0][0] = 10000000.0;
1335                cut1[0] = 0; 
1336                cut2[0] = 0;
1337                crossscore[count+1][count+1] = 10000000.0;
1338                cut1[count+1] = len1;
1339                cut2[count+1] = len2;
1340                count += 2;
1341                count0 = count;
1342       
1343                blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
1344
1345//              if( count-count0 )
1346//                      fprintf( stderr, "%d unused anchors\n", count0-count );
1347
1348                if( !kobetsubunkatsu && fftkeika )
1349                        fprintf( stderr, "%d anchors found\n", count );
1350                if( fftkeika )
1351                {
1352                        if( count0 > count )
1353                        {
1354#if 0
1355                                fprintf( stderr, "\7 REPEAT!? \n" );
1356#else
1357                                fprintf( stderr, "REPEAT!? \n" ); 
1358#endif
1359                                if( fftRepeatStop ) exit( 1 );
1360                        }
1361#if KEIKA
1362                        else fprintf( stderr, "done\n" );
1363#endif
1364                }
1365        }
1366
1367#if 0
1368        fftfp = fopen( "fft", "a" );
1369        fprintf( fftfp, "RESULT after sort:\n" );
1370        for( l=0; l<count; l++ )
1371        {
1372                fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
1373                fprintf( fftfp, "%d\n", segment2[l].center );
1374        }
1375        fclose( fftfp );
1376#endif
1377
1378#if 0
1379        fprintf( stderr, "RESULT after blckalign:\n" );
1380        for( l=0; l<count+1; l++ )
1381        {
1382                fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
1383        }
1384#endif
1385
1386#if 0
1387        fprintf( trap_g, "Devided to %d segments\n", count-1 );
1388        fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
1389#endif
1390
1391        totallen = 0;
1392        for( j=0; j<clus1; j++ ) result1[j][0] = 0;
1393        for( j=0; j<clus2; j++ ) result2[j][0] = 0;
1394        totalscore = 0.0;
1395        *fftlog = -1;
1396        for( i=0; i<count-1; i++ )
1397        {
1398                *fftlog += 1;
1399                if( i == 0 ) headgp = outgap; else headgp = 1;
1400                if( i == count-2 ) tailgp = outgap; else tailgp = 1;
1401
1402
1403                if( cut1[i] ) // chuui
1404                {
1405//                      getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
1406//                      getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
1407                        getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
1408                        getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
1409                }
1410                else
1411                {
1412                        for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
1413                        for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
1414                }
1415                if( cut1[i+1] != len1 ) // chuui
1416                {       
1417                        getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
1418                        getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
1419                }       
1420                else   
1421                {       
1422                        for( j=0; j<clus1; j++ ) egap1[j] = 'o';
1423                        for( j=0; j<clus2; j++ ) egap2[j] = 'o';
1424                }
1425#if 0
1426                {
1427                        fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
1428                        for( j=0; j<clus1; j++ )
1429                                fprintf( stderr, "%c", sgap1[j] );
1430                        fprintf( stderr, "=kyokkaigap1-start\n" );
1431                }
1432                {
1433                        fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
1434                        for( j=0; j<clus2; j++ )
1435                                fprintf( stderr, "%c", sgap2[j] );
1436                        fprintf( stderr, "=kyokkaigap2-start\n" );
1437                }
1438                {
1439                        fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
1440                        for( j=0; j<clus1; j++ )
1441                                fprintf( stderr, "%c", egap1[j] );
1442                        fprintf( stderr, "=kyokkaigap1-end\n" );
1443                }
1444                {
1445                        fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
1446                        for( j=0; j<clus2; j++ )
1447                                fprintf( stderr, "%c", egap2[j] );
1448                        fprintf( stderr, "=kyokkaigap2-end\n" );
1449                }
1450#endif
1451
1452#if DEBUG
1453                fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
1454#else
1455#if KEIKA
1456                fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
1457#endif
1458#endif
1459                for( j=0; j<clus1; j++ )
1460                {
1461                        strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
1462                        tmpres1[j][cut1[i+1]-cut1[i]] = 0;
1463                }
1464                if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr $B$K8F$P$l$?$H$-(B fftkeika=1
1465//              if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
1466                for( j=0; j<clus2; j++ )
1467                {
1468                        strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
1469                        tmpres2[j][cut2[i+1]-cut2[i]] = 0;
1470                }
1471                if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr $B$K8F$P$l$?$H$-(B fftkeika=1
1472//              if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );
1473
1474                if( constraint )
1475                {
1476                        fprintf( stderr, "Not supported\n" );
1477                        exit( 1 );
1478                }
1479#if 0
1480                fprintf( stderr, "i=%d, before alignment", i );
1481                fprintf( stderr, "%4d\n", totallen );
1482                fprintf( stderr, "\n\n" );
1483                for( j=0; j<clus1; j++ )
1484                {
1485                        fprintf( stderr, "%s\n", tmpres1[j] );
1486                }
1487                fprintf( stderr, "-------\n" );
1488                for( j=0; j<clus2; j++ )
1489                {
1490                        fprintf( stderr, "%s\n", tmpres2[j] );
1491                }
1492#endif
1493
1494#if 0
1495                fprintf( stdout, "writing input\n" );
1496                for( j=0; j<clus1; j++ )
1497                {
1498                        fprintf( stdout, ">%d of GROUP1\n", j );
1499                        fprintf( stdout, "%s\n", tmpres1[j] );
1500                }
1501                for( j=0; j<clus2; j++ )
1502                {
1503                        fprintf( stdout, ">%d of GROUP2\n", j );
1504                        fprintf( stdout, "%s\n", tmpres2[j] );
1505                }
1506                fflush( stdout );
1507#endif
1508                switch( alg )
1509                {
1510                        case( 'a' ):
1511                                totalscore += Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
1512                                break;
1513                        case( 'M' ):
1514                                        totalscore += MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1515                                break;
1516                        case( 'A' ):
1517                                if( clus1 == 1 && clus2 == 1 )
1518                                {
1519                                        totalscore += G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
1520                                }
1521                                else
1522                                        totalscore += A__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1523                                break;
1524                        case( 'H' ):
1525                                if( clus1 == 1 && clus2 == 1 )
1526                                {
1527                                        totalscore += G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
1528                                }
1529                                else
1530                                        totalscore += H__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2 );
1531                                break;
1532                        case( 'Q' ):
1533                                if( clus1 == 1 && clus2 == 1 )
1534                                {
1535                                        totalscore += G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
1536                                }
1537                                else
1538                                        totalscore += Q__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2 );
1539                                break;
1540                        default:
1541                                fprintf( stderr, "alg = %c\n", alg );
1542                                ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
1543                                break;
1544                }
1545
1546#ifdef enablemultithread
1547                if( chudanres && *chudanres )
1548                {
1549//                      fprintf( stderr, "\n\n## CHUUDAN!!! at Falign_localhom\n" );
1550                        return( -1.0 );
1551                }
1552#endif
1553
1554                nlen = strlen( tmpres1[0] );
1555                if( totallen + nlen > alloclen )
1556                {
1557                        fprintf( stderr, "totallen=%d +  nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
1558                        ErrorExit( "LENGTH OVER in Falign\n " );
1559                }
1560                for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
1561                for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
1562                totallen += nlen;
1563#if 0
1564                fprintf( stderr, "$#####$$$$ i=%d", i );
1565                fprintf( stderr, "%4d\n", totallen );
1566                fprintf( stderr, "\n\n" );
1567                for( j=0; j<clus1; j++ )
1568                {
1569                        fprintf( stderr, "%s\n", tmpres1[j] );
1570                }
1571                fprintf( stderr, "-------\n" );
1572                for( j=0; j<clus2; j++ )
1573                {
1574                        fprintf( stderr, "%s\n", tmpres2[j] );
1575                }
1576#endif
1577        }
1578#if KEIKA
1579        fprintf( stderr, "DP ... done   \n" );
1580#endif
1581
1582        for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
1583        for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
1584#if 0
1585        for( j=0; j<clus1; j++ )
1586        {
1587                fprintf( stderr, "in Falign, %s\n", result1[j] );
1588        }
1589        fprintf( stderr, "- - - - - - - - - - -\n" );
1590        for( j=0; j<clus2; j++ )
1591        {
1592                fprintf( stderr, "in Falign, %s\n", result2[j] );
1593        }
1594#endif
1595        return( totalscore );
1596}
1597
1598
1599
1600
1601
1602
1603
1604
1605/*
1606sakujo wo kentou (2010/10/05)
1607*/
1608float Falign_udpari_long( char  **seq1, char  **seq2, 
1609                          double *eff1, double *eff2, 
1610                          int    clus1, int    clus2,
1611                          int alloclen, int *fftlog )
1612{
1613        int i, j, k, l, m, maxk;
1614        int nlen, nlen2, nlen4;
1615        static TLS int prevalloclen = 0;
1616        static TLS int crossscoresize = 0;
1617        static TLS char **tmpseq1 = NULL;
1618        static TLS char **tmpseq2 = NULL;
1619        static TLS char **tmpptr1 = NULL;
1620        static TLS char **tmpptr2 = NULL;
1621        static TLS char **tmpres1 = NULL;
1622        static TLS char **tmpres2 = NULL;
1623        static TLS char **result1 = NULL;
1624        static TLS char **result2 = NULL;
1625#if RND
1626        static TLS char **rndseq1 = NULL;
1627        static TLS char **rndseq2 = NULL;
1628#endif
1629        static TLS Fukusosuu **seqVector1 = NULL;
1630        static TLS Fukusosuu **seqVector2 = NULL;
1631        static TLS Fukusosuu **naiseki = NULL;   
1632        static TLS Fukusosuu *naisekiNoWa = NULL; 
1633        static TLS double *soukan = NULL;
1634        static TLS double **crossscore = NULL;
1635        int nlentmp;
1636        static TLS int *kouho = NULL;
1637        static TLS Segment *segment = NULL;
1638        static TLS Segment *segment1 = NULL;
1639        static TLS Segment *segment2 = NULL;
1640        static TLS Segment **sortedseg1 = NULL;
1641        static TLS Segment **sortedseg2 = NULL;
1642        static TLS int *cut1 = NULL;
1643        static TLS int *cut2 = NULL;
1644        static TLS char *sgap1, *egap1, *sgap2, *egap2;
1645        static TLS int localalloclen = 0;
1646        int lag;
1647        int tmpint;
1648        int count, count0;
1649        int len1, len2;
1650        int totallen;
1651        float totalscore;
1652        int nkouho = 0;
1653        int headgp, tailgp;
1654//      float dumfl = 0.0;
1655
1656        if( seq1 == NULL )
1657        {
1658                if( result1 ) 
1659                {
1660//                      fprintf( stderr, "### Freeing localarrays in Falign\n" );
1661                        localalloclen = 0;
1662                        mymergesort( 0, 0, NULL );
1663                        alignableReagion( 0, 0, NULL, NULL, NULL, NULL, NULL );
1664                        fft( 0, NULL, 1 );
1665                        A__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
1666                        G__align11( NULL, NULL, 0, 0, 0 );
1667                        blockAlign2( NULL, NULL, NULL, NULL, NULL, NULL );
1668                        if( crossscore ) FreeDoubleMtx( crossscore );
1669                        FreeCharMtx( result1 ); result1 = NULL;
1670                        FreeCharMtx( result2 );
1671                        FreeCharMtx( tmpres1 );
1672                        FreeCharMtx( tmpres2 );
1673                        FreeCharMtx( tmpseq1 );
1674                        FreeCharMtx( tmpseq2 );
1675                        free( sgap1 );
1676                        free( egap1 );
1677                        free( sgap2 );
1678                        free( egap2 );
1679                        free( kouho );
1680                        free( cut1 );
1681                        free( cut2 );
1682                        free( tmpptr1 );
1683                        free( tmpptr2 );
1684                        free( segment );
1685                        free( segment1 );
1686                        free( segment2 );
1687                        free( sortedseg1 );
1688                        free( sortedseg2 );
1689                        if( !kobetsubunkatsu )
1690                        {
1691                                FreeFukusosuuMtx ( seqVector1 );
1692                                FreeFukusosuuMtx ( seqVector2 );
1693                                FreeFukusosuuVec( naisekiNoWa );
1694                                FreeFukusosuuMtx( naiseki );
1695                                FreeDoubleVec( soukan );
1696                        }
1697                }
1698                else
1699                {
1700//                      fprintf( stderr, "Did not allocate localarrays in Falign\n" );
1701                }
1702
1703                return( 0.0 );
1704        }
1705
1706
1707
1708        len1 = strlen( seq1[0] );
1709        len2 = strlen( seq2[0] );
1710        nlentmp = MAX( len1, len2 );
1711
1712        nlen = 1;
1713        while( nlentmp >= nlen ) nlen <<= 1;
1714#if 0
1715        fprintf( stderr, "###   nlen    = %d\n", nlen );
1716#endif
1717
1718        nlen2 = nlen/2; nlen4 = nlen2 / 2;
1719
1720#if 0
1721        fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
1722        fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
1723#endif
1724
1725        if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
1726        {
1727                if( prevalloclen )
1728                {
1729                        FreeCharMtx( result1 );
1730                        FreeCharMtx( result2 );
1731                        FreeCharMtx( tmpres1 );
1732                        FreeCharMtx( tmpres2 );
1733                }
1734//              fprintf( stderr, "\n\n\nreallocating ...\n" );
1735                result1 = AllocateCharMtx( njob, alloclen );
1736                result2 = AllocateCharMtx( njob, alloclen );
1737                tmpres1 = AllocateCharMtx( njob, alloclen );
1738                tmpres2 = AllocateCharMtx( njob, alloclen );
1739                prevalloclen = alloclen;
1740        }
1741
1742        if( !localalloclen )
1743        {
1744                sgap1 = AllocateCharVec( njob );
1745                egap1 = AllocateCharVec( njob );
1746                sgap2 = AllocateCharVec( njob );
1747                egap2 = AllocateCharVec( njob );
1748                kouho = AllocateIntVec( NKOUHO_LONG );
1749                cut1 = AllocateIntVec( MAXSEG );
1750                cut2 = AllocateIntVec( MAXSEG );
1751                tmpptr1 = AllocateCharMtx( njob, 0 );
1752                tmpptr2 = AllocateCharMtx( njob, 0 );
1753                segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1754                segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1755                segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1756                sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
1757                sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
1758                if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
1759                        ErrorExit( "Allocation error\n" );
1760
1761                if     ( scoremtx == -1 ) n20or4or2 = 1;
1762                else if( fftscore )       n20or4or2 = 1;
1763                else                      n20or4or2 = 20;
1764        }
1765        if( localalloclen < nlen )
1766        {
1767                if( localalloclen )
1768                {
1769#if 1
1770                        if( !kobetsubunkatsu )
1771                        {
1772                                FreeFukusosuuMtx ( seqVector1 );
1773                                FreeFukusosuuMtx ( seqVector2 );
1774                                FreeFukusosuuVec( naisekiNoWa );
1775                                FreeFukusosuuMtx( naiseki );
1776                                FreeDoubleVec( soukan );
1777                        }
1778                        FreeCharMtx( tmpseq1 );
1779                        FreeCharMtx( tmpseq2 );
1780#endif
1781#if RND
1782                        FreeCharMtx( rndseq1 );
1783                        FreeCharMtx( rndseq2 );
1784#endif
1785                }
1786
1787
1788                tmpseq1 = AllocateCharMtx( njob, nlen );
1789                tmpseq2 = AllocateCharMtx( njob, nlen );
1790                if( !kobetsubunkatsu )
1791                {
1792                        naisekiNoWa = AllocateFukusosuuVec( nlen );
1793                        naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
1794                        seqVector1 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
1795                        seqVector2 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
1796                        soukan = AllocateDoubleVec( nlen+1 );
1797                }
1798#if RND
1799                rndseq1 = AllocateCharMtx( njob, nlen );
1800                rndseq2 = AllocateCharMtx( njob, nlen );
1801                for( i=0; i<njob; i++ )
1802                {
1803                        generateRndSeq( rndseq1[i], nlen );
1804                        generateRndSeq( rndseq2[i], nlen );
1805                }
1806#endif
1807                localalloclen = nlen;
1808        }
1809       
1810        for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
1811        for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
1812
1813#if 0
1814fftfp = fopen( "input_of_Falign", "w" );
1815fprintf( fftfp, "nlen = %d\n", nlen );
1816fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
1817for( i=0; i<clus1; i++ )
1818        fprintf( fftfp, "%s\n", seq1[i] );
1819fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
1820for( i=0; i<clus2; i++ )
1821        fprintf( fftfp, "%s\n", seq2[i] );
1822fclose( fftfp );
1823system( "less input_of_Falign < /dev/tty > /dev/tty" );
1824#endif
1825        if( !kobetsubunkatsu )
1826        {
1827                fprintf( stderr,  " FFT ... " );
1828
1829                for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
1830                if( scoremtx == -1 )
1831                {
1832                        for( i=0; i<clus1; i++ )
1833                                seq_vec_4( seqVector1[0], eff1[i], tmpseq1[i] );
1834                }
1835                else if( fftscore )
1836                {
1837                        for( i=0; i<clus1; i++ )
1838                        {
1839#if 0
1840                                seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
1841                                seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
1842#else
1843                                seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
1844#endif
1845                        }
1846                }
1847                else
1848                {
1849                        for( i=0; i<clus1; i++ )
1850                                seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
1851                }
1852#if RND
1853                for( i=0; i<clus1; i++ )
1854                {
1855                        vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
1856                }
1857#endif
1858#if 0
1859fftfp = fopen( "seqVec", "w" );
1860fprintf( fftfp, "before transform\n" );
1861for( k=0; k<n20or4or2; k++ )
1862{
1863   fprintf( fftfp, "nlen=%d\n", nlen );
1864   fprintf( fftfp, "%c\n", amino[k] );
1865   for( l=0; l<nlen; l++ )
1866   fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
1867}
1868fclose( fftfp );
1869system( "less seqVec < /dev/tty > /dev/tty" );
1870#endif
1871
1872                for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
1873                if( scoremtx == -1 )
1874                {
1875                        for( i=0; i<clus2; i++ )
1876                                seq_vec_4( seqVector2[0], eff2[i], tmpseq2[i] );
1877                }
1878                else if( fftscore )
1879                {
1880                        for( i=0; i<clus2; i++ )
1881                        {
1882#if 0
1883                                seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
1884                                seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
1885#else
1886                                seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
1887#endif
1888                        }
1889                }
1890                else
1891                {
1892                        for( i=0; i<clus2; i++ )
1893                                seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
1894                }
1895#if RND
1896                for( i=0; i<clus2; i++ )
1897                {
1898                        vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
1899                }
1900#endif
1901
1902#if 0
1903fftfp = fopen( "seqVec2", "w" );
1904fprintf( fftfp, "before fft\n" );
1905for( k=0; k<n20or4or2; k++ )
1906{
1907   fprintf( fftfp, "%c\n", amino[k] );
1908   for( l=0; l<nlen; l++ )
1909   fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1910}
1911fclose( fftfp );
1912system( "less seqVec2 < /dev/tty > /dev/tty" );
1913#endif
1914
1915                for( j=0; j<n20or4or2; j++ )
1916                {
1917                        fft( nlen, seqVector2[j], 0 );
1918                        fft( nlen, seqVector1[j], 0 );
1919                }
1920#if 0
1921fftfp = fopen( "seqVec2", "w" );
1922fprintf( fftfp, "#after fft\n" );
1923for( k=0; k<n20or4or2; k++ )
1924{
1925   fprintf( fftfp, "#%c\n", amino[k] );
1926   for( l=0; l<nlen; l++ )
1927           fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1928}
1929fclose( fftfp );
1930system( "less seqVec2 < /dev/tty > /dev/tty" );
1931#endif
1932
1933                for( k=0; k<n20or4or2; k++ ) 
1934                {
1935                        for( l=0; l<nlen; l++ ) 
1936                                calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
1937                }
1938                for( l=0; l<nlen; l++ ) 
1939                {
1940                        naisekiNoWa[l].R = 0.0;
1941                        naisekiNoWa[l].I = 0.0;
1942                        for( k=0; k<n20or4or2; k++ ) 
1943                        {
1944                                naisekiNoWa[l].R += naiseki[k][l].R;
1945                                naisekiNoWa[l].I += naiseki[k][l].I;
1946                        }
1947                }
1948       
1949#if 0
1950        fftfp = fopen( "naisekiNoWa", "w" );
1951        fprintf( fftfp, "#Before fft\n" );
1952        for( l=0; l<nlen; l++ )
1953                fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
1954        fclose( fftfp );
1955        system( "less naisekiNoWa < /dev/tty > /dev/tty " );
1956#endif
1957
1958                fft( -nlen, naisekiNoWa, 0 );
1959       
1960                for( m=0; m<=nlen2; m++ ) 
1961                        soukan[m] = naisekiNoWa[nlen2-m].R;
1962                for( m=nlen2+1; m<nlen; m++ ) 
1963                        soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
1964
1965#if 0
1966        fftfp = fopen( "naisekiNoWa", "w" );
1967        fprintf( fftfp, "#After fft\n" );
1968        for( l=0; l<nlen; l++ )
1969                fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R );
1970        fclose( fftfp );
1971        fftfp = fopen( "list.plot", "w"  );
1972        fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
1973        fclose( fftfp );
1974        system( "/usr/bin/gnuplot list.plot &" );
1975#endif
1976#if 0
1977        fprintf( stderr, "soukan\n" );
1978        for( l=0; l<nlen; l++ )
1979                fprintf( stderr, "%d  %f\n", l-nlen2, soukan[l] );
1980#if 0
1981        fftfp = fopen( "list.plot", "w"  );
1982        fprintf( fftfp, "plot 'frt'\n pause +1" );
1983        fclose( fftfp );
1984        system( "/usr/bin/gnuplot list.plot" );
1985#endif
1986#endif
1987
1988
1989                nkouho = getKouho( kouho, NKOUHO_LONG, soukan, nlen );
1990
1991#if 0
1992                for( i=0; i<nkouho; i++ )
1993                {
1994                        fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
1995                }
1996#endif
1997        }
1998
1999#if KEIKA
2000        fprintf( stderr, "Searching anchors ... " );
2001#endif
2002        count = 0;
2003
2004
2005
2006#define CAND 0
2007#if CAND
2008        fftfp = fopen( "cand", "w" );
2009        fclose( fftfp );
2010#endif
2011        if( kobetsubunkatsu )
2012        {
2013                maxk = 1;
2014                kouho[0] = 0;
2015        }
2016        else
2017        {
2018                maxk = nkouho;
2019        }
2020
2021        for( k=0; k<maxk; k++ ) 
2022        {
2023                lag = kouho[k];
2024                if( lag <= -len1 || len2 <= lag ) continue;
2025//              fprintf( stderr, "k=%d, lag=%d\n", k, lag );
2026                zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
2027#if CAND
2028                fftfp = fopen( "cand", "a" );
2029                fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
2030                fprintf( fftfp, "%s\n", tmpptr1[0] );
2031                fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
2032                fprintf( fftfp, "%s\n", tmpptr2[0] );
2033                fprintf( fftfp, ">\n", k+1, lag );
2034                fclose( fftfp );
2035#endif
2036
2037//              fprintf( stderr, "lag = %d\n", lag );
2038                tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
2039//              fprintf( stderr, "lag = %d, %d found\n", lag, tmpint );
2040
2041//              if( lag == -50 ) exit( 1 );
2042               
2043                if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
2044
2045//              fprintf( stderr, "##### k=%d / %d\n", k, maxk );
2046//              if( tmpint == 0 ) break; // 060430 iinoka ? // 090530 yameta
2047                while( tmpint-- > 0 )
2048                {
2049#if 0
2050                        if( segment[count].end - segment[count].start < fftWinSize )
2051                        {
2052                                count++;
2053                                continue;
2054                        }
2055#endif
2056                        if( lag > 0 )
2057                        {
2058                                segment1[count].start  = segment[count].start ;
2059                                segment1[count].end    = segment[count].end   ;
2060                                segment1[count].center = segment[count].center;
2061                                segment1[count].score  = segment[count].score;
2062
2063                                segment2[count].start  = segment[count].start  + lag;
2064                                segment2[count].end    = segment[count].end    + lag;
2065                                segment2[count].center = segment[count].center + lag;
2066                                segment2[count].score  = segment[count].score       ;
2067                        }
2068                        else
2069                        {
2070                                segment1[count].start  = segment[count].start  - lag;
2071                                segment1[count].end    = segment[count].end    - lag;
2072                                segment1[count].center = segment[count].center - lag;
2073                                segment1[count].score  = segment[count].score       ;
2074
2075                                segment2[count].start  = segment[count].start ;
2076                                segment2[count].end    = segment[count].end   ;
2077                                segment2[count].center = segment[count].center;
2078                                segment2[count].score  = segment[count].score ;
2079                        }
2080#if 0
2081                        fprintf( stderr, "##### k=%d / %d\n", k, maxk );
2082                        fprintf( stderr, "anchor %d, score = %f\n", count, segment1[count].score );
2083                        fprintf( stderr, "in 1 %d\n", segment1[count].center );
2084                        fprintf( stderr, "in 2 %d\n", segment2[count].center );
2085#endif
2086                        segment1[count].pair = &segment2[count];
2087                        segment2[count].pair = &segment1[count];
2088                        count++;
2089#if 0
2090                        fprintf( stderr, "count=%d\n", count );
2091#endif
2092                }
2093        }
2094#if 1
2095        if( !kobetsubunkatsu )
2096                fprintf( stderr, "done. (%d anchors) ", count );
2097#endif
2098        if( !count && fftNoAnchStop )
2099                ErrorExit( "Cannot detect anchor!" );
2100#if 0
2101        fprintf( stderr, "RESULT before sort:\n" );
2102        for( l=0; l<count+1; l++ )
2103        {
2104                fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
2105                fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
2106        }
2107#endif
2108
2109        for( i=0; i<count; i++ )
2110        {
2111                sortedseg1[i] = &segment1[i];
2112                sortedseg2[i] = &segment2[i];
2113        }
2114#if 0
2115        tmpsort( count, sortedseg1 );
2116        tmpsort( count, sortedseg2 );
2117        qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
2118        qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
2119#else
2120        mymergesort( 0, count-1, sortedseg1 ); 
2121        mymergesort( 0, count-1, sortedseg2 ); 
2122#endif
2123        for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
2124        for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
2125
2126
2127
2128        if( kobetsubunkatsu )
2129        {
2130                for( i=0; i<count; i++ )
2131            {
2132                        cut1[i+1] = sortedseg1[i]->center;
2133                        cut2[i+1] = sortedseg2[i]->center;
2134                }
2135                cut1[0] = 0;
2136                cut2[0] = 0;
2137                cut1[count+1] = len1;
2138                cut2[count+1] = len2;
2139                count += 2;
2140        }
2141
2142        else
2143        {
2144                if( count < 5000 )
2145                {
2146                        if( crossscoresize < count+2 )
2147                        {
2148                                crossscoresize = count+2;
2149#if 1
2150                                if( fftkeika ) fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
2151#endif
2152                                if( crossscore ) FreeDoubleMtx( crossscore );
2153                                crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
2154                        }
2155                        for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
2156                                crossscore[i][j] = 0.0;
2157                        for( i=0; i<count; i++ )
2158                        {
2159                                crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
2160                                cut1[i+1] = sortedseg1[i]->center;
2161                                cut2[i+1] = sortedseg2[i]->center;
2162                        }
2163       
2164#if 0
2165                        fprintf( stderr, "AFTER SORT\n" );
2166                        for( i=0; i<count+1; i++ ) fprintf( stderr, "%d, %d\n", cut1[i], cut2[i] );
2167                        fprintf( stderr, "crossscore = \n" );
2168                        for( i=0; i<count+1; i++ )
2169                        {
2170                                for( j=0; j<count+1; j++ )
2171                                        fprintf( stderr, "%.0f ", crossscore[i][j] );
2172                                fprintf( stderr, "\n" );
2173                        }
2174#endif
2175
2176                        crossscore[0][0] = 10000000.0;
2177                        cut1[0] = 0; 
2178                        cut2[0] = 0;
2179                        crossscore[count+1][count+1] = 10000000.0;
2180                        cut1[count+1] = len1;
2181                        cut2[count+1] = len2;
2182                        count += 2;
2183                        count0 = count;
2184               
2185//                      fprintf( stderr, "\n\n\ncalling blockAlign2\n\n\n\n" );
2186                        blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
2187       
2188//                      if( count-count0 )
2189//                              fprintf( stderr, "%d unused anchors\n", count0-count );
2190       
2191                        if( !kobetsubunkatsu && fftkeika )
2192                                fprintf( stderr, "%d anchors found\n", count );
2193                        if( fftkeika )
2194                        {
2195                                if( count0 > count )
2196                                {
2197#if 0
2198                                        fprintf( stderr, "\7 REPEAT!? \n" );
2199#else
2200                                        fprintf( stderr, "REPEAT!? \n" ); 
2201#endif
2202                                        if( fftRepeatStop ) exit( 1 );
2203                                }
2204#if KEIKA
2205                                else fprintf( stderr, "done\n" );
2206#endif
2207                        }
2208                }
2209
2210
2211                else
2212                {
2213                        fprintf( stderr, "\nMany anchors were found. The upper-level DP is skipped.\n\n" );
2214
2215                        cut1[0] = 0; 
2216                        cut2[0] = 0;
2217                        count0 = 0;
2218                        for( i=0; i<count; i++ )
2219                        {
2220//                              fprintf( stderr, "i=%d, %d-%d ?\n", i, sortedseg1[i]->center, sortedseg1[i]->pair->center );
2221                                if( sortedseg1[i]->center > cut1[count0]
2222                                 && sortedseg1[i]->pair->center > cut2[count0] )
2223                                {
2224                                        count0++;
2225                                        cut1[count0] = sortedseg1[i]->center;
2226                                        cut2[count0] = sortedseg1[i]->pair->center;
2227                                }
2228                                else
2229                                {
2230                                        if( i && sortedseg1[i]->score > sortedseg1[i-1]->score )
2231                                        {
2232                                                if( sortedseg1[i]->center > cut1[count0-1]
2233                                                 && sortedseg1[i]->pair->center > cut2[count0-1] )
2234                                                {
2235                                                        cut1[count0] = sortedseg1[i]->center;
2236                                                        cut2[count0] = sortedseg1[i]->pair->center;
2237                                                }
2238                                                else
2239                                                {
2240//                                                      count0--;
2241                                                }
2242                                        }
2243                                }
2244                        }
2245//                      if( count-count0 )
2246//                              fprintf( stderr, "%d anchors unused\n", count-count0 );
2247                        cut1[count0+1] = len1;
2248                        cut2[count0+1] = len2;
2249                        count = count0 + 2;
2250                        count0 = count;
2251       
2252                }
2253        }
2254
2255//      exit( 0 );
2256
2257#if 0
2258        fftfp = fopen( "fft", "a" );
2259        fprintf( fftfp, "RESULT after sort:\n" );
2260        for( l=0; l<count; l++ )
2261        {
2262                fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
2263                fprintf( fftfp, "%d\n", segment2[l].center );
2264        }
2265        fclose( fftfp );
2266#endif
2267
2268#if 0
2269        fprintf( stderr, "RESULT after blckalign:\n" );
2270        for( l=0; l<count+1; l++ )
2271        {
2272                fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
2273        }
2274#endif
2275
2276#if 0
2277        fprintf( trap_g, "Devided to %d segments\n", count-1 );
2278        fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
2279#endif
2280
2281        totallen = 0;
2282        for( j=0; j<clus1; j++ ) result1[j][0] = 0;
2283        for( j=0; j<clus2; j++ ) result2[j][0] = 0;
2284        totalscore = 0.0;
2285        *fftlog = -1;
2286        for( i=0; i<count-1; i++ )
2287        {
2288                *fftlog += 1;
2289                if( i == 0 ) headgp = outgap; else headgp = 1;
2290                if( i == count-2 ) tailgp = outgap; else tailgp = 1;
2291
2292                if( cut1[i] )
2293                {
2294//                      getkyokaigap( sgap1, seq1, cut1[i]-1, clus1 );
2295//                      getkyokaigap( sgap2, seq2, cut2[i]-1, clus2 );
2296                        getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
2297                        getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
2298                }
2299                else
2300                {
2301                        for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
2302                        for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
2303                }
2304                if( cut1[i+1] != len1 )
2305                {       
2306                        getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
2307                        getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
2308                }       
2309                else   
2310                {       
2311                        for( j=0; j<clus1; j++ ) egap1[j] = 'o';
2312                        for( j=0; j<clus2; j++ ) egap2[j] = 'o';
2313                }
2314#if DEBUG
2315                fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
2316#else
2317#if 1
2318                fprintf( stderr, "DP %05d / %05d \b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b", i+1, count-1 );
2319#endif
2320#endif
2321                for( j=0; j<clus1; j++ )
2322                {
2323                        strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
2324                        tmpres1[j][cut1[i+1]-cut1[i]] = 0;
2325                }
2326                if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr $B$K8F$P$l$?$H$-(B fftkeika=1
2327//              if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
2328                for( j=0; j<clus2; j++ )
2329                {
2330//                      fprintf( stderr, "### cut2[i+1]-cut2[i] = %d\n", cut2[i+1]-cut2[i] );
2331                        if( cut2[i+1]-cut2[i] <= 0 )
2332                                fprintf( stderr, "### cut2[i+1]=%d, cut2[i]=%d\n", cut2[i+1], cut2[i] );
2333                        strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
2334                        tmpres2[j][cut2[i+1]-cut2[i]] = 0;
2335                }
2336                if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr $B$K8F$P$l$?$H$-(B fftkeika=1
2337//              if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );
2338
2339                if( constraint )
2340                {
2341                        fprintf( stderr, "Not supported\n" );
2342                        exit( 1 );
2343                }
2344#if 0
2345                fprintf( stderr, "i=%d, before alignment", i );
2346                fprintf( stderr, "%4d\n", totallen );
2347                fprintf( stderr, "\n\n" );
2348                for( j=0; j<clus1; j++ )
2349                {
2350                        fprintf( stderr, "%s\n", tmpres1[j] );
2351                }
2352                fprintf( stderr, "-------\n" );
2353                for( j=0; j<clus2; j++ )
2354                {
2355                        fprintf( stderr, "%s\n", tmpres2[j] );
2356                }
2357#endif
2358
2359#if 0
2360                fprintf( stdout, "writing input\n" );
2361                for( j=0; j<clus1; j++ )
2362                {
2363                        fprintf( stdout, ">%d of GROUP1\n", j );
2364                        fprintf( stdout, "%s\n", tmpres1[j] );
2365                }
2366                for( j=0; j<clus2; j++ )
2367                {
2368                        fprintf( stdout, ">%d of GROUP2\n", j );
2369                        fprintf( stdout, "%s\n", tmpres2[j] );
2370                }
2371                fflush( stdout );
2372#endif
2373                switch( alg )
2374                {
2375                        case( 'M' ):
2376                                        totalscore += MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, NULL, 0, NULL, headgp, tailgp );
2377                                break;
2378                        default:
2379                                fprintf( stderr, "alg = %c\n", alg );
2380                                ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
2381                                break;
2382                }
2383
2384                nlen = strlen( tmpres1[0] );
2385                if( totallen + nlen > alloclen )
2386                {
2387                        fprintf( stderr, "totallen=%d +  nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
2388                        ErrorExit( "LENGTH OVER in Falign\n " );
2389                }
2390                for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
2391                for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
2392                totallen += nlen;
2393#if 0
2394                fprintf( stderr, "i=%d", i );
2395                fprintf( stderr, "%4d\n", totallen );
2396                fprintf( stderr, "\n\n" );
2397                for( j=0; j<clus1; j++ )
2398                {
2399                        fprintf( stderr, "%s\n", tmpres1[j] );
2400                }
2401                fprintf( stderr, "-------\n" );
2402                for( j=0; j<clus2; j++ )
2403                {
2404                        fprintf( stderr, "%s\n", tmpres2[j] );
2405                }
2406#endif
2407        }
2408#if KEIKA
2409        fprintf( stderr, "DP ... done   \n" );
2410#endif
2411
2412        for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
2413        for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
2414#if 0
2415        for( j=0; j<clus1; j++ )
2416        {
2417                fprintf( stderr, "%s\n", result1[j] );
2418        }
2419        fprintf( stderr, "- - - - - - - - - - -\n" );
2420        for( j=0; j<clus2; j++ )
2421        {
2422                fprintf( stderr, "%s\n", result2[j] );
2423        }
2424#endif
2425        return( totalscore );
2426}
Note: See TracBrowser for help on using the repository browser.