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