source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/core/genGalign11.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: 9.1 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}
18#else
19static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
20{
21        int j;
22
23        for( j=0; j<lgth2; j++ )
24                match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
25}
26#endif
27
28static float genGtracking( float *lasthorizontalw, float *lastverticalw, 
29                                                char **seq1, char **seq2, 
30                        char **mseq1, char **mseq2, 
31                        float **cpmx1, float **cpmx2, 
32                        int **ijpi, int **ijpj )
33{
34        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
35//      char gap[] = "-";
36        char *gap;
37        gap = newgapstr;
38        lgth1 = strlen( seq1[0] );
39        lgth2 = strlen( seq2[0] );
40
41
42#if 0
43        for( i=0; i<lgth1; i++ )
44        {
45                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
46        }
47#endif
48 
49    for( i=0; i<lgth1+1; i++ ) 
50    {
51        ijpi[i][0] = -1;
52                ijpj[i][0] = -1; // ???
53    }
54    for( j=0; j<lgth2+1; j++ ) 
55    {
56                ijpi[0][j] = -1; // ???
57        ijpj[0][j] = -1;
58    }
59
60
61        mseq1[0] += lgth1+lgth2;
62        *mseq1[0] = 0;
63        mseq2[0] += lgth1+lgth2;
64        *mseq2[0] = 0;
65        iin = lgth1; jin = lgth2;
66        limk = lgth1+lgth2 + 1;
67        for( k=0; k<limk; k++ ) 
68        {
69                ifi = ( ijpi[iin][jin] );
70                jfi = ( ijpj[iin][jin] );
71                l = iin - ifi;
72                while( --l ) 
73                {
74                        *--mseq1[0] = seq1[0][ifi+l];
75                        *--mseq2[0] = *gap;
76                        k++;
77                }
78                l= jin - jfi;
79                while( --l )
80                {
81                        *--mseq1[0] = *gap;
82                        *--mseq2[0] = seq2[0][jfi+l];
83                        k++;
84                }
85                if( iin <= 0 || jin <= 0 ) break;
86                *--mseq1[0] = seq1[0][ifi];
87                *--mseq2[0] = seq2[0][jfi];
88
89//              fprintf( stdout, "mseq1 = %s\n", mseq1[0] );
90//              fprintf( stdout, "mseq2 = %s\n", mseq2[0] );
91
92                k++;
93                iin = ifi; jin = jfi;
94        }
95        return( 0.0 );
96}
97
98
99float genG__align11( char **seq1, char **seq2, int alloclen )
100/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
101{
102//      int k;
103        register int i, j;
104        int lasti;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
105        int lgth1, lgth2;
106        int resultlen;
107        float wm;   /* int ?????? */
108        float g;
109        float *currentw, *previousw;
110        float fpenalty = (float)penalty;
111        float fpenalty_OP = (float)penalty_OP;
112#if USE_PENALTY_EX
113        float fpenalty_ex = (float)penalty_ex;
114#endif
115#if 1
116        float *wtmp;
117        int *ijpipt;
118        int *ijpjpt;
119        float *mjpt, *Mjpt, *prept, *curpt;
120        int *mpjpt, *Mpjpt;
121#endif
122        static float mi, *m;
123        static float Mi, *largeM;
124        static int **ijpi;
125        static int **ijpj;
126        static int mpi, *mp;
127        static int Mpi, *Mp;
128        static float *w1, *w2;
129        static float *match;
130        static float *initverticalw;    /* kufuu sureba iranai */
131        static float *lastverticalw;    /* kufuu sureba iranai */
132        static char **mseq1;
133        static char **mseq2;
134        static char **mseq;
135        static float **cpmx1;
136        static float **cpmx2;
137        static int **intwork;
138        static float **floatwork;
139        static int orlgth1 = 0, orlgth2 = 0;
140        float tbk;
141        int tbki, tbkj;
142
143        wm = 0.0;
144
145        if( orlgth1 == 0 )
146        {
147                mseq1 = AllocateCharMtx( njob, 0 );
148                mseq2 = AllocateCharMtx( njob, 0 );
149        }
150
151
152        lgth1 = strlen( seq1[0] );
153        lgth2 = strlen( seq2[0] );
154
155
156
157        if( lgth1 <= 0 || lgth2 <= 0 )
158        {
159                fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
160        }
161
162        if( lgth1 > orlgth1 || lgth2 > orlgth2 )
163        {
164                int ll1, ll2;
165
166                if( orlgth1 > 0 && orlgth2 > 0 )
167                {
168                        FreeFloatVec( w1 );
169                        FreeFloatVec( w2 );
170                        FreeFloatVec( match );
171                        FreeFloatVec( initverticalw );
172                        FreeFloatVec( lastverticalw );
173
174                        FreeFloatVec( m );
175                        FreeIntVec( mp );
176                        FreeFloatVec( largeM );
177                        FreeIntVec( Mp );
178
179                        FreeCharMtx( mseq );
180
181
182                        FreeFloatMtx( cpmx1 );
183                        FreeFloatMtx( cpmx2 );
184
185                        FreeFloatMtx( floatwork );
186                        FreeIntMtx( intwork );
187                }
188
189                ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
190                ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
191
192#if DEBUG
193                fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
194#endif
195
196                w1 = AllocateFloatVec( ll2+2 );
197                w2 = AllocateFloatVec( ll2+2 );
198                match = AllocateFloatVec( ll2+2 );
199
200                initverticalw = AllocateFloatVec( ll1+2 );
201                lastverticalw = AllocateFloatVec( ll1+2 );
202
203                m = AllocateFloatVec( ll2+2 );
204                mp = AllocateIntVec( ll2+2 );
205                largeM = AllocateFloatVec( ll2+2 );
206                Mp = AllocateIntVec( ll2+2 );
207
208                mseq = AllocateCharMtx( njob, ll1+ll2 );
209
210                cpmx1 = AllocateFloatMtx( 26, ll1+2 );
211                cpmx2 = AllocateFloatMtx( 26, ll2+2 );
212
213                floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 ); 
214                intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 ); 
215
216#if DEBUG
217                fprintf( stderr, "succeeded\n" );
218#endif
219
220                orlgth1 = ll1 - 100;
221                orlgth2 = ll2 - 100;
222        }
223
224
225        mseq1[0] = mseq[0];
226        mseq2[0] = mseq[1];
227
228
229        if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
230        {
231                int ll1, ll2;
232
233                if( commonAlloc1 && commonAlloc2 )
234                {
235                        FreeIntMtx( commonIP );
236                        FreeIntMtx( commonJP );
237                }
238
239                ll1 = MAX( orlgth1, commonAlloc1 );
240                ll2 = MAX( orlgth2, commonAlloc2 );
241
242#if DEBUG
243                fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
244#endif
245
246                commonIP = AllocateIntMtx( ll1+10, ll2+10 );
247                commonJP = AllocateIntMtx( ll1+10, ll2+10 );
248
249#if DEBUG
250                fprintf( stderr, "succeeded\n\n" );
251#endif
252
253                commonAlloc1 = ll1;
254                commonAlloc2 = ll2;
255        }
256        ijpi = commonIP;
257        ijpj = commonJP;
258
259
260#if 0
261        for( i=0; i<lgth1; i++ )
262                fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
263#endif
264
265        currentw = w1;
266        previousw = w2;
267
268
269        match_calc( initverticalw, seq2, seq1, 0, lgth1 );
270
271
272        match_calc( currentw, seq1, seq2, 0, lgth2 );
273
274        if( outgap == 1 )
275        {
276                for( i=1; i<lgth1+1; i++ )
277                {
278                        initverticalw[i] += fpenalty;
279                }
280                for( j=1; j<lgth2+1; j++ )
281                {
282                        currentw[j] += fpenalty;
283                }
284        }
285
286        for( j=1; j<lgth2+1; ++j ) 
287        {
288                m[j] = currentw[j-1]; mp[j] = 0;
289                largeM[j] = currentw[j-1]; Mp[j] = 0;
290        }
291
292        if( lgth2 == 0 )
293                lastverticalw[0] = 0.0;               // lgth2==0 no toki error
294        else
295                lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
296
297        if( outgap ) lasti = lgth1+1; else lasti = lgth1;
298
299#if XXXXXXX
300fprintf( stderr, "currentw = \n" );
301for( i=0; i<lgth1+1; i++ )
302{
303        fprintf( stderr, "%5.2f ", currentw[i] );
304}
305fprintf( stderr, "\n" );
306fprintf( stderr, "initverticalw = \n" );
307for( i=0; i<lgth2+1; i++ )
308{
309        fprintf( stderr, "%5.2f ", initverticalw[i] );
310}
311fprintf( stderr, "\n" );
312#endif
313
314        for( i=1; i<lasti; i++ )
315        {
316                wtmp = previousw; 
317                previousw = currentw;
318                currentw = wtmp;
319
320                previousw[0] = initverticalw[i-1];
321
322                match_calc( currentw, seq1, seq2, i, lgth2 );
323#if XXXXXXX
324fprintf( stderr, "\n" );
325fprintf( stderr, "i=%d\n", i );
326fprintf( stderr, "currentw = \n" );
327for( j=0; j<lgth2; j++ )
328{
329        fprintf( stderr, "%5.2f ", currentw[j] );
330}
331fprintf( stderr, "\n" );
332#endif
333#if XXXXXXX
334fprintf( stderr, "\n" );
335fprintf( stderr, "i=%d\n", i );
336fprintf( stderr, "currentw = \n" );
337for( j=0; j<lgth2; j++ )
338{
339        fprintf( stderr, "%5.2f ", currentw[j] );
340}
341fprintf( stderr, "\n" );
342#endif
343                currentw[0] = initverticalw[i];
344
345                mi = previousw[0]; mpi = 0;
346                Mi = previousw[0]; Mpi = 0;
347
348                ijpipt = ijpi[i] + 1;
349                ijpjpt = ijpj[i] + 1;
350                mjpt = m + 1;
351                Mjpt = largeM + 1;
352                prept = previousw;
353                curpt = currentw + 1;
354                mpjpt = mp + 1;
355                Mpjpt = Mp + 1;
356                tbk = -9999999.9;
357                tbki = 0;
358                tbkj = 0;
359                for( j=1; j<lgth2+1; j++ )
360                {
361                        wm = *prept;
362                        *ijpipt = i-1;
363                        *ijpjpt = j-1;
364
365#if 0
366                        fprintf( stderr, "%5.0f->", wm );
367#endif
368#if 0
369                        fprintf( stderr, "%5.0f?", g );
370#endif
371                        if( (g=mi+fpenalty) > wm )
372                        {
373                                wm = g;
374//                              *ijpipt = i - 1; // iranai
375                                *ijpjpt = mpi;
376                        }
377                        if( (g=*prept) >= mi )
378                        {
379                                mi = g;
380                                mpi = j-1;
381                        }
382#if USE_PENALTY_EX
383                        mi += fpenalty_ex;
384#endif
385
386#if 0 
387                        fprintf( stderr, "%5.0f?", g );
388#endif
389                        if( (g=*mjpt + fpenalty) > wm )
390                        {
391                                wm = g;
392                                *ijpipt = *mpjpt;
393                                *ijpjpt = j - 1; //IRU!
394                        }
395                        if( (g=*prept) >= *mjpt )
396                        {
397                                *mjpt = g;
398                                *mpjpt = i-1;
399                        }
400#if USE_PENALTY_EX
401                        m[j] += fpenalty_ex;
402#endif
403
404#if 1
405                        g =  tbk + fpenalty_OP;
406                        if( g > wm )
407                        {
408                                wm = g;
409                                *ijpipt = tbki;
410                                *ijpjpt = tbkj;
411//                              fprintf( stderr, "hit! i%d, j%d, ijpi = %d, ijpj = %d\n", i, j, *ijpipt, *ijpjpt );
412                        }
413                        if( Mi > tbk )
414                        {
415                                tbk = Mi; //error desu.
416                                tbki = i-1;
417                                tbkj = Mpi;
418                        }
419                        if( *Mjpt > tbk )
420                        {
421                                tbk = *Mjpt;
422                                tbki = *Mpjpt;
423                                tbkj = j-1;
424                        }
425
426                        if( *prept > *Mjpt )
427                        {
428                                *Mjpt = *prept;
429                                *Mpjpt = i-1;
430                        }
431                        if( *prept > Mi )
432                        {
433                                Mi = *prept;
434                                Mpi = j-1;
435                        }
436
437#endif
438
439
440#if 0
441                        fprintf( stderr, "%5.0f ", wm );
442#endif
443                        *curpt++ += wm;
444                        ijpipt++;
445                        ijpjpt++;
446                        mjpt++;
447                        Mjpt++;
448                        prept++;
449                        mpjpt++;
450                        Mpjpt++;
451                }
452                lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
453        }
454#if 0
455        for( i=0; i<lgth1; i++ )
456        {
457                for( j=0; j<lgth2; j++ )
458                {
459                        fprintf( stdout, "i,j=%d,%d - ijpi,ijpj=%d,%d\n", i, j, ijpi[i][j], ijpj[i][j] );
460                }
461        }
462        fflush( stdout );
463#endif
464
465        genGtracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijpi, ijpj );
466
467        resultlen = strlen( mseq1[0] );
468        if( alloclen < resultlen || resultlen > N )
469        {
470                fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
471                ErrorExit( "LENGTH OVER!\n" );
472        }
473
474
475        strcpy( seq1[0], mseq1[0] );
476        strcpy( seq2[0], mseq2[0] );
477#if 0
478        fprintf( stderr, "\n" );
479        fprintf( stderr, ">\n%s\n", mseq1[0] );
480        fprintf( stderr, ">\n%s\n", mseq2[0] );
481        fprintf( stderr, "wm = %f\n", wm );
482#endif
483
484        return( wm );
485}
486
Note: See TracBrowser for help on using the repository browser.