source: branches/items/GDE/MAFFT/mafft-7.055-with-extensions/core/Lalignmm.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: 53.6 KB
Line 
1#include "mltaln.h"
2#include "dp.h"
3
4#define MEMSAVE 1
5
6#define DEBUG 0
7#define USE_PENALTY_EX  0
8#define STOREWM 1
9
10#define DPTANNI 10
11
12#define LOCAL 0
13
14static int reccycle = 0;
15
16static float localthr;
17
18static void match_ribosum( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
19{
20        int j, k, l;
21        float scarr[38];
22        float **cpmxpd = floatwork;
23        int **cpmxpdn = intwork;
24        int count = 0;
25        float *matchpt;
26        float **cpmxpdpt;
27        int **cpmxpdnpt;
28        int cpkd;
29
30        if( initialize )
31        {
32                for( j=0; j<lgth2; j++ )
33                {
34                        count = 0;
35                        for( l=0; l<37; l++ )
36                        {
37                                if( cpmx2[j][l] )
38                                {
39                                        cpmxpd[j][count] = cpmx2[j][l];
40                                        cpmxpdn[j][count] = l;
41                                        count++;
42                                }
43                        }
44                        cpmxpdn[j][count] = -1;
45                }
46        }
47
48        for( l=0; l<37; l++ )
49        {
50                scarr[l] = 0.0;
51                for( k=0; k<37; k++ )
52                {
53                        scarr[l] += ribosumdis[k][l] * cpmx1[i1][k];
54                }
55        }
56#if 0 /* €³€ì€ò»È€Š€È€­€Ïfloatwork€Î¥¢¥í¥±¡Œ¥È€òµÕ€Ë€¹€ë */
57        {
58                float *fpt, **fptpt, *fpt2;
59                int *ipt, **iptpt;
60                fpt2 = match;
61                iptpt = cpmxpdn;
62                fptpt = cpmxpd;
63                while( lgth2-- )
64                {
65                        *fpt2 = 0.0;
66                        ipt=*iptpt,fpt=*fptpt;
67                        while( *ipt > -1 )
68                                *fpt2 += scarr[*ipt++] * *fpt++;
69                        fpt2++,iptpt++,fptpt++;
70                }
71        }
72        for( j=0; j<lgth2; j++ )
73        {
74                match[j] = 0.0;
75                for( k=0; cpmxpdn[j][k]>-1; k++ )
76                        match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
77        }
78#else
79        matchpt = match;
80        cpmxpdnpt = cpmxpdn;
81        cpmxpdpt = cpmxpd;
82        while( lgth2-- )
83        {
84                *matchpt = 0.0;
85                for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
86                        *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
87                matchpt++;
88                cpmxpdnpt++;
89                cpmxpdpt++;
90        }
91#endif
92}
93
94static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
95{
96        int j, k, l;
97        float scarr[26];
98        float **cpmxpd = floatwork;
99        int **cpmxpdn = intwork;
100        int count = 0;
101        float *matchpt;
102        float **cpmxpdpt;
103        int **cpmxpdnpt;
104        int cpkd;
105
106        if( initialize )
107        {
108                for( j=0; j<lgth2; j++ )
109                {
110                        count = 0;
111                        for( l=0; l<26; l++ )
112                        {
113                                if( cpmx2[j][l] )
114                                {
115                                        cpmxpd[j][count] = cpmx2[j][l];
116                                        cpmxpdn[j][count] = l;
117                                        count++;
118                                }
119                        }
120                        cpmxpdn[j][count] = -1;
121                }
122        }
123
124        for( l=0; l<26; l++ )
125        {
126                scarr[l] = 0.0;
127                for( k=0; k<26; k++ )
128                {
129                        scarr[l] += (n_dis[k][l]-RNAthr) * cpmx1[i1][k];
130                }
131        }
132#if 0 /* €³€ì€ò»È€Š€È€­€Ïfloatwork€Î¥¢¥í¥±¡Œ¥È€òµÕ€Ë€¹€ë */
133        {
134                float *fpt, **fptpt, *fpt2;
135                int *ipt, **iptpt;
136                fpt2 = match;
137                iptpt = cpmxpdn;
138                fptpt = cpmxpd;
139                while( lgth2-- )
140                {
141                        *fpt2 = 0.0;
142                        ipt=*iptpt,fpt=*fptpt;
143                        while( *ipt > -1 )
144                                *fpt2 += scarr[*ipt++] * *fpt++;
145                        fpt2++,iptpt++,fptpt++;
146                }
147        }
148        for( j=0; j<lgth2; j++ )
149        {
150                match[j] = 0.0;
151                for( k=0; cpmxpdn[j][k]>-1; k++ )
152                        match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
153        }
154#else
155        matchpt = match;
156        cpmxpdnpt = cpmxpdn;
157        cpmxpdpt = cpmxpd;
158        while( lgth2-- )
159        {
160                *matchpt = 0.0;
161                for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
162                        *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
163                matchpt++;
164                cpmxpdnpt++;
165                cpmxpdpt++;
166        }
167#endif
168}
169
170#if 0
171static void match_add( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
172{
173        int j, k, l;
174        float scarr[26];
175        float **cpmxpd = floatwork;
176        int **cpmxpdn = intwork;
177        int count = 0;
178        float *matchpt;
179        float **cpmxpdpt;
180        int **cpmxpdnpt;
181        int cpkd;
182
183
184        if( initialize )
185        {
186                for( j=0; j<lgth2; j++ )
187                {
188                        count = 0;
189                        for( l=0; l<26; l++ )
190                        {
191                                if( cpmx2[j][l] )
192                                {
193                                        cpmxpd[j][count] = cpmx2[j][l];
194                                        cpmxpdn[j][count] = l;
195                                        count++;
196                                }
197                        }
198                        cpmxpdn[j][count] = -1;
199                }
200        }
201
202        for( l=0; l<26; l++ )
203        {
204                scarr[l] = 0.0;
205                for( k=0; k<26; k++ )
206                {
207                        scarr[l] += n_dis[k][l] * cpmx1[i1][k];
208                }
209        }
210#if 0 /* €³€ì€ò»È€Š€È€­€Ïfloatwork€Î¥¢¥í¥±¡Œ¥È€òµÕ€Ë€¹€ë */
211        {
212                float *fpt, **fptpt, *fpt2;
213                int *ipt, **iptpt;
214                fpt2 = match;
215                iptpt = cpmxpdn;
216                fptpt = cpmxpd;
217                while( lgth2-- )
218                {
219                        *fpt2 = 0.0;
220                        ipt=*iptpt,fpt=*fptpt;
221                        while( *ipt > -1 )
222                                *fpt2 += scarr[*ipt++] * *fpt++;
223                        fpt2++,iptpt++,fptpt++;
224                }
225        }
226        for( j=0; j<lgth2; j++ )
227        {
228                match[j] = 0.0;
229                for( k=0; cpmxpdn[j][k]>-1; k++ )
230                        match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
231        }
232#else
233        matchpt = match;
234        cpmxpdnpt = cpmxpdn;
235        cpmxpdpt = cpmxpd;
236        while( lgth2-- )
237        {
238//              *matchpt = 0.0; // add dakara
239                for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
240                        *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
241                matchpt++;
242                cpmxpdnpt++;
243                cpmxpdpt++;
244        }
245#endif
246}
247#endif
248
249#if 0
250static float Atracking(
251                                                char **seq1, char **seq2,
252                        char **mseq1, char **mseq2,
253                        int **ijp, int icyc, int jcyc,
254                                                int ist, int ien, int jst, int jen )
255{
256        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, klim;
257        char *gaptable1, *gt1bk;
258        char *gaptable2, *gt2bk;
259        lgth1 = ien-ist+1;
260        lgth2 = jen-jst+1;
261
262        gt1bk = AllocateCharVec( lgth1+lgth2+1 );
263        gt2bk = AllocateCharVec( lgth1+lgth2+1 );
264
265#if 0
266        for( i=0; i<lgth1; i++ )
267        {
268                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
269        }
270#endif
271
272
273//      fprintf( stderr, "in Atracking, lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
274 
275    for( i=0; i<lgth1+1; i++ )
276    {
277        ijp[i][0] = i + 1;
278    }
279    for( j=0; j<lgth2+1; j++ )
280    {
281        ijp[0][j] = -( j + 1 );
282    }
283
284
285        gaptable1 = gt1bk + lgth1+lgth2;
286        *gaptable1 = 0;
287        gaptable2 = gt2bk + lgth1+lgth2;
288        *gaptable2 = 0;
289
290//      if( lgth2 == 1 ) fprintf( stderr, "in Atracking, mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
291
292        iin = lgth1; jin = lgth2;
293        klim = lgth1+lgth2;
294        for( k=0; k<=klim; k++ )
295        {
296                if( ijp[iin][jin] < 0 )
297                {
298                        ifi = iin-1; jfi = jin+ijp[iin][jin];
299                }
300                else if( ijp[iin][jin] > 0 )
301                {
302                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
303                }
304                else
305                {
306                        ifi = iin-1; jfi = jin-1;
307                }
308                l = iin - ifi;
309                while( --l )
310                {
311                        *--gaptable1 = 'o';
312                        *--gaptable2 = '-';
313                        k++;
314                }
315                l= jin - jfi;
316                while( --l )
317                {
318                        *--gaptable1 = '-';
319                        *--gaptable2 = 'o';
320                        k++;
321                }
322                if( iin <= 0 || jin <= 0 ) break;
323                *--gaptable1 = 'o';
324                *--gaptable2 = 'o';
325                k++;
326                iin = ifi; jin = jfi;
327
328        }
329
330        for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i]+ist, gaptable1 );
331        for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j]+jst, gaptable2 );
332
333        free( gt1bk );
334        free( gt2bk );
335
336//      fprintf( stderr, "in Atracking (owari), mseq1 = %s\n", mseq1[0] );
337//      fprintf( stderr, "in Atracking (owari), mseq2 = %s\n", mseq2[0] );
338        return( 0.0 );
339}
340#endif
341
342
343static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, float **cpmx1, float **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, int depth, float **gapinfo, float **map )
344/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
345{
346        float value = 0.0;
347        register int i, j;
348        char **aseq1, **aseq2;
349        int ll1, ll2;
350        int lasti, lastj, imid, jmid = 0;
351        float wm = 0.0;   /* int ?????? */
352        float g;
353        float *currentw, *previousw;
354#if USE_PENALTY_EX
355        float fpenalty_ex = (float)penalty_ex;
356#endif
357//      float fpenalty = (float)penalty;
358        float *wtmp;
359//      short *ijppt;
360        int *mpjpt;
361//      short **ijp;
362        int *mp;
363        int mpi;
364        float *mjpt, *prept, *curpt;
365        float mi;
366        float *m;
367        float *w1, *w2;
368//      float *match;
369        float *initverticalw;    /* kufuu sureba iranai */
370        float *lastverticalw;    /* kufuu sureba iranai */
371        int **intwork;
372        float **floatwork;
373//      short **shortmtx;
374#if STOREWM
375        float **WMMTX;
376        float **WMMTX2;
377#endif
378        float *midw;
379        float *midm;
380        float *midn;
381        int lgth1, lgth2;
382        float maxwm = 0.0;
383        int *jumpforwi;
384        int *jumpforwj;
385        int *jumpbacki;
386        int *jumpbackj;
387        int *jumpdummi; //muda
388        int *jumpdummj; //muda
389        int jumpi, jumpj = 0;
390        char *gaps;
391        int ijpi, ijpj;
392        float *ogcp1;
393        float *fgcp1;
394        float *ogcp2;
395        float *fgcp2;
396        float firstm;
397        int firstmp;
398#if 0
399        static char ttt1[50000];
400        static char ttt2[50000];
401#endif
402
403        localthr = -offset + 500; // 0?
404
405        ogcp1 = gapinfo[0] + ist;
406        fgcp1 = gapinfo[1] + ist;
407        ogcp2 = gapinfo[2] + jst;
408        fgcp2 = gapinfo[3] + jst;
409
410        depth++;
411        reccycle++;
412
413        lgth1 = ien-ist+1;
414        lgth2 = jen-jst+1;
415
416//      if( lgth1 < 5 )
417//              fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
418//      if( lgth2 < 5 )
419//              fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
420//
421
422#if 0
423        fprintf( stderr, "==== MSalign (depth=%d, reccycle=%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, reccycle, ist, ien, jst, jen );
424        strncpy( ttt1, seq1[0]+ist, lgth1 );
425        strncpy( ttt2, seq2[0]+jst, lgth2 );
426        ttt1[lgth1] = 0;
427        ttt2[lgth2] = 0;
428        fprintf( stderr, "seq1 = %s\n", ttt1 );
429        fprintf( stderr, "seq2 = %s\n", ttt2 );
430#endif
431        if( lgth2 <= 0 ) // lgth1 <= 0 ha?
432        {
433//              fprintf( stderr, "\n\n==== jimei\n\n" );
434//              exit( 1 );
435                for( i=0; i<icyc; i++ ) 
436                {
437                        strncpy( mseq1[i], seq1[i]+ist, lgth1 );
438                        mseq1[i][lgth1] = 0;
439                }
440                for( i=0; i<jcyc; i++ ) 
441                {
442                        mseq2[i][0] = 0;
443                        for( j=0; j<lgth1; j++ )
444                                strcat( mseq2[i], "-" );
445                }
446
447//              fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
448//              fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
449
450                return( 0.0 );
451        }
452
453#if MEMSAVE
454        aseq1 = AllocateCharMtx( icyc, 0 );
455        aseq2 = AllocateCharMtx( jcyc, 0 );
456        for( i=0; i<icyc; i++ ) aseq1[i] = mseq1[i];
457        for( i=0; i<jcyc; i++ ) aseq2[i] = mseq2[i];
458#else
459        aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
460        aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
461#endif
462
463//  if( lgth1 < DPTANNI && lgth2 < DPTANNI ) // & dato lgth ==1 no kanousei ga arunode yokunai
464//    if( lgth1 < DPTANNI ) // kore mo lgth2 ga mijikasugiru kanousei ari
465    if( lgth1 < DPTANNI || lgth2 < DPTANNI ) // zettai ni anzen ka?
466        {
467//              fprintf( stderr, "==== Going to _tanni\n" );
468
469//              value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2, gapinfo );       
470
471
472#if MEMSAVE
473                free( aseq1 );
474                free( aseq2 );
475#else
476                for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
477                for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
478
479                FreeCharMtx( aseq1 );
480                FreeCharMtx( aseq2 );
481#endif
482
483//              fprintf( stderr, "value = %f\n", value );
484
485                return( value );
486        }
487//      fprintf( stderr, "Trying to divide the mtx\n" );
488
489        ll1 = ( (int)(lgth1) ) + 100;
490        ll2 = ( (int)(lgth2) ) + 100;
491
492//      fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
493
494        w1 = AllocateFloatVec( ll2+2 );
495        w2 = AllocateFloatVec( ll2+2 );
496//      match = AllocateFloatVec( ll2+2 );
497        midw = AllocateFloatVec( ll2+2 );
498        midn = AllocateFloatVec( ll2+2 );
499        midm = AllocateFloatVec( ll2+2 );
500        jumpbacki = AllocateIntVec( ll2+2 );
501        jumpbackj = AllocateIntVec( ll2+2 );
502        jumpforwi = AllocateIntVec( ll2+2 );
503        jumpforwj = AllocateIntVec( ll2+2 );
504        jumpdummi = AllocateIntVec( ll2+2 );
505        jumpdummj = AllocateIntVec( ll2+2 );
506
507        initverticalw = AllocateFloatVec( ll1+2 );
508        lastverticalw = AllocateFloatVec( ll1+2 );
509
510        m = AllocateFloatVec( ll2+2 );
511        mp = AllocateIntVec( ll2+2 );
512        gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
513
514        floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 ); 
515        intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 ); 
516
517#if DEBUG
518        fprintf( stderr, "succeeded\n" );
519#endif
520
521#if STOREWM
522        WMMTX = AllocateFloatMtx( ll1, ll2 );
523        WMMTX2 = AllocateFloatMtx( ll1, ll2 );
524#endif
525#if 0
526        shortmtx = AllocateShortMtx( ll1, ll2 );
527
528#if DEBUG
529        fprintf( stderr, "succeeded\n\n" );
530#endif
531
532        ijp = shortmtx;
533#endif
534
535        currentw = w1;
536        previousw = w2;
537
538        match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
539
540        match_ribosum( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
541
542        for( i=1; i<lgth1+1; i++ )
543        {
544                initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] );
545        }
546        for( j=1; j<lgth2+1; j++ )
547        {
548                currentw[j] += ( ogcp2[0] + fgcp2[j-1] );
549        }
550
551#if STOREWM
552        WMMTX[0][0] = initverticalw[0];
553        for( i=1; i<lgth1+1; i++ )
554        {
555                WMMTX[i][0] = initverticalw[i];
556        }
557        for( j=1; j<lgth2+1; j++ )
558        {
559                WMMTX[0][j] = currentw[j];
560        }
561#endif
562
563
564        for( j=1; j<lgth2+1; ++j ) 
565        {
566                m[j] = currentw[j-1] + ogcp1[1];
567//              m[j] = currentw[j-1];
568                mp[j] = 0;
569        }
570
571        lastverticalw[0] = currentw[lgth2-1];
572
573        imid = lgth1 * 0.5;
574
575        jumpi = 0; // atode kawaru.
576        lasti = lgth1+1;
577#if STOREWM
578        for( i=1; i<lasti; i++ )
579#else
580        for( i=1; i<=imid; i++ )
581#endif
582        {
583                wtmp = previousw; 
584                previousw = currentw;
585                currentw = wtmp;
586
587                previousw[0] = initverticalw[i-1];
588
589                match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
590                currentw[0] = initverticalw[i];
591
592                m[0] = ogcp1[i];
593#if STOREM
594                WMMTX2[i][0] = m[0];
595#endif
596                if( i == imid ) midm[0] = m[0];
597
598                mi = previousw[0] + ogcp2[1]; 
599//              mi = previousw[0];
600                mpi = 0;
601
602
603//              ijppt = ijp[i] + 1;
604                mjpt = m + 1;
605                prept = previousw;
606                curpt = currentw + 1;
607                mpjpt = mp + 1;
608
609
610                lastj = lgth2+1;
611                for( j=1; j<lastj; j++ )
612                {
613
614                        wm = *prept;
615
616#if 0
617                        fprintf( stderr, "%5.0f->", wm );
618#endif
619                        g = mi + fgcp2[j-1];
620//                      g = mi + fpenalty;
621#if 0
622                        fprintf( stderr, "%5.0f?", g );
623#endif
624                        if( g > wm )
625                        {
626                                wm = g;
627//                              *ijppt = -( j - mpi );
628                        }
629                        g = *prept + ogcp2[j];
630//                      g = *prept;
631                        if( g >= mi )
632                        {
633                                mi = g;
634                                mpi = j-1;
635                        }
636#if USE_PENALTY_EX
637                        mi += fpenalty_ex;
638#endif
639
640                        g = *mjpt + fgcp1[i-1];
641//                      g = *mjpt + fpenalty;
642#if 0 
643                        fprintf( stderr, "%5.0f?", g );
644#endif
645                        if( g > wm )
646                        {
647                                wm = g;
648//                              *ijppt = +( i - *mpjpt );
649                        }
650
651
652                        g = *prept + ogcp1[i];
653//                      g = *prept;
654                        if( g >= *mjpt )
655                        {
656                                *mjpt = g;
657                                *mpjpt = i-1;
658                        }
659#if USE_PENALTY_EX
660                        m[j] += fpenalty_ex;
661#endif
662#if LOCAL
663            if( wm < localthr )
664            {       
665//                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
666                                wm = 0;
667            }       
668#endif
669
670#if 0
671                        fprintf( stderr, "%5.0f ", wm );
672#endif
673                        *curpt += wm;
674
675
676#if STOREWM
677                        WMMTX[i][j] = *curpt;
678                        WMMTX2[i][j] = *mjpt;
679#endif
680
681                        if( i == imid ) //muda
682                        {       
683                                jumpbackj[j] = *mpjpt; // muda atode matomeru
684                                jumpbacki[j] = mpi; // muda atode matomeru
685//                              fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
686//                              fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
687                                midw[j] = *curpt;
688                                midm[j] = *mjpt;
689                                midn[j] = mi;
690                        }
691
692//                      fprintf( stderr, "m[%d] = %f\n", j, m[j] );
693                        mjpt++;
694                        prept++;
695                        mpjpt++;
696                        curpt++;
697
698                }
699                lastverticalw[i] = currentw[lgth2-1];
700
701#if STOREWM
702                WMMTX2[i][lgth2] = m[lgth2-1];
703#endif
704
705#if 0  // ue
706                if( i == imid )
707                {
708                        for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
709                        for( j=0; j<lgth2; j++ ) midm[j] = m[j];
710                }
711#endif
712        }
713//      for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
714//      for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
715
716#if 0
717    for( i=0; i<lgth1; i++ )
718    {
719        for( j=0; j<lgth2; j++ )
720        {
721            fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
722        }
723        fprintf( stderr, "\n" );
724    }
725        fprintf( stderr, "\n" );
726        fprintf( stderr, "WMMTX2 = \n" );
727    for( i=0; i<lgth1; i++ )
728    {
729        for( j=0; j<lgth2; j++ )
730        {
731            fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
732        }
733        fprintf( stderr, "\n" );
734    }
735        fprintf( stderr, "\n" );
736#endif
737
738// gyakudp
739
740        match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
741        match_ribosum( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
742
743        for( i=0; i<lgth1-1; i++ )
744        {
745                initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
746//              initverticalw[i] += fpenalty;
747        }
748        for( j=0; j<lgth2-1; j++ )
749        {
750                currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
751//              currentw[j] += fpenalty;
752        }
753
754#if STOREWM
755        for( i=0; i<lgth1-1; i++ )
756        {
757                WMMTX[i][lgth2-1] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
758//              fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
759        }
760        for( j=0; j<lgth2-1; j++ )
761        {
762                WMMTX[lgth1-1][j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
763//              fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
764        }
765#endif
766
767
768
769
770
771
772        for( j=lgth2-1; j>0; --j )
773        {
774                m[j-1] = currentw[j] + fgcp2[lgth2-2];
775//              m[j-1] = currentw[j];
776                mp[j] = lgth1-1;
777        }
778
779//      for( j=0; j<lgth2; j++ ) m[j] = 0.0;
780        // m[lgth2-1] ha irunoka?
781
782
783//      for( i=lgth1-2; i>=imid; i-- )
784        firstm = -9999999.9;
785        firstmp = lgth1-1;
786        for( i=lgth1-2; i>-1; i-- )
787        {
788                wtmp = previousw;
789                previousw = currentw;
790                currentw = wtmp;
791                previousw[lgth2-1] = initverticalw[i+1];
792//              match_calc( currentw, seq1, seq2, i, lgth2 );
793                match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
794
795                currentw[lgth2-1] = initverticalw[i];
796
797//              m[lgth2] = fgcp1[i];
798//              WMMTX2[i][lgth2] += m[lgth2];
799//              fprintf( stderr, "m[] = %f\n", m[lgth2] );
800
801                mi = previousw[lgth2-1] + fgcp2[lgth2-2];
802//              mi = previousw[lgth2-1];
803                mpi = lgth2 - 1;
804
805                mjpt = m + lgth2 - 2;
806                prept = previousw + lgth2 - 1;
807                curpt = currentw + lgth2 - 2;
808                mpjpt = mp + lgth2 - 2;
809
810
811                for( j=lgth2-2; j>-1; j-- )
812                {
813                        wm = *prept;
814                        ijpi = i+1;
815                        ijpj = j+1;
816
817                        g = mi + ogcp2[j+1];
818//                      g = mi + fpenalty;
819                        if( g > wm )
820                        {
821                                wm = g;
822                                ijpj = mpi;
823                                ijpi = i+1;
824                        }
825
826                        g = *prept + fgcp2[j];
827//                      g = *prept;
828                        if( g >= mi )
829                        {
830//                              fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
831                                mi = g;
832                                mpi = j + 1;
833                        }
834
835#if USE_PENALTY_EX
836                        mi += fpenalty_ex;
837#endif
838
839//                      fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
840                        g = *mjpt + ogcp1[i+1];
841//                      g = *mjpt + fpenalty;
842                        if( g > wm )
843                        {
844                                wm = g;
845                                ijpi = *mpjpt;
846                                ijpj = j+1;
847                        }
848
849//                      if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
850                        g = *prept + fgcp1[i];
851//                      g = *prept;
852                        if( g >= *mjpt )
853                        {
854                                *mjpt = g;
855                                *mpjpt = i + 1;
856                        }
857
858#if USE_PENALTY_EX
859                        m[j] += fpenalty_ex;
860#endif
861
862                        if( i == jumpi || i == imid - 1 )
863                        {
864                                jumpforwi[j] = ijpi; //muda
865                                jumpforwj[j] = ijpj; //muda
866//                              fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
867//                              fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
868                        }
869                        if( i == imid ) // muda
870                        {
871                                midw[j] += wm;
872//                              midm[j+1] += *mjpt + fpenalty; //??????
873                                midm[j+1] += *mjpt; //??????
874                        }
875                        if( i == imid - 1 )
876                        {
877//                              midn[j] += mi + fpenalty;  //????
878                                midn[j] += mi;  //????
879                        }
880#if LOCAL
881            if( wm < localthr )
882            {       
883//                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
884                                wm = 0;
885            }       
886#endif
887
888#if STOREWM
889                        WMMTX[i][j] += wm;
890//                      WMMTX2[i][j+1] += *mjpt + fpenalty;
891                        WMMTX2[i][j] += *curpt;
892#endif
893                        *curpt += wm;
894
895                        mjpt--;
896                        prept--;
897                        mpjpt--;
898                        curpt--;
899                }
900//              fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
901                g = *prept + fgcp1[i];
902                if( firstm < g ) 
903                {
904                        firstm = g;
905                        firstmp = i + 1;
906                }
907#if STOREWM
908//              WMMTX2[i][j+1] += firstm;
909#endif
910                if( i == imid ) midm[j+1] += firstm;
911
912                if( i == imid - 1 )     
913                {
914                        maxwm = midw[1];
915                        jmid = 0;
916//                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
917                        for( j=2; j<lgth2-1; j++ )
918                        {
919                                wm = midw[j];
920                                if( wm > maxwm )
921                                {
922                                        jmid = j;
923                                        maxwm = wm;
924                                }
925//                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
926                        }
927                        for( j=0; j<lgth2+1; j++ )
928                        {
929                                wm = midm[j];
930                                if( wm > maxwm )
931                                {
932                                        jmid = j;
933                                        maxwm = wm;
934                                }
935//                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
936                        }
937
938//                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
939
940
941//                      fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
942                        wm = midw[jmid];
943                        jumpi = imid-1;
944                        jumpj = jmid-1;
945                        if( jmid > 0 && midn[jmid-1] > wm ) //060413
946                        {
947                                jumpi = imid-1;
948                                jumpj = jumpbacki[jmid];
949                                wm = midn[jmid-1];
950//                              fprintf( stderr, "rejump (n)\n" );
951                        }
952                        if( midm[jmid] > wm )
953                        {
954                                jumpi = jumpbackj[jmid];
955                                jumpj = jmid-1;
956                                wm = midm[jmid];
957//                              fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
958                        }
959
960
961//                      fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
962//                      fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
963#if 0
964                        fprintf( stderr, "imid = %d\n", imid );
965                        fprintf( stderr, "midn = \n" );
966                        for( j=0; j<lgth2; j++ )
967                        {
968                                fprintf( stderr, "% 7.1f ", midn[j] );
969                        }
970                        fprintf( stderr, "\n" );
971                        fprintf( stderr, "midw = \n" );
972                        for( j=0; j<lgth2; j++ )
973                        {
974                                fprintf( stderr, "% 7.1f ", midw[j] );
975                        }
976                        fprintf( stderr, "\n" );
977                        fprintf( stderr, "midm = \n" );
978                        for( j=0; j<lgth2; j++ )
979                        {
980                                fprintf( stderr, "% 7.1f ", midm[j] );
981                        }
982                        fprintf( stderr, "\n" );
983#endif
984//                      fprintf( stderr, "maxwm = %f\n", maxwm );
985                }
986                if( i == jumpi ) //saki?
987                {
988//                      fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
989//                      fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
990                        if( jmid == 0 )
991                        {
992//                              fprintf( stderr, "CHUI2!\n" );
993                                jumpj = 0; jmid = 1;
994                                jumpi = firstmp - 1;
995                                imid = firstmp;
996                        }
997
998#if 0
999                        else if( jmid == lgth2 )
1000                        {
1001                                fprintf( stderr, "CHUI1!\n" );
1002                                jumpi=0; jumpj=0;
1003                                imid=jumpforwi[0]; jmid=lgth2-1;
1004                        }
1005#else // 060414
1006                        else if( jmid >= lgth2 ) 
1007                        {
1008//                              fprintf( stderr, "CHUI1!\n" );
1009                                jumpi=imid-1; jmid=lgth2;
1010                                jumpj = lgth2-1;
1011                        }
1012#endif
1013                        else
1014                        {
1015                                imid = jumpforwi[jumpj];
1016                                jmid = jumpforwj[jumpj];
1017                        }
1018#if 0
1019                        fprintf( stderr, "jumpi -> %d\n", jumpi );
1020                        fprintf( stderr, "jumpj -> %d\n", jumpj );
1021                        fprintf( stderr, "imid -> %d\n", imid );
1022                        fprintf( stderr, "jmid -> %d\n", jmid );
1023#endif
1024
1025#if STOREWM
1026//                      break;
1027#else
1028                        break;
1029#endif
1030                }
1031        }
1032#if 0
1033                jumpi=0; jumpj=0;
1034                imid=lgth1-1; jmid=lgth2-1;
1035        }
1036#endif
1037
1038//      fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1039//      fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1040
1041//      for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1042//      for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1043//      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1044//      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1045
1046
1047
1048        for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
1049                map[i][j] = WMMTX[i][j] / maxwm;
1050//              map[i][j] = WMMTX2[i][j] / maxwm;
1051
1052#if STOREWM
1053
1054#if 0
1055        for( i=0; i<lgth1; i++ )
1056        {
1057                float maxpairscore = -9999.9;
1058                float tmpscore;
1059
1060                for( j=0; j<lgth2; j++ )
1061                {
1062                        if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1063                        {
1064                                map12[i].pos = j;
1065                                map12[i].score = tmpscore;
1066                                maxpairscore = tmpscore;
1067                        }
1068                }
1069
1070                for( k=0; k<lgth1; k++ )
1071                {
1072                        if( i == k ) continue;
1073                        if( map12[i].score <= WMMTX[k][map12[i].pos] )
1074                                break;
1075                }
1076                if( k != lgth1 )
1077                {
1078                        map12[i].pos = -1;
1079                        map12[i].score = -1.0;
1080                }
1081                fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, map12[i].pos, map12[i].score, seq1[0][i], seq2[0][map12[i].pos] );
1082        }
1083        for( j=0; j<lgth2; j++ )
1084        {
1085                float maxpairscore = -9999.9;
1086                float tmpscore;
1087
1088                for( i=0; i<lgth1; i++ )
1089                {
1090                        if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1091                        {
1092                                map21[j].pos = i;
1093                                map21[j].score = tmpscore;
1094                                maxpairscore = tmpscore;
1095                        }
1096                }
1097
1098                for( k=0; k<lgth2; k++ )
1099                {
1100                        if( j == k ) continue;
1101                        if( map21[j].score <= WMMTX[map21[j].pos][k] )
1102                                break;
1103                }
1104                if( k != lgth2 )
1105                {
1106                        map21[j].pos = -1;
1107                        map21[j].score = -1.0;
1108                }
1109                fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", j, map21[j].pos, map21[j].score, seq2[0][j], seq1[0][map21[j].pos] );
1110        }
1111
1112        for( i=0; i<lgth1; i++ )
1113        {
1114                if( map12[i].pos != -1 ) if( map21[map12[i].pos].pos != i )
1115                        fprintf( stderr, "ERROR i=%d, but map12[i].pos=%d and map21[map12[i].pos]=%d\n", i, map12[i].pos, map21[map12[i].pos].pos );
1116        }
1117#endif
1118
1119#if 0
1120        fprintf( stderr, "WMMTX = \n" );
1121    for( i=0; i<lgth1; i++ )
1122    {
1123        fprintf( stderr, "%d ", i );
1124        for( j=0; j<lgth2; j++ )
1125        {
1126            fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
1127        }
1128        fprintf( stderr, "\n" );
1129    }
1130//      fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
1131    for( i=0; i<lgth1; i++ )
1132    {
1133        fprintf( stderr, "%d ", i );
1134        for( j=0; j<lgth2+1; j++ )
1135        {
1136            fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
1137        }
1138        fprintf( stderr, "\n" );
1139    }
1140#endif
1141
1142#if 0
1143    fprintf( stdout, "#WMMTX = \n" );
1144    for( i=0; i<lgth1; i++ )
1145    {
1146//        fprintf( stdout, "%d ", i );
1147        for( j=0; j<lgth2; j++ )
1148        {
1149//          if( WMMTX[i][j] > amino_dis['a']['g'] -1 )
1150                fprintf( stdout, "%d %d %8.1f", i, j, WMMTX[i][j] );
1151                        if( WMMTX[i][j] == maxwm )
1152                fprintf( stdout, "selected \n" );
1153                        else
1154                fprintf( stdout, "\n" );
1155        }
1156        fprintf( stdout, "\n" );
1157    }
1158#endif
1159
1160#if 0
1161
1162        fprintf( stderr, "jumpbacki = \n" );
1163        for( j=0; j<lgth2; j++ )
1164        {
1165                fprintf( stderr, "% 7d ", jumpbacki[j] );
1166        }
1167        fprintf( stderr, "\n" );
1168        fprintf( stderr, "jumpbackj = \n" );
1169        for( j=0; j<lgth2; j++ )
1170        {
1171                fprintf( stderr, "% 7d ", jumpbackj[j] );
1172        }
1173        fprintf( stderr, "\n" );
1174        fprintf( stderr, "jumpforwi = \n" );
1175        for( j=0; j<lgth2; j++ )
1176        {
1177                fprintf( stderr, "% 7d ", jumpforwi[j] );
1178        }
1179        fprintf( stderr, "\n" );
1180        fprintf( stderr, "jumpforwj = \n" );
1181        for( j=0; j<lgth2; j++ )
1182        {
1183                fprintf( stderr, "% 7d ", jumpforwj[j] );
1184        }
1185        fprintf( stderr, "\n" );
1186#endif
1187
1188
1189#endif
1190
1191
1192//      Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1193
1194#if 0 // irukamo
1195        resultlen = strlen( mseq1[0] );
1196        if( alloclen < resultlen || resultlen > N )
1197        {
1198                fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1199                ErrorExit( "LENGTH OVER!\n" );
1200        }
1201#endif
1202
1203
1204
1205#if 0
1206        fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
1207        fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
1208
1209        fprintf( stderr, "imid = %d\n", imid );
1210        fprintf( stderr, "jmid = %d\n", jmid );
1211#endif
1212
1213
1214        FreeFloatVec( w1 );
1215        FreeFloatVec( w2 );
1216        FreeFloatVec( initverticalw );
1217        FreeFloatVec( lastverticalw );
1218        FreeFloatVec( midw );
1219        FreeFloatVec( midm );
1220        FreeFloatVec( midn );
1221
1222        FreeIntVec( jumpbacki );
1223        FreeIntVec( jumpbackj );
1224        FreeIntVec( jumpforwi );
1225        FreeIntVec( jumpforwj );
1226        FreeIntVec( jumpdummi );
1227        FreeIntVec( jumpdummj );
1228
1229        FreeFloatVec( m );
1230        FreeIntVec( mp );
1231
1232        FreeFloatMtx( floatwork );
1233        FreeIntMtx( intwork );
1234
1235#if STOREWM
1236        FreeFloatMtx( WMMTX );
1237        FreeFloatMtx( WMMTX2 );
1238#endif
1239
1240        return( value );
1241
1242}
1243static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, float **cpmx1, float **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, int depth, float **gapinfo, float **map )
1244/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1245{
1246        float value = 0.0;
1247        register int i, j;
1248        char **aseq1, **aseq2;
1249        int ll1, ll2;
1250        int lasti, lastj, imid, jmid=0;
1251        float wm = 0.0;   /* int ?????? */
1252        float g;
1253        float *currentw, *previousw;
1254#if USE_PENALTY_EX
1255        float fpenalty_ex = (float)RNApenalty_ex;
1256#endif
1257//      float fpenalty = (float)penalty;
1258        float *wtmp;
1259//      short *ijppt;
1260        int *mpjpt;
1261//      short **ijp;
1262        int *mp;
1263        int mpi;
1264        float *mjpt, *prept, *curpt;
1265        float mi;
1266        float *m;
1267        float *w1, *w2;
1268//      float *match;
1269        float *initverticalw;    /* kufuu sureba iranai */
1270        float *lastverticalw;    /* kufuu sureba iranai */
1271        int **intwork;
1272        float **floatwork;
1273//      short **shortmtx;
1274#if STOREWM
1275        float **WMMTX;
1276        float **WMMTX2;
1277#endif
1278        float *midw;
1279        float *midm;
1280        float *midn;
1281        int lgth1, lgth2;
1282        float maxwm = 0.0;
1283        int *jumpforwi;
1284        int *jumpforwj;
1285        int *jumpbacki;
1286        int *jumpbackj;
1287        int *jumpdummi; //muda
1288        int *jumpdummj; //muda
1289        int jumpi, jumpj = 0;
1290        char *gaps;
1291        int ijpi, ijpj;
1292        float *ogcp1;
1293        float *fgcp1;
1294        float *ogcp2;
1295        float *fgcp2;
1296        float firstm;
1297        int firstmp;
1298#if 0
1299        static char ttt1[50000];
1300        static char ttt2[50000];
1301#endif
1302
1303        localthr = -offset + 500; // 0?
1304
1305        ogcp1 = gapinfo[0] + ist;
1306        fgcp1 = gapinfo[1] + ist;
1307        ogcp2 = gapinfo[2] + jst;
1308        fgcp2 = gapinfo[3] + jst;
1309
1310        depth++;
1311        reccycle++;
1312
1313        lgth1 = ien-ist+1;
1314        lgth2 = jen-jst+1;
1315
1316//      if( lgth1 < 5 )
1317//              fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
1318//      if( lgth2 < 5 )
1319//              fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
1320//
1321
1322#if 0
1323        fprintf( stderr, "==== MSalign (depth=%d, reccycle=%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, reccycle, ist, ien, jst, jen );
1324        strncpy( ttt1, seq1[0]+ist, lgth1 );
1325        strncpy( ttt2, seq2[0]+jst, lgth2 );
1326        ttt1[lgth1] = 0;
1327        ttt2[lgth2] = 0;
1328        fprintf( stderr, "seq1 = %s\n", ttt1 );
1329        fprintf( stderr, "seq2 = %s\n", ttt2 );
1330#endif
1331        if( lgth2 <= 0 ) // lgth1 <= 0 ha?
1332        {
1333//              fprintf( stderr, "\n\n==== jimei\n\n" );
1334//              exit( 1 );
1335                for( i=0; i<icyc; i++ ) 
1336                {
1337                        strncpy( mseq1[i], seq1[i]+ist, lgth1 );
1338                        mseq1[i][lgth1] = 0;
1339                }
1340                for( i=0; i<jcyc; i++ ) 
1341                {
1342                        mseq2[i][0] = 0;
1343                        for( j=0; j<lgth1; j++ )
1344                                strcat( mseq2[i], "-" );
1345                }
1346
1347//              fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
1348//              fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
1349
1350                return( 0.0 );
1351        }
1352
1353#if MEMSAVE
1354        aseq1 = AllocateCharMtx( icyc, 0 );
1355        aseq2 = AllocateCharMtx( jcyc, 0 );
1356        for( i=0; i<icyc; i++ ) aseq1[i] = mseq1[i];
1357        for( i=0; i<jcyc; i++ ) aseq2[i] = mseq2[i];
1358#else
1359        aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
1360        aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
1361#endif
1362
1363//  if( lgth1 < DPTANNI && lgth2 < DPTANNI ) // & dato lgth ==1 no kanousei ga arunode yokunai
1364//    if( lgth1 < DPTANNI ) // kore mo lgth2 ga mijikasugiru kanousei ari
1365    if( lgth1 < DPTANNI || lgth2 < DPTANNI ) // zettai ni anzen ka?
1366        {
1367//              fprintf( stderr, "==== Going to _tanni\n" );
1368
1369//              value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2, gapinfo );       
1370
1371
1372#if MEMSAVE
1373                free( aseq1 );
1374                free( aseq2 );
1375#else
1376                for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
1377                for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
1378
1379                FreeCharMtx( aseq1 );
1380                FreeCharMtx( aseq2 );
1381#endif
1382
1383//              fprintf( stderr, "value = %f\n", value );
1384
1385                return( value );
1386        }
1387//      fprintf( stderr, "Trying to divide the mtx\n" );
1388
1389        ll1 = ( (int)(lgth1) ) + 100;
1390        ll2 = ( (int)(lgth2) ) + 100;
1391
1392//      fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
1393
1394        w1 = AllocateFloatVec( ll2+2 );
1395        w2 = AllocateFloatVec( ll2+2 );
1396//      match = AllocateFloatVec( ll2+2 );
1397        midw = AllocateFloatVec( ll2+2 );
1398        midn = AllocateFloatVec( ll2+2 );
1399        midm = AllocateFloatVec( ll2+2 );
1400        jumpbacki = AllocateIntVec( ll2+2 );
1401        jumpbackj = AllocateIntVec( ll2+2 );
1402        jumpforwi = AllocateIntVec( ll2+2 );
1403        jumpforwj = AllocateIntVec( ll2+2 );
1404        jumpdummi = AllocateIntVec( ll2+2 );
1405        jumpdummj = AllocateIntVec( ll2+2 );
1406
1407        initverticalw = AllocateFloatVec( ll1+2 );
1408        lastverticalw = AllocateFloatVec( ll1+2 );
1409
1410        m = AllocateFloatVec( ll2+2 );
1411        mp = AllocateIntVec( ll2+2 );
1412        gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
1413
1414        floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 ); 
1415        intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 ); 
1416
1417#if DEBUG
1418        fprintf( stderr, "succeeded\n" );
1419#endif
1420
1421#if STOREWM
1422        WMMTX = AllocateFloatMtx( ll1, ll2 );
1423        WMMTX2 = AllocateFloatMtx( ll1, ll2 );
1424#endif
1425#if 0
1426        shortmtx = AllocateShortMtx( ll1, ll2 );
1427
1428#if DEBUG
1429        fprintf( stderr, "succeeded\n\n" );
1430#endif
1431
1432        ijp = shortmtx;
1433#endif
1434
1435        currentw = w1;
1436        previousw = w2;
1437
1438        match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
1439
1440        match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
1441
1442        for( i=1; i<lgth1+1; i++ )
1443        {
1444                initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] );
1445        }
1446        for( j=1; j<lgth2+1; j++ )
1447        {
1448                currentw[j] += ( ogcp2[0] + fgcp2[j-1] );
1449        }
1450
1451#if STOREWM
1452        WMMTX[0][0] = initverticalw[0];
1453        for( i=1; i<lgth1+1; i++ )
1454        {
1455                WMMTX[i][0] = initverticalw[i];
1456        }
1457        for( j=1; j<lgth2+1; j++ )
1458        {
1459                WMMTX[0][j] = currentw[j];
1460        }
1461#endif
1462
1463
1464        for( j=1; j<lgth2+1; ++j ) 
1465        {
1466                m[j] = currentw[j-1] + ogcp1[1];
1467//              m[j] = currentw[j-1];
1468                mp[j] = 0;
1469        }
1470
1471        lastverticalw[0] = currentw[lgth2-1];
1472
1473        imid = lgth1 * 0.5;
1474
1475        jumpi = 0; // atode kawaru.
1476        lasti = lgth1+1;
1477#if STOREWM
1478        for( i=1; i<lasti; i++ )
1479#else
1480        for( i=1; i<=imid; i++ )
1481#endif
1482        {
1483                wtmp = previousw; 
1484                previousw = currentw;
1485                currentw = wtmp;
1486
1487                previousw[0] = initverticalw[i-1];
1488
1489                match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
1490                currentw[0] = initverticalw[i];
1491
1492                m[0] = ogcp1[i];
1493#if STOREM
1494                WMMTX2[i][0] = m[0];
1495#endif
1496                if( i == imid ) midm[0] = m[0];
1497
1498                mi = previousw[0] + ogcp2[1]; 
1499//              mi = previousw[0];
1500                mpi = 0;
1501
1502
1503//              ijppt = ijp[i] + 1;
1504                mjpt = m + 1;
1505                prept = previousw;
1506                curpt = currentw + 1;
1507                mpjpt = mp + 1;
1508
1509
1510                lastj = lgth2+1;
1511                for( j=1; j<lastj; j++ )
1512                {
1513
1514                        wm = *prept;
1515
1516#if 0
1517                        fprintf( stderr, "%5.0f->", wm );
1518#endif
1519                        g = mi + fgcp2[j-1];
1520//                      g = mi + fpenalty;
1521#if 0
1522                        fprintf( stderr, "%5.0f?", g );
1523#endif
1524                        if( g > wm )
1525                        {
1526                                wm = g;
1527//                              *ijppt = -( j - mpi );
1528                        }
1529                        g = *prept + ogcp2[j];
1530//                      g = *prept;
1531                        if( g >= mi )
1532                        {
1533                                mi = g;
1534                                mpi = j-1;
1535                        }
1536#if USE_PENALTY_EX
1537                        mi += fpenalty_ex;
1538#endif
1539
1540                        g = *mjpt + fgcp1[i-1];
1541//                      g = *mjpt + fpenalty;
1542#if 0 
1543                        fprintf( stderr, "%5.0f?", g );
1544#endif
1545                        if( g > wm )
1546                        {
1547                                wm = g;
1548//                              *ijppt = +( i - *mpjpt );
1549                        }
1550
1551
1552                        g = *prept + ogcp1[i];
1553//                      g = *prept;
1554                        if( g >= *mjpt )
1555                        {
1556                                *mjpt = g;
1557                                *mpjpt = i-1;
1558                        }
1559#if USE_PENALTY_EX
1560                        m[j] += fpenalty_ex;
1561#endif
1562#if LOCAL
1563            if( wm < localthr )
1564            {       
1565//                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
1566                                wm = 0;
1567            }       
1568#endif
1569
1570#if 0
1571                        fprintf( stderr, "%5.0f ", wm );
1572#endif
1573                        *curpt += wm;
1574
1575
1576#if STOREWM
1577                        WMMTX[i][j] = *curpt;
1578                        WMMTX2[i][j] = *mjpt;
1579#endif
1580
1581                        if( i == imid ) //muda
1582                        {       
1583                                jumpbackj[j] = *mpjpt; // muda atode matomeru
1584                                jumpbacki[j] = mpi; // muda atode matomeru
1585//                              fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
1586//                              fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
1587                                midw[j] = *curpt;
1588                                midm[j] = *mjpt;
1589                                midn[j] = mi;
1590                        }
1591
1592//                      fprintf( stderr, "m[%d] = %f\n", j, m[j] );
1593                        mjpt++;
1594                        prept++;
1595                        mpjpt++;
1596                        curpt++;
1597
1598                }
1599                lastverticalw[i] = currentw[lgth2-1];
1600
1601#if STOREWM
1602                WMMTX2[i][lgth2] = m[lgth2-1];
1603#endif
1604
1605#if 0  // ue
1606                if( i == imid )
1607                {
1608                        for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
1609                        for( j=0; j<lgth2; j++ ) midm[j] = m[j];
1610                }
1611#endif
1612        }
1613//      for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
1614//      for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
1615
1616#if 0
1617    for( i=0; i<lgth1; i++ )
1618    {
1619        for( j=0; j<lgth2; j++ )
1620        {
1621            fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
1622        }
1623        fprintf( stderr, "\n" );
1624    }
1625        fprintf( stderr, "\n" );
1626        fprintf( stderr, "WMMTX2 = \n" );
1627    for( i=0; i<lgth1; i++ )
1628    {
1629        for( j=0; j<lgth2; j++ )
1630        {
1631            fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
1632        }
1633        fprintf( stderr, "\n" );
1634    }
1635        fprintf( stderr, "\n" );
1636#endif
1637
1638// gyakudp
1639
1640        match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
1641        match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
1642
1643        for( i=0; i<lgth1-1; i++ )
1644        {
1645                initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
1646//              initverticalw[i] += fpenalty;
1647        }
1648        for( j=0; j<lgth2-1; j++ )
1649        {
1650                currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
1651//              currentw[j] += fpenalty;
1652        }
1653
1654#if STOREWM
1655        for( i=0; i<lgth1-1; i++ )
1656        {
1657                WMMTX[i][lgth2-1] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
1658//              fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
1659        }
1660        for( j=0; j<lgth2-1; j++ )
1661        {
1662                WMMTX[lgth1-1][j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
1663//              fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
1664        }
1665#endif
1666
1667
1668
1669
1670
1671
1672        for( j=lgth2-1; j>0; --j )
1673        {
1674                m[j-1] = currentw[j] + fgcp2[lgth2-2];
1675//              m[j-1] = currentw[j];
1676                mp[j] = lgth1-1;
1677        }
1678
1679//      for( j=0; j<lgth2; j++ ) m[j] = 0.0;
1680        // m[lgth2-1] ha irunoka?
1681
1682
1683//      for( i=lgth1-2; i>=imid; i-- )
1684        firstm = -9999999.9;
1685        firstmp = lgth1-1;
1686        for( i=lgth1-2; i>-1; i-- )
1687        {
1688                wtmp = previousw;
1689                previousw = currentw;
1690                currentw = wtmp;
1691                previousw[lgth2-1] = initverticalw[i+1];
1692//              match_calc( currentw, seq1, seq2, i, lgth2 );
1693                match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
1694
1695                currentw[lgth2-1] = initverticalw[i];
1696
1697//              m[lgth2] = fgcp1[i];
1698//              WMMTX2[i][lgth2] += m[lgth2];
1699//              fprintf( stderr, "m[] = %f\n", m[lgth2] );
1700
1701                mi = previousw[lgth2-1] + fgcp2[lgth2-2];
1702//              mi = previousw[lgth2-1];
1703                mpi = lgth2 - 1;
1704
1705                mjpt = m + lgth2 - 2;
1706                prept = previousw + lgth2 - 1;
1707                curpt = currentw + lgth2 - 2;
1708                mpjpt = mp + lgth2 - 2;
1709
1710
1711                for( j=lgth2-2; j>-1; j-- )
1712                {
1713                        wm = *prept;
1714                        ijpi = i+1;
1715                        ijpj = j+1;
1716
1717                        g = mi + ogcp2[j+1];
1718//                      g = mi + fpenalty;
1719                        if( g > wm )
1720                        {
1721                                wm = g;
1722                                ijpj = mpi;
1723                                ijpi = i+1;
1724                        }
1725
1726                        g = *prept + fgcp2[j];
1727//                      g = *prept;
1728                        if( g >= mi )
1729                        {
1730//                              fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
1731                                mi = g;
1732                                mpi = j + 1;
1733                        }
1734
1735#if USE_PENALTY_EX
1736                        mi += fpenalty_ex;
1737#endif
1738
1739//                      fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
1740                        g = *mjpt + ogcp1[i+1];
1741//                      g = *mjpt + fpenalty;
1742                        if( g > wm )
1743                        {
1744                                wm = g;
1745                                ijpi = *mpjpt;
1746                                ijpj = j+1;
1747                        }
1748
1749//                      if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
1750                        g = *prept + fgcp1[i];
1751//                      g = *prept;
1752                        if( g >= *mjpt )
1753                        {
1754                                *mjpt = g;
1755                                *mpjpt = i + 1;
1756                        }
1757
1758#if USE_PENALTY_EX
1759                        m[j] += fpenalty_ex;
1760#endif
1761
1762                        if( i == jumpi || i == imid - 1 )
1763                        {
1764                                jumpforwi[j] = ijpi; //muda
1765                                jumpforwj[j] = ijpj; //muda
1766//                              fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
1767//                              fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
1768                        }
1769                        if( i == imid ) // muda
1770                        {
1771                                midw[j] += wm;
1772//                              midm[j+1] += *mjpt + fpenalty; //??????
1773                                midm[j+1] += *mjpt; //??????
1774                        }
1775                        if( i == imid - 1 )
1776                        {
1777//                              midn[j] += mi + fpenalty;  //????
1778                                midn[j] += mi;  //????
1779                        }
1780#if LOCAL
1781            if( wm < localthr )
1782            {       
1783//                              fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
1784                                wm = 0;
1785            }       
1786#endif
1787
1788#if STOREWM
1789                        WMMTX[i][j] += wm;
1790//                      WMMTX2[i][j+1] += *mjpt + fpenalty;
1791                        WMMTX2[i][j] += *curpt;
1792#endif
1793                        *curpt += wm;
1794
1795                        mjpt--;
1796                        prept--;
1797                        mpjpt--;
1798                        curpt--;
1799                }
1800//              fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
1801                g = *prept + fgcp1[i];
1802                if( firstm < g ) 
1803                {
1804                        firstm = g;
1805                        firstmp = i + 1;
1806                }
1807#if STOREWM
1808//              WMMTX2[i][j+1] += firstm;
1809#endif
1810                if( i == imid ) midm[j+1] += firstm;
1811
1812                if( i == imid - 1 )     
1813                {
1814                        maxwm = midw[1];
1815                        jmid = 0;
1816//                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1817                        for( j=2; j<lgth2-1; j++ )
1818                        {
1819                                wm = midw[j];
1820                                if( wm > maxwm )
1821                                {
1822                                        jmid = j;
1823                                        maxwm = wm;
1824                                }
1825//                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1826                        }
1827                        for( j=0; j<lgth2+1; j++ )
1828                        {
1829                                wm = midm[j];
1830                                if( wm > maxwm )
1831                                {
1832                                        jmid = j;
1833                                        maxwm = wm;
1834                                }
1835//                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1836                        }
1837
1838//                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1839
1840
1841//                      fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
1842                        wm = midw[jmid];
1843                        jumpi = imid-1;
1844                        jumpj = jmid-1;
1845                        if( jmid > 0 && midn[jmid-1] > wm ) //060413
1846                        {
1847                                jumpi = imid-1;
1848                                jumpj = jumpbacki[jmid];
1849                                wm = midn[jmid-1];
1850//                              fprintf( stderr, "rejump (n)\n" );
1851                        }
1852                        if( midm[jmid] > wm )
1853                        {
1854                                jumpi = jumpbackj[jmid];
1855                                jumpj = jmid-1;
1856                                wm = midm[jmid];
1857//                              fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
1858                        }
1859
1860
1861//                      fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
1862//                      fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
1863#if 0
1864                        fprintf( stderr, "imid = %d\n", imid );
1865                        fprintf( stderr, "midn = \n" );
1866                        for( j=0; j<lgth2; j++ )
1867                        {
1868                                fprintf( stderr, "% 7.1f ", midn[j] );
1869                        }
1870                        fprintf( stderr, "\n" );
1871                        fprintf( stderr, "midw = \n" );
1872                        for( j=0; j<lgth2; j++ )
1873                        {
1874                                fprintf( stderr, "% 7.1f ", midw[j] );
1875                        }
1876                        fprintf( stderr, "\n" );
1877                        fprintf( stderr, "midm = \n" );
1878                        for( j=0; j<lgth2; j++ )
1879                        {
1880                                fprintf( stderr, "% 7.1f ", midm[j] );
1881                        }
1882                        fprintf( stderr, "\n" );
1883#endif
1884//                      fprintf( stderr, "maxwm = %f\n", maxwm );
1885                }
1886                if( i == jumpi ) //saki?
1887                {
1888//                      fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
1889//                      fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
1890                        if( jmid == 0 )
1891                        {
1892//                              fprintf( stderr, "CHUI2!\n" );
1893                                jumpj = 0; jmid = 1;
1894                                jumpi = firstmp - 1;
1895                                imid = firstmp;
1896                        }
1897
1898#if 0
1899                        else if( jmid == lgth2 )
1900                        {
1901                                fprintf( stderr, "CHUI1!\n" );
1902                                jumpi=0; jumpj=0;
1903                                imid=jumpforwi[0]; jmid=lgth2-1;
1904                        }
1905#else // 060414
1906                        else if( jmid >= lgth2 ) 
1907                        {
1908//                              fprintf( stderr, "CHUI1!\n" );
1909                                jumpi=imid-1; jmid=lgth2;
1910                                jumpj = lgth2-1;
1911                        }
1912#endif
1913                        else
1914                        {
1915                                imid = jumpforwi[jumpj];
1916                                jmid = jumpforwj[jumpj];
1917                        }
1918#if 0
1919                        fprintf( stderr, "jumpi -> %d\n", jumpi );
1920                        fprintf( stderr, "jumpj -> %d\n", jumpj );
1921                        fprintf( stderr, "imid -> %d\n", imid );
1922                        fprintf( stderr, "jmid -> %d\n", jmid );
1923#endif
1924
1925#if STOREWM
1926//                      break;
1927#else
1928                        break;
1929#endif
1930                }
1931        }
1932#if 0
1933                jumpi=0; jumpj=0;
1934                imid=lgth1-1; jmid=lgth2-1;
1935        }
1936#endif
1937
1938//      fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1939//      fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1940
1941//      for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1942//      for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1943//      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1944//      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1945
1946
1947
1948        for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
1949                map[i][j] = WMMTX[i][j] / maxwm;
1950//              map[i][j] = WMMTX2[i][j] / maxwm;
1951
1952#if STOREWM
1953
1954#if 0
1955        for( i=0; i<lgth1; i++ )
1956        {
1957                float maxpairscore = -9999.9;
1958                float tmpscore;
1959
1960                for( j=0; j<lgth2; j++ )
1961                {
1962                        if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1963                        {
1964                                map12[i].pos = j;
1965                                map12[i].score = tmpscore;
1966                                maxpairscore = tmpscore;
1967                        }
1968                }
1969
1970                for( k=0; k<lgth1; k++ )
1971                {
1972                        if( i == k ) continue;
1973                        if( map12[i].score <= WMMTX[k][map12[i].pos] )
1974                                break;
1975                }
1976                if( k != lgth1 )
1977                {
1978                        map12[i].pos = -1;
1979                        map12[i].score = -1.0;
1980                }
1981                fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, map12[i].pos, map12[i].score, seq1[0][i], seq2[0][map12[i].pos] );
1982        }
1983        for( j=0; j<lgth2; j++ )
1984        {
1985                float maxpairscore = -9999.9;
1986                float tmpscore;
1987
1988                for( i=0; i<lgth1; i++ )
1989                {
1990                        if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1991                        {
1992                                map21[j].pos = i;
1993                                map21[j].score = tmpscore;
1994                                maxpairscore = tmpscore;
1995                        }
1996                }
1997
1998                for( k=0; k<lgth2; k++ )
1999                {
2000                        if( j == k ) continue;
2001                        if( map21[j].score <= WMMTX[map21[j].pos][k] )
2002                                break;
2003                }
2004                if( k != lgth2 )
2005                {
2006                        map21[j].pos = -1;
2007                        map21[j].score = -1.0;
2008                }
2009                fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", j, map21[j].pos, map21[j].score, seq2[0][j], seq1[0][map21[j].pos] );
2010        }
2011
2012        for( i=0; i<lgth1; i++ )
2013        {
2014                if( map12[i].pos != -1 ) if( map21[map12[i].pos].pos != i )
2015                        fprintf( stderr, "ERROR i=%d, but map12[i].pos=%d and map21[map12[i].pos]=%d\n", i, map12[i].pos, map21[map12[i].pos].pos );
2016        }
2017#endif
2018
2019#if 0
2020        fprintf( stderr, "WMMTX = \n" );
2021    for( i=0; i<lgth1; i++ )
2022    {
2023        fprintf( stderr, "%d ", i );
2024        for( j=0; j<lgth2; j++ )
2025        {
2026            fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
2027        }
2028        fprintf( stderr, "\n" );
2029    }
2030//      fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
2031    for( i=0; i<lgth1; i++ )
2032    {
2033        fprintf( stderr, "%d ", i );
2034        for( j=0; j<lgth2+1; j++ )
2035        {
2036            fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
2037        }
2038        fprintf( stderr, "\n" );
2039    }
2040        exit( 1 );
2041#endif
2042
2043#if 0
2044    fprintf( stdout, "#WMMTX = \n" );
2045    for( i=0; i<lgth1; i++ )
2046    {
2047//        fprintf( stdout, "%d ", i );
2048        for( j=0; j<lgth2; j++ )
2049        {
2050//          if( WMMTX[i][j] > amino_dis['a']['g'] -1 )
2051                fprintf( stdout, "%d %d %8.1f", i, j, WMMTX[i][j] );
2052                        if( WMMTX[i][j] == maxwm )
2053                fprintf( stdout, "selected \n" );
2054                        else
2055                fprintf( stdout, "\n" );
2056        }
2057        fprintf( stdout, "\n" );
2058    }
2059        exit( 1 );
2060#endif
2061
2062#if 0
2063
2064        fprintf( stderr, "jumpbacki = \n" );
2065        for( j=0; j<lgth2; j++ )
2066        {
2067                fprintf( stderr, "% 7d ", jumpbacki[j] );
2068        }
2069        fprintf( stderr, "\n" );
2070        fprintf( stderr, "jumpbackj = \n" );
2071        for( j=0; j<lgth2; j++ )
2072        {
2073                fprintf( stderr, "% 7d ", jumpbackj[j] );
2074        }
2075        fprintf( stderr, "\n" );
2076        fprintf( stderr, "jumpforwi = \n" );
2077        for( j=0; j<lgth2; j++ )
2078        {
2079                fprintf( stderr, "% 7d ", jumpforwi[j] );
2080        }
2081        fprintf( stderr, "\n" );
2082        fprintf( stderr, "jumpforwj = \n" );
2083        for( j=0; j<lgth2; j++ )
2084        {
2085                fprintf( stderr, "% 7d ", jumpforwj[j] );
2086        }
2087        fprintf( stderr, "\n" );
2088#endif
2089
2090
2091#endif
2092
2093
2094//      Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
2095
2096#if 0 // irukamo
2097        resultlen = strlen( mseq1[0] );
2098        if( alloclen < resultlen || resultlen > N )
2099        {
2100                fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
2101                ErrorExit( "LENGTH OVER!\n" );
2102        }
2103#endif
2104
2105
2106
2107#if 0
2108        fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
2109        fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
2110
2111        fprintf( stderr, "imid = %d\n", imid );
2112        fprintf( stderr, "jmid = %d\n", jmid );
2113#endif
2114
2115
2116        FreeFloatVec( w1 );
2117        FreeFloatVec( w2 );
2118        FreeFloatVec( initverticalw );
2119        FreeFloatVec( lastverticalw );
2120        FreeFloatVec( midw );
2121        FreeFloatVec( midm );
2122        FreeFloatVec( midn );
2123
2124        FreeIntVec( jumpbacki );
2125        FreeIntVec( jumpbackj );
2126        FreeIntVec( jumpforwi );
2127        FreeIntVec( jumpforwj );
2128        FreeIntVec( jumpdummi );
2129        FreeIntVec( jumpdummj );
2130
2131        FreeFloatVec( m );
2132        FreeIntVec( mp );
2133
2134        FreeFloatMtx( floatwork );
2135        FreeIntMtx( intwork );
2136
2137#if STOREWM
2138        FreeFloatMtx( WMMTX );
2139        FreeFloatMtx( WMMTX2 );
2140#endif
2141
2142
2143        free( gaps );
2144#if MEMSAVE
2145        free( aseq1 );
2146        free( aseq2 );
2147#else
2148        FreeCharMtx( aseq1 );
2149        FreeCharMtx( aseq2 );
2150#endif
2151       
2152        return( value );
2153}
2154
2155
2156
2157float Lalignmm_hmout( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, float **map )
2158/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
2159{
2160//      int k;
2161        int i, j;
2162        int ll1, ll2;
2163        int lgth1, lgth2;
2164        float wm = 0.0;   /* int ?????? */
2165        char **mseq1;
2166        char **mseq2;
2167//      char **mseq;
2168        float *ogcp1;
2169        float *ogcp2;
2170        float *fgcp1;
2171        float *fgcp2;
2172        float **cpmx1;
2173        float **cpmx2;
2174        float **gapinfo;
2175//      float fpenalty;
2176        float fpenalty = (float)RNApenalty;
2177        int nglen1, nglen2;
2178
2179
2180
2181
2182
2183#if 0
2184        fprintf( stderr, "eff in SA+++align\n" );
2185        for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
2186#endif
2187
2188        nglen1 = seqlen( seq1[0] );
2189        nglen2 = seqlen( seq2[0] );
2190
2191        lgth1 = strlen( seq1[0] );
2192        lgth2 = strlen( seq2[0] );
2193
2194        ll1 = ( (int)(lgth1) ) + 100;
2195        ll2 = ( (int)(lgth2) ) + 100;
2196
2197        mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
2198        mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
2199
2200        gapinfo = AllocateFloatMtx( 4, 0 );
2201        ogcp1 = AllocateFloatVec( ll1+2 );
2202        ogcp2 = AllocateFloatVec( ll2+2 );
2203        fgcp1 = AllocateFloatVec( ll1+2 );
2204        fgcp2 = AllocateFloatVec( ll2+2 );
2205
2206
2207        cpmx1 = AllocateFloatMtx( ll1+2, 27 );
2208        cpmx2 = AllocateFloatMtx( ll2+2, 27 );
2209
2210        for( i=0; i<icyc; i++ ) 
2211        {
2212                if( strlen( seq1[i] ) != lgth1 )
2213                {
2214                        fprintf( stderr, "i = %d / %d\n", i, icyc );
2215                        fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2216                        exit( 1 );
2217                }
2218        }
2219        for( j=0; j<jcyc; j++ )
2220        {
2221                if( strlen( seq2[j] ) != lgth2 )
2222                {
2223                        fprintf( stderr, "j = %d / %d\n", j, jcyc );
2224                        fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2225                        exit( 1 );
2226                }
2227        }
2228
2229        MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
2230        MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
2231
2232
2233#if 1
2234
2235        if( sgap1 )
2236        {
2237                new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
2238                new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
2239                new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap2 );
2240                new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
2241        }
2242        else
2243        {
2244                st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
2245                st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
2246                st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
2247                st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
2248        }
2249
2250#if 1
2251        for( i=0; i<lgth1; i++ ) 
2252        {
2253                ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
2254                fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
2255//              fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
2256        }
2257        for( i=0; i<lgth2; i++ ) 
2258        {
2259                ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
2260                fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
2261//              fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
2262        }
2263#else
2264        for( i=0; i<lgth1; i++ ) 
2265        {
2266                ogcp1[i] = 0.5 * fpenalty;
2267                fgcp1[i] = 0.5 * fpenalty;
2268        }
2269        for( i=0; i<lgth2; i++ ) 
2270        {
2271                ogcp2[i] = 0.5 * fpenalty;
2272                fgcp2[i] = 0.5 * fpenalty;
2273        }
2274#endif
2275
2276        gapinfo[0] = ogcp1;
2277        gapinfo[1] = fgcp1;
2278        gapinfo[2] = ogcp2;
2279        gapinfo[3] = fgcp2;
2280#endif
2281
2282#if 0
2283        fprintf( stdout, "in MSalignmm.c\n" );
2284        for( i=0; i<icyc; i++ )
2285        {
2286                fprintf( stdout, ">%d of GROUP1\n", i );
2287                fprintf( stdout, "%s\n", seq1[i] );
2288        }
2289        for( i=0; i<jcyc; i++ )
2290        {
2291                fprintf( stdout, ">%d of GROUP2\n", i );
2292                fprintf( stdout, "%s\n", seq2[i] );
2293        }
2294        fflush( stdout );
2295#endif
2296
2297        wm = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo, map );
2298#if DEBUG
2299                fprintf( stderr, " seq1[0] = %s\n", seq1[0] );
2300                fprintf( stderr, " seq2[0] = %s\n", seq2[0] );
2301                fprintf( stderr, "mseq1[0] = %s\n", mseq1[0] );
2302                fprintf( stderr, "mseq2[0] = %s\n", mseq2[0] );
2303#endif
2304
2305//      fprintf( stderr, "wm = %f\n", wm );
2306
2307#if 0
2308
2309        for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2310        for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
2311
2312        if( seqlen( seq1[0] ) != nglen1 )
2313        {
2314                fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
2315                fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2316                exit( 1 );
2317        }
2318        if( seqlen( seq2[0] ) != nglen2 )
2319        {
2320                fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
2321                exit( 1 );
2322        }
2323#endif
2324
2325        FreeFloatVec( ogcp1 );
2326        FreeFloatVec( ogcp2 );
2327        FreeFloatVec( fgcp1 );
2328        FreeFloatVec( fgcp2 );
2329        FreeFloatMtx( cpmx1 );
2330        FreeFloatMtx( cpmx2 );
2331        free( (void *)gapinfo );
2332
2333        FreeCharMtx( mseq1 );
2334        FreeCharMtx( mseq2 );
2335
2336        lgth1 = strlen( seq1[0] );
2337        lgth2 = strlen( seq2[0] );
2338        for( i=0; i<icyc; i++ ) 
2339        {
2340                if( strlen( seq1[i] ) != lgth1 )
2341                {
2342                        fprintf( stderr, "i = %d / %d\n", i, icyc );
2343                        fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2344                        exit( 1 );
2345                }
2346        }
2347        for( j=0; j<jcyc; j++ )
2348        {
2349                if( strlen( seq2[j] ) != lgth2 )
2350                {
2351                        fprintf( stderr, "j = %d / %d\n", j, jcyc );
2352                        fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2353                        exit( 1 );
2354                }
2355        }
2356
2357        return( wm );
2358}
2359
2360float Lalign2m2m_hmout( char **seq1, char **seq2, char **seq1r, char **seq2r, char *dir1, char *dir2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, float **map )
2361/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
2362{
2363//      int k;
2364        int i, j;
2365        int ll1, ll2;
2366        int lgth1, lgth2;
2367        float wm = 0.0;   /* int ?????? */
2368        char **mseq1;
2369        char **mseq2;
2370        float *ogcp1;
2371        float *ogcp2;
2372        float *fgcp1;
2373        float *fgcp2;
2374        float **cpmx1;
2375        float **cpmx2;
2376        float **gapinfo;
2377        float fpenalty = (float)penalty;
2378        int nglen1, nglen2;
2379
2380#if 0
2381        fprintf( stderr, "eff in SA+++align\n" );
2382        for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
2383#endif
2384
2385        nglen1 = seqlen( seq1[0] );
2386        nglen2 = seqlen( seq2[0] );
2387
2388        lgth1 = strlen( seq1[0] );
2389        lgth2 = strlen( seq2[0] );
2390
2391        ll1 = ( (int)(lgth1) ) + 100;
2392        ll2 = ( (int)(lgth2) ) + 100;
2393
2394        mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
2395        mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
2396
2397        gapinfo = AllocateFloatMtx( 4, 0 );
2398        ogcp1 = AllocateFloatVec( ll1+2 );
2399        ogcp2 = AllocateFloatVec( ll2+2 );
2400        fgcp1 = AllocateFloatVec( ll1+2 );
2401        fgcp2 = AllocateFloatVec( ll2+2 );
2402
2403
2404        cpmx1 = AllocateFloatMtx( ll1+2, 39 );
2405        cpmx2 = AllocateFloatMtx( ll2+2, 39 );
2406
2407        for( i=0; i<icyc; i++ ) 
2408        {
2409                if( strlen( seq1[i] ) != lgth1 )
2410                {
2411                        fprintf( stderr, "i = %d / %d\n", i, icyc );
2412                        fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2413                        exit( 1 );
2414                }
2415        }
2416        for( j=0; j<jcyc; j++ )
2417        {
2418                if( strlen( seq2[j] ) != lgth2 )
2419                {
2420                        fprintf( stderr, "j = %d / %d\n", j, jcyc );
2421                        fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2422                        exit( 1 );
2423                }
2424        }
2425
2426#if 0
2427        MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
2428        MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
2429#else
2430        cpmx_ribosum( seq1, seq1r, dir1, cpmx1, eff1, lgth1, icyc );
2431        cpmx_ribosum( seq2, seq2r, dir2, cpmx2, eff2, lgth2, jcyc );
2432#endif
2433
2434
2435#if 1
2436
2437        if( sgap1 )
2438        {
2439                new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
2440                new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
2441                new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap2 );
2442                new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
2443        }
2444        else
2445        {
2446                st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
2447                st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
2448                st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
2449                st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
2450        }
2451
2452#if 1
2453        for( i=0; i<lgth1; i++ ) 
2454        {
2455                ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
2456                fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
2457//              fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
2458        }
2459        for( i=0; i<lgth2; i++ ) 
2460        {
2461                ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
2462                fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
2463//              fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
2464        }
2465#else
2466        for( i=0; i<lgth1; i++ ) 
2467        {
2468                ogcp1[i] = 0.5 * fpenalty;
2469                fgcp1[i] = 0.5 * fpenalty;
2470        }
2471        for( i=0; i<lgth2; i++ ) 
2472        {
2473                ogcp2[i] = 0.5 * fpenalty;
2474                fgcp2[i] = 0.5 * fpenalty;
2475        }
2476#endif
2477
2478        gapinfo[0] = ogcp1;
2479        gapinfo[1] = fgcp1;
2480        gapinfo[2] = ogcp2;
2481        gapinfo[3] = fgcp2;
2482#endif
2483
2484#if 0
2485        fprintf( stdout, "in MSalignmm.c\n" );
2486        for( i=0; i<icyc; i++ )
2487        {
2488                fprintf( stdout, ">%d of GROUP1\n", i );
2489                fprintf( stdout, "%s\n", seq1[i] );
2490        }
2491        for( i=0; i<jcyc; i++ )
2492        {
2493                fprintf( stdout, ">%d of GROUP2\n", i );
2494                fprintf( stdout, "%s\n", seq2[i] );
2495        }
2496        fflush( stdout );
2497#endif
2498
2499        wm = MSalign2m2m_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo, map );
2500#if DEBUG
2501                fprintf( stderr, " seq1[0] = %s\n", seq1[0] );
2502                fprintf( stderr, " seq2[0] = %s\n", seq2[0] );
2503                fprintf( stderr, "mseq1[0] = %s\n", mseq1[0] );
2504                fprintf( stderr, "mseq2[0] = %s\n", mseq2[0] );
2505#endif
2506
2507//      fprintf( stderr, "wm = %f\n", wm );
2508
2509#if 0
2510
2511        for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2512        for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
2513
2514        if( seqlen( seq1[0] ) != nglen1 )
2515        {
2516                fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
2517                fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2518                exit( 1 );
2519        }
2520        if( seqlen( seq2[0] ) != nglen2 )
2521        {
2522                fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
2523                exit( 1 );
2524        }
2525#endif
2526
2527        FreeFloatVec( ogcp1 );
2528        FreeFloatVec( ogcp2 );
2529        FreeFloatVec( fgcp1 );
2530        FreeFloatVec( fgcp2 );
2531        FreeFloatMtx( cpmx1 );
2532        FreeFloatMtx( cpmx2 );
2533        free( (void *)gapinfo );
2534
2535        FreeCharMtx( mseq1 );
2536        FreeCharMtx( mseq2 );
2537
2538        lgth1 = strlen( seq1[0] );
2539        lgth2 = strlen( seq2[0] );
2540        for( i=0; i<icyc; i++ ) 
2541        {
2542                if( strlen( seq1[i] ) != lgth1 )
2543                {
2544                        fprintf( stderr, "i = %d / %d\n", i, icyc );
2545                        fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2546                        exit( 1 );
2547                }
2548        }
2549        for( j=0; j<jcyc; j++ )
2550        {
2551                if( strlen( seq2[j] ) != lgth2 )
2552                {
2553                        fprintf( stderr, "j = %d / %d\n", j, jcyc );
2554                        fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2555                        exit( 1 );
2556                }
2557        }
2558
2559        return( wm );
2560}
Note: See TracBrowser for help on using the repository browser.