source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/core/Lalign11.c

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

updated mafft version. Added extensions (no svn ignore, yet)
File size: 18.2 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
9
10static TLS int localstop; // 060910
11
12#if 1
13static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 ) 
14{
15        char *seq2 = s2[0];
16        int *intptr;
17
18        intptr = amino_dis[(int)s1[0][i1]];
19        while( lgth2-- )
20                *match++ = intptr[(int)*seq2++];
21}
22#else
23static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
24{
25        int j;
26
27        for( j=0; j<lgth2; j++ )
28                match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
29}
30#endif
31
32#if 0
33static void match_calc_bk( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
34{
35        int j, k, l;
36        float scarr[26];
37        float **cpmxpd = floatwork;
38        int **cpmxpdn = intwork;
39        int count = 0;
40
41        if( initialize )
42        {
43                for( j=0; j<lgth2; j++ )
44                {
45                        count = 0;
46                        for( l=0; l<26; l++ )
47                        {
48                                if( cpmx2[l][j] )
49                                {
50                                        cpmxpd[count][j] = cpmx2[l][j];
51                                        cpmxpdn[count][j] = l;
52                                        count++;
53                                }
54                        }
55                        cpmxpdn[count][j] = -1;
56                }
57        }
58
59        for( l=0; l<26; l++ )
60        {
61                scarr[l] = 0.0;
62                for( k=0; k<26; k++ )
63                        scarr[l] += n_dis[k][l] * cpmx1[k][i1];
64        }
65#if 0 /* €³€ì€ò»È€Š€È€­€Ïfloatwork€Î¥¢¥í¥±¡Œ¥È€òµÕ€Ë€¹€ë */
66        {
67                float *fpt, **fptpt, *fpt2;
68                int *ipt, **iptpt;
69                fpt2 = match;
70                iptpt = cpmxpdn;
71                fptpt = cpmxpd;
72                while( lgth2-- )
73                {
74                        *fpt2 = 0.0;
75                        ipt=*iptpt,fpt=*fptpt;
76                        while( *ipt > -1 )
77                                *fpt2 += scarr[*ipt++] * *fpt++;
78                        fpt2++,iptpt++,fptpt++;
79                }
80        }
81#else
82        for( j=0; j<lgth2; j++ )
83        {
84                match[j] = 0.0;
85                for( k=0; cpmxpdn[k][j]>-1; k++ )
86                        match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
87        }
88#endif
89}
90#endif
91
92static float Ltracking( float *lasthorizontalw, float *lastverticalw, 
93                                                char **seq1, char **seq2, 
94                        char **mseq1, char **mseq2, 
95                        int **ijp, int *off1pt, int *off2pt, int endi, int endj )
96{
97        int i, j, l, iin, jin, lgth1, lgth2, k, limk;
98        int ifi=0, jfi=0; // by D.Mathog, a guess
99//      char gap[] = "-";
100        char *gap;
101        gap = newgapstr;
102        lgth1 = strlen( seq1[0] );
103        lgth2 = strlen( seq2[0] );
104
105#if 0
106        for( i=0; i<lgth1; i++ )
107        {
108                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
109        }
110#endif
111 
112    for( i=0; i<lgth1+1; i++ ) 
113    {
114        ijp[i][0] = localstop;
115    }
116    for( j=0; j<lgth2+1; j++ ) 
117    {
118        ijp[0][j] = localstop;
119    }
120
121        mseq1[0] += lgth1+lgth2;
122        *mseq1[0] = 0;
123        mseq2[0] += lgth1+lgth2;
124        *mseq2[0] = 0;
125        iin = endi; jin = endj;
126        limk = lgth1+lgth2;
127        for( k=0; k<=limk; k++ ) 
128        {
129                if( ijp[iin][jin] < 0 ) 
130                {
131                        ifi = iin-1; jfi = jin+ijp[iin][jin];
132                }
133                else if( ijp[iin][jin] > 0 )
134                {
135                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
136                }
137                else
138                {
139                        ifi = iin-1; jfi = jin-1;
140                }
141                l = iin - ifi;
142                while( --l ) 
143                {
144                        *--mseq1[0] = seq1[0][ifi+l];
145                        *--mseq2[0] = *gap;
146                        k++;
147                }
148                l= jin - jfi;
149                while( --l )
150                {
151                        *--mseq1[0] = *gap;
152                        *--mseq2[0] = seq2[0][jfi+l];
153                        k++;
154                }
155
156                if( iin <= 0 || jin <= 0 ) break;
157                *--mseq1[0] = seq1[0][ifi];
158                *--mseq2[0] = seq2[0][jfi];
159                if( ijp[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 L__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        int i, j;
178        int lasti, lastj;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
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 *ijppt;
187        float *mjpt, *prept, *curpt;
188        int *mpjpt;
189#endif
190        static TLS float mi, *m;
191        static TLS int **ijp;
192        static TLS int mpi, *mp;
193        static TLS float *w1, *w2;
194        static TLS float *match;
195        static TLS float *initverticalw;    /* kufuu sureba iranai */
196        static TLS float *lastverticalw;    /* kufuu sureba iranai */
197        static TLS char **mseq1;
198        static TLS char **mseq2;
199        static TLS char **mseq;
200//      static TLS int **intwork;
201//      static TLS float **floatwork;
202        static TLS int orlgth1 = 0, orlgth2 = 0;
203        float maxwm;
204        int endali = 0, endalj = 0; // by D.Mathog, a guess
205//      int endali, endalj;
206        float localthr = -offset;
207        float localthr2 = -offset;
208//      float localthr = 100;
209//      float localthr2 = 100;
210        float fpenalty = (float)penalty;
211        float fpenalty_ex = (float)penalty_ex;
212
213        if( seq1 == NULL )
214        {
215                if( orlgth1 > 0 && orlgth2 > 0 )
216                {
217                        orlgth1 = 0;
218                        orlgth2 = 0;
219                        free( mseq1 );
220                        free( mseq2 );
221                        FreeFloatVec( w1 );
222                        FreeFloatVec( w2 );
223                        FreeFloatVec( match );
224                        FreeFloatVec( initverticalw );
225                        FreeFloatVec( lastverticalw );
226
227                        FreeFloatVec( m );
228                        FreeIntVec( mp );
229
230                        FreeCharMtx( mseq );
231
232                }
233                return( 0.0 );
234        }
235
236
237        if( orlgth1 == 0 )
238        {
239                mseq1 = AllocateCharMtx( njob, 0 );
240                mseq2 = AllocateCharMtx( njob, 0 );
241        }
242
243
244        lgth1 = strlen( seq1[0] );
245        lgth2 = strlen( seq2[0] );
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        match_calc( initverticalw, seq2, seq1, 0, lgth1 );
343
344        match_calc( currentw, seq1, seq2, 0, lgth2 );
345
346
347        lasti = lgth2+1;
348        for( j=1; j<lasti; ++j ) 
349        {
350                m[j] = currentw[j-1]; mp[j] = 0;
351#if 0
352                if( m[j] < localthr ) m[j] = localthr2;
353#endif
354        }
355
356        lastverticalw[0] = currentw[lgth2-1];
357
358        lasti = lgth1+1;
359
360#if 0
361fprintf( stderr, "currentw = \n" );
362for( i=0; i<lgth1+1; i++ )
363{
364        fprintf( stderr, "%5.2f ", currentw[i] );
365}
366fprintf( stderr, "\n" );
367fprintf( stderr, "initverticalw = \n" );
368for( i=0; i<lgth2+1; i++ )
369{
370        fprintf( stderr, "%5.2f ", initverticalw[i] );
371}
372fprintf( stderr, "\n" );
373#endif
374#if DEBUG2
375        fprintf( stderr, "\n" );
376        fprintf( stderr, "       " );
377        for( j=0; j<lgth2; j++ )
378                fprintf( stderr, "%c     ", seq2[0][j] );
379        fprintf( stderr, "\n" );
380#endif
381
382        localstop = lgth1+lgth2+1;
383        maxwm = -999999999.9;
384#if DEBUG2
385        fprintf( stderr, "\n" );
386        fprintf( stderr, "%c   ", seq1[0][0] );
387
388        for( j=0; j<lgth2+1; j++ )
389                fprintf( stderr, "%5.0f ", currentw[j] );
390        fprintf( stderr, "\n" );
391#endif
392
393        for( i=1; i<lasti; i++ )
394        {
395                wtmp = previousw; 
396                previousw = currentw;
397                currentw = wtmp;
398
399                previousw[0] = initverticalw[i-1];
400
401                match_calc( currentw, seq1, seq2, i, lgth2 );
402#if DEBUG2
403                fprintf( stderr, "%c   ", seq1[0][i] );
404                fprintf( stderr, "%5.0f ", currentw[0] );
405#endif
406
407#if XXXXXXX
408fprintf( stderr, "\n" );
409fprintf( stderr, "i=%d\n", i );
410fprintf( stderr, "currentw = \n" );
411for( j=0; j<lgth2; j++ )
412{
413        fprintf( stderr, "%5.2f ", currentw[j] );
414}
415fprintf( stderr, "\n" );
416#endif
417#if XXXXXXX
418fprintf( stderr, "\n" );
419fprintf( stderr, "i=%d\n", i );
420fprintf( stderr, "currentw = \n" );
421for( j=0; j<lgth2; j++ )
422{
423        fprintf( stderr, "%5.2f ", currentw[j] );
424}
425fprintf( stderr, "\n" );
426#endif
427                currentw[0] = initverticalw[i];
428
429                mi = previousw[0]; mpi = 0;
430
431#if 0
432                if( mi < localthr ) mi = localthr2;
433#endif
434
435                ijppt = ijp[i] + 1;
436                mjpt = m + 1;
437                prept = previousw;
438                curpt = currentw + 1;
439                mpjpt = mp + 1;
440                lastj = lgth2+1;
441                for( j=1; j<lastj; j++ )
442                {
443                        wm = *prept;
444                        *ijppt = 0;
445
446#if 0
447                        fprintf( stderr, "%5.0f->", wm );
448#endif
449#if 0
450                        fprintf( stderr, "%5.0f?", g );
451#endif
452                        if( (g=mi+fpenalty) > wm )
453                        {
454                                wm = g;
455                                *ijppt = -( j - mpi );
456                        }
457                        if( *prept > mi )
458                        {
459                                mi = *prept;
460                                mpi = j-1;
461                        }
462
463#if USE_PENALTY_EX
464                        mi += fpenalty_ex;
465#endif
466
467#if 0 
468                        fprintf( stderr, "%5.0f?", g );
469#endif
470                        if( (g=*mjpt+fpenalty) > wm )
471                        {
472                                wm = g;
473                                *ijppt = +( i - *mpjpt );
474                        }
475                        if( *prept > *mjpt )
476                        {
477                                *mjpt = *prept;
478                                *mpjpt = i-1;
479                        }
480#if USE_PENALTY_EX
481                        *mjpt += fpenalty_ex;
482#endif
483
484                        if( maxwm < wm )
485                        {
486                                maxwm = wm;
487                                endali = i;
488                                endalj = j;
489                        }
490#if 1
491                        if( wm < localthr )
492                        {
493//                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
494                                *ijppt = localstop;
495                                wm = localthr2;
496                        }
497#endif
498#if 0
499                        fprintf( stderr, "%5.0f ", *curpt );
500#endif
501#if DEBUG2
502                        fprintf( stderr, "%5.0f ", wm );
503//                      fprintf( stderr, "%c-%c *ijppt = %d, localstop = %d\n", seq1[0][i], seq2[0][j], *ijppt, localstop );
504#endif
505
506                        *curpt++ += wm;
507                        ijppt++;
508                        mjpt++;
509                        prept++;
510                        mpjpt++;
511                }
512#if DEBUG2
513                fprintf( stderr, "\n" );
514#endif
515
516                lastverticalw[i] = currentw[lgth2-1];
517        }
518
519
520#if 0
521        fprintf( stderr, "maxwm = %f\n", maxwm );
522        fprintf( stderr, "endali = %d\n", endali );
523        fprintf( stderr, "endalj = %d\n", endalj );
524#endif
525
526        if( ijp[endali][endalj] == localstop )
527        {
528                strcpy( seq1[0], "" );
529                strcpy( seq2[0], "" );
530                *off1pt = *off2pt = 0;
531                fprintf( stderr, "maxwm <- 0.0 \n" );
532                return( 0.0 );
533        }
534               
535        Ltracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, off1pt, off2pt, endali, endalj );
536
537
538        resultlen = strlen( mseq1[0] );
539        if( alloclen < resultlen || resultlen > N )
540        {
541                fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
542                ErrorExit( "LENGTH OVER!\n" );
543        }
544
545
546        strcpy( seq1[0], mseq1[0] );
547        strcpy( seq2[0], mseq2[0] );
548
549#if 0
550        fprintf( stderr, "wm=%f\n", wm );
551        fprintf( stderr, ">\n%s\n", mseq1[0] );
552        fprintf( stderr, ">\n%s\n", mseq2[0] );
553
554        fprintf( stderr, "maxwm = %f\n", maxwm );
555        fprintf( stderr, "   wm = %f\n",    wm );
556#endif
557
558        return( maxwm );
559}
560
561
562float L__align11_noalign( char **seq1, char **seq2 )
563/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
564{
565//      int k;
566        int i, j;
567        int lasti, lastj;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
568        int lgth1, lgth2;
569//      int resultlen;
570        float wm = 0.0;   /* int ?????? */
571        float g;
572        float *currentw, *previousw;
573#if 1
574        float *wtmp;
575//      int *ijppt;
576        float *mjpt, *prept, *curpt;
577//      int *mpjpt;
578#endif
579        static TLS float mi, *m;
580//      static TLS int **ijp;
581//      static TLS int mpi, *mp;
582        static TLS float *w1, *w2;
583        static TLS float *match;
584        static TLS float *initverticalw;    /* kufuu sureba iranai */
585        static TLS float *lastverticalw;    /* kufuu sureba iranai */
586//      static TLS char **mseq1;
587//      static TLS char **mseq2;
588//      static TLS char **mseq;
589//      static TLS int **intwork;
590//      static TLS float **floatwork;
591        static TLS int orlgth1 = 0, orlgth2 = 0;
592        float maxwm;
593//      int endali = 0, endalj = 0; // by D.Mathog, a guess
594//      int endali, endalj;
595        float localthr = -offset;
596        float localthr2 = -offset;
597//      float localthr = 100;
598//      float localthr2 = 100;
599        float fpenalty = (float)penalty;
600        float fpenalty_ex = (float)penalty_ex;
601
602        if( seq1 == NULL )
603        {
604                if( orlgth1 > 0 && orlgth2 > 0 )
605                {
606                        orlgth1 = 0;
607                        orlgth2 = 0;
608//                      free( mseq1 );
609//                      free( mseq2 );
610                        FreeFloatVec( w1 );
611                        FreeFloatVec( w2 );
612                        FreeFloatVec( match );
613                        FreeFloatVec( initverticalw );
614                        FreeFloatVec( lastverticalw );
615
616                        FreeFloatVec( m );
617//                      FreeIntVec( mp );
618
619//                      FreeCharMtx( mseq );
620
621                }
622                return( 0.0 );
623        }
624
625
626//      if( orlgth1 == 0 )
627//      {
628//              mseq1 = AllocateCharMtx( njob, 0 );
629//              mseq2 = AllocateCharMtx( njob, 0 );
630//      }
631
632
633        lgth1 = strlen( seq1[0] );
634        lgth2 = strlen( seq2[0] );
635
636        if( lgth1 > orlgth1 || lgth2 > orlgth2 )
637        {
638                int ll1, ll2;
639
640                if( orlgth1 > 0 && orlgth2 > 0 )
641                {
642                        FreeFloatVec( w1 );
643                        FreeFloatVec( w2 );
644                        FreeFloatVec( match );
645                        FreeFloatVec( initverticalw );
646                        FreeFloatVec( lastverticalw );
647
648                        FreeFloatVec( m );
649//                      FreeIntVec( mp );
650
651//                      FreeCharMtx( mseq );
652
653
654
655//                      FreeFloatMtx( floatwork );
656//                      FreeIntMtx( intwork );
657                }
658
659                ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
660                ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
661
662#if DEBUG
663                fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
664#endif
665
666                w1 = AllocateFloatVec( ll2+2 );
667                w2 = AllocateFloatVec( ll2+2 );
668                match = AllocateFloatVec( ll2+2 );
669
670                initverticalw = AllocateFloatVec( ll1+2 );
671                lastverticalw = AllocateFloatVec( ll1+2 );
672
673                m = AllocateFloatVec( ll2+2 );
674//              mp = AllocateIntVec( ll2+2 );
675
676//              mseq = AllocateCharMtx( njob, ll1+ll2 );
677
678
679//              floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
680//              intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
681
682#if DEBUG
683                fprintf( stderr, "succeeded\n" );
684#endif
685
686                orlgth1 = ll1 - 100;
687                orlgth2 = ll2 - 100;
688        }
689
690
691//      mseq1[0] = mseq[0];
692//      mseq2[0] = mseq[1];
693
694
695//      if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
696//      {
697//              int ll1, ll2;
698//
699//              if( commonAlloc1 && commonAlloc2 )
700//              {
701//                      FreeIntMtx( commonIP );
702//              }
703//
704//              ll1 = MAX( orlgth1, commonAlloc1 );
705//              ll2 = MAX( orlgth2, commonAlloc2 );
706
707#if DEBUG
708//              fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
709#endif
710
711//              commonIP = AllocateIntMtx( ll1+10, ll2+10 );
712
713#if DEBUG
714//              fprintf( stderr, "succeeded\n\n" );
715#endif
716
717//              commonAlloc1 = ll1;
718//              commonAlloc2 = ll2;
719//      }
720//      ijp = commonIP;
721
722
723#if 0
724        for( i=0; i<lgth1; i++ )
725                fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
726#endif
727
728        currentw = w1;
729        previousw = w2;
730
731        match_calc( initverticalw, seq2, seq1, 0, lgth1 );
732
733        match_calc( currentw, seq1, seq2, 0, lgth2 );
734
735
736        lasti = lgth2+1;
737        for( j=1; j<lasti; ++j ) 
738        {
739                m[j] = currentw[j-1]; 
740//              mp[j] = 0;
741#if 0
742                if( m[j] < localthr ) m[j] = localthr2;
743#endif
744        }
745
746        lastverticalw[0] = currentw[lgth2-1];
747
748        lasti = lgth1+1;
749
750#if 0
751fprintf( stderr, "currentw = \n" );
752for( i=0; i<lgth1+1; i++ )
753{
754        fprintf( stderr, "%5.2f ", currentw[i] );
755}
756fprintf( stderr, "\n" );
757fprintf( stderr, "initverticalw = \n" );
758for( i=0; i<lgth2+1; i++ )
759{
760        fprintf( stderr, "%5.2f ", initverticalw[i] );
761}
762fprintf( stderr, "\n" );
763#endif
764#if DEBUG2
765        fprintf( stderr, "\n" );
766        fprintf( stderr, "       " );
767        for( j=0; j<lgth2; j++ )
768                fprintf( stderr, "%c     ", seq2[0][j] );
769        fprintf( stderr, "\n" );
770#endif
771
772        localstop = lgth1+lgth2+1;
773        maxwm = -999999999.9;
774#if DEBUG2
775        fprintf( stderr, "\n" );
776        fprintf( stderr, "%c   ", seq1[0][0] );
777
778        for( j=0; j<lgth2+1; j++ )
779                fprintf( stderr, "%5.0f ", currentw[j] );
780        fprintf( stderr, "\n" );
781#endif
782
783        for( i=1; i<lasti; i++ )
784        {
785                wtmp = previousw; 
786                previousw = currentw;
787                currentw = wtmp;
788
789                previousw[0] = initverticalw[i-1];
790
791                match_calc( currentw, seq1, seq2, i, lgth2 );
792#if DEBUG2
793                fprintf( stderr, "%c   ", seq1[0][i] );
794                fprintf( stderr, "%5.0f ", currentw[0] );
795#endif
796
797#if XXXXXXX
798fprintf( stderr, "\n" );
799fprintf( stderr, "i=%d\n", i );
800fprintf( stderr, "currentw = \n" );
801for( j=0; j<lgth2; j++ )
802{
803        fprintf( stderr, "%5.2f ", currentw[j] );
804}
805fprintf( stderr, "\n" );
806#endif
807#if XXXXXXX
808fprintf( stderr, "\n" );
809fprintf( stderr, "i=%d\n", i );
810fprintf( stderr, "currentw = \n" );
811for( j=0; j<lgth2; j++ )
812{
813        fprintf( stderr, "%5.2f ", currentw[j] );
814}
815fprintf( stderr, "\n" );
816#endif
817                currentw[0] = initverticalw[i];
818
819                mi = previousw[0]; 
820//              mpi = 0;
821
822#if 0
823                if( mi < localthr ) mi = localthr2;
824#endif
825
826//              ijppt = ijp[i] + 1;
827                mjpt = m + 1;
828                prept = previousw;
829                curpt = currentw + 1;
830//              mpjpt = mp + 1;
831                lastj = lgth2+1;
832                for( j=1; j<lastj; j++ )
833                {
834                        wm = *prept;
835//                      *ijppt = 0;
836
837#if 0
838                        fprintf( stderr, "%5.0f->", wm );
839#endif
840#if 0
841                        fprintf( stderr, "%5.0f?", g );
842#endif
843                        if( (g=mi+fpenalty) > wm )
844                        {
845                                wm = g;
846//                              *ijppt = -( j - mpi );
847                        }
848                        if( *prept > mi )
849                        {
850                                mi = *prept;
851//                              mpi = j-1;
852                        }
853
854#if USE_PENALTY_EX
855                        mi += fpenalty_ex;
856#endif
857
858#if 0 
859                        fprintf( stderr, "%5.0f?", g );
860#endif
861                        if( (g=*mjpt+fpenalty) > wm )
862                        {
863                                wm = g;
864//                              *ijppt = +( i - *mpjpt );
865                        }
866                        if( *prept > *mjpt )
867                        {
868                                *mjpt = *prept;
869//                              *mpjpt = i-1;
870                        }
871#if USE_PENALTY_EX
872                        *mjpt += fpenalty_ex;
873#endif
874
875                        if( maxwm < wm )
876                        {
877                                maxwm = wm;
878//                              endali = i;
879//                              endalj = j;
880                        }
881#if 1
882                        if( wm < localthr )
883                        {
884//                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
885//                              *ijppt = localstop;
886                                wm = localthr2;
887                        }
888#endif
889#if 0
890                        fprintf( stderr, "%5.0f ", *curpt );
891#endif
892#if DEBUG2
893                        fprintf( stderr, "%5.0f ", wm );
894//                      fprintf( stderr, "%c-%c *ijppt = %d, localstop = %d\n", seq1[0][i], seq2[0][j], *ijppt, localstop );
895#endif
896
897                        *curpt++ += wm;
898//                      ijppt++;
899                        mjpt++;
900                        prept++;
901//                      mpjpt++;
902                }
903#if DEBUG2
904                fprintf( stderr, "\n" );
905#endif
906
907                lastverticalw[i] = currentw[lgth2-1];
908        }
909
910
911#if 0
912        fprintf( stderr, "maxwm = %f\n", maxwm );
913        fprintf( stderr, "endali = %d\n", endali );
914        fprintf( stderr, "endalj = %d\n", endalj );
915#endif
916
917
918#if 0 // IRUKAMO!!!!
919        if( ijp[endali][endalj] == localstop )
920        {
921                strcpy( seq1[0], "" );
922                strcpy( seq2[0], "" );
923                *off1pt = *off2pt = 0;
924                fprintf( stderr, "maxwm <- 0.0 \n" );
925                return( 0.0 );
926        }
927#else
928        if( maxwm < localthr )
929        {
930                fprintf( stderr, "maxwm <- 0.0 \n" );
931                return( 0.0 );
932        }
933#endif
934               
935//      Ltracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, off1pt, off2pt, endali, endalj );
936
937
938//      resultlen = strlen( mseq1[0] );
939//      if( alloclen < resultlen || resultlen > N )
940//      {
941//              fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
942//              ErrorExit( "LENGTH OVER!\n" );
943//      }
944
945
946//      strcpy( seq1[0], mseq1[0] );
947//      strcpy( seq2[0], mseq2[0] );
948
949#if 0
950        fprintf( stderr, "wm=%f\n", wm );
951        fprintf( stderr, ">\n%s\n", mseq1[0] );
952        fprintf( stderr, ">\n%s\n", mseq2[0] );
953
954        fprintf( stderr, "maxwm = %f\n", maxwm );
955        fprintf( stderr, "   wm = %f\n",    wm );
956#endif
957
958        return( maxwm );
959}
Note: See TracBrowser for help on using the repository browser.