source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/core/Galign11.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: 14.4 KB
Line 
1#include "mltaln.h"
2#include "dp.h"
3
4#define DEBUG 0
5#define XXXXXXX    0
6#define USE_PENALTY_EX  1
7
8
9#if 1
10static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 ) 
11{
12        char *seq2 = s2[0];
13        int *intptr = amino_dis[(int)s1[0][i1]];
14
15        while( lgth2-- )
16                *match++ = intptr[(int)*seq2++];
17}
18static void match_calc_mtx( int mtx[0x80][0x80], float *match, char **s1, char **s2, int i1, int lgth2 ) 
19{
20        char *seq2 = s2[0];
21        int *intptr = mtx[(int)s1[0][i1]];
22
23        while( lgth2-- )
24                *match++ = intptr[(int)*seq2++];
25}
26#else
27static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
28{
29        int j;
30
31        for( j=0; j<lgth2; j++ )
32                match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
33}
34#endif
35
36static float Atracking( float *lasthorizontalw, float *lastverticalw, 
37                                                char **seq1, char **seq2, 
38                        char **mseq1, char **mseq2, 
39                        int **ijp,
40                                                int tailgp )
41{
42        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
43//      char gap[] = "-";
44        char *gap;
45        gap = newgapstr;
46        lgth1 = strlen( seq1[0] );
47        lgth2 = strlen( seq2[0] );
48        float wm;
49
50
51#if 0
52        for( i=0; i<lgth1; i++ )
53        {
54                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
55        }
56#endif
57 
58    for( i=0; i<lgth1+1; i++ ) 
59    {
60        ijp[i][0] = i + 1;
61    }
62    for( j=0; j<lgth2+1; j++ ) 
63    {
64        ijp[0][j] = -( j + 1 );
65    }
66
67        if( tailgp == 1 )
68                ;
69        else
70        {
71                wm = lastverticalw[0];
72                for( i=0; i<lgth1; i++ )
73                {
74                        if( lastverticalw[i] >= wm )
75                        {
76                                wm = lastverticalw[i];
77                                iin = i; jin = lgth2-1;
78                                ijp[lgth1][lgth2] = +( lgth1 - i );
79                        }
80                }
81                for( j=0; j<lgth2; j++ )
82                {
83                        if( lasthorizontalw[j] >= wm )
84                        {
85                                wm = lasthorizontalw[j];
86                                iin = lgth1-1; jin = j;
87                                ijp[lgth1][lgth2] = -( lgth2 - j );
88                        }
89                }
90        }
91
92
93
94        mseq1[0] += lgth1+lgth2;
95        *mseq1[0] = 0;
96        mseq2[0] += lgth1+lgth2;
97        *mseq2[0] = 0;
98
99
100        iin = lgth1; jin = lgth2;
101        limk = lgth1+lgth2 + 1;
102        for( k=0; k<limk; k++ ) 
103        {
104                if( ijp[iin][jin] < 0 ) 
105                {
106                        ifi = iin-1; jfi = jin+ijp[iin][jin];
107                }
108                else if( ijp[iin][jin] > 0 )
109                {
110                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
111                }
112                else
113                {
114                        ifi = iin-1; jfi = jin-1;
115                }
116                l = iin - ifi;
117                while( --l ) 
118                {
119                        *--mseq1[0] = seq1[0][ifi+l];
120                        *--mseq2[0] = *gap;
121                        k++;
122                }
123                l= jin - jfi;
124                while( --l )
125                {
126                        *--mseq1[0] = *gap;
127                        *--mseq2[0] = seq2[0][jfi+l];
128                        k++;
129                }
130                if( iin <= 0 || jin <= 0 ) break;
131                *--mseq1[0] = seq1[0][ifi];
132                *--mseq2[0] = seq2[0][jfi];
133                k++;
134                iin = ifi; jin = jfi;
135        }
136        return( 0.0 );
137}
138
139
140float G__align11( char **seq1, char **seq2, int alloclen, int headgp, int tailgp )
141/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
142{
143//      int k;
144        register int i, j;
145        int lasti;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
146        int lgth1, lgth2;
147        int resultlen;
148        float wm;   /* int ?????? */
149        float g;
150        float *currentw, *previousw;
151        float fpenalty = (float)penalty;
152#if USE_PENALTY_EX
153        float fpenalty_ex = (float)penalty_ex;
154#endif
155#if 1
156        float *wtmp;
157        int *ijppt;
158        float *mjpt, *prept, *curpt;
159        int *mpjpt;
160#endif
161        static TLS float mi = 0.0;
162        static TLS float *m = NULL;
163        static TLS int **ijp = NULL;
164        static TLS int mpi = 0;
165        static TLS int *mp = NULL;
166        static TLS float *w1 = NULL;
167        static TLS float *w2 = NULL;
168        static TLS float *match = NULL;
169        static TLS float *initverticalw = NULL;    /* kufuu sureba iranai */
170        static TLS float *lastverticalw = NULL;    /* kufuu sureba iranai */
171        static TLS char **mseq1 = NULL;
172        static TLS char **mseq2 = NULL;
173        static TLS char **mseq = NULL;
174        static TLS int **intwork = NULL;
175        static TLS float **floatwork = NULL;
176        static TLS int orlgth1 = 0, orlgth2 = 0;
177
178        if( seq1 == NULL )
179        {
180                if( orlgth1 > 0 && orlgth2 > 0 )
181                {
182                        orlgth1 = 0;
183                        orlgth2 = 0;
184                        if( mseq1 ) free( mseq1 ); mseq1 = NULL;
185                        if( mseq2 ) free( mseq2 ); mseq2 = NULL;
186                        if( w1 ) FreeFloatVec( w1 ); w1 = NULL;
187                        if( w2 ) FreeFloatVec( w2 ); w2 = NULL;
188                        if( match ) FreeFloatVec( match ); match = NULL;
189                        if( initverticalw ) FreeFloatVec( initverticalw ); initverticalw = NULL;
190                        if( lastverticalw ) FreeFloatVec( lastverticalw ); lastverticalw = NULL;
191
192                        if( m ) FreeFloatVec( m ); m = NULL;
193                        if( mp ) FreeIntVec( mp ); mp = NULL;
194
195                        if( mseq ) FreeCharMtx( mseq ); mseq = NULL;
196
197
198
199                        if( floatwork ) FreeFloatMtx( floatwork ); floatwork = NULL;
200                        if( intwork ) FreeIntMtx( intwork ); intwork = NULL;
201                }
202                return( 0.0 );
203        }
204
205        lgth1 = strlen( seq1[0] );
206        lgth2 = strlen( seq2[0] );
207
208#if 0
209        if( lgth1 <= 0 || lgth2 <= 0 )
210        {
211                fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
212        }
213#endif
214
215#if 1
216        if( lgth1 == 0 && lgth2 == 0 )
217                return( 0.0 );
218
219        if( lgth1 == 0 )
220        {
221                seq1[0][lgth2] = 0;
222                while( lgth2 ) seq1[0][--lgth2] = *newgapstr;
223//              fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
224                return( 0.0 );
225        }
226
227        if( lgth2 == 0 )
228        {
229                seq2[0][lgth1] = 0;
230                while( lgth1 ) seq2[0][--lgth1] = *newgapstr;
231//              fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
232                return( 0.0 );
233        }
234#endif
235
236
237        wm = 0.0;
238
239        if( orlgth1 == 0 )
240        {
241                mseq1 = AllocateCharMtx( njob, 0 );
242                mseq2 = AllocateCharMtx( njob, 0 );
243        }
244
245
246
247        if( lgth1 > orlgth1 || lgth2 > orlgth2 )
248        {
249                int ll1, ll2;
250
251                if( orlgth1 > 0 && orlgth2 > 0 )
252                {
253                        FreeFloatVec( w1 );
254                        FreeFloatVec( w2 );
255                        FreeFloatVec( match );
256                        FreeFloatVec( initverticalw );
257                        FreeFloatVec( lastverticalw );
258
259                        FreeFloatVec( m );
260                        FreeIntVec( mp );
261
262                        FreeCharMtx( mseq );
263
264
265
266                        FreeFloatMtx( floatwork );
267                        FreeIntMtx( intwork );
268                }
269
270                ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
271                ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
272
273#if DEBUG
274                fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
275#endif
276
277                w1 = AllocateFloatVec( ll2+2 );
278                w2 = AllocateFloatVec( ll2+2 );
279                match = AllocateFloatVec( ll2+2 );
280
281                initverticalw = AllocateFloatVec( ll1+2 );
282                lastverticalw = AllocateFloatVec( ll1+2 );
283
284                m = AllocateFloatVec( ll2+2 );
285                mp = AllocateIntVec( ll2+2 );
286
287                mseq = AllocateCharMtx( njob, ll1+ll2 );
288
289
290                floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 ); 
291                intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 ); 
292
293#if DEBUG
294                fprintf( stderr, "succeeded\n" );
295#endif
296
297                orlgth1 = ll1 - 100;
298                orlgth2 = ll2 - 100;
299        }
300
301
302        mseq1[0] = mseq[0];
303        mseq2[0] = mseq[1];
304
305
306        if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
307        {
308                int ll1, ll2;
309
310                if( commonAlloc1 && commonAlloc2 )
311                {
312                        FreeIntMtx( commonIP );
313                }
314
315                ll1 = MAX( orlgth1, commonAlloc1 );
316                ll2 = MAX( orlgth2, commonAlloc2 );
317
318#if DEBUG
319                fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
320#endif
321
322                commonIP = AllocateIntMtx( ll1+10, ll2+10 );
323
324#if DEBUG
325                fprintf( stderr, "succeeded\n\n" );
326#endif
327
328                commonAlloc1 = ll1;
329                commonAlloc2 = ll2;
330        }
331        ijp = commonIP;
332
333
334#if 0
335        for( i=0; i<lgth1; i++ )
336                fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
337#endif
338
339        currentw = w1;
340        previousw = w2;
341
342
343        match_calc( initverticalw, seq2, seq1, 0, lgth1 );
344
345
346        match_calc( currentw, seq1, seq2, 0, lgth2 );
347
348        if( headgp == 1 )
349        {
350                for( i=1; i<lgth1+1; i++ )
351                {
352                        initverticalw[i] += fpenalty;
353                }
354                for( j=1; j<lgth2+1; j++ )
355                {
356                        currentw[j] += fpenalty;
357                }
358        }
359
360        for( j=1; j<lgth2+1; ++j ) 
361        {
362                m[j] = currentw[j-1]; mp[j] = 0;
363        }
364
365        if( lgth2 == 0 )
366                lastverticalw[0] = 0.0;               // lgth2==0 no toki error
367        else
368                lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
369
370        if( tailgp ) lasti = lgth1+1; else lasti = lgth1;
371
372#if XXXXXXX
373fprintf( stderr, "currentw = \n" );
374for( i=0; i<lgth1+1; i++ )
375{
376        fprintf( stderr, "%5.2f ", currentw[i] );
377}
378fprintf( stderr, "\n" );
379fprintf( stderr, "initverticalw = \n" );
380for( i=0; i<lgth2+1; i++ )
381{
382        fprintf( stderr, "%5.2f ", initverticalw[i] );
383}
384fprintf( stderr, "\n" );
385#endif
386
387        for( i=1; i<lasti; i++ )
388        {
389                wtmp = previousw; 
390                previousw = currentw;
391                currentw = wtmp;
392
393                previousw[0] = initverticalw[i-1];
394
395                match_calc( currentw, seq1, seq2, i, lgth2 );
396#if XXXXXXX
397fprintf( stderr, "\n" );
398fprintf( stderr, "i=%d\n", i );
399fprintf( stderr, "currentw = \n" );
400for( j=0; j<lgth2; j++ )
401{
402        fprintf( stderr, "%5.2f ", currentw[j] );
403}
404fprintf( stderr, "\n" );
405#endif
406#if XXXXXXX
407fprintf( stderr, "\n" );
408fprintf( stderr, "i=%d\n", i );
409fprintf( stderr, "currentw = \n" );
410for( j=0; j<lgth2; j++ )
411{
412        fprintf( stderr, "%5.2f ", currentw[j] );
413}
414fprintf( stderr, "\n" );
415#endif
416                currentw[0] = initverticalw[i];
417
418                mi = previousw[0]; mpi = 0;
419
420                ijppt = ijp[i] + 1;
421                mjpt = m + 1;
422                prept = previousw;
423                curpt = currentw + 1;
424                mpjpt = mp + 1;
425                for( j=1; j<lgth2+1; j++ )
426                {
427                        wm = *prept;
428                        *ijppt = 0;
429
430#if 0
431                        fprintf( stderr, "%5.0f->", wm );
432#endif
433#if 0
434                        fprintf( stderr, "%5.0f?", g );
435#endif
436                        if( (g=mi+fpenalty) > wm )
437                        {
438                                wm = g;
439                                *ijppt = -( j - mpi );
440                        }
441                        if( (g=*prept) >= mi )
442                        {
443                                mi = g;
444                                mpi = j-1;
445                        }
446#if USE_PENALTY_EX
447                        mi += fpenalty_ex;
448#endif
449
450#if 0 
451                        fprintf( stderr, "%5.0f?", g );
452#endif
453                        if( (g=*mjpt + fpenalty) > wm )
454                        {
455                                wm = g;
456                                *ijppt = +( i - *mpjpt );
457                        }
458                        if( (g=*prept) >= *mjpt )
459                        {
460                                *mjpt = g;
461                                *mpjpt = i-1;
462                        }
463#if USE_PENALTY_EX
464                        m[j] += fpenalty_ex;
465#endif
466
467#if 0
468                        fprintf( stderr, "%5.0f ", wm );
469#endif
470                        *curpt++ += wm;
471                        ijppt++;
472                        mjpt++;
473                        prept++;
474                        mpjpt++;
475                }
476                lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
477        }
478
479        Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, tailgp );
480
481        resultlen = strlen( mseq1[0] );
482        if( alloclen < resultlen || resultlen > N )
483        {
484                fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
485                ErrorExit( "LENGTH OVER!\n" );
486        }
487
488
489        strcpy( seq1[0], mseq1[0] );
490        strcpy( seq2[0], mseq2[0] );
491#if 0
492        fprintf( stderr, "\n" );
493        fprintf( stderr, ">\n%s\n", mseq1[0] );
494        fprintf( stderr, ">\n%s\n", mseq2[0] );
495        fprintf( stderr, "wm = %f\n", wm );
496#endif
497
498        return( wm );
499}
500
501float G__align11_noalign( int scoremtx[0x80][0x80], int penal, int penal_ex, char **seq1, char **seq2, int alloclen )
502/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
503{
504//      int k;
505        register int i, j;
506        int lasti;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
507        int lgth1, lgth2;
508//      int resultlen;
509        float wm;   /* int ?????? */
510        float g;
511        float *currentw, *previousw;
512        float fpenalty = (float)penal;
513#if USE_PENALTY_EX
514        float fpenalty_ex = (float)penal_ex;
515#endif
516#if 1
517        float *wtmp;
518        float *mjpt, *prept, *curpt;
519//      int *mpjpt;
520#endif
521        static TLS float mi, *m;
522        static TLS float *w1, *w2;
523        static TLS float *match;
524        static TLS float *initverticalw;    /* kufuu sureba iranai */
525        static TLS float *lastverticalw;    /* kufuu sureba iranai */
526        static TLS int **intwork;
527        static TLS float **floatwork;
528        static TLS int orlgth1 = 0, orlgth2 = 0;
529
530        if( seq1 == NULL )
531        {
532                if( orlgth1 > 0 && orlgth2 > 0 )
533                {
534                        orlgth1 = 0;
535                        orlgth2 = 0;
536                        FreeFloatVec( w1 );
537                        FreeFloatVec( w2 );
538                        FreeFloatVec( match );
539                        FreeFloatVec( initverticalw );
540                        FreeFloatVec( lastverticalw );
541                        free( m );
542
543
544                        FreeFloatMtx( floatwork );
545                        FreeIntMtx( intwork );
546                }
547                return( 0.0 );
548        }
549
550        wm = 0.0;
551
552
553
554        lgth1 = strlen( seq1[0] );
555        lgth2 = strlen( seq2[0] );
556
557
558
559#if 0
560        if( lgth1 <= 0 || lgth2 <= 0 )
561        {
562                fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
563        }
564#endif
565
566        if( lgth1 > orlgth1 || lgth2 > orlgth2 )
567        {
568                int ll1, ll2;
569
570                if( orlgth1 > 0 && orlgth2 > 0 )
571                {
572                        FreeFloatVec( w1 );
573                        FreeFloatVec( w2 );
574                        FreeFloatVec( match );
575                        FreeFloatVec( initverticalw );
576                        FreeFloatVec( lastverticalw );
577
578                        FreeFloatVec( m );
579
580
581
582
583                        FreeFloatMtx( floatwork );
584                        FreeIntMtx( intwork );
585                }
586
587                ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
588                ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
589
590#if DEBUG
591                fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
592#endif
593
594                w1 = AllocateFloatVec( ll2+2 );
595                w2 = AllocateFloatVec( ll2+2 );
596                match = AllocateFloatVec( ll2+2 );
597
598                initverticalw = AllocateFloatVec( ll1+2 );
599                lastverticalw = AllocateFloatVec( ll1+2 );
600
601                m = AllocateFloatVec( ll2+2 );
602
603
604
605                floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 ); 
606                intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 ); 
607
608#if DEBUG
609                fprintf( stderr, "succeeded\n" );
610#endif
611
612                orlgth1 = ll1 - 100;
613                orlgth2 = ll2 - 100;
614        }
615
616
617
618
619
620
621#if 0
622        for( i=0; i<lgth1; i++ )
623                fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
624#endif
625
626        currentw = w1;
627        previousw = w2;
628
629
630        match_calc_mtx( scoremtx, initverticalw, seq2, seq1, 0, lgth1 );
631
632
633        match_calc_mtx( scoremtx, currentw, seq1, seq2, 0, lgth2 );
634
635        if( 1 ) // tsuneni outgap-1
636        {
637                for( i=1; i<lgth1+1; i++ )
638                {
639                        initverticalw[i] += fpenalty;
640                }
641                for( j=1; j<lgth2+1; j++ )
642                {
643                        currentw[j] += fpenalty;
644                }
645        }
646
647        for( j=1; j<lgth2+1; ++j ) 
648        {
649                m[j] = currentw[j-1];
650        }
651
652        if( lgth2 == 0 )
653                lastverticalw[0] = 0.0;               // lgth2==0 no toki error
654        else
655                lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
656
657        if( 1 ) lasti = lgth1+1; else lasti = lgth1; // tsuneni outgap-1
658
659#if XXXXXXX
660fprintf( stderr, "currentw = \n" );
661for( i=0; i<lgth1+1; i++ )
662{
663        fprintf( stderr, "%5.2f ", currentw[i] );
664}
665fprintf( stderr, "\n" );
666fprintf( stderr, "initverticalw = \n" );
667for( i=0; i<lgth2+1; i++ )
668{
669        fprintf( stderr, "%5.2f ", initverticalw[i] );
670}
671fprintf( stderr, "\n" );
672#endif
673
674        for( i=1; i<lasti; i++ )
675        {
676                wtmp = previousw; 
677                previousw = currentw;
678                currentw = wtmp;
679
680                previousw[0] = initverticalw[i-1];
681
682                match_calc_mtx( scoremtx, currentw, seq1, seq2, i, lgth2 );
683#if XXXXXXX
684fprintf( stderr, "\n" );
685fprintf( stderr, "i=%d\n", i );
686fprintf( stderr, "currentw = \n" );
687for( j=0; j<lgth2; j++ )
688{
689        fprintf( stderr, "%5.2f ", currentw[j] );
690}
691fprintf( stderr, "\n" );
692#endif
693#if XXXXXXX
694fprintf( stderr, "\n" );
695fprintf( stderr, "i=%d\n", i );
696fprintf( stderr, "currentw = \n" );
697for( j=0; j<lgth2; j++ )
698{
699        fprintf( stderr, "%5.2f ", currentw[j] );
700}
701fprintf( stderr, "\n" );
702#endif
703                currentw[0] = initverticalw[i];
704
705                mi = previousw[0];
706
707                mjpt = m + 1;
708                prept = previousw;
709                curpt = currentw + 1;
710                for( j=1; j<lgth2+1; j++ )
711                {
712                        wm = *prept;
713
714#if 0
715                        fprintf( stderr, "%5.0f->", wm );
716#endif
717#if 0
718                        fprintf( stderr, "%5.0f?", g );
719#endif
720                        if( (g=mi+fpenalty) > wm )
721                        {
722                                wm = g;
723                        }
724                        if( (g=*prept) >= mi )
725                        {
726                                mi = g;
727                        }
728#if USE_PENALTY_EX
729                        mi += fpenalty_ex;
730#endif
731
732#if 0 
733                        fprintf( stderr, "%5.0f?", g );
734#endif
735                        if( (g=*mjpt + fpenalty) > wm )
736                        {
737                                wm = g;
738                        }
739                        if( (g=*prept) >= *mjpt )
740                        {
741                                *mjpt = g;
742                        }
743#if USE_PENALTY_EX
744                        m[j] += fpenalty_ex;
745#endif
746
747#if 0
748                        fprintf( stderr, "%5.0f ", wm );
749#endif
750                        *curpt++ += wm;
751                        mjpt++;
752                        prept++;
753                }
754                lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
755        }
756
757#if 0
758        fprintf( stderr, "\n" );
759        fprintf( stderr, ">\n%s\n", mseq1[0] );
760        fprintf( stderr, ">\n%s\n", mseq2[0] );
761        fprintf( stderr, "wm = %f\n", wm );
762#endif
763
764        return( wm );
765}
Note: See TracBrowser for help on using the repository browser.