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

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

updated mafft version. Added extensions (no svn ignore, yet)
File size: 54.0 KB
Line 
1#include "mltaln.h"
2#include "dp.h"
3
4#define MACHIGAI 0
5#define OUTGAP0TRY 1
6#define DEBUG 0
7#define XXXXXXX    0
8#define USE_PENALTY_EX  0
9#define FASTMATCHCALC 1
10
11
12
13static TLS float **impmtx = NULL;
14static TLS int impalloclen = 0;
15#if 1 // tditeration to naiveQscore_imp de tsukawareru.
16float imp_match_out_scQ( int i1, int j1 )
17{
18//      fprintf( stderr, "imp+match = %f\n", impmtx[i1][j1] * fastathreshold );
19//      fprintf( stderr, "val = %f\n", impmtx[i1][j1] );
20        return( impmtx[i1][j1] );
21}
22#endif
23
24static void imp_match_out_veadQ_gapmap( float *imp, int i1, int lgth2, int *gapmap2 )
25{
26#if FASTMATCHCALC
27        float *pt = impmtx[i1];
28        int *gapmappt = gapmap2;
29        while( lgth2-- )
30                *imp++ += pt[*gapmappt++];
31#else
32        int j;
33        float *pt = impmtx[i1];
34        for( j=0; j<lgth2; j++ )
35                *imp++ += pt[gapmap2[j]];
36#endif
37}
38
39
40static void imp_match_out_vead_tateQ_gapmap( float *imp, int j1, int lgth1, int *gapmap1 )
41{
42#if FASTMATCHCALC
43        int *gapmappt = gapmap1;
44        while( lgth1-- )
45                *imp++ += impmtx[*gapmappt++][j1];
46#else
47        int i;
48        for( i=0; i<lgth1; i++ )
49                *imp++ += impmtx[gapmap1[i]][j1];
50#endif
51}
52
53
54static void imp_match_out_veadQ( float *imp, int i1, int lgth2 )
55{
56#if FASTMATCHCALC
57        float *pt = impmtx[i1];
58        while( lgth2-- )
59                *imp++ += *pt++;
60#else
61        int j;
62        float *pt = impmtx[i1];
63        for( j=0; j<lgth2; j++ )
64                *imp++ += pt[j];
65#endif
66}
67static void imp_match_out_vead_tateQ( float *imp, int j1, int lgth1 )
68{
69        int i;
70        for( i=0; i<lgth1; i++ )
71                *imp++ += impmtx[i][j1];
72}
73
74void imp_rnaQ( int nseq1, int nseq2, char **seq1, char **seq2, double *eff1, double *eff2, RNApair ***grouprna1, RNApair ***grouprna2, int *gapmap1, int *gapmap2, RNApair *additionalpair )
75{
76        foldrna( nseq1, nseq2, seq1, seq2, eff1, eff2, grouprna1, grouprna2, impmtx, gapmap1, gapmap2, additionalpair );
77}
78
79#if 1 // tbfast.c kara yobareru.
80void imp_match_init_strictQ( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom, int forscore )
81{
82        int i, j, k1, k2, tmpint, start1, start2, end1, end2;
83        float effij;
84        double effijx;
85        char *pt, *pt1, *pt2;
86        static TLS char *nocount1 = NULL;
87        static TLS char *nocount2 = NULL;
88        LocalHom *tmpptr;
89
90        if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
91        {
92                if( impmtx ) FreeFloatMtx( impmtx );
93                if( nocount1 ) free( nocount1 );
94                if( nocount2 ) free( nocount2 );
95                impalloclen = MAX( lgth1, lgth2 ) + 2;
96                impmtx = AllocateFloatMtx( impalloclen, impalloclen );
97                nocount1 = AllocateCharVec( impalloclen );
98                nocount2 = AllocateCharVec( impalloclen );
99        }
100
101        for( i=0; i<lgth1; i++ )
102        {
103                for( j=0; j<clus1; j++ )
104                        if( seq1[j][i] == '-' ) break;
105                if( j != clus1 ) nocount1[i] = 1; 
106                else                     nocount1[i] = 0;
107        }
108        for( i=0; i<lgth2; i++ )
109        {
110                for( j=0; j<clus2; j++ )
111                        if( seq2[j][i] == '-' ) break;
112                if( j != clus2 ) nocount2[i] = 1;
113                else                     nocount2[i] = 0;
114        }
115
116#if 0
117fprintf( stderr, "nocount2 =\n" );
118for( i = 0; i<impalloclen; i++ )
119{
120        fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
121}
122#endif
123
124
125
126#if 0
127        fprintf( stderr, "eff1 in _init_strict = \n" );
128        for( i=0; i<clus1; i++ )
129                fprintf( stderr, "eff1[] = %f\n", eff1[i] );
130        for( i=0; i<clus2; i++ )
131                fprintf( stderr, "eff2[] = %f\n", eff2[i] );
132#endif
133
134        for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
135                impmtx[i][j] = 0.0;
136        effijx =  fastathreshold;
137        for( i=0; i<clus1; i++ )
138        {
139                for( j=0; j<clus2; j++ )
140                {
141                        effij = (float)( eff1[i] * eff2[j] * effijx );
142                        tmpptr = localhom[i][j];
143                        while( tmpptr )
144                        {
145//                              fprintf( stderr, "start1 = %d\n", tmpptr->start1 );
146//                              fprintf( stderr, "end1   = %d\n", tmpptr->end1   );
147//                              fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] );
148//                              fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] );
149                                pt = seq1[i];
150                                tmpint = -1;
151                                while( *pt != 0 )
152                                {
153                                        if( *pt++ != '-' ) tmpint++;
154                                        if( tmpint == tmpptr->start1 ) break;
155                                }
156                                start1 = pt - seq1[i] - 1;
157       
158                                if( tmpptr->start1 == tmpptr->end1 ) end1 = start1;
159                                else
160                                {
161#if MACHIGAI
162                                        while( *pt != 0 )
163                                        {
164//                                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
165                                                if( tmpint == tmpptr->end1 ) break;
166                                                if( *pt++ != '-' ) tmpint++;
167                                        }
168                                        end1 = pt - seq1[i] - 0;
169#else
170                                        while( *pt != 0 )
171                                        {
172//                                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
173                                                if( *pt++ != '-' ) tmpint++;
174                                                if( tmpint == tmpptr->end1 ) break;
175                                        }
176                                        end1 = pt - seq1[i] - 1;
177#endif
178                                }
179       
180                                pt = seq2[j];
181                                tmpint = -1;
182                                while( *pt != 0 )
183                                {
184                                        if( *pt++ != '-' ) tmpint++;
185                                        if( tmpint == tmpptr->start2 ) break;
186                                }
187                                start2 = pt - seq2[j] - 1;
188                                if( tmpptr->start2 == tmpptr->end2 ) end2 = start2;
189                                else
190                                {
191#if MACHIGAI
192                                        while( *pt != 0 )
193                                        {
194                                                if( tmpint == tmpptr->end2 ) break;
195                                                if( *pt++ != '-' ) tmpint++;
196                                        }
197                                        end2 = pt - seq2[j] - 0;
198#else
199                                        while( *pt != 0 )
200                                        {
201                                                if( *pt++ != '-' ) tmpint++;
202                                                if( tmpint == tmpptr->end2 ) break;
203                                        }
204                                        end2 = pt - seq2[j] - 1;
205#endif
206                                }
207//                              fprintf( stderr, "start1 = %d (%c), end1 = %d (%c), start2 = %d (%c), end2 = %d (%c)\n", start1, seq1[i][start1], end1, seq1[i][end1], start2, seq2[j][start2], end2, seq2[j][end2] );
208//                              fprintf( stderr, "step 0\n" );
209                                if( end1 - start1 != end2 - start2 )
210                                {
211//                                      fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
212                                }
213
214#if 1
215                                k1 = start1; k2 = start2;
216                                pt1 = seq1[i] + k1;
217                                pt2 = seq2[j] + k2;
218                                while( *pt1 && *pt2 )
219                                {
220                                        if( *pt1 != '-' && *pt2 != '-' )
221                                        {
222// œÅ€ß€òÆóœÅ€Ë€«€±€Ê€€€è€Š€ËÃí°Õ€·€Æ²Œ€µ€€¡£
223//                                              impmtx[k1][k2] += tmpptr->wimportance * fastathreshold;
224//                                              impmtx[k1][k2] += tmpptr->importance * effij;
225                                                impmtx[k1][k2] += tmpptr->fimportance * effij;
226//                                              fprintf( stderr, "#### impmtx[k1][k2] = %f, tmpptr->fimportance=%f, effij=%f\n", impmtx[k1][k2], tmpptr->fimportance, effij );
227//                                              fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
228//                                              fprintf( stderr, "%d (%c) - %d (%c)  - %f\n", k1, *pt1, k2, *pt2, tmpptr->fimportance * effij );
229                                                k1++; k2++;
230                                                pt1++; pt2++;
231                                        }
232                                        else if( *pt1 != '-' && *pt2 == '-' )
233                                        {
234//                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
235                                                k2++; pt2++;
236                                        }
237                                        else if( *pt1 == '-' && *pt2 != '-' )
238                                        {
239//                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
240                                                k1++; pt1++;
241                                        }
242                                        else if( *pt1 == '-' && *pt2 == '-' )
243                                        {
244//                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
245                                                k1++; pt1++;
246                                                k2++; pt2++;
247                                        }
248                                        if( k1 > end1 || k2 > end2 ) break;
249                                }
250#else
251                                while( k1 <= end1 && k2 <= end2 )
252                                {
253                                        fprintf( stderr, "k1,k2=%d,%d - ", k1, k2 );
254                                        if( !nocount1[k1] && !nocount2[k2] )
255                                        {
256                                                impmtx[k1][k2] += tmpptr->wimportance * eff1[i] * eff2[j]  * fastathreshold;
257                                                fprintf( stderr, "marked\n" );
258                                        }
259                                        else
260                                                fprintf( stderr, "no count\n" );
261                                        k1++; k2++;
262                                }
263#endif
264                                tmpptr = tmpptr->next;
265                        }
266                }
267        }
268#if 0
269        if( clus1 == 1 && clus2 == 6 )
270        {
271                fprintf( stderr, "\n" );
272                fprintf( stderr, "seq1[0] =  %s\n", seq1[0] );
273                fprintf( stderr, "seq2[0] =  %s\n", seq2[0] );
274                fprintf( stderr, "impmtx = \n" );
275                for( k2=0; k2<lgth2; k2++ )
276                        fprintf( stderr, "%6.3f ", (double)k2 );
277                fprintf( stderr, "\n" );
278                for( k1=0; k1<lgth1; k1++ )
279                {
280                        fprintf( stderr, "%d ", k1 );
281                        for( k2=0; k2<3; k2++ )
282                                fprintf( stderr, "%2.1f ", impmtx[k1][k2] );
283                        fprintf( stderr, "\n" );
284                }
285                exit( 1 );
286        }
287#endif
288}
289#endif
290
291
292static void clearvec( float *match, int lgth )
293{
294        while( lgth-- ) *match++ = 0.0;
295}
296
297static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
298{
299#if FASTMATCHCALC
300        int j, l;
301        float scarr[26];
302        float **cpmxpd = floatwork;
303        int **cpmxpdn = intwork;
304        float *matchpt, *cpmxpdpt, **cpmxpdptpt;
305        int *cpmxpdnpt, **cpmxpdnptpt;
306        if( initialize )
307        {
308                int count = 0;
309                for( j=0; j<lgth2; j++ )
310                {
311                        count = 0;
312                        for( l=0; l<26; l++ )
313                        {
314                                if( cpmx2[l][j] )
315                                {
316                                        cpmxpd[j][count] = cpmx2[l][j];
317                                        cpmxpdn[j][count] = l;
318                                        count++;
319                                }
320                        }
321                        cpmxpdn[j][count] = -1;
322                }
323        }
324
325        {
326                for( l=0; l<26; l++ )
327                {
328                        scarr[l] = 0.0;
329                        for( j=0; j<26; j++ )
330                                scarr[l] += n_dis_consweight_multi[j][l] * cpmx1[j][i1];
331//                      scarr[l] *= consweight_multi;
332                }
333                matchpt = match;
334                cpmxpdnptpt = cpmxpdn;
335                cpmxpdptpt = cpmxpd;
336                while( lgth2-- )
337                {
338                        *matchpt = 0.0;
339                        cpmxpdnpt = *cpmxpdnptpt++;
340                        cpmxpdpt = *cpmxpdptpt++;
341                        while( *cpmxpdnpt>-1 )
342                                *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;
343                        matchpt++;
344                } 
345        }
346#else
347        int j, k, l;
348        float scarr[26];
349        float **cpmxpd = floatwork;
350        int **cpmxpdn = intwork;
351// simple
352        if( initialize )
353        {
354                int count = 0;
355                for( j=0; j<lgth2; j++ )
356                {
357                        count = 0;
358                        for( l=0; l<26; l++ )
359                        {
360                                if( cpmx2[l][j] )
361                                {
362                                        cpmxpd[count][j] = cpmx2[l][j];
363                                        cpmxpdn[count][j] = l;
364                                        count++;
365                                }
366                        }
367                        cpmxpdn[count][j] = -1;
368                }
369        }
370        for( l=0; l<26; l++ )
371        {
372                scarr[l] = 0.0;
373                for( k=0; k<26; k++ )
374                        scarr[l] += n_dis_consweight_multi[k][l] * cpmx1[k][i1];
375//              scarr[l] *= consweight_multi;
376        }
377        for( j=0; j<lgth2; j++ )
378        {
379                match[j] = 0.0;
380                for( k=0; cpmxpdn[k][j]>-1; k++ )
381                        match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
382        } 
383#endif
384}
385
386static void Atracking_localhom_gapmap( float *impwmpt, float *lasthorizontalw, float *lastverticalw, 
387                                                char **seq1, char **seq2, 
388                        char **mseq1, char **mseq2, 
389                        float **cpmx1, float **cpmx2, 
390                        int **ijp, int icyc, int jcyc,
391                                                int *gapmap1, int *gapmap2 )
392{
393        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
394//      char gap[] = "-";
395        char *gap;
396        float wm;
397        gap = newgapstr;
398        lgth1 = strlen( seq1[0] );
399        lgth2 = strlen( seq2[0] );
400
401#if 0
402        for( i=0; i<lgth1; i++ )
403        {
404                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
405        }
406#endif
407 
408        if( outgap == 1 )
409                ;
410        else
411        {
412                wm = lastverticalw[0];
413                for( i=0; i<lgth1; i++ )
414                {
415                        if( lastverticalw[i] >= wm )
416                        {
417                                wm = lastverticalw[i];
418                                iin = i; jin = lgth2-1;
419                                ijp[lgth1][lgth2] = +( lgth1 - i );
420                        }
421                }
422                for( j=0; j<lgth2; j++ )
423                {
424                        if( lasthorizontalw[j] >= wm )
425                        {
426                                wm = lasthorizontalw[j];
427                                iin = lgth1-1; jin = j;
428                                ijp[lgth1][lgth2] = -( lgth2 - j );
429                        }
430                }
431        }
432
433    for( i=0; i<lgth1+1; i++ ) 
434    {
435        ijp[i][0] = i + 1;
436    }
437    for( j=0; j<lgth2+1; j++ ) 
438    {
439        ijp[0][j] = -( j + 1 );
440    }
441
442        for( i=0; i<icyc; i++ )
443        {
444                mseq1[i] += lgth1+lgth2;
445                *mseq1[i] = 0;
446        }
447        for( j=0; j<jcyc; j++ )
448        {
449                mseq2[j] += lgth1+lgth2;
450                *mseq2[j] = 0;
451        }
452        iin = lgth1; jin = lgth2;
453        *impwmpt = 0.0;
454        for( k=0; k<=lgth1+lgth2; k++ ) 
455        {
456                if( ijp[iin][jin] < 0 ) 
457                {
458                        ifi = iin-1; jfi = jin+ijp[iin][jin];
459                }
460                else if( ijp[iin][jin] > 0 )
461                {
462                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
463                }
464                else
465                {
466                        ifi = iin-1; jfi = jin-1;
467                }
468                l = iin - ifi;
469                while( --l ) 
470                {
471                        for( i=0; i<icyc; i++ )
472                                *--mseq1[i] = seq1[i][ifi+l];
473                        for( j=0; j<jcyc; j++ ) 
474                                *--mseq2[j] = *gap;
475                        k++;
476                }
477                l= jin - jfi;
478                while( --l )
479                {
480                        for( i=0; i<icyc; i++ ) 
481                                *--mseq1[i] = *gap;
482                        for( j=0; j<jcyc; j++ ) 
483                                *--mseq2[j] = seq2[j][jfi+l];
484                        k++;
485                }
486                if( iin != lgth1 && jin != lgth2 ) // ??
487                {
488                        *impwmpt += imp_match_out_scQ( gapmap1[iin], gapmap2[jin] );
489//                      fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
490                }
491                if( iin <= 0 || jin <= 0 ) break;
492                for( i=0; i<icyc; i++ ) 
493                        *--mseq1[i] = seq1[i][ifi];
494                for( j=0; j<jcyc; j++ ) 
495                        *--mseq2[j] = seq2[j][jfi];
496                k++;
497                iin = ifi; jin = jfi;
498        }
499}
500
501#if 0
502static void Atracking_localhom_gapmap_bk( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
503                                                char **seq1, char **seq2,
504                        char **mseq1, char **mseq2,
505                        float **cpmx1, float **cpmx2,
506                        int **ijp, int icyc, int jcyc,
507                                                int *gapmap1, int *gapmap2 )
508{
509        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
510        float wm;
511        char *gaptable1, *gt1bk;
512        char *gaptable2, *gt2bk;
513        lgth1 = strlen( seq1[0] );
514        lgth2 = strlen( seq2[0] );
515        gt1bk = AllocateCharVec( lgth1+lgth2+1 );
516        gt2bk = AllocateCharVec( lgth1+lgth2+1 );
517
518#if 0
519        for( i=0; i<lgth1; i++ )
520        {
521                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
522        }
523#endif
524 
525        if( outgap == 1 )
526                ;
527        else
528        {
529                wm = lastverticalw[0];
530                for( i=0; i<lgth1; i++ )
531                {
532                        if( lastverticalw[i] >= wm )
533                        {
534                                wm = lastverticalw[i];
535                                iin = i; jin = lgth2-1;
536                                ijp[lgth1][lgth2] = +( lgth1 - i );
537                        }
538                }
539                for( j=0; j<lgth2; j++ )
540                {
541                        if( lasthorizontalw[j] >= wm )
542                        {
543                                wm = lasthorizontalw[j];
544                                iin = lgth1-1; jin = j;
545                                ijp[lgth1][lgth2] = -( lgth2 - j );
546                        }
547                }
548        }
549
550    for( i=0; i<lgth1+1; i++ )
551    {
552        ijp[i][0] = i + 1;
553    }
554    for( j=0; j<lgth2+1; j++ )
555    {
556        ijp[0][j] = -( j + 1 );
557    }
558
559        gaptable1 = gt1bk + lgth1+lgth2;
560        *gaptable1 = 0;
561        gaptable2 = gt2bk + lgth1+lgth2;
562        *gaptable2 = 0;
563
564        iin = lgth1; jin = lgth2;
565        *impwmpt = 0.0;
566        for( k=0; k<=lgth1+lgth2; k++ )
567        {
568                if( ijp[iin][jin] < 0 )
569                {
570                        ifi = iin-1; jfi = jin+ijp[iin][jin];
571                }
572                else if( ijp[iin][jin] > 0 )
573                {
574                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
575                }
576                else
577                {
578                        ifi = iin-1; jfi = jin-1;
579                }
580                l = iin - ifi;
581                while( --l )
582                {
583                        *--gaptable1 = 'o';
584                        *--gaptable2 = '-';
585                        k++;
586                }
587                l= jin - jfi;
588                while( --l )
589                {
590                        *--gaptable1 = '-';
591                        *--gaptable2 = 'o';
592                        k++;
593                }
594                if( iin == lgth1 || jin == lgth2 )
595                        ;
596                else
597                {
598                        *impwmpt += imp_match_out_scQ( gapmap1[iin], gapmap2[jin] );
599
600//              fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
601                }
602                if( iin <= 0 || jin <= 0 ) break;
603                *--gaptable1 = '-';
604                *--gaptable2 = '-';
605                k++;
606                iin = ifi; jin = jfi;
607        }
608        for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );
609        for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );
610
611        fprintf( stderr, "mseq1[0] = %s\n", mseq1[0] );
612        fprintf( stderr, "mseq2[0] = %s\n", mseq2[0] );
613
614
615        free( gt1bk );
616        free( gt2bk );
617}
618
619#endif
620
621static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw, 
622                                                char **seq1, char **seq2, 
623                        char **mseq1, char **mseq2, 
624                        float **cpmx1, float **cpmx2, 
625                        int **ijp, int icyc, int jcyc )
626{
627        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
628        float wm;
629        char *gaptable1, *gt1bk;
630        char *gaptable2, *gt2bk;
631        lgth1 = strlen( seq1[0] );
632        lgth2 = strlen( seq2[0] );
633        gt1bk = AllocateCharVec( lgth1+lgth2+1 );
634        gt2bk = AllocateCharVec( lgth1+lgth2+1 );
635
636#if 0
637        for( i=0; i<lgth1; i++ )
638        {
639                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
640        }
641#endif
642 
643        if( outgap == 1 )
644                ;
645        else
646        {
647                wm = lastverticalw[0];
648                for( i=0; i<lgth1; i++ )
649                {
650                        if( lastverticalw[i] >= wm )
651                        {
652                                wm = lastverticalw[i];
653                                iin = i; jin = lgth2-1;
654                                ijp[lgth1][lgth2] = +( lgth1 - i );
655                        }
656                }
657                for( j=0; j<lgth2; j++ )
658                {
659                        if( lasthorizontalw[j] >= wm )
660                        {
661                                wm = lasthorizontalw[j];
662                                iin = lgth1-1; jin = j;
663                                ijp[lgth1][lgth2] = -( lgth2 - j );
664                        }
665                }
666        }
667
668    for( i=0; i<lgth1+1; i++ ) 
669    {
670        ijp[i][0] = i + 1;
671    }
672    for( j=0; j<lgth2+1; j++ ) 
673    {
674        ijp[0][j] = -( j + 1 );
675    }
676
677        gaptable1 = gt1bk + lgth1+lgth2;
678        *gaptable1 = 0;
679        gaptable2 = gt2bk + lgth1+lgth2;
680        *gaptable2 = 0;
681
682        iin = lgth1; jin = lgth2;
683        *impwmpt = 0.0;
684        for( k=0; k<=lgth1+lgth2; k++ ) 
685        {
686                if( ijp[iin][jin] < 0 ) 
687                {
688                        ifi = iin-1; jfi = jin+ijp[iin][jin];
689                }
690                else if( ijp[iin][jin] > 0 )
691                {
692                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
693                }
694                else
695                {
696                        ifi = iin-1; jfi = jin-1;
697                }
698                l = iin - ifi;
699                while( --l ) 
700                {
701                        *--gaptable1 = 'o';
702                        *--gaptable2 = '-';
703                        k++;
704                }
705                l= jin - jfi;
706                while( --l )
707                {
708                        *--gaptable1 = '-';
709                        *--gaptable2 = 'o';
710                        k++;
711                }
712                if( iin == lgth1 || jin == lgth2 )
713                        ;
714                else
715                {
716                        *impwmpt += imp_match_out_scQ( iin, jin );
717
718//              fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
719                }
720                if( iin <= 0 || jin <= 0 ) break;
721                *--gaptable1 = 'o';
722                *--gaptable2 = 'o';
723                k++;
724                iin = ifi; jin = jfi;
725        }
726
727        for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );
728        for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );
729
730        free( gt1bk );
731        free( gt2bk );
732}
733
734
735static float Atracking( float *lasthorizontalw, float *lastverticalw, 
736                                                char **seq1, char **seq2, 
737                        char **mseq1, char **mseq2, 
738                        float **cpmx1, float **cpmx2, 
739                        int **ijp, int icyc, int jcyc )
740{
741        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
742        float wm;
743        char *gaptable1, *gt1bk;
744        char *gaptable2, *gt2bk;
745        lgth1 = strlen( seq1[0] );
746        lgth2 = strlen( seq2[0] );
747
748        gt1bk = AllocateCharVec( lgth1+lgth2+1 );
749        gt2bk = AllocateCharVec( lgth1+lgth2+1 );
750
751#if 0
752        for( i=0; i<lgth1; i++ )
753        {
754                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
755        }
756#endif
757 
758        if( outgap == 1 )
759                ;
760        else
761        {
762                wm = lastverticalw[0];
763                for( i=0; i<lgth1; i++ )
764                {
765                        if( lastverticalw[i] >= wm )
766                        {
767                                wm = lastverticalw[i];
768                                iin = i; jin = lgth2-1;
769                                ijp[lgth1][lgth2] = +( lgth1 - i );
770                        }
771                }
772                for( j=0; j<lgth2; j++ )
773                {
774                        if( lasthorizontalw[j] >= wm )
775                        {
776                                wm = lasthorizontalw[j];
777                                iin = lgth1-1; jin = j;
778                                ijp[lgth1][lgth2] = -( lgth2 - j );
779                        }
780                }
781        }
782
783    for( i=0; i<lgth1+1; i++ ) 
784    {
785        ijp[i][0] = i + 1;
786    }
787    for( j=0; j<lgth2+1; j++ ) 
788    {
789        ijp[0][j] = -( j + 1 );
790    }
791
792        gaptable1 = gt1bk + lgth1+lgth2;
793        *gaptable1 = 0;
794        gaptable2 = gt2bk + lgth1+lgth2;
795        *gaptable2 = 0;
796
797        iin = lgth1; jin = lgth2;
798        for( k=0; k<=lgth1+lgth2; k++ ) 
799        {
800                if( ijp[iin][jin] < 0 ) 
801                {
802                        ifi = iin-1; jfi = jin+ijp[iin][jin];
803                }
804                else if( ijp[iin][jin] > 0 )
805                {
806                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
807                }
808                else
809                {
810                        ifi = iin-1; jfi = jin-1;
811                }
812                l = iin - ifi;
813                while( --l ) 
814                {
815                        *--gaptable1 = 'o';
816                        *--gaptable2 = '-';
817                        k++;
818                }
819                l= jin - jfi;
820                while( --l )
821                {
822                        *--gaptable1 = '-';
823                        *--gaptable2 = 'o';
824                        k++;
825                }
826                if( iin <= 0 || jin <= 0 ) break;
827                *--gaptable1 = 'o';
828                *--gaptable2 = 'o';
829                k++;
830                iin = ifi; jin = jfi;
831        }
832
833        for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );
834        for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );
835
836        free( gt1bk );
837        free( gt2bk );
838
839        return( 0.0 );
840}
841
842float Q__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, char *sgap1, char *sgap2, char *egap1, char *egap2 )
843/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
844{
845//      int k;
846        register int i, j;
847        int lasti, lastj;      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
848        int lgth1, lgth2;
849        int resultlen;
850        float wm = 0.0;   /* int ?????? */
851        float g;
852        float *currentw, *previousw;
853//      float fpenalty = (float)penalty;
854#if USE_PENALTY_EX
855        float fpenalty_ex = (float)penalty_ex;
856        fprintf( stderr, "fpenalty_ex = %f\n", fpenalty_ex );
857#endif
858#if 1
859        float *wtmp;
860        int *ijppt;
861        float *mjpt, *prept, *curpt;
862        int *mpjpt;
863#endif
864        static TLS float mi, *m;
865        static TLS int **ijp;
866        static TLS int mpi, *mp;
867        static TLS float *w1, *w2;
868        static TLS float *match;
869        static TLS float *initverticalw;    /* kufuu sureba iranai */
870        static TLS float *lastverticalw;    /* kufuu sureba iranai */
871        static TLS char **mseq1;
872        static TLS char **mseq2;
873        static TLS char **mseq;
874        static TLS float *digf1;
875        static TLS float *digf2;
876        static TLS float *diaf1;
877        static TLS float *diaf2;
878        static TLS float *gapz1;
879        static TLS float *gapz2;
880        static TLS float *gapf1;
881        static TLS float *gapf2;
882        static TLS float *ogcp1g;
883        static TLS float *ogcp2g;
884        static TLS float *fgcp1g;
885        static TLS float *fgcp2g;
886        static TLS float *og_h_dg_n1_p;
887        static TLS float *og_h_dg_n2_p;
888        static TLS float *fg_h_dg_n1_p;
889        static TLS float *fg_h_dg_n2_p;
890        static TLS float *og_t_fg_h_dg_n1_p;
891        static TLS float *og_t_fg_h_dg_n2_p;
892        static TLS float *fg_t_og_h_dg_n1_p;
893        static TLS float *fg_t_og_h_dg_n2_p;
894        static TLS float *gapz_n1;
895        static TLS float *gapz_n2;
896        static TLS float **cpmx1;
897        static TLS float **cpmx2;
898        static TLS int **intwork;
899        static TLS float **floatwork;
900        static TLS int orlgth1 = 0, orlgth2 = 0;
901        float tmppenal;
902        float *fg_t_og_h_dg_n2_p_pt;
903        float *og_t_fg_h_dg_n2_p_pt;
904        float *og_h_dg_n2_p_pt;
905        float *fg_h_dg_n2_p_pt;
906        float *gapz_n2_pt0;
907        float *gapz_n2_pt1;
908        float *fgcp2pt;
909        float *ogcp2pt;
910        float fg_t_og_h_dg_n1_p_va;
911        float og_t_fg_h_dg_n1_p_va;
912        float og_h_dg_n1_p_va;
913        float fg_h_dg_n1_p_va;
914        float gapz_n1_va0;
915        float gapz_n1_va1;
916        float fgcp1va;
917        float ogcp1va;
918        float kyokaipenal;
919#if 1
920        float fpenalty = (float)penalty;
921#else
922        float fpenalty;
923        if( RNAscoremtx != 'r' ) fpenalty = (float)penalty;
924        else fpenalty = (float)penalty * 10;
925#endif
926
927#if 0
928        fprintf( stderr, "####  seq1[0] = %s\n", seq1[0] );
929        fprintf( stderr, "####  seq2[0] = %s\n", seq2[0] );
930#endif
931
932
933
934        if( orlgth1 == 0 )
935        {
936                mseq1 = AllocateCharMtx( njob, 0 );
937                mseq2 = AllocateCharMtx( njob, 0 );
938        }
939
940
941        lgth1 = strlen( seq1[0] );
942        lgth2 = strlen( seq2[0] );
943#if 0
944        if( lgth1 == 0 || lgth2 == 0 )
945        {
946                fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
947        }
948#endif
949
950        if( lgth1 > orlgth1 || lgth2 > orlgth2 )
951        {
952                int ll1, ll2;
953
954                if( orlgth1 > 0 && orlgth2 > 0 )
955                {
956                        FreeFloatVec( w1 );
957                        FreeFloatVec( w2 );
958                        FreeFloatVec( match );
959                        FreeFloatVec( initverticalw );
960                        FreeFloatVec( lastverticalw );
961
962                        FreeFloatVec( m );
963                        FreeIntVec( mp );
964
965                        FreeCharMtx( mseq );
966
967                        FreeFloatVec( digf1 );
968                        FreeFloatVec( digf2 );
969                        FreeFloatVec( diaf1 );
970                        FreeFloatVec( diaf2 );
971                        FreeFloatVec( gapz1 );
972                        FreeFloatVec( gapz2 );
973                        FreeFloatVec( gapf1 );
974                        FreeFloatVec( gapf2 );
975                        FreeFloatVec( ogcp1g );
976                        FreeFloatVec( ogcp2g );
977                        FreeFloatVec( fgcp1g );
978                        FreeFloatVec( fgcp2g );
979                        FreeFloatVec( og_h_dg_n1_p );
980                        FreeFloatVec( og_h_dg_n2_p );
981                        FreeFloatVec( fg_h_dg_n1_p );
982                        FreeFloatVec( fg_h_dg_n2_p );
983                        FreeFloatVec( og_t_fg_h_dg_n1_p );
984                        FreeFloatVec( og_t_fg_h_dg_n2_p );
985                        FreeFloatVec( fg_t_og_h_dg_n1_p );
986                        FreeFloatVec( fg_t_og_h_dg_n2_p );
987                        FreeFloatVec( gapz_n1 );
988                        FreeFloatVec( gapz_n2 );
989
990                        FreeFloatMtx( cpmx1 );
991                        FreeFloatMtx( cpmx2 );
992
993                        FreeFloatMtx( floatwork );
994                        FreeIntMtx( intwork );
995                }
996
997                ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
998                ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
999
1000#if DEBUG
1001                fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1002#endif
1003
1004                w1 = AllocateFloatVec( ll2+2 );
1005                w2 = AllocateFloatVec( ll2+2 );
1006                match = AllocateFloatVec( ll2+2 );
1007
1008                initverticalw = AllocateFloatVec( ll1+2 );
1009                lastverticalw = AllocateFloatVec( ll1+2 );
1010
1011                m = AllocateFloatVec( ll2+2 );
1012                mp = AllocateIntVec( ll2+2 );
1013
1014                mseq = AllocateCharMtx( njob, ll1+ll2 );
1015
1016                digf1 = AllocateFloatVec( ll1+2 );
1017                digf2 = AllocateFloatVec( ll2+2 );
1018                diaf1 = AllocateFloatVec( ll1+2 );
1019                diaf2 = AllocateFloatVec( ll2+2 );
1020                gapz1 = AllocateFloatVec( ll1+2 );
1021                gapz2 = AllocateFloatVec( ll2+2 );
1022                gapf1 = AllocateFloatVec( ll1+2 );
1023                gapf2 = AllocateFloatVec( ll2+2 );
1024                ogcp1g = AllocateFloatVec( ll1+2 );
1025                ogcp2g = AllocateFloatVec( ll2+2 );
1026                fgcp1g = AllocateFloatVec( ll1+2 );
1027                fgcp2g = AllocateFloatVec( ll2+2 );
1028                og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
1029                og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
1030                fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
1031                fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
1032                og_t_fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
1033                og_t_fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
1034                fg_t_og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
1035                fg_t_og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
1036                gapz_n1 = AllocateFloatVec( ll1+2 );
1037                gapz_n2 = AllocateFloatVec( ll2+2 );
1038
1039                cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1040                cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1041
1042#if FASTMATCHCALC
1043                floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 ); 
1044                intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 ); 
1045#else
1046                floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 ); 
1047                intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 ); 
1048#endif
1049
1050#if DEBUG
1051                fprintf( stderr, "succeeded\n" );
1052#endif
1053
1054                orlgth1 = ll1 - 100;
1055                orlgth2 = ll2 - 100;
1056        }
1057
1058
1059        for( i=0; i<icyc; i++ )
1060        {
1061                mseq1[i] = mseq[i];
1062                seq1[i][lgth1] = 0;
1063        }
1064        for( j=0; j<jcyc; j++ )
1065        {
1066                mseq2[j] = mseq[icyc+j];
1067                seq2[j][lgth2] = 0;
1068        }
1069
1070
1071        if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1072        {
1073                int ll1, ll2;
1074
1075                if( commonAlloc1 && commonAlloc2 )
1076                {
1077                        FreeIntMtx( commonIP );
1078                }
1079
1080                ll1 = MAX( orlgth1, commonAlloc1 );
1081                ll2 = MAX( orlgth2, commonAlloc2 );
1082
1083#if DEBUG
1084                fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1085#endif
1086
1087                commonIP = AllocateIntMtx( ll1+10, ll2+10 );
1088
1089#if DEBUG
1090                fprintf( stderr, "succeeded\n\n" );
1091#endif
1092
1093                commonAlloc1 = ll1;
1094                commonAlloc2 = ll2;
1095        }
1096        ijp = commonIP;
1097
1098#if 0
1099        {
1100                float t = 0.0;
1101                for( i=0; i<icyc; i++ )
1102                        t += eff1[i];
1103        fprintf( stderr, "## totaleff = %f\n", t );
1104        }
1105#endif
1106
1107        cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1108        cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1109
1110        if( sgap1 )
1111        {
1112                new_OpeningGapCount_zure( ogcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
1113                new_OpeningGapCount_zure( ogcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
1114                new_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
1115                new_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
1116                getdigapfreq_part( digf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
1117                getdigapfreq_part( digf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
1118                getdiaminofreq_part( diaf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
1119                getdiaminofreq_part( diaf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
1120                getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
1121                getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
1122                getgapfreq_zure_part( gapz1, icyc, seq1, eff1, lgth1, sgap1 );
1123                getgapfreq_zure_part( gapz2, jcyc, seq2, eff2, lgth2, sgap1 );
1124        }
1125        else
1126        {
1127                st_OpeningGapCount( ogcp1g, icyc, seq1, eff1, lgth1 );
1128                st_OpeningGapCount( ogcp2g, jcyc, seq2, eff2, lgth2 );
1129                st_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1 );
1130                st_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2 );
1131                getdigapfreq_st( digf1, icyc, seq1, eff1, lgth1 );
1132                getdigapfreq_st( digf2, jcyc, seq2, eff2, lgth2 );
1133                getdiaminofreq_x( diaf1, icyc, seq1, eff1, lgth1 );
1134                getdiaminofreq_x( diaf2, jcyc, seq2, eff2, lgth2 );
1135                getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
1136                getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
1137                getgapfreq_zure( gapz1, icyc, seq1, eff1, lgth1 );
1138                getgapfreq_zure( gapz2, jcyc, seq2, eff2, lgth2 );
1139        }
1140
1141#if 1
1142        lastj = lgth2+2;
1143        for( i=0; i<lastj; i++ )
1144        {
1145                og_h_dg_n2_p[i] = ( 1.0-ogcp2g[i]-digf2[i] ) * fpenalty * 0.5;
1146                fg_h_dg_n2_p[i] = ( 1.0-fgcp2g[i]-digf2[i] ) * fpenalty * 0.5;
1147                og_t_fg_h_dg_n2_p[i] = (1.0-ogcp2g[i]+fgcp2g[i]-digf2[i]) * 0.5 * fpenalty;
1148                fg_t_og_h_dg_n2_p[i] = (1.0-fgcp2g[i]+ogcp2g[i]-digf2[i]) * 0.5 * fpenalty;
1149                gapz_n2[i] = (1.0-gapz2[i]);
1150        }
1151        lastj = lgth1+2;
1152        for( i=0; i<lastj; i++ )
1153        {
1154                og_h_dg_n1_p[i] = ( 1.0-ogcp1g[i]-digf1[i] ) * fpenalty * 0.5;
1155                fg_h_dg_n1_p[i] = ( 1.0-fgcp1g[i]-digf1[i] ) * fpenalty * 0.5;
1156                og_t_fg_h_dg_n1_p[i] = (1.0-ogcp1g[i]+fgcp1g[i]-digf1[i]) * 0.5 * fpenalty;
1157                fg_t_og_h_dg_n1_p[i] = (1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) * 0.5 * fpenalty;
1158                gapz_n1[i] = (1.0-gapz1[i]);
1159        }
1160#endif
1161
1162
1163
1164#if 0
1165        for( i=0; i<lgth1; i++ )
1166                fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1167#endif
1168
1169        currentw = w1;
1170        previousw = w2;
1171
1172        if( RNAscoremtx != 'r' )
1173                match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1174        else
1175                clearvec( initverticalw, lgth1 );
1176        if( localhom )
1177                imp_match_out_vead_tateQ( initverticalw, 0, lgth1 ); // 060306
1178
1179        if( RNAscoremtx != 'r' )
1180                match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1181        else
1182                clearvec( currentw, lgth2 );
1183        if( localhom )
1184                imp_match_out_veadQ( currentw, 0, lgth2 ); // 060306
1185
1186
1187#if 0 // -> tbfast.c
1188        if( localhom )
1189                imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1190
1191#endif
1192
1193
1194
1195        kyokaipenal = 0.0;
1196        if( outgap == 1 )
1197        {
1198                g = 0.0;
1199
1200                g += ogcp1g[0] * og_h_dg_n2_p[0];
1201//              g += ogcp1g[0] * ( 1.0-ogcp2g[0]-digf2[0] ) * fpenalty * 0.5;
1202//              if( g ) fprintf( stderr, "init-match penal1=%f, %c-%c\n", g, seq1[0][0], seq2[0][0] );
1203
1204                g += ogcp2g[0] * og_h_dg_n1_p[0];
1205//              g += ogcp2g[0] * ( 1.0-ogcp1g[0]-digf1[0] ) * fpenalty * 0.5;
1206//              if( g ) fprintf( stderr, "init-match penal2=%f, %c-%c\n", g, seq1[0][0], seq2[0][0] );
1207
1208                g += fgcp1g[0] * fg_h_dg_n2_p[0];
1209//              g += fgcp1g[0] * ( 1.0-fgcp2g[0]-digf2[0] ) * fpenalty * 0.5;
1210//              if( g ) fprintf( stderr, "match penal1=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1211
1212                g += fgcp2g[0] * fg_h_dg_n1_p[0];
1213//              g += fgcp2g[0] * ( 1.0-fgcp1g[0]-digf1[0] ) * fpenalty * 0.5;
1214//              if( g ) fprintf( stderr, "match penal2=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1215
1216                kyokaipenal = g;
1217                initverticalw[0] += g;
1218                currentw[0] += g;
1219
1220                for( i=1; i<lgth1+1; i++ )
1221                {
1222                        tmppenal = gapz_n2[0]*og_t_fg_h_dg_n1_p[0];
1223//                      tmppenal = ( (1.0-gapz2[0])*(1.0-ogcp1g[0]+fgcp1g[0]-digf1[0]) ) * 0.5 * fpenalty; // mada
1224                        initverticalw[i] += tmppenal;
1225
1226                        tmppenal = gapz_n2[1]*fg_t_og_h_dg_n1_p[i];
1227//                      tmppenal = ( (1.0-gapz2[1])*(1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada
1228                        initverticalw[i] += tmppenal;
1229
1230                }
1231                for( j=1; j<lgth2+1; j++ )
1232                {
1233                        tmppenal = gapz_n1[0]*og_t_fg_h_dg_n2_p[0];
1234//                      tmppenal = ( (1.0-gapz1[0])*(1.0-ogcp2g[0]+fgcp2g[0]-digf2[0]) ) * 0.5 * fpenalty; // mada
1235                        currentw[j] += tmppenal;
1236
1237                        tmppenal = gapz_n1[1]*fg_t_og_h_dg_n2_p[j];
1238//                      tmppenal = ( (1.0-gapz1[1])*(1.0-fgcp2g[j]+ogcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada
1239                        currentw[j] += tmppenal;
1240                }
1241        }
1242#if OUTGAP0TRY
1243        else
1244        {
1245                for( j=1; j<lgth2+1; j++ )
1246                        currentw[j] -= offset * j / 2.0;
1247                for( i=1; i<lgth1+1; i++ )
1248                        initverticalw[i] -= offset * i / 2.0;
1249        }
1250#endif
1251
1252        m[0] = 0.0; // iranai
1253        for( j=1; j<lgth2+1; ++j ) 
1254        {
1255                mp[j] = 0;
1256                m[j] = currentw[j-1] + 10000 * fpenalty; //iinoka?
1257        }
1258        if( lgth2 == 0 )
1259                lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari
1260        else
1261                lastverticalw[0] = currentw[lgth2-1];
1262
1263        if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1264
1265#if XXXXXXX
1266fprintf( stderr, "currentw = \n" );
1267for( i=0; i<lgth1+1; i++ )
1268{
1269        fprintf( stderr, "%5.2f ", currentw[i] );
1270}
1271fprintf( stderr, "\n" );
1272fprintf( stderr, "initverticalw = \n" );
1273for( i=0; i<lgth2+1; i++ )
1274{
1275        fprintf( stderr, "%5.2f ", initverticalw[i] );
1276}
1277fprintf( stderr, "\n" );
1278fprintf( stderr, "fcgp\n" );
1279for( i=0; i<lgth1; i++ ) 
1280        fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1281for( i=0; i<lgth2; i++ ) 
1282        fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1283#endif
1284
1285        for( i=1; i<lasti; i++ )
1286        {
1287                wtmp = previousw; 
1288                previousw = currentw;
1289                currentw = wtmp;
1290
1291                previousw[0] = initverticalw[i-1];
1292
1293                if( RNAscoremtx != 'r' )
1294                        match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1295                else
1296                        clearvec( currentw, lgth2 );
1297#if XXXXXXX
1298fprintf( stderr, "\n" );
1299fprintf( stderr, "i=%d\n", i );
1300fprintf( stderr, "currentw = \n" );
1301for( j=0; j<lgth2; j++ )
1302{
1303        fprintf( stderr, "%5.2f ", currentw[j] );
1304}
1305fprintf( stderr, "\n" );
1306#endif
1307                if( localhom )
1308                {
1309//                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1310#if  0
1311                        imp_match_out_veadQ( currentw, i, lgth2 );
1312#else
1313                        imp_match_out_veadQ( currentw, i, lgth2 );
1314#endif
1315                }
1316#if XXXXXXX
1317fprintf( stderr, "\n" );
1318fprintf( stderr, "i=%d\n", i );
1319fprintf( stderr, "currentw = \n" );
1320for( j=0; j<lgth2; j++ )
1321{
1322        fprintf( stderr, "%5.2f ", currentw[j] );
1323}
1324fprintf( stderr, "\n" );
1325#endif
1326                currentw[0] = initverticalw[i];
1327
1328                mpi = 0;
1329                mi = previousw[0] + 10000 * fpenalty;
1330
1331                ijppt = ijp[i] + 1;
1332                mjpt = m + 1;
1333                prept = previousw;
1334                curpt = currentw + 1;
1335                mpjpt = mp + 1;
1336                fg_t_og_h_dg_n2_p_pt = fg_t_og_h_dg_n2_p + 1;
1337                og_t_fg_h_dg_n2_p_pt = og_t_fg_h_dg_n2_p + 1;
1338                og_h_dg_n2_p_pt = og_h_dg_n2_p + 1;
1339                fg_h_dg_n2_p_pt = fg_h_dg_n2_p + 1;
1340                gapz_n2_pt0 = gapz_n2 + 1;
1341                gapz_n2_pt1 = gapz_n2 + 2;
1342                fgcp2pt = fgcp2g + 1;
1343                ogcp2pt = ogcp2g + 1;
1344
1345                fg_t_og_h_dg_n1_p_va = fg_t_og_h_dg_n1_p[i];
1346                og_t_fg_h_dg_n1_p_va = og_t_fg_h_dg_n1_p[i];
1347                og_h_dg_n1_p_va = og_h_dg_n1_p[i];
1348                fg_h_dg_n1_p_va = fg_h_dg_n1_p[i];
1349                gapz_n1_va0 = gapz_n1[i];
1350                gapz_n1_va1 = gapz_n1[i+1];
1351                fgcp1va = fgcp1g[i];
1352                ogcp1va = ogcp1g[i];
1353
1354                lastj = lgth2+1;
1355                for( j=1; j<lastj; j++ )
1356                {
1357                        wm = *prept;
1358
1359                        g = ogcp1va * *og_h_dg_n2_p_pt;
1360//                      g = ogcp1g[i] * og_h_dg_n2_p[j];
1361//                      g = ogcp1g[i] * ( 1.0-ogcp2g[j]-digf2[j] ) * fpenalty * 0.5;
1362//                      if( g && i==j ) fprintf( stderr, "match penal1=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1363                        wm += g;
1364
1365                        g = *ogcp2pt * og_h_dg_n1_p_va;
1366//                      g = ogcp2g[j] * og_h_dg_n1_p[i];
1367//                      g = ogcp2g[j] * ( 1.0-ogcp1g[i]-digf1[i] ) * fpenalty * 0.5;
1368//                      if( g && i==j ) fprintf( stderr, "match penal2=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1369                        wm += g;
1370
1371                        g = fgcp1va * *fg_h_dg_n2_p_pt;
1372//                      g = fgcp1g[i] * fg_h_dg_n2_p[j];
1373//                      g = fgcp1g[i] * ( 1.0-fgcp2g[j]-digf2[j] ) * fpenalty * 0.5;
1374//                      if( g && i==j ) fprintf( stderr, "match penal3=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1375                        wm += g;
1376
1377                        g = *fgcp2pt * fg_h_dg_n1_p_va;
1378//                      g = fgcp2g[j] * fg_h_dg_n1_p[i];
1379//                      g = fgcp2g[j] * ( 1.0-fgcp1g[i]-digf1[i] ) * fpenalty * 0.5;
1380//                      if( g && i==j ) fprintf( stderr, "match penal4=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1381                        wm += g;
1382
1383                        *ijppt = 0;
1384
1385#if 0
1386                        fprintf( stderr, "%5.0f->", wm );
1387#endif
1388#if 0
1389                        fprintf( stderr, "%5.0f?", g );
1390#endif
1391                        tmppenal = gapz_n1_va1 * *fg_t_og_h_dg_n2_p_pt;
1392//                      tmppenal = gapz_n1[i+1] * fg_t_og_h_dg_n2_p[j];
1393//                      tmppenal = ( (1.0-gapz1[i+1])*(1.0-fgcp2g[j]+ogcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada
1394                        if( (g=mi+tmppenal) > wm )
1395                        {
1396//                              fprintf( stderr, "jump i start=%f (j=%d, fgcp2g[j]=%f, digf2[j]=%f, diaf2[j]=%f), %c-%c\n", g-mi, j, fgcp2g[j], digf2[j], diaf2[j], seq1[0][i], seq2[0][j] );
1397                                wm = g;
1398                                *ijppt = -( j - mpi );
1399                        }
1400                        tmppenal = gapz_n1_va0 * *og_t_fg_h_dg_n2_p_pt;
1401//                      tmppenal = gapz_n1[i] * og_t_fg_h_dg_n2_p[j];
1402//                      tmppenal = ( (1.0-gapz1[i])*(1.0-ogcp2g[j]+fgcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada
1403                        if( (g=*prept+tmppenal) >= mi )
1404                        {
1405//                              fprintf( stderr, "jump i end=%f, %c-%c\n", g-*prept, seq1[0][i-1], seq2[0][j-1] );
1406                                mi = g;
1407                                mpi = j-1;
1408                        }
1409
1410#if USE_PENALTY_EX
1411            mi += fpenalty_ex;
1412#endif
1413
1414#if 0 
1415                        fprintf( stderr, "%5.0f?", g );
1416#endif
1417                        tmppenal = *gapz_n2_pt1 * fg_t_og_h_dg_n1_p_va;
1418//                      tmppenal = gapz_n2[j+1] * fg_t_og_h_dg_n1_p[i];
1419//                      tmppenal = ( (1.0-gapz2[j+1])*(1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada
1420                        if( (g=*mjpt+tmppenal) > wm )
1421                        {
1422                                wm = g;
1423                                *ijppt = +( i - *mpjpt );
1424                        }
1425                        tmppenal = *gapz_n2_pt0 * og_t_fg_h_dg_n1_p_va;
1426//                      tmppenal = gapz_n2[j] * og_t_fg_h_dg_n1_p[i];
1427//                      tmppenal = ( (1.0-gapz2[j])*(1.0-ogcp1g[i]+fgcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada
1428                        if( (g=*prept+tmppenal) >= *mjpt )
1429                        {
1430                                *mjpt = g;
1431                                *mpjpt = i-1;
1432                        }
1433#if 0
1434                        fprintf( stderr, "%5.0f ", wm );
1435#endif
1436
1437#if USE_PENALTY_EX
1438            m[j] += fpenalty_ex;
1439#endif
1440
1441
1442
1443
1444
1445                        *curpt++ += wm;
1446                        ijppt++;
1447                        mjpt++;
1448                        prept++;
1449                        mpjpt++;
1450                        fg_t_og_h_dg_n2_p_pt++;
1451                        og_t_fg_h_dg_n2_p_pt++;
1452                        og_h_dg_n2_p_pt++;
1453                        fg_h_dg_n2_p_pt++;
1454                        gapz_n2_pt0++;
1455                        gapz_n2_pt1++;
1456                        fgcp2pt++;
1457                        ogcp2pt++;
1458                }
1459                lastverticalw[i] = currentw[lgth2-1];
1460        }
1461
1462//      fprintf( stderr, "wm = %f\n", wm );
1463
1464#if OUTGAP0TRY
1465        if( !outgap )
1466        {
1467                for( j=1; j<lgth2+1; j++ )
1468                        currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1469                for( i=1; i<lgth1+1; i++ )
1470                        lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1471        }
1472#endif
1473               
1474        /*
1475        fprintf( stderr, "\n" );
1476        for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1477        fprintf( stderr, "#####\n" );
1478        for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1479        fprintf( stderr, "====>" );
1480        for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1481        for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1482        */
1483        if( localhom )
1484        {
1485                Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1486        }
1487        else
1488                Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1489
1490//      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1491
1492        resultlen = strlen( mseq1[0] );
1493        if( alloclen < resultlen || resultlen > N )
1494        {
1495                fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1496                ErrorExit( "LENGTH OVER!\n" );
1497        }
1498
1499
1500        for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1501        for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1502        /*
1503        fprintf( stderr, "\n" );
1504        for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1505        fprintf( stderr, "#####\n" );
1506        for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1507        */
1508
1509//      fprintf( stderr, "wm = %f\n", wm );
1510
1511
1512        return( wm );
1513}
1514
1515float Q__align_gapmap( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, char *sgap1, char *sgap2, char *egap1, char *egap2, int *gapmap1, int *gapmap2 )
1516/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1517{
1518//      int k;
1519        register int i, j;
1520        int lasti, lastj;      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
1521        int lgth1, lgth2;
1522        int resultlen;
1523        float wm = 0.0;   /* int ?????? */
1524        float g;
1525        float *currentw, *previousw;
1526//      float fpenalty = (float)penalty;
1527#if USE_PENALTY_EX
1528        float fpenalty_ex = (float)penalty_ex;
1529        fprintf( stderr, "fpenalty_ex = %f\n", fpenalty_ex );
1530#endif
1531#if 1
1532        float *wtmp;
1533        int *ijppt;
1534        float *mjpt, *prept, *curpt;
1535        int *mpjpt;
1536#endif
1537        static TLS float mi, *m;
1538        static TLS int **ijp;
1539        static TLS int mpi, *mp;
1540        static TLS float *w1, *w2;
1541        static TLS float *match;
1542        static TLS float *initverticalw;    /* kufuu sureba iranai */
1543        static TLS float *lastverticalw;    /* kufuu sureba iranai */
1544        static TLS char **mseq1;
1545        static TLS char **mseq2;
1546        static TLS char **mseq;
1547        static TLS float *digf1;
1548        static TLS float *digf2;
1549        static TLS float *diaf1;
1550        static TLS float *diaf2;
1551        static TLS float *gapz1;
1552        static TLS float *gapz2;
1553        static TLS float *gapf1;
1554        static TLS float *gapf2;
1555        static TLS float *ogcp1g;
1556        static TLS float *ogcp2g;
1557        static TLS float *fgcp1g;
1558        static TLS float *fgcp2g;
1559        static TLS float *og_h_dg_n1_p;
1560        static TLS float *og_h_dg_n2_p;
1561        static TLS float *fg_h_dg_n1_p;
1562        static TLS float *fg_h_dg_n2_p;
1563        static TLS float *og_t_fg_h_dg_n1_p;
1564        static TLS float *og_t_fg_h_dg_n2_p;
1565        static TLS float *fg_t_og_h_dg_n1_p;
1566        static TLS float *fg_t_og_h_dg_n2_p;
1567        static TLS float *gapz_n1;
1568        static TLS float *gapz_n2;
1569        static TLS float **cpmx1;
1570        static TLS float **cpmx2;
1571        static TLS int **intwork;
1572        static TLS float **floatwork;
1573        static TLS int orlgth1 = 0, orlgth2 = 0;
1574        float tmppenal;
1575        float *fg_t_og_h_dg_n2_p_pt;
1576        float *og_t_fg_h_dg_n2_p_pt;
1577        float *og_h_dg_n2_p_pt;
1578        float *fg_h_dg_n2_p_pt;
1579        float *gapz_n2_pt0;
1580        float *gapz_n2_pt1;
1581        float *fgcp2pt;
1582        float *ogcp2pt;
1583        float fg_t_og_h_dg_n1_p_va;
1584        float og_t_fg_h_dg_n1_p_va;
1585        float og_h_dg_n1_p_va;
1586        float fg_h_dg_n1_p_va;
1587        float gapz_n1_va0;
1588        float gapz_n1_va1;
1589        float fgcp1va;
1590        float ogcp1va;
1591        float kyokaipenal;
1592        float fpenalty = (float)penalty;
1593
1594#if 0
1595        fprintf( stderr, "####  seq1[0] = %s\n", seq1[0] );
1596        fprintf( stderr, "####  seq2[0] = %s\n", seq2[0] );
1597#endif
1598
1599
1600
1601        if( orlgth1 == 0 )
1602        {
1603                mseq1 = AllocateCharMtx( njob, 0 );
1604                mseq2 = AllocateCharMtx( njob, 0 );
1605        }
1606
1607
1608        lgth1 = strlen( seq1[0] );
1609        lgth2 = strlen( seq2[0] );
1610#if 0
1611        if( lgth1 == 0 || lgth2 == 0 )
1612        {
1613                fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
1614        }
1615#endif
1616
1617        if( lgth1 > orlgth1 || lgth2 > orlgth2 )
1618        {
1619                int ll1, ll2;
1620
1621                if( orlgth1 > 0 && orlgth2 > 0 )
1622                {
1623                        FreeFloatVec( w1 );
1624                        FreeFloatVec( w2 );
1625                        FreeFloatVec( match );
1626                        FreeFloatVec( initverticalw );
1627                        FreeFloatVec( lastverticalw );
1628
1629                        FreeFloatVec( m );
1630                        FreeIntVec( mp );
1631
1632                        FreeCharMtx( mseq );
1633
1634                        FreeFloatVec( digf1 );
1635                        FreeFloatVec( digf2 );
1636                        FreeFloatVec( diaf1 );
1637                        FreeFloatVec( diaf2 );
1638                        FreeFloatVec( gapz1 );
1639                        FreeFloatVec( gapz2 );
1640                        FreeFloatVec( gapf1 );
1641                        FreeFloatVec( gapf2 );
1642                        FreeFloatVec( ogcp1g );
1643                        FreeFloatVec( ogcp2g );
1644                        FreeFloatVec( fgcp1g );
1645                        FreeFloatVec( fgcp2g );
1646                        FreeFloatVec( og_h_dg_n1_p );
1647                        FreeFloatVec( og_h_dg_n2_p );
1648                        FreeFloatVec( fg_h_dg_n1_p );
1649                        FreeFloatVec( fg_h_dg_n2_p );
1650                        FreeFloatVec( og_t_fg_h_dg_n1_p );
1651                        FreeFloatVec( og_t_fg_h_dg_n2_p );
1652                        FreeFloatVec( fg_t_og_h_dg_n1_p );
1653                        FreeFloatVec( fg_t_og_h_dg_n2_p );
1654                        FreeFloatVec( gapz_n1 );
1655                        FreeFloatVec( gapz_n2 );
1656
1657                        FreeFloatMtx( cpmx1 );
1658                        FreeFloatMtx( cpmx2 );
1659
1660                        FreeFloatMtx( floatwork );
1661                        FreeIntMtx( intwork );
1662                }
1663
1664                ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
1665                ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
1666
1667#if DEBUG
1668                fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1669#endif
1670
1671                w1 = AllocateFloatVec( ll2+2 );
1672                w2 = AllocateFloatVec( ll2+2 );
1673                match = AllocateFloatVec( ll2+2 );
1674
1675                initverticalw = AllocateFloatVec( ll1+2 );
1676                lastverticalw = AllocateFloatVec( ll1+2 );
1677
1678                m = AllocateFloatVec( ll2+2 );
1679                mp = AllocateIntVec( ll2+2 );
1680
1681                mseq = AllocateCharMtx( njob, ll1+ll2 );
1682
1683                digf1 = AllocateFloatVec( ll1+2 );
1684                digf2 = AllocateFloatVec( ll2+2 );
1685                diaf1 = AllocateFloatVec( ll1+2 );
1686                diaf2 = AllocateFloatVec( ll2+2 );
1687                gapz1 = AllocateFloatVec( ll1+2 );
1688                gapz2 = AllocateFloatVec( ll2+2 );
1689                gapf1 = AllocateFloatVec( ll1+2 );
1690                gapf2 = AllocateFloatVec( ll2+2 );
1691                ogcp1g = AllocateFloatVec( ll1+2 );
1692                ogcp2g = AllocateFloatVec( ll2+2 );
1693                fgcp1g = AllocateFloatVec( ll1+2 );
1694                fgcp2g = AllocateFloatVec( ll2+2 );
1695                og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
1696                og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
1697                fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
1698                fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
1699                og_t_fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
1700                og_t_fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
1701                fg_t_og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
1702                fg_t_og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
1703                gapz_n1 = AllocateFloatVec( ll1+2 );
1704                gapz_n2 = AllocateFloatVec( ll2+2 );
1705
1706                cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1707                cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1708
1709#if FASTMATCHCALC
1710                floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 ); 
1711                intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 ); 
1712#else
1713                floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 ); 
1714                intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 ); 
1715#endif
1716
1717#if DEBUG
1718                fprintf( stderr, "succeeded\n" );
1719#endif
1720
1721                orlgth1 = ll1 - 100;
1722                orlgth2 = ll2 - 100;
1723        }
1724
1725
1726        for( i=0; i<icyc; i++ )
1727        {
1728                mseq1[i] = mseq[i];
1729                seq1[i][lgth1] = 0;
1730        }
1731        for( j=0; j<jcyc; j++ )
1732        {
1733                mseq2[j] = mseq[icyc+j];
1734                seq2[j][lgth2] = 0;
1735        }
1736
1737
1738        if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1739        {
1740                int ll1, ll2;
1741
1742                if( commonAlloc1 && commonAlloc2 )
1743                {
1744                        FreeIntMtx( commonIP );
1745                }
1746
1747                ll1 = MAX( orlgth1, commonAlloc1 );
1748                ll2 = MAX( orlgth2, commonAlloc2 );
1749
1750#if DEBUG
1751                fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1752#endif
1753
1754                commonIP = AllocateIntMtx( ll1+10, ll2+10 );
1755
1756#if DEBUG
1757                fprintf( stderr, "succeeded\n\n" );
1758#endif
1759
1760                commonAlloc1 = ll1;
1761                commonAlloc2 = ll2;
1762        }
1763        ijp = commonIP;
1764
1765#if 0
1766        {
1767                float t = 0.0;
1768                for( i=0; i<icyc; i++ )
1769                        t += eff1[i];
1770        fprintf( stderr, "## totaleff = %f\n", t );
1771        }
1772#endif
1773
1774        cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1775        cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1776
1777        if( sgap1 )
1778        {
1779                new_OpeningGapCount_zure( ogcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
1780                new_OpeningGapCount_zure( ogcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
1781                new_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
1782                new_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
1783                getdigapfreq_part( digf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
1784                getdigapfreq_part( digf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
1785                getdiaminofreq_part( diaf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
1786                getdiaminofreq_part( diaf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
1787                getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
1788                getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
1789                getgapfreq_zure_part( gapz1, icyc, seq1, eff1, lgth1, sgap1 );
1790                getgapfreq_zure_part( gapz2, jcyc, seq2, eff2, lgth2, sgap1 );
1791        }
1792        else
1793        {
1794                st_OpeningGapCount( ogcp1g, icyc, seq1, eff1, lgth1 );
1795                st_OpeningGapCount( ogcp2g, jcyc, seq2, eff2, lgth2 );
1796                st_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1 );
1797                st_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2 );
1798                getdigapfreq_st( digf1, icyc, seq1, eff1, lgth1 );
1799                getdigapfreq_st( digf2, jcyc, seq2, eff2, lgth2 );
1800                getdiaminofreq_x( diaf1, icyc, seq1, eff1, lgth1 );
1801                getdiaminofreq_x( diaf2, jcyc, seq2, eff2, lgth2 );
1802                getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
1803                getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
1804                getgapfreq_zure( gapz1, icyc, seq1, eff1, lgth1 );
1805                getgapfreq_zure( gapz2, jcyc, seq2, eff2, lgth2 );
1806        }
1807
1808#if 1
1809        lastj = lgth2+2;
1810        for( i=0; i<lastj; i++ )
1811        {
1812                og_h_dg_n2_p[i] = ( 1.0-ogcp2g[i]-digf2[i] ) * fpenalty * 0.5;
1813                fg_h_dg_n2_p[i] = ( 1.0-fgcp2g[i]-digf2[i] ) * fpenalty * 0.5;
1814                og_t_fg_h_dg_n2_p[i] = (1.0-ogcp2g[i]+fgcp2g[i]-digf2[i]) * 0.5 * fpenalty;
1815                fg_t_og_h_dg_n2_p[i] = (1.0-fgcp2g[i]+ogcp2g[i]-digf2[i]) * 0.5 * fpenalty;
1816                gapz_n2[i] = (1.0-gapz2[i]);
1817        }
1818        lastj = lgth1+2;
1819        for( i=0; i<lastj; i++ )
1820        {
1821                og_h_dg_n1_p[i] = ( 1.0-ogcp1g[i]-digf1[i] ) * fpenalty * 0.5;
1822                fg_h_dg_n1_p[i] = ( 1.0-fgcp1g[i]-digf1[i] ) * fpenalty * 0.5;
1823                og_t_fg_h_dg_n1_p[i] = (1.0-ogcp1g[i]+fgcp1g[i]-digf1[i]) * 0.5 * fpenalty;
1824                fg_t_og_h_dg_n1_p[i] = (1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) * 0.5 * fpenalty;
1825                gapz_n1[i] = (1.0-gapz1[i]);
1826        }
1827#endif
1828
1829
1830#if 0
1831        for( i=0; i<lgth1; i++ )
1832                fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1833#endif
1834
1835        currentw = w1;
1836        previousw = w2;
1837
1838        if( RNAscoremtx != 'r' )
1839                match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1840        else
1841                clearvec( initverticalw, lgth1 );
1842        if( localhom )
1843                imp_match_out_vead_tateQ_gapmap( initverticalw, gapmap2[0], lgth1, gapmap1 ); // 060306
1844
1845        if( RNAscoremtx != 'r' )
1846                match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1847        else
1848                clearvec( currentw, lgth2 );
1849        if( localhom )
1850                imp_match_out_veadQ_gapmap( currentw, gapmap1[0], lgth2, gapmap2 ); // 060306
1851
1852
1853
1854#if 0 // -> tbfast.c
1855        if( localhom )
1856                imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1857
1858#endif
1859
1860
1861
1862        kyokaipenal = 0.0;
1863        if( outgap == 1 )
1864        {
1865                g = 0.0;
1866
1867                g += ogcp1g[0] * og_h_dg_n2_p[0];
1868//              g += ogcp1g[0] * ( 1.0-ogcp2g[0]-digf2[0] ) * fpenalty * 0.5;
1869//              if( g ) fprintf( stderr, "init-match penal1=%f, %c-%c\n", g, seq1[0][0], seq2[0][0] );
1870
1871                g += ogcp2g[0] * og_h_dg_n1_p[0];
1872//              g += ogcp2g[0] * ( 1.0-ogcp1g[0]-digf1[0] ) * fpenalty * 0.5;
1873//              if( g ) fprintf( stderr, "init-match penal2=%f, %c-%c\n", g, seq1[0][0], seq2[0][0] );
1874
1875                g += fgcp1g[0] * fg_h_dg_n2_p[0];
1876//              g += fgcp1g[0] * ( 1.0-fgcp2g[0]-digf2[0] ) * fpenalty * 0.5;
1877//              if( g ) fprintf( stderr, "match penal1=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1878
1879                g += fgcp2g[0] * fg_h_dg_n1_p[0];
1880//              g += fgcp2g[0] * ( 1.0-fgcp1g[0]-digf1[0] ) * fpenalty * 0.5;
1881//              if( g ) fprintf( stderr, "match penal2=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
1882
1883                kyokaipenal = g;
1884                initverticalw[0] += g;
1885                currentw[0] += g;
1886
1887                for( i=1; i<lgth1+1; i++ )
1888                {
1889                        tmppenal = gapz_n2[0]*og_t_fg_h_dg_n1_p[0];
1890//                      tmppenal = ( (1.0-gapz2[0])*(1.0-ogcp1g[0]+fgcp1g[0]-digf1[0]) ) * 0.5 * fpenalty; // mada
1891                        initverticalw[i] += tmppenal;
1892
1893                        tmppenal = gapz_n2[1]*fg_t_og_h_dg_n1_p[i];
1894//                      tmppenal = ( (1.0-gapz2[1])*(1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada
1895                        initverticalw[i] += tmppenal;
1896
1897                }
1898                for( j=1; j<lgth2+1; j++ )
1899                {
1900                        tmppenal = gapz_n1[0]*og_t_fg_h_dg_n2_p[0];
1901//                      tmppenal = ( (1.0-gapz1[0])*(1.0-ogcp2g[0]+fgcp2g[0]-digf2[0]) ) * 0.5 * fpenalty; // mada
1902                        currentw[j] += tmppenal;
1903
1904                        tmppenal = gapz_n1[1]*fg_t_og_h_dg_n2_p[j];
1905//                      tmppenal = ( (1.0-gapz1[1])*(1.0-fgcp2g[j]+ogcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada
1906                        currentw[j] += tmppenal;
1907                }
1908        }
1909#if OUTGAP0TRY
1910        else
1911        {
1912                for( j=1; j<lgth2+1; j++ )
1913                        currentw[j] -= offset * j / 2.0;
1914                for( i=1; i<lgth1+1; i++ )
1915                        initverticalw[i] -= offset * i / 2.0;
1916        }
1917#endif
1918
1919        m[0] = 0.0; // iranai
1920        for( j=1; j<lgth2+1; ++j ) 
1921        {
1922                mp[j] = 0;
1923                m[j] = currentw[j-1] + 10000 * fpenalty; //iinoka?
1924        }
1925        if( lgth2 == 0 )
1926                lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari
1927        else
1928                lastverticalw[0] = currentw[lgth2-1];
1929
1930        if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1931
1932#if XXXXXXX
1933fprintf( stderr, "currentw = \n" );
1934for( i=0; i<lgth1+1; i++ )
1935{
1936        fprintf( stderr, "%5.2f ", currentw[i] );
1937}
1938fprintf( stderr, "\n" );
1939fprintf( stderr, "initverticalw = \n" );
1940for( i=0; i<lgth2+1; i++ )
1941{
1942        fprintf( stderr, "%5.2f ", initverticalw[i] );
1943}
1944fprintf( stderr, "\n" );
1945fprintf( stderr, "fcgp\n" );
1946for( i=0; i<lgth1; i++ ) 
1947        fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1948for( i=0; i<lgth2; i++ ) 
1949        fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1950#endif
1951
1952        for( i=1; i<lasti; i++ )
1953        {
1954                wtmp = previousw; 
1955                previousw = currentw;
1956                currentw = wtmp;
1957
1958                previousw[0] = initverticalw[i-1];
1959
1960                if( RNAscoremtx != 'r' )
1961                        match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1962                else
1963                        clearvec( currentw, lgth2 );
1964#if XXXXXXX
1965fprintf( stderr, "\n" );
1966fprintf( stderr, "i=%d\n", i );
1967fprintf( stderr, "currentw = \n" );
1968for( j=0; j<lgth2; j++ )
1969{
1970        fprintf( stderr, "%5.2f ", currentw[j] );
1971}
1972fprintf( stderr, "\n" );
1973#endif
1974                if( localhom )
1975                {
1976//                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1977#if  0
1978                        imp_match_out_veadQ( currentw, i, lgth2 );
1979#else
1980                        imp_match_out_veadQ_gapmap( currentw, gapmap1[i], lgth2, gapmap2 );
1981#endif
1982                }
1983#if XXXXXXX
1984fprintf( stderr, "\n" );
1985fprintf( stderr, "i=%d\n", i );
1986fprintf( stderr, "currentw = \n" );
1987for( j=0; j<lgth2; j++ )
1988{
1989        fprintf( stderr, "%5.2f ", currentw[j] );
1990}
1991fprintf( stderr, "\n" );
1992#endif
1993                currentw[0] = initverticalw[i];
1994
1995                mpi = 0;
1996                mi = previousw[0] + 10000 * fpenalty;
1997
1998                ijppt = ijp[i] + 1;
1999                mjpt = m + 1;
2000                prept = previousw;
2001                curpt = currentw + 1;
2002                mpjpt = mp + 1;
2003                fg_t_og_h_dg_n2_p_pt = fg_t_og_h_dg_n2_p + 1;
2004                og_t_fg_h_dg_n2_p_pt = og_t_fg_h_dg_n2_p + 1;
2005                og_h_dg_n2_p_pt = og_h_dg_n2_p + 1;
2006                fg_h_dg_n2_p_pt = fg_h_dg_n2_p + 1;
2007                gapz_n2_pt0 = gapz_n2 + 1;
2008                gapz_n2_pt1 = gapz_n2 + 2;
2009                fgcp2pt = fgcp2g + 1;
2010                ogcp2pt = ogcp2g + 1;
2011
2012                fg_t_og_h_dg_n1_p_va = fg_t_og_h_dg_n1_p[i];
2013                og_t_fg_h_dg_n1_p_va = og_t_fg_h_dg_n1_p[i];
2014                og_h_dg_n1_p_va = og_h_dg_n1_p[i];
2015                fg_h_dg_n1_p_va = fg_h_dg_n1_p[i];
2016                gapz_n1_va0 = gapz_n1[i];
2017                gapz_n1_va1 = gapz_n1[i+1];
2018                fgcp1va = fgcp1g[i];
2019                ogcp1va = ogcp1g[i];
2020
2021                lastj = lgth2+1;
2022                for( j=1; j<lastj; j++ )
2023                {
2024                        wm = *prept;
2025
2026                        g = ogcp1va * *og_h_dg_n2_p_pt;
2027//                      g = ogcp1g[i] * og_h_dg_n2_p[j];
2028//                      g = ogcp1g[i] * ( 1.0-ogcp2g[j]-digf2[j] ) * fpenalty * 0.5;
2029//                      if( g && i==j ) fprintf( stderr, "match penal1=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
2030                        wm += g;
2031
2032                        g = *ogcp2pt * og_h_dg_n1_p_va;
2033//                      g = ogcp2g[j] * og_h_dg_n1_p[i];
2034//                      g = ogcp2g[j] * ( 1.0-ogcp1g[i]-digf1[i] ) * fpenalty * 0.5;
2035//                      if( g && i==j ) fprintf( stderr, "match penal2=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
2036                        wm += g;
2037
2038                        g = fgcp1va * *fg_h_dg_n2_p_pt;
2039//                      g = fgcp1g[i] * fg_h_dg_n2_p[j];
2040//                      g = fgcp1g[i] * ( 1.0-fgcp2g[j]-digf2[j] ) * fpenalty * 0.5;
2041//                      if( g && i==j ) fprintf( stderr, "match penal3=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
2042                        wm += g;
2043
2044                        g = *fgcp2pt * fg_h_dg_n1_p_va;
2045//                      g = fgcp2g[j] * fg_h_dg_n1_p[i];
2046//                      g = fgcp2g[j] * ( 1.0-fgcp1g[i]-digf1[i] ) * fpenalty * 0.5;
2047//                      if( g && i==j ) fprintf( stderr, "match penal4=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );
2048                        wm += g;
2049
2050                        *ijppt = 0;
2051
2052#if 0
2053                        fprintf( stderr, "%5.0f->", wm );
2054#endif
2055#if 0
2056                        fprintf( stderr, "%5.0f?", g );
2057#endif
2058                        tmppenal = gapz_n1_va1 * *fg_t_og_h_dg_n2_p_pt;
2059//                      tmppenal = gapz_n1[i+1] * fg_t_og_h_dg_n2_p[j];
2060//                      tmppenal = ( (1.0-gapz1[i+1])*(1.0-fgcp2g[j]+ogcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada
2061                        if( (g=mi+tmppenal) > wm )
2062                        {
2063//                              fprintf( stderr, "jump i start=%f (j=%d, fgcp2g[j]=%f, digf2[j]=%f, diaf2[j]=%f), %c-%c\n", g-mi, j, fgcp2g[j], digf2[j], diaf2[j], seq1[0][i], seq2[0][j] );
2064                                wm = g;
2065                                *ijppt = -( j - mpi );
2066                        }
2067                        tmppenal = gapz_n1_va0 * *og_t_fg_h_dg_n2_p_pt;
2068//                      tmppenal = gapz_n1[i] * og_t_fg_h_dg_n2_p[j];
2069//                      tmppenal = ( (1.0-gapz1[i])*(1.0-ogcp2g[j]+fgcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada
2070                        if( (g=*prept+tmppenal) >= mi )
2071                        {
2072//                              fprintf( stderr, "jump i end=%f, %c-%c\n", g-*prept, seq1[0][i-1], seq2[0][j-1] );
2073                                mi = g;
2074                                mpi = j-1;
2075                        }
2076
2077#if USE_PENALTY_EX
2078            mi += fpenalty_ex;
2079#endif
2080
2081#if 0 
2082                        fprintf( stderr, "%5.0f?", g );
2083#endif
2084                        tmppenal = *gapz_n2_pt1 * fg_t_og_h_dg_n1_p_va;
2085//                      tmppenal = gapz_n2[j+1] * fg_t_og_h_dg_n1_p[i];
2086//                      tmppenal = ( (1.0-gapz2[j+1])*(1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada
2087                        if( (g=*mjpt+tmppenal) > wm )
2088                        {
2089                                wm = g;
2090                                *ijppt = +( i - *mpjpt );
2091                        }
2092                        tmppenal = *gapz_n2_pt0 * og_t_fg_h_dg_n1_p_va;
2093//                      tmppenal = gapz_n2[j] * og_t_fg_h_dg_n1_p[i];
2094//                      tmppenal = ( (1.0-gapz2[j])*(1.0-ogcp1g[i]+fgcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada
2095                        if( (g=*prept+tmppenal) >= *mjpt )
2096                        {
2097                                *mjpt = g;
2098                                *mpjpt = i-1;
2099                        }
2100#if 0
2101                        fprintf( stderr, "%5.0f ", wm );
2102#endif
2103
2104#if USE_PENALTY_EX
2105            m[j] += fpenalty_ex;
2106#endif
2107
2108
2109
2110
2111
2112                        *curpt++ += wm;
2113                        ijppt++;
2114                        mjpt++;
2115                        prept++;
2116                        mpjpt++;
2117                        fg_t_og_h_dg_n2_p_pt++;
2118                        og_t_fg_h_dg_n2_p_pt++;
2119                        og_h_dg_n2_p_pt++;
2120                        fg_h_dg_n2_p_pt++;
2121                        gapz_n2_pt0++;
2122                        gapz_n2_pt1++;
2123                        fgcp2pt++;
2124                        ogcp2pt++;
2125                }
2126                lastverticalw[i] = currentw[lgth2-1];
2127        }
2128
2129//      fprintf( stderr, "wm = %f\n", wm );
2130
2131#if OUTGAP0TRY
2132        if( !outgap )
2133        {
2134                for( j=1; j<lgth2+1; j++ )
2135                        currentw[j] -= offset * ( lgth2 - j ) / 2.0;
2136                for( i=1; i<lgth1+1; i++ )
2137                        lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
2138        }
2139#endif
2140               
2141        /*
2142        fprintf( stderr, "\n" );
2143        for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
2144        fprintf( stderr, "#####\n" );
2145        for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
2146        fprintf( stderr, "====>" );
2147        for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
2148        for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
2149        */
2150        if( localhom )
2151        {
2152                Atracking_localhom_gapmap( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, gapmap1, gapmap2 );
2153        }
2154        else
2155                Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
2156
2157//      fprintf( stderr, "### impmatch = %f\n", *impmatch );
2158
2159        resultlen = strlen( mseq1[0] );
2160        if( alloclen < resultlen || resultlen > N )
2161        {
2162                fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
2163                ErrorExit( "LENGTH OVER!\n" );
2164        }
2165
2166
2167        for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2168        for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
2169        /*
2170        fprintf( stderr, "\n" );
2171        for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
2172        fprintf( stderr, "#####\n" );
2173        for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
2174        */
2175
2176//      fprintf( stderr, "wm = %f\n", wm );
2177
2178
2179        return( wm );
2180}
2181
Note: See TracBrowser for help on using the repository browser.