source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/core/MSalignmm.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: 39.4 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 0
9
10#define DPTANNI 100
11
12
13static TLS int reccycle = 0;
14
15
16static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
17{
18        int j, k, l;
19        float scarr[26];
20        float **cpmxpd = floatwork;
21        int **cpmxpdn = intwork;
22        int count = 0;
23        float *matchpt;
24        float **cpmxpdpt;
25        int **cpmxpdnpt;
26        int cpkd;
27
28        if( initialize )
29        {
30                for( j=0; j<lgth2; j++ )
31                {
32                        count = 0;
33                        for( l=0; l<26; l++ )
34                        {
35                                if( cpmx2[j][l] )
36                                {
37                                        cpmxpd[j][count] = cpmx2[j][l];
38                                        cpmxpdn[j][count] = l;
39                                        count++;
40                                }
41                        }
42                        cpmxpdn[j][count] = -1;
43                }
44        }
45
46        for( l=0; l<26; l++ )
47        {
48                scarr[l] = 0.0;
49                for( k=0; k<26; k++ )
50                {
51                        scarr[l] += n_dis[k][l] * cpmx1[i1][k];
52                }
53        }
54#if 0 /* €³€ì€ò»È€Š€È€­€Ïfloatwork€Î¥¢¥í¥±¡Œ¥È€òµÕ€Ë€¹€ë */
55        {
56                float *fpt, **fptpt, *fpt2;
57                int *ipt, **iptpt;
58                fpt2 = match;
59                iptpt = cpmxpdn;
60                fptpt = cpmxpd;
61                while( lgth2-- )
62                {
63                        *fpt2 = 0.0;
64                        ipt=*iptpt,fpt=*fptpt;
65                        while( *ipt > -1 )
66                                *fpt2 += scarr[*ipt++] * *fpt++;
67                        fpt2++,iptpt++,fptpt++;
68                }
69        }
70        for( j=0; j<lgth2; j++ )
71        {
72                match[j] = 0.0;
73                for( k=0; cpmxpdn[j][k]>-1; k++ )
74                        match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
75        }
76#else
77        matchpt = match;
78        cpmxpdnpt = cpmxpdn;
79        cpmxpdpt = cpmxpd;
80        while( lgth2-- )
81        {
82                *matchpt = 0.0;
83                for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
84                        *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
85                matchpt++;
86                cpmxpdnpt++;
87                cpmxpdpt++;
88        }
89#endif
90}
91
92static float Atracking( float *lasthorizontalw, float *lastverticalw,
93                                                char **seq1, char **seq2, 
94                        char **mseq1, char **mseq2, 
95                        int **ijp, int icyc, int jcyc,
96                                                int ist, int ien, int jst, int jen, 
97                                                int fulllen1, int fulllen2, int tailgp )
98{
99        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, klim;
100        char *gaptable1, *gt1bk;
101        char *gaptable2, *gt2bk;
102        float wm;
103        lgth1 = ien-ist+1;
104        lgth2 = jen-jst+1;
105
106        gt1bk = AllocateCharVec( lgth1+lgth2+1 );
107        gt2bk = AllocateCharVec( lgth1+lgth2+1 );
108
109#if 0
110        for( i=0; i<lgth1; i++ )
111        {
112                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
113        }
114#endif
115
116
117//      fprintf( stderr, "in Atracking, tailgp=%d, ien=%d, fulllen1=%d, jen=%d, fulllen2=%d\n", tailgp, ien, fulllen1, jen, fulllen2 );
118
119        if( tailgp == 1 )
120                ;
121        else if( ien == fulllen1-1 || jen == fulllen2-1 )
122        {
123//              fprintf( stderr, "searching lasthorizontalw and lasthorizontalw\n" );
124                wm = lastverticalw[0];
125                for( i=0; i<lgth1; i++ )
126                {
127                        if( lastverticalw[i] >= wm )
128                        {
129                                wm = lastverticalw[i];
130                                iin = i; jin = lgth2-1;
131                                ijp[lgth1][lgth2] = +( lgth1 - i );
132                        }
133                }
134                for( j=0; j<lgth2; j++ )
135                {
136                        if( lasthorizontalw[j] >= wm )
137                        {
138                                wm = lasthorizontalw[j];
139                                iin = lgth1-1; jin = j;
140                                ijp[lgth1][lgth2] = -( lgth2 - j );
141                        }
142                }
143        }
144#if 0
145        else if( jen == fulllen2-1 )
146        {
147                fprintf( stderr, "searching lastverticalw\n" );
148                wm = lastverticalw[0];
149                for( i=0; i<lgth1; i++ )
150                {
151                        if( lastverticalw[i] >= wm )
152                        {
153                                wm = lastverticalw[i];
154                                iin = i; jin = lgth2-1;
155                                ijp[lgth1][lgth2] = +( lgth1 - i );
156                        }
157                }
158        }
159        else if( ien == fulllen1-1 )
160        {
161                fprintf( stderr, "searching lasthorizontalw\n" );
162                wm = lasthorizontalw[0];
163                for( j=0; j<lgth2; j++ )
164                {
165                        if( lasthorizontalw[j] >= wm )
166                        {
167                                wm = lasthorizontalw[j];
168                                iin = lgth1-1; jin = j;
169                                ijp[lgth1][lgth2] = -( lgth2 - j );
170                        }
171                }
172        }
173#endif
174 
175    for( i=0; i<lgth1+1; i++ ) 
176    {
177        ijp[i][0] = i + 1;
178    }
179    for( j=0; j<lgth2+1; j++ ) 
180    {
181        ijp[0][j] = -( j + 1 );
182    }
183
184
185        gaptable1 = gt1bk + lgth1+lgth2;
186        *gaptable1 = 0;
187        gaptable2 = gt2bk + lgth1+lgth2;
188        *gaptable2 = 0;
189
190//      if( lgth2 == 1 ) fprintf( stderr, "in Atracking, mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
191
192        iin = lgth1; jin = lgth2;
193        klim = lgth1+lgth2;
194        for( k=0; k<=klim; k++ ) 
195        {
196                if( ijp[iin][jin] < 0 ) 
197                {
198                        ifi = iin-1; jfi = jin+ijp[iin][jin];
199                }
200                else if( ijp[iin][jin] > 0 )
201                {
202                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
203                }
204                else
205                {
206                        ifi = iin-1; jfi = jin-1;
207                }
208                l = iin - ifi;
209                while( --l ) 
210                {
211                        *--gaptable1 = 'o';
212                        *--gaptable2 = '-';
213                        k++;
214                }
215                l= jin - jfi;
216                while( --l )
217                {
218                        *--gaptable1 = '-';
219                        *--gaptable2 = 'o';
220                        k++;
221                }
222                if( iin <= 0 || jin <= 0 ) break;
223                *--gaptable1 = 'o';
224                *--gaptable2 = 'o';
225                k++;
226                iin = ifi; jin = jfi;
227
228        }
229
230        for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i]+ist, gaptable1 );
231        for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j]+jst, gaptable2 );
232
233        free( gt1bk );
234        free( gt2bk );
235
236//      fprintf( stderr, "in Atracking (owari), mseq1 = %s\n", mseq1[0] );
237//      fprintf( stderr, "in Atracking (owari), mseq2 = %s\n", mseq2[0] );
238        return( 0.0 );
239}
240
241static float MSalignmm_tanni( 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, int fulllen1, int fulllen2, char **mseq1, char **mseq2, float **gapinfo, int headgp, int tailgp )
242/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
243{
244//      int k;
245        register int i, j;
246        int ll1, ll2;
247        int lasti, lastj;
248        float wm = 0.0;   /* int ?????? */
249        float g;
250        float *currentw, *previousw;
251#if USE_PENALTY_EX
252        float fpenalty_ex = (float)penalty_ex;
253#endif
254#if 1
255        float *wtmp;
256        int *ijppt;
257        float *mjpt, *prept, *curpt;
258        int *mpjpt;
259#endif
260        float mi, *m;
261        int **ijp;
262        int mpi, *mp;
263        float *w1, *w2;
264        float *initverticalw;    /* kufuu sureba iranai */
265        float *lastverticalw;    /* kufuu sureba iranai */
266        int **intwork;
267        float **floatwork;
268        int **intmtx;
269        int lgth1, lgth2;
270        float *ogcp1;
271        float *fgcp1;
272        float *ogcp2;
273        float *fgcp2;
274//      char **aseq1;
275//      char **aseq2;
276//      char **aseq1bk, **aseq2bk;
277
278
279        ogcp1 = gapinfo[0] + ist;
280        fgcp1 = gapinfo[1] + ist;
281        ogcp2 = gapinfo[2] + jst;
282        fgcp2 = gapinfo[3] + jst;
283
284#if STOREWM
285        char ttt1[10000], ttt2[10000];
286#endif
287
288
289        lgth1 = ien-ist+1;
290        lgth2 = jen-jst+1;
291
292#if STOREWM
293        strncpy( ttt1, seq1[0]+ist, lgth1 ); ttt1[lgth1] = 0;
294        strncpy( ttt2, seq2[0]+jst, lgth2 ); ttt2[lgth2] = 0;
295
296        fprintf( stderr, "in _tanni ist,ien = %d,%d, lgth1=%d\n", ist, ien, lgth1 );
297        fprintf( stderr, "in _tanni jst,jen = %d,%d, lgth2=%d\n", jst, jen, lgth2 );
298        fprintf( stderr, "ttt1 = %s\n", ttt1 );
299        fprintf( stderr, "ttt2 = %s\n", ttt2 );
300#endif
301
302#if 0
303        fprintf( stderr, "in _tanni ist,ien = %d,%d, fulllen1=%d\n", ist, ien, fulllen1 );
304        fprintf( stderr, "in _tanni jst,jen = %d,%d, fulllen2=%d\n", jst, jen, fulllen2 );
305        fprintf( stderr, "in _tanni seq1[0] = %-*.*s\n", ien-ist+1, ien-ist+1, seq1[0]+ist );
306        fprintf( stderr, "in _tanni seq2[0] = %-*.*s\n", jen-jst+1, jen-jst+1, seq2[0]+jst );
307#endif
308
309
310        ll1 = ( (int)(lgth1) ) + 100;
311        ll2 = ( (int)(lgth2) ) + 100;
312
313//      aseq1 = AllocateCharMtx( icyc, 0 );
314//      aseq2 = AllocateCharMtx( jcyc, 0 );
315//      aseq1bk = AllocateCharMtx( icyc, lgth1+lgth2+100 );
316//      aseq2bk = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
317//      for( i=0; i<icyc; i++ ) aseq1[i] = aseq1bk[i];
318//      for( i=0; i<jcyc; i++ ) aseq2[i] = aseq2bk[i];
319
320        w1 = AllocateFloatVec( ll2+2 );
321        w2 = AllocateFloatVec( ll2+2 );
322
323        initverticalw = AllocateFloatVec( ll1+2 );
324        lastverticalw = AllocateFloatVec( ll1+2 );
325
326        m = AllocateFloatVec( ll2+2 );
327        mp = AllocateIntVec( ll2+2 );
328
329        floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 27 ); 
330        intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 ); 
331
332
333        intmtx = AllocateIntMtx( ll1+1, ll2+1 );
334
335        ijp = intmtx;
336
337        currentw = w1;
338        previousw = w2;
339
340        match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
341
342        match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
343
344        if( headgp || ist != 0 )
345        {
346                for( i=1; i<lgth1+1; i++ )
347                {
348                        initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
349                }
350        }
351        if( headgp || jst != 0 )
352        {
353                for( j=1; j<lgth2+1; j++ )
354                {
355                        currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
356                }
357        }
358
359        for( j=1; j<lgth2+1; ++j ) 
360        {
361                m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;;
362        }
363
364        lastverticalw[0] = currentw[lgth2-1];
365
366
367
368        if( tailgp || jen != fulllen2-1 ) lasti = lgth1+1; else lasti = lgth1;
369//      if( 1 ) lasti = lgth1+1; else lasti = lgth1;
370        for( i=1; i<lasti; i++ )
371        {
372                wtmp = previousw; 
373                previousw = currentw;
374                currentw = wtmp;
375
376                previousw[0] = initverticalw[i-1];
377
378                match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
379                currentw[0] = initverticalw[i];
380
381                mi = previousw[0] + ogcp2[1]; 
382                mpi = 0;
383
384                ijppt = ijp[i] + 1;
385                mjpt = m + 1;
386                prept = previousw;
387                curpt = currentw + 1;
388                mpjpt = mp + 1;
389//              if( tailgp && jen != fulllen2-1 ) lastj = lgth2+1; else lastj = lgth2;
390                lastj = lgth2+1; 
391                for( j=1; j<lastj; j++ )
392                {
393                        wm = *prept;
394                        *ijppt = 0;
395
396#if 0
397                        fprintf( stderr, "%5.0f->", wm );
398#endif
399                        g = mi + fgcp2[j-1];
400#if 0
401                        fprintf( stderr, "%5.0f?", g );
402#endif
403                        if( g > wm )
404                        {
405                                wm = g;
406                                *ijppt = -( j - mpi );
407                        }
408                        g = *prept + ogcp2[j];
409                        if( g >= mi )
410                        {
411                                mi = g;
412                                mpi = j-1;
413                        }
414#if USE_PENALTY_EX
415                        mi += fpenalty_ex;
416#endif
417
418                        g = *mjpt + fgcp1[i-1];
419#if 0 
420                        fprintf( stderr, "%5.0f?", g );
421#endif
422                        if( g > wm )
423                        {
424                                wm = g;
425                                *ijppt = +( i - *mpjpt );
426                        }
427                        g = *prept + ogcp1[i];
428                        if( g >= *mjpt )
429                        {
430                                *mjpt = g;
431                                *mpjpt = i-1;
432                        }
433#if USE_PENALTY_EX
434                        m[j] += fpenalty_ex;
435#endif
436
437#if 0
438                        fprintf( stderr, "%5.0f ", wm );
439#endif
440                        *curpt += wm;
441
442
443                        ijppt++;
444                        mjpt++;
445                        prept++;
446                        mpjpt++;
447                        curpt++;
448                }
449                lastverticalw[i] = currentw[lgth2-1];
450        }
451
452//      fprintf( stderr, "wm = %f\n", wm );
453
454        Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, ist, ien, jst, jen, fulllen1, fulllen2, tailgp );
455#if 0
456        fprintf( stderr, "res in _tanni mseq1[0] = %s\n", mseq1[0] );
457        fprintf( stderr, "res in _tanni mseq2[0] = %s\n", mseq2[0] );
458#endif
459
460//      for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
461//      for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
462
463//      fprintf( stderr, "in _tanni, aseq1 = %s\n", aseq1[0] );
464//      fprintf( stderr, "in _tanni, mseq1 = %s\n", mseq1[0] );
465
466        FreeFloatVec( w1 );
467        FreeFloatVec( w2 );
468        FreeFloatVec( initverticalw );
469        FreeFloatVec( lastverticalw );
470
471        FreeFloatVec( m );
472        FreeIntVec( mp );
473
474
475        FreeFloatMtx( floatwork );
476        FreeIntMtx( intwork );
477
478        FreeIntMtx( intmtx );
479
480
481//      FreeCharMtx( aseq1bk );
482//      FreeCharMtx( aseq2bk );
483
484//      free( aseq1 );
485//      free( aseq2 );
486
487        return( wm );
488
489}
490
491static void freearrays_rec1(
492        float *w1, float *w2, float *initverticalw, float *lastverticalw,
493        float *midw, float *midm, float *midn,
494        int *jumpbacki, int *jumpbackj, int *jumpforwi, int *jumpforwj, int *jumpdummi, int *jumpdummj,
495        float *m, int *mp,
496        float **floatwork, int **intwork
497#if STOREWM
498        , float **WMMTX, float **WMMTX2
499#endif
500)
501{
502
503        FreeFloatVec( w1 );
504        FreeFloatVec( w2 );
505        FreeFloatVec( initverticalw );
506        FreeFloatVec( lastverticalw );
507        FreeFloatVec( midw );
508        FreeFloatVec( midm );
509        FreeFloatVec( midn );
510
511        FreeIntVec( jumpbacki );
512        FreeIntVec( jumpbackj );
513        FreeIntVec( jumpforwi );
514        FreeIntVec( jumpforwj );
515        FreeIntVec( jumpdummi );
516        FreeIntVec( jumpdummj );
517
518        FreeFloatVec( m );
519        FreeIntVec( mp );
520
521        FreeFloatMtx( floatwork );
522        FreeIntMtx( intwork );
523
524#if STOREWM
525        FreeFloatMtx( WMMTX );
526        FreeFloatMtx( WMMTX2 );
527#endif
528}
529
530static void freearrays_rec2( char *gaps, char **aseq1, char **aseq2 )
531{
532        free( gaps );
533#if MEMSAVE
534        free( aseq1 );
535        free( aseq2 );
536#else
537        FreeCharMtx( aseq1 );
538        FreeCharMtx( aseq2 );
539#endif
540}
541
542static 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, int fulllen1, int fulllen2, char **mseq1, char **mseq2, int depth, float **gapinfo, int *chudanpt, int chudanref, int *chudanres, int headgp, int tailgp )
543/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
544{
545//      int k;
546        int alnlen;
547        float value = 0.0;
548        register int i, j;
549        char **aseq1, **aseq2;
550        int ll1, ll2, l, len;
551        int lasti, lastj, imid;
552        int jmid = 0; // by D.Mathog, a guess
553        float wm = 0.0;   /* int ?????? */
554        float g;
555        float *currentw, *previousw;
556#if USE_PENALTY_EX
557        float fpenalty_ex = (float)penalty_ex;
558#endif
559        float *wtmp;
560//      short *ijppt;
561        int *mpjpt;
562//      short **ijp;
563        int *mp;
564        int mpi;
565        float *mjpt, *prept, *curpt;
566        float mi;
567        float *m;
568        float *w1, *w2;
569//      float *match;
570        float *initverticalw;    /* kufuu sureba iranai */
571        float *lastverticalw;    /* kufuu sureba iranai */
572        int **intwork;
573        float **floatwork;
574//      short **shortmtx;
575#if STOREWM
576        float **WMMTX;
577        float **WMMTX2;
578#endif
579        float *midw;
580        float *midm;
581        float *midn;
582        int lgth1, lgth2;
583        float maxwm;
584        int *jumpforwi;
585        int *jumpforwj;
586        int *jumpbacki;
587        int *jumpbackj;
588        int *jumpdummi; //muda
589        int *jumpdummj = NULL; // by D.Mathog, a guess
590        int jumpi, jumpj = 0; // by D.Mathog, a guess
591        char *gaps;
592        int ijpi, ijpj;
593        float *ogcp1;
594        float *fgcp1;
595        float *ogcp2;
596        float *fgcp2;
597        float firstm;
598        int firstmp;
599#if STOREWM
600        static TLS char ttt1[50000];
601        static TLS char ttt2[50000];
602#endif
603
604#if 0
605        int nglen1, nglen2;
606        nglen1 = seqlen( seq1[0] );
607        nglen2 = seqlen( seq2[0] );
608#endif
609
610//      fprintf( stderr, "fulllen1 = %d, fulllen2 = %d, headgp = %d, tailgp = %d\n", fulllen1, fulllen2, headgp, tailgp );
611
612        ogcp1 = gapinfo[0] + ist;
613        fgcp1 = gapinfo[1] + ist;
614        ogcp2 = gapinfo[2] + jst;
615        fgcp2 = gapinfo[3] + jst;
616
617        depth++;
618        reccycle++;
619
620        lgth1 = ien-ist+1;
621        lgth2 = jen-jst+1;
622
623//      if( lgth1 < 5 )
624//              fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
625//      if( lgth2 < 5 )
626//              fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
627//
628
629
630#if STOREWM
631        fprintf( stderr, "==== MSalign (depth=%d, reccycle=%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, reccycle, ist, ien, jst, jen );
632        strncpy( ttt1, seq1[0]+ist, lgth1 );
633        strncpy( ttt2, seq2[0]+jst, lgth2 );
634        ttt1[lgth1] = 0;
635        ttt2[lgth2] = 0;
636        fprintf( stderr, "seq1 = %s\n", ttt1 );
637        fprintf( stderr, "seq2 = %s\n", ttt2 );
638#endif
639        if( lgth2 <= 0 ) // lgth1 <= 0 ha?
640        {
641//              fprintf( stderr, "\n\n==== jimei\n\n" );
642//              exit( 1 );
643                for( i=0; i<icyc; i++ ) 
644                {
645                        strncpy( mseq1[i], seq1[i]+ist, lgth1 );
646                        mseq1[i][lgth1] = 0;
647                }
648                for( i=0; i<jcyc; i++ ) 
649                {
650                        mseq2[i][0] = 0;
651                        for( j=0; j<lgth1; j++ )
652//                              strcat( mseq2[i], "-" );
653                                strcat( mseq2[i], newgapstr );
654                }
655
656//              fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
657//              fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
658
659                return( 0.0 );
660        }
661
662#if MEMSAVE
663        aseq1 = AllocateCharMtx( icyc, 0 );
664        aseq2 = AllocateCharMtx( jcyc, 0 );
665        for( i=0; i<icyc; i++ ) aseq1[i] = mseq1[i];
666        for( i=0; i<jcyc; i++ ) aseq2[i] = mseq2[i];
667#else
668        aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
669        aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
670#endif
671
672//      fprintf( stderr, "####(s) seq1[0] (%d) = \n%-*.*s\n (a%d-a%d)\n", depth, ien-ist+1, ien-ist+1, seq1[0]+ist, ist, ien );
673//      fprintf( stderr, "####(s) seq2[0] (%d) = \n%-*.*s\n (b%d-b%d)\n", depth, jen-jst+1, jen-jst+1, seq2[0]+jst, jst, jen );
674
675//  if( lgth1 < DPTANNI && lgth2 < DPTANNI ) // & dato lgth ==1 no kanousei ga arunode yokunai
676//    if( lgth1 < DPTANNI ) // kore mo lgth2 ga mijikasugiru kanousei ari
677    if( lgth1 < DPTANNI || lgth2 < DPTANNI ) // zettai ni anzen ka?
678        {
679//              fprintf( stderr, "==== Going to _tanni\n" );
680
681                value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, fulllen1, fulllen2, aseq1, aseq2, gapinfo, headgp, tailgp );   
682
683
684#if MEMSAVE
685                free( aseq1 );
686                free( aseq2 );
687#else
688                for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
689                for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
690
691                FreeCharMtx( aseq1 );
692                FreeCharMtx( aseq2 );
693#endif
694
695//              fprintf( stderr, "value = %f\n", value );
696
697                return( value );
698        }
699//      fprintf( stderr, "Trying to divide the mtx\n" );
700
701        ll1 = ( (int)(lgth1) ) + 100;
702        ll2 = ( (int)(lgth2) ) + 100;
703
704//      fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
705
706        w1 = AllocateFloatVec( ll2+2 );
707        w2 = AllocateFloatVec( ll2+2 );
708//      match = AllocateFloatVec( ll2+2 );
709        midw = AllocateFloatVec( ll2+2 );
710        midn = AllocateFloatVec( ll2+2 );
711        midm = AllocateFloatVec( ll2+2 );
712        jumpbacki = AllocateIntVec( ll2+2 );
713        jumpbackj = AllocateIntVec( ll2+2 );
714        jumpforwi = AllocateIntVec( ll2+2 );
715        jumpforwj = AllocateIntVec( ll2+2 );
716        jumpdummi = AllocateIntVec( ll2+2 );
717        jumpdummj = AllocateIntVec( ll2+2 );
718
719        initverticalw = AllocateFloatVec( ll1+2 );
720        lastverticalw = AllocateFloatVec( ll1+2 );
721
722        m = AllocateFloatVec( ll2+2 );
723        mp = AllocateIntVec( ll2+2 );
724        gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
725
726        floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 ); 
727        intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 ); 
728
729#if DEBUG
730        fprintf( stderr, "succeeded\n" );
731#endif
732
733#if STOREWM
734        WMMTX = AllocateFloatMtx( ll1, ll2 );
735        WMMTX2 = AllocateFloatMtx( ll1, ll2 );
736#endif
737#if 0
738        shortmtx = AllocateShortMtx( ll1, ll2 );
739
740#if DEBUG
741        fprintf( stderr, "succeeded\n\n" );
742#endif
743
744        ijp = shortmtx;
745#endif
746
747        currentw = w1;
748        previousw = w2;
749
750        match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
751
752        match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
753
754        for( i=1; i<lgth1+1; i++ )
755        {
756                initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
757        }
758        for( j=1; j<lgth2+1; j++ )
759        {
760                currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
761        }
762
763#if STOREWM
764        WMMTX[0][0] = initverticalw[0];
765        for( i=1; i<lgth1+1; i++ )
766        {
767                WMMTX[i][0] = initverticalw[i];
768        }
769        for( j=1; j<lgth2+1; j++ )
770        {
771                WMMTX[0][j] = currentw[j];
772        }
773#endif
774
775
776        for( j=1; j<lgth2+1; ++j ) 
777        {
778                m[j] = currentw[j-1] + ogcp1[1];
779//              m[j] = currentw[j-1];
780                mp[j] = 0;
781        }
782
783        lastverticalw[0] = currentw[lgth2-1];
784
785        imid = lgth1 * 0.5;
786
787        jumpi = 0; // atode kawaru.
788        lasti = lgth1+1;
789#if STOREWM
790        for( i=1; i<lasti; i++ )
791#else
792        for( i=1; i<=imid; i++ )
793#endif
794        {
795#ifdef enablemultithread
796//              fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
797                if( chudanpt && *chudanpt != chudanref ) 
798                {
799//                      fprintf( stderr, "\n\n## CHUUDAN!!! zenhan\n" );
800                        *chudanres = 1;
801                        freearrays_rec1
802                        (
803                                w1, w2, initverticalw, lastverticalw, midw, midm, midn,
804                                jumpbacki, jumpbackj, jumpforwi, jumpforwj, jumpdummi, jumpdummj,
805                                m, mp,
806                                floatwork, intwork
807#if STOREWM
808                                , WMMTX, WMMTX2
809#endif
810                        );
811                        freearrays_rec2( gaps, aseq1, aseq2 );
812                        return( -1.0 );
813                }
814#endif
815                wtmp = previousw; 
816                previousw = currentw;
817                currentw = wtmp;
818
819                previousw[0] = initverticalw[i-1];
820
821                match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
822                currentw[0] = initverticalw[i];
823
824                m[0] = ogcp1[i];
825#if STOREM
826                WMMTX2[i][0] = m[0];
827#endif
828                if( i == imid ) midm[0] = m[0];
829
830                mi = previousw[0] + ogcp2[1]; 
831//              mi = previousw[0];
832                mpi = 0;
833
834
835//              ijppt = ijp[i] + 1;
836                mjpt = m + 1;
837                prept = previousw;
838                curpt = currentw + 1;
839                mpjpt = mp + 1;
840
841
842                lastj = lgth2+1;
843                for( j=1; j<lastj; j++ )
844                {
845
846                        wm = *prept;
847
848#if 0
849                        fprintf( stderr, "%5.0f->", wm );
850#endif
851                        g = mi + fgcp2[j-1];
852//                      g = mi + fpenalty;
853#if 0
854                        fprintf( stderr, "%5.0f?", g );
855#endif
856                        if( g > wm )
857                        {
858                                wm = g;
859//                              *ijppt = -( j - mpi );
860                        }
861                        g = *prept + ogcp2[j];
862//                      g = *prept;
863                        if( g >= mi )
864                        {
865                                mi = g;
866                                mpi = j-1;
867                        }
868#if USE_PENALTY_EX
869                        mi += fpenalty_ex;
870#endif
871
872                        g = *mjpt + fgcp1[i-1];
873//                      g = *mjpt + fpenalty;
874#if 0 
875                        fprintf( stderr, "%5.0f?", g );
876#endif
877                        if( g > wm )
878                        {
879                                wm = g;
880//                              *ijppt = +( i - *mpjpt );
881                        }
882
883
884                        g = *prept + ogcp1[i];
885//                      g = *prept;
886                        if( g >= *mjpt )
887                        {
888                                *mjpt = g;
889                                *mpjpt = i-1;
890                        }
891#if USE_PENALTY_EX
892                        m[j] += fpenalty_ex;
893#endif
894
895#if 0
896                        fprintf( stderr, "%5.0f ", wm );
897#endif
898                        *curpt += wm;
899
900#if STOREWM
901                        WMMTX[i][j] = *curpt;
902                        WMMTX2[i][j] = *mjpt;
903#endif
904
905                        if( i == imid ) //muda
906                        {       
907                                jumpbackj[j] = *mpjpt; // muda atode matomeru
908                                jumpbacki[j] = mpi; // muda atode matomeru
909//                              fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
910//                              fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
911                                midw[j] = *curpt;
912                                midm[j] = *mjpt;
913                                midn[j] = mi;
914                        }
915
916//                      fprintf( stderr, "m[%d] = %f\n", j, m[j] );
917                        mjpt++;
918                        prept++;
919                        mpjpt++;
920                        curpt++;
921
922                }
923                lastverticalw[i] = currentw[lgth2-1];
924
925#if STOREWM
926                WMMTX2[i][lgth2] = m[lgth2-1];
927#endif
928
929#if 0  // ue
930                if( i == imid )
931                {
932                        for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
933                        for( j=0; j<lgth2; j++ ) midm[j] = m[j];
934                }
935#endif
936        }
937//      for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
938//      for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
939
940#if 0
941    for( i=0; i<lgth1; i++ )
942    {
943        for( j=0; j<lgth2; j++ )
944        {
945            fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
946        }
947        fprintf( stderr, "\n" );
948    }
949        fprintf( stderr, "\n" );
950        fprintf( stderr, "WMMTX2 = \n" );
951    for( i=0; i<lgth1; i++ )
952    {
953        for( j=0; j<lgth2; j++ )
954        {
955            fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
956        }
957        fprintf( stderr, "\n" );
958    }
959        fprintf( stderr, "\n" );
960#endif
961
962// gyakudp
963
964        match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
965        match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
966
967        for( i=0; i<lgth1-1; i++ )
968        {
969                initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
970//              initverticalw[i] += fpenalty;
971        }
972        for( j=0; j<lgth2-1; j++ )
973        {
974                currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
975//              currentw[j] += fpenalty;
976        }
977
978#if STOREWM
979        for( i=0; i<lgth1-1; i++ )
980        {
981                WMMTX[i][lgth2-1] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
982                fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
983        }
984        for( j=0; j<lgth2-1; j++ )
985        {
986                WMMTX[lgth1-1][j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
987                fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
988        }
989#endif
990
991
992
993
994
995
996        for( j=lgth2-1; j>0; --j )
997        {
998                m[j-1] = currentw[j] + fgcp2[lgth2-2];
999//              m[j-1] = currentw[j];
1000                mp[j] = lgth1-1;
1001        }
1002
1003//      for( j=0; j<lgth2; j++ ) m[j] = 0.0;
1004        // m[lgth2-1] ha irunoka?
1005
1006
1007//      for( i=lgth1-2; i>=imid; i-- )
1008        firstm = -9999999.9;
1009//      firstmp = lgth1-1;
1010        firstmp = lgth1;
1011        for( i=lgth1-2; i>-1; i-- )
1012        {
1013#ifdef enablemultithread
1014//              fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
1015                if( chudanpt && *chudanpt != chudanref ) 
1016                {
1017//                      fprintf( stderr, "\n\n## CHUUDAN!!! kouhan\n" );
1018                        *chudanres = 1;
1019                        freearrays_rec1
1020                        (
1021                                w1, w2, initverticalw, lastverticalw, midw, midm, midn,
1022                                jumpbacki, jumpbackj, jumpforwi, jumpforwj, jumpdummi, jumpdummj,
1023                                m, mp,
1024                                floatwork, intwork
1025#if STOREWM
1026                                , WMMTX, WMMTX2
1027#endif
1028                        );
1029                        freearrays_rec2( gaps, aseq1, aseq2 );
1030                        return( -1.0 );
1031                }
1032#endif
1033                wtmp = previousw;
1034                previousw = currentw;
1035                currentw = wtmp;
1036                previousw[lgth2-1] = initverticalw[i+1];
1037//              match_calc( currentw, seq1, seq2, i, lgth2 );
1038                match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
1039
1040                currentw[lgth2-1] = initverticalw[i];
1041
1042//              m[lgth2] = fgcp1[i];
1043//              WMMTX2[i][lgth2] += m[lgth2];
1044//              fprintf( stderr, "m[] = %f\n", m[lgth2] );
1045
1046                mi = previousw[lgth2-1] + fgcp2[lgth2-2];
1047//              mi = previousw[lgth2-1];
1048                mpi = lgth2 - 1;
1049
1050                mjpt = m + lgth2 - 2;
1051                prept = previousw + lgth2 - 1;
1052                curpt = currentw + lgth2 - 2;
1053                mpjpt = mp + lgth2 - 2;
1054
1055
1056                for( j=lgth2-2; j>-1; j-- )
1057                {
1058                        wm = *prept;
1059                        ijpi = i+1;
1060                        ijpj = j+1;
1061
1062                        g = mi + ogcp2[j+1];
1063//                      g = mi + fpenalty;
1064                        if( g > wm )
1065                        {
1066                                wm = g;
1067                                ijpj = mpi;
1068                                ijpi = i+1;
1069                        }
1070
1071                        g = *prept + fgcp2[j];
1072//                      g = *prept;
1073                        if( g >= mi )
1074                        {
1075//                              fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
1076                                mi = g;
1077                                mpi = j + 1;
1078                        }
1079
1080#if USE_PENALTY_EX
1081                        mi += fpenalty_ex;
1082#endif
1083
1084//                      fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
1085                        g = *mjpt + ogcp1[i+1];
1086//                      g = *mjpt + fpenalty;
1087                        if( g > wm )
1088                        {
1089                                wm = g;
1090                                ijpi = *mpjpt;
1091                                ijpj = j+1;
1092                        }
1093
1094//                      if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
1095                        g = *prept + fgcp1[i];
1096//                      g = *prept;
1097                        if( g >= *mjpt )
1098                        {
1099                                *mjpt = g;
1100                                *mpjpt = i + 1;
1101                        }
1102
1103#if USE_PENALTY_EX
1104                        m[j] += fpenalty_ex;
1105#endif
1106
1107                        if( i == jumpi || i == imid - 1 )
1108                        {
1109                                jumpforwi[j] = ijpi; //muda
1110                                jumpforwj[j] = ijpj; //muda
1111//                              fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
1112//                              fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
1113                        }
1114                        if( i == imid ) // muda
1115                        {
1116                                midw[j] += wm;
1117//                              midm[j+1] += *mjpt + fpenalty; //??????
1118                                midm[j+1] += *mjpt; //??????
1119                        }
1120                        if( i == imid - 1 )
1121                        {
1122//                              midn[j] += mi + fpenalty;  //????
1123                                midn[j] += mi;  //????
1124                        }
1125#if STOREWM
1126                        WMMTX[i][j] += wm;
1127//                      WMMTX2[i][j+1] += *mjpt + fpenalty;
1128                        WMMTX2[i][j+1] += *mjpt;
1129#endif
1130                        *curpt += wm;
1131
1132                        mjpt--;
1133                        prept--;
1134                        mpjpt--;
1135                        curpt--;
1136                }
1137//              fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
1138                g = *prept + fgcp1[i];
1139                if( firstm < g ) 
1140                {
1141                        firstm = g;
1142                        firstmp = i + 1;
1143                }
1144#if STOREWM
1145                WMMTX2[i][j+1] += firstm;
1146#endif
1147                if( i == imid ) midm[j+1] += firstm;
1148
1149
1150                if( i == imid - 1 )     
1151                {
1152                        maxwm = midw[1];
1153                        jmid = 0;
1154//                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1155                        for( j=2; j<lgth2-1; j++ )
1156                        {
1157                                wm = midw[j];
1158                                if( wm > maxwm )
1159                                {
1160                                        jmid = j;
1161                                        maxwm = wm;
1162                                }
1163//                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1164                        }
1165                        for( j=0; j<lgth2+1; j++ )
1166                        {
1167                                wm = midm[j];
1168                                if( wm > maxwm )
1169                                {
1170                                        jmid = j;
1171                                        maxwm = wm;
1172                                }
1173//                              if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1174                        }
1175
1176//                      if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1177
1178
1179//                      fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
1180                        wm = midw[jmid];
1181                        jumpi = imid-1;
1182                        jumpj = jmid-1;
1183                        if( jmid > 0 && midn[jmid-1] > wm ) //060413
1184                        {
1185                                jumpi = imid-1;
1186                                jumpj = jumpbacki[jmid];
1187                                wm = midn[jmid-1];
1188//                              fprintf( stderr, "rejump (n)\n" );
1189                        }
1190                        if( midm[jmid] > wm )
1191                        {
1192                                jumpi = jumpbackj[jmid];
1193                                jumpj = jmid-1;
1194                                wm = midm[jmid];
1195//                              fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
1196                        }
1197
1198
1199//                      fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
1200//                      fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
1201#if STOREWM
1202                        fprintf( stderr, "imid = %d\n", imid );
1203                        fprintf( stderr, "midn = \n" );
1204                        for( j=0; j<lgth2; j++ )
1205                        {
1206                                fprintf( stderr, "% 7.1f ", midn[j] );
1207                        }
1208                        fprintf( stderr, "\n" );
1209                        fprintf( stderr, "midw = \n" );
1210                        for( j=0; j<lgth2; j++ )
1211                        {
1212                                fprintf( stderr, "% 7.1f ", midw[j] );
1213                        }
1214                        fprintf( stderr, "\n" );
1215                        fprintf( stderr, "midm = \n" );
1216                        for( j=0; j<lgth2; j++ )
1217                        {
1218                                fprintf( stderr, "% 7.1f ", midm[j] );
1219                        }
1220                        fprintf( stderr, "\n" );
1221#endif
1222//                      fprintf( stderr, "maxwm = %f\n", maxwm );
1223                }
1224                if( i == jumpi ) //saki?
1225                {
1226//                      fprintf( stderr, "#### FIRST i=%d, jumpi<imid=%d<%d, ist=%d, ien=%d, firstmp=%d, lgth1=%d\n", i, jumpi, imid, ist, ien, firstmp, lgth1 );
1227//                      fprintf( stderr, "#### mark 1\n" );
1228//                      fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
1229//                      fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
1230
1231                        if( jmid == 0 )
1232                        {
1233//                              fprintf( stderr, "#### CHUI2!\n" );
1234                                jumpj = 0; jmid = 1;
1235#if 0 // v6.823 made
1236                                jumpi = firstmp-1;
1237                                imid = firstmp;
1238#endif
1239#if 0
1240                                jumpi = 0;
1241                                imid = 1;
1242#else
1243//                              if( 1 || firstmp > 100 ) // naze 100
1244                                if( imid < firstmp-1 ) // naze 100
1245                                {
1246                                        jumpi = firstmp;
1247                                        imid = firstmp+1;
1248                                }
1249#if 0
1250                                else
1251                                {
1252                                        jumpi = 0;
1253                                        imid = 1;
1254                                }
1255#endif
1256#endif
1257                        }
1258
1259#if 0
1260                        else if( jmid == lgth2 )
1261                        {
1262                                fprintf( stderr, "CHUI1!\n" );
1263                                jumpi=0; jumpj=0;
1264                                imid=jumpforwi[0]; jmid=lgth2-1;
1265                        }
1266#else // 060414
1267                        else if( jmid >= lgth2 ) 
1268                        {
1269//                              fprintf( stderr, "CHUI1!\n" );
1270                                jumpi=imid-1; jmid=lgth2;
1271                                jumpj = lgth2-1;
1272                        }
1273#endif
1274                        else
1275                        {
1276//                              fprintf( stderr, "#### CHUI3!\n" );
1277                                imid = jumpforwi[jumpj];
1278                                jmid = jumpforwj[jumpj];
1279                                if( imid == jumpi ) jumpi = imid-1;
1280                        }
1281#if 0
1282                        fprintf( stderr, "jumpi -> %d\n", jumpi );
1283                        fprintf( stderr, "jumpj -> %d\n", jumpj );
1284                        fprintf( stderr, "imid -> %d\n", imid );
1285                        fprintf( stderr, "jmid -> %d\n", jmid );
1286#endif
1287//                      fprintf( stderr, "#### FINAL i=%d, jumpi<imid=%d<%d, ist=%d, ien=%d, firstmp=%d\n", i, jumpi, imid, ist, ien, firstmp );
1288
1289#if STOREWM
1290                        break;
1291#else
1292                        break;
1293#endif
1294                }
1295        }
1296
1297//      fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1298//      fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1299
1300//      for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1301//      for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1302//      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1303//      for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1304
1305
1306#if STOREWM
1307        fprintf( stderr, "WMMTX = \n" );
1308    for( i=0; i<lgth1; i++ )
1309    {
1310        fprintf( stderr, "%d ", i );
1311        for( j=0; j<lgth2; j++ )
1312        {
1313            fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
1314        }
1315        fprintf( stderr, "\n" );
1316    }
1317//      fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
1318    for( i=0; i<lgth1; i++ )
1319    {
1320        fprintf( stderr, "%d ", i );
1321        for( j=0; j<lgth2+1; j++ )
1322        {
1323            fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
1324        }
1325        fprintf( stderr, "\n" );
1326    }
1327
1328        fprintf( stderr, "jumpbacki = \n" );
1329        for( j=0; j<lgth2; j++ )
1330        {
1331                fprintf( stderr, "% 7d ", jumpbacki[j] );
1332        }
1333        fprintf( stderr, "\n" );
1334        fprintf( stderr, "jumpbackj = \n" );
1335        for( j=0; j<lgth2; j++ )
1336        {
1337                fprintf( stderr, "% 7d ", jumpbackj[j] );
1338        }
1339        fprintf( stderr, "\n" );
1340        fprintf( stderr, "jumpforwi = \n" );
1341        for( j=0; j<lgth2; j++ )
1342        {
1343                fprintf( stderr, "% 7d ", jumpforwi[j] );
1344        }
1345        fprintf( stderr, "\n" );
1346        fprintf( stderr, "jumpforwj = \n" );
1347        for( j=0; j<lgth2; j++ )
1348        {
1349                fprintf( stderr, "% 7d ", jumpforwj[j] );
1350        }
1351        fprintf( stderr, "\n" );
1352
1353
1354#endif
1355
1356
1357//      Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1358
1359#if 0 // irukamo
1360        resultlen = strlen( mseq1[0] );
1361        if( alloclen < resultlen || resultlen > N )
1362        {
1363                fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1364                ErrorExit( "LENGTH OVER!\n" );
1365        }
1366#endif
1367
1368
1369
1370#if 0
1371        fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
1372        fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
1373
1374        fprintf( stderr, "imid = %d\n", imid );
1375        fprintf( stderr, "jmid = %d\n", jmid );
1376#endif
1377
1378        freearrays_rec1
1379        (
1380                w1, w2, initverticalw, lastverticalw, midw, midm, midn,
1381                jumpbacki, jumpbackj, jumpforwi, jumpforwj, jumpdummi, jumpdummj,
1382                m, mp,
1383                floatwork, intwork
1384#if STOREWM
1385                , WMMTX, WMMTX2
1386#endif
1387        );
1388
1389
1390//      fprintf( stderr, "==== calling myself (first)\n" );
1391
1392        value = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ist+jumpi, jst, jst+jumpj, alloclen, fulllen1, fulllen2, aseq1, aseq2, depth, gapinfo, NULL, 0, NULL, headgp, tailgp ); // chudan mada   
1393#if 0
1394                fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
1395                fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
1396#endif
1397#if MEMSAVE
1398#else
1399        for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
1400        for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
1401#endif
1402
1403//      fprintf( stderr, "====(f) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
1404//      fprintf( stderr, "====(f) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
1405
1406        len = strlen( mseq1[0] );
1407//      fprintf( stderr, "len = %d\n", len );
1408        l = jmid - jumpj - 1;
1409//      fprintf( stderr, "l=%d\n", l );
1410        if( l > 0 )
1411        {
1412//              for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
1413                for( i=0; i<l; i++ ) gaps[i] = *newgapstr; gaps[i] = 0;
1414                for( i=0; i<icyc; i++ ) 
1415                {
1416                        strcat( mseq1[i], gaps );
1417                        mseq1[i][len+l] = 0;
1418                }
1419                for( j=0; j<jcyc; j++ )
1420                {
1421                        strncat( mseq2[j], seq2[j]+jst+jumpj+1, l );
1422                        mseq2[j][len+l] = 0;
1423                }
1424//              fprintf( stderr, "penalizing (2) .. %f(%d), %f(%d)\n", ogcp2[jumpj+1], jumpj+1, fgcp2[jmid-1], jmid-1 );
1425                value +=  ( ogcp2[jumpj+1] + fgcp2[jmid-1] );
1426//              value += fpenalty;
1427        }
1428        len = strlen( mseq1[0] );
1429        l = imid - jumpi - 1;
1430//      fprintf( stderr, "l=%d\n", l );
1431        if( l > 0 )
1432        {
1433//              for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
1434                for( i=0; i<l; i++ ) gaps[i] = *newgapstr; gaps[i] = 0;
1435                for( i=0; i<icyc; i++ )
1436                {
1437                        strncat( mseq1[i], seq1[i]+ist+jumpi+1, l );
1438                        mseq1[i][len+l] = 0;
1439                }
1440                for( j=0; j<jcyc; j++ ) 
1441                {
1442                        strcat( mseq2[j], gaps );
1443                        mseq2[j][len+l] = 0;
1444                }
1445
1446//              for( i=0; i<lgth1; i++ ) fprintf( stderr, "ogcp1[%d] = %f\n", i, ogcp1[i] );
1447//              for( i=0; i<lgth1; i++ ) fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
1448
1449
1450//              fprintf( stderr, "penalizing (1) .. ogcp1[%d] = %f, fgcp1[%d] = %f\n", jumpi+1, ogcp1[jumpi+1], imid-1, fgcp1[imid-1] );
1451                value += ( ogcp1[jumpi+1] + fgcp1[imid-1] );
1452//              value += fpenalty;
1453        }
1454#if 0
1455        for( i=0; i<icyc; i++ ) fprintf( stderr, "after gapfill mseq1[%d]=%s\n", i, mseq1[i] );
1456        for( i=0; i<jcyc; i++ ) fprintf( stderr, "after gapfill mseq2[%d]=%s\n", i, mseq2[i] );
1457#endif
1458
1459//      fprintf( stderr, "==== calling myself (second)\n" );
1460
1461#if MEMSAVE
1462        alnlen = strlen( aseq1[0] );
1463        for( i=0; i<icyc; i++ ) aseq1[i] += alnlen;
1464        for( i=0; i<jcyc; i++ ) aseq2[i] += alnlen;
1465#endif
1466
1467        value += MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist+imid, ien, jst+jmid, jen, alloclen, fulllen1, fulllen2, aseq1, aseq2, depth, gapinfo, NULL, 0, NULL, headgp, tailgp ); // chudan mada
1468#if 0
1469                fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
1470                fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
1471#endif
1472
1473
1474
1475#if DEBUG
1476        if( value - maxwm > 1 || maxwm - value > 1 )
1477        {
1478                fprintf( stderr, "WARNING value  = %f, but maxwm = %f\n", value, maxwm );
1479                for( i=0; i<icyc; i++ )
1480                {
1481                        fprintf( stderr, ">1-%d\n%s\n", i, mseq1[i] );
1482                        fprintf( stderr, "%s\n", aseq1[i] );
1483                }
1484                for( i=0; i<jcyc; i++ )
1485                {
1486                        fprintf( stderr, ">2-%d\n%s\n", i, mseq2[i] );
1487                        fprintf( stderr, "%s\n", aseq2[i] );
1488                }
1489
1490//              exit( 1 );
1491        }
1492        else
1493        {
1494                fprintf( stderr, "value = %.0f, maxwm = %.0f -> ok\n", value, maxwm );
1495        }
1496#endif
1497
1498#if MEMSAVE
1499#else
1500        for( i=0; i<icyc; i++ ) strcat( mseq1[i], aseq1[i] );
1501        for( i=0; i<jcyc; i++ ) strcat( mseq2[i], aseq2[i] );
1502#endif
1503
1504//      fprintf( stderr, "====(s) aseq1[0] (%d) = \n%s\n (a%d-a%d)\n", depth, mseq1[0], ist, ien );
1505//      fprintf( stderr, "====(s) aseq2[0] (%d) = \n%s\n (b%d-b%d)\n", depth, mseq2[0], jst, jen );
1506
1507        freearrays_rec2( gaps, aseq1, aseq2 );
1508
1509#if 0
1510        if( seqlen( seq1[0] ) != nglen1 )
1511        {
1512                fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
1513                fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
1514                exit( 1 );
1515        }
1516        else
1517                fprintf( stderr, "nglen1 is ok in _rec\n" );
1518        if( seqlen( seq2[0] ) != nglen2 )
1519        {
1520                fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
1521                exit( 1 );
1522        }
1523        else
1524                fprintf( stderr, "nglen2 is ok in _rec\n" );
1525#endif
1526
1527        return( value );
1528}
1529
1530static void freearrays( 
1531        float *ogcp1, 
1532        float *ogcp2, 
1533        float *fgcp1, 
1534        float *fgcp2, 
1535        float **cpmx1,
1536        float **cpmx2, 
1537        float **gapinfo, 
1538        char **mseq1, 
1539        char **mseq2
1540)
1541{
1542        FreeFloatVec( ogcp1 );
1543        FreeFloatVec( ogcp2 );
1544        FreeFloatVec( fgcp1 );
1545        FreeFloatVec( fgcp2 );
1546        FreeFloatMtx( cpmx1 );
1547        FreeFloatMtx( cpmx2 );
1548        free( (void *)gapinfo );
1549
1550        FreeCharMtx( mseq1 );
1551        FreeCharMtx( mseq2 );
1552}
1553
1554float MSalignmm( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, int *chudanpt, int chudanref, int *chudanres, int headgp, int tailgp )
1555/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1556{
1557//      int k;
1558        int i, j;
1559        int ll1, ll2;
1560        int lgth1, lgth2;
1561        float wm = 0.0;   /* int ?????? */
1562        char **mseq1;
1563        char **mseq2;
1564        float *ogcp1;
1565        float *ogcp2;
1566        float *fgcp1;
1567        float *fgcp2;
1568        float **cpmx1;
1569        float **cpmx2;
1570        float **gapinfo;
1571        float fpenalty = (float)penalty;
1572        int nglen1, nglen2;
1573
1574#if 0
1575        fprintf( stderr, "eff in SA+++align\n" );
1576        for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1577#endif
1578
1579        nglen1 = seqlen( seq1[0] );
1580        nglen2 = seqlen( seq2[0] );
1581
1582#if 0
1583        fprintf( stderr, "\n" );
1584        for( i=0; i<icyc; i++ ) fprintf( stderr, "seq1[%d] at root = %s\n", i, seq1[i] );
1585        for( j=0; j<jcyc; j++ ) fprintf( stderr, "seq2[%d] at root = %s\n", j, seq2[j] );
1586        fprintf( stderr, "\n" );
1587#endif
1588
1589        lgth1 = strlen( seq1[0] );
1590        lgth2 = strlen( seq2[0] );
1591
1592        ll1 = ( (int)(lgth1) ) + 100;
1593        ll2 = ( (int)(lgth2) ) + 100;
1594
1595        mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
1596        mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
1597
1598        gapinfo = AllocateFloatMtx( 4, 0 );
1599        ogcp1 = AllocateFloatVec( ll1+2 );
1600        ogcp2 = AllocateFloatVec( ll2+2 );
1601        fgcp1 = AllocateFloatVec( ll1+2 );
1602        fgcp2 = AllocateFloatVec( ll2+2 );
1603
1604
1605        cpmx1 = AllocateFloatMtx( ll1+2, 27 );
1606        cpmx2 = AllocateFloatMtx( ll2+2, 27 );
1607
1608        for( i=0; i<icyc; i++ ) 
1609        {
1610                if( strlen( seq1[i] ) != lgth1 )
1611                {
1612                        fprintf( stderr, "i = %d / %d\n", i, icyc );
1613                        fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
1614                        exit( 1 );
1615                }
1616        }
1617        for( j=0; j<jcyc; j++ )
1618        {
1619                if( strlen( seq2[j] ) != lgth2 )
1620                {
1621                        fprintf( stderr, "j = %d / %d\n", j, jcyc );
1622                        fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
1623                        exit( 1 );
1624                }
1625        }
1626
1627        MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1628        MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1629
1630
1631#if 1
1632
1633        if( sgap1 )
1634        {
1635                new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
1636                new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
1637                new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap2 );
1638                new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
1639        }
1640        else
1641        {
1642                st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1643                st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1644                st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1645                st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1646        }
1647
1648#if 1
1649        for( i=0; i<lgth1; i++ ) 
1650        {
1651                ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
1652                fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
1653//              fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
1654        }
1655        for( i=0; i<lgth2; i++ ) 
1656        {
1657                ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
1658                fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
1659//              fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
1660        }
1661#else
1662        for( i=0; i<lgth1; i++ ) 
1663        {
1664                ogcp1[i] = 0.5 * fpenalty;
1665                fgcp1[i] = 0.5 * fpenalty;
1666        }
1667        for( i=0; i<lgth2; i++ ) 
1668        {
1669                ogcp2[i] = 0.5 * fpenalty;
1670                fgcp2[i] = 0.5 * fpenalty;
1671        }
1672#endif
1673
1674        gapinfo[0] = ogcp1;
1675        gapinfo[1] = fgcp1;
1676        gapinfo[2] = ogcp2;
1677        gapinfo[3] = fgcp2;
1678#endif
1679
1680#if 0
1681        fprintf( stdout, "in MSalignmm.c\n" );
1682        for( i=0; i<icyc; i++ )
1683        {
1684                fprintf( stdout, ">%d of GROUP1\n", i );
1685                fprintf( stdout, "%s\n", seq1[i] );
1686        }
1687        for( i=0; i<jcyc; i++ )
1688        {
1689                fprintf( stdout, ">%d of GROUP2\n", i );
1690                fprintf( stdout, "%s\n", seq2[i] );
1691        }
1692        fflush( stdout );
1693#endif
1694
1695        wm = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, lgth1, lgth2, mseq1, mseq2, 0, gapinfo, chudanpt, chudanref, chudanres, headgp, tailgp );
1696#ifdef enablemultithread
1697        if( chudanres && *chudanres ) 
1698        {
1699//              fprintf( stderr, "\n\n## CHUUDAN!!! relay\n" );
1700                *chudanres = 1;
1701                freearrays( ogcp1, ogcp2, fgcp1, fgcp2, cpmx1, cpmx2, gapinfo, mseq1, mseq2 );
1702                return( -1.0 );
1703        }
1704#endif
1705
1706#if 0
1707                fprintf( stderr, "\n" );
1708                fprintf( stderr, " seq1[0] = %s\n", seq1[0] );
1709                fprintf( stderr, " seq2[0] = %s\n", seq2[0] );
1710                fprintf( stderr, "mseq1[0] = %s\n", mseq1[0] );
1711                fprintf( stderr, "mseq2[0] = %s\n", mseq2[0] );
1712                fprintf( stderr, "\n" );
1713#endif
1714
1715//      fprintf( stderr, "wm = %f\n", wm );
1716
1717
1718        for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1719        for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
1720
1721        if( seqlen( seq1[0] ) != nglen1 )
1722        {
1723                fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
1724                fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
1725                exit( 1 );
1726        }
1727        if( seqlen( seq2[0] ) != nglen2 )
1728        {
1729                fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
1730                exit( 1 );
1731        }
1732
1733
1734        freearrays( ogcp1, ogcp2, fgcp1, fgcp2, cpmx1, cpmx2, gapinfo, mseq1, mseq2 );
1735
1736        lgth1 = strlen( seq1[0] );
1737        lgth2 = strlen( seq2[0] );
1738        for( i=0; i<icyc; i++ ) 
1739        {
1740                if( strlen( seq1[i] ) != lgth1 )
1741                {
1742                        fprintf( stderr, "i = %d / %d\n", i, icyc );
1743                        fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
1744                        exit( 1 );
1745                }
1746        }
1747        for( j=0; j<jcyc; j++ )
1748        {
1749                if( strlen( seq2[j] ) != lgth2 )
1750                {
1751                        fprintf( stderr, "j = %d / %d\n", j, jcyc );
1752                        fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
1753                        exit( 1 );
1754                }
1755        }
1756
1757#if 0
1758        fprintf( stderr, "\n" );
1759        for( i=0; i<icyc; i++ ) fprintf( stderr, " seq1[i] = %s\n", seq1[i] );
1760        for( j=0; j<jcyc; j++ ) fprintf( stderr, " seq2[j] = %s\n", seq2[j] );
1761        fprintf( stderr, "\n" );
1762#endif
1763
1764        return( wm );
1765}
Note: See TracBrowser for help on using the repository browser.