source: branches/stable/GDE/MAFFT/mafft-7.055-with-extensions/core/partQalignmm.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: 30.1 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
12static int impalloclen = 0;
13static float **impmtx = NULL;
14float part_imp_match_out_scQ( int i1, int j1 )
15{
16//      fprintf( stderr, "impalloclen = %d\n", impalloclen );
17//      fprintf( stderr, "i1,j1=%d,%d -> impmtx=%f\n", i1, j1, impmtx[i1][j1] );
18        return( impmtx[i1][j1] );
19#if 0
20        if( i1 == l1 || j1 == l2 ) return( 0.0 );
21        return( impmtx[i1+start1][j1+start2] );
22#endif
23}
24static void part_imp_match_out_vead_gapmapQ( float *imp, int i1, int lgth2, int start2, int *gapmap2 )
25{
26#if FASTMACHCALC
27        float *pt = imp;
28        int *gapmappt = gapmap2;
29        while( lgth2-- )
30                *pt++ += impmtx[i1][start2+*gapmappt++];
31#else
32        int j;
33        for( j=0; j<lgth2; j++ )
34        {
35                imp[j] += impmtx[i1][start2+gapmap2[j]];
36        }
37#endif
38}
39
40static void part_imp_match_out_vead_tate_gapmapQ( float *imp, int j1, int lgth1, int start1, int *gapmap1 )
41{
42#if FASTMACHCALC
43        float *pt = imp;
44        int *gapmappt = gapmap1;
45        while( lgth1-- )
46                *pt++ = impmtx[start1+*gapmappt++][j1];
47#else
48        int i;
49        for( i=0; i<lgth1; i++ )
50        {
51                imp[i] += impmtx[start1+gapmap1[i]][j1];
52        }
53#endif
54}
55
56void part_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 )
57{
58        int i, j, k1, k2, tmpint, start1, start2, end1, end2;
59        double effij, effijx; 
60        char *pt, *pt1, *pt2;
61        LocalHom *tmpptr;
62
63        if( impalloclen <= lgth1 + 2 || impalloclen <= lgth2 + 2 )
64        {
65                if( impmtx ) FreeFloatMtx( impmtx );
66                impalloclen = MAX( lgth1, lgth2 ) + 2;
67                impmtx = AllocateFloatMtx( impalloclen+100, impalloclen+100 );
68        }
69
70
71#if 0
72        fprintf( stderr, "eff1 in _init_strict = \n" );
73        for( i=0; i<clus1; i++ )
74                fprintf( stderr, "eff1[] = %f\n", eff1[i] );
75        for( i=0; i<clus2; i++ )
76                fprintf( stderr, "eff2[] = %f\n", eff2[i] );
77#endif
78
79        for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
80                impmtx[i][j] = 0.0;
81        effijx = 1.0 * fastathreshold;
82        for( i=0; i<clus1; i++ )
83        {
84                for( j=0; j<clus2; j++ )
85                {
86                        effij = eff1[i] * eff2[j] * effijx;
87                        tmpptr = localhom[i][j];
88                        while( tmpptr )
89                        {
90//                              fprintf( stderr, "start1 = %d\n", tmpptr->start1 );
91//                              fprintf( stderr, "end1   = %d\n", tmpptr->end1   );
92//                              fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] );
93//                              fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] );
94                                pt = seq1[i];
95                                tmpint = -1;
96                                while( *pt != 0 )
97                                {
98                                        if( *pt++ != '-' ) tmpint++;
99                                        if( tmpint == tmpptr->start1 ) break;
100                                }
101                                start1 = (int)( pt - seq1[i] ) - 1;
102       
103                                if( tmpptr->start1 == tmpptr->end1 ) end1 = start1;
104                                else
105                                {
106#if MACHIGAI
107                                        while( *pt != 0 )
108                                        {
109                                                if( tmpint == tmpptr->end1 ) break;
110                                                if( *pt++ != '-' ) tmpint++;
111                                        }
112                                        end1 = (int)( pt - seq1[i] ) - 1;
113#else
114                                        while( *pt != 0 )
115                                        {
116//                                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
117                                                if( *pt++ != '-' ) tmpint++;
118                                                if( tmpint == tmpptr->end1 ) break;
119                                        }
120                                        end1 = (int)( pt - seq1[i] ) - 1;
121#endif
122                                }
123       
124                                pt = seq2[j];
125                                tmpint = -1;
126                                while( *pt != 0 )
127                                {
128                                        if( *pt++ != '-' ) tmpint++;
129                                        if( tmpint == tmpptr->start2 ) break;
130                                }
131                                start2 = (int)( pt - seq2[j] ) - 1;
132                                if( tmpptr->start2 == tmpptr->end2 ) end2 = start2;
133                                else
134                                {
135#if MACHIGAI
136                                        while( *pt != 0 )
137                                        {
138                                                if( tmpint == tmpptr->end2 ) break;
139                                                if( *pt++ != '-' ) tmpint++;
140                                        }
141                                        end2 = (int)( pt - seq2[j] ) - 1;
142#else
143                                        while( *pt != 0 )
144                                        {
145                                                if( *pt++ != '-' ) tmpint++;
146                                                if( tmpint == tmpptr->end2 ) break;
147                                        }
148                                        end2 = (int)( pt - seq2[j] ) - 1;
149#endif
150                                }
151//                              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] );
152//                              fprintf( stderr, "step 0\n" );
153                                if( end1 - start1 != end2 - start2 )
154                                {
155//                                      fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
156                                }
157
158                                k1 = start1; k2 = start2;
159                                pt1 = seq1[i] + k1;
160                                pt2 = seq2[j] + k2;
161                                while( *pt1 && *pt2 )
162                                {
163                                        if( *pt1 != '-' && *pt2 != '-' )
164                                        {
165// œÅ€ß€òÆóœÅ€Ë€«€±€Ê€€€è€Š€ËÃí°Õ€·€Æ²Œ€µ€€¡£
166//                                              impmtx[k1][k2] += tmpptr->wimportance * fastathreshold;
167//                                              impmtx[k1][k2] += tmpptr->importance * effij;
168                                                impmtx[k1][k2] += tmpptr->fimportance * effij;
169//                                              fprintf( stderr, "k1=%d, k2=%d, impalloclen=%d\n", k1, k2, impalloclen );
170//                                              fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
171                                                k1++; k2++;
172                                                pt1++; pt2++;
173                                        }
174                                        else if( *pt1 != '-' && *pt2 == '-' )
175                                        {
176//                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
177                                                k2++; pt2++;
178                                        }
179                                        else if( *pt1 == '-' && *pt2 != '-' )
180                                        {
181//                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
182                                                k1++; pt1++;
183                                        }
184                                        else if( *pt1 == '-' && *pt2 == '-' )
185                                        {
186//                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
187                                                k1++; pt1++;
188                                                k2++; pt2++;
189                                        }
190                                        if( k1 > end1 || k2 > end2 ) break;
191                                }
192                                tmpptr = tmpptr->next;
193                        }
194                }
195        }
196#if 0
197        fprintf( stderr, "impmtx = \n" );
198        for( k2=0; k2<lgth2; k2++ )
199                fprintf( stderr, "%6.3f ", (double)k2 );
200        fprintf( stderr, "\n" );
201        for( k1=0; k1<lgth1; k1++ )
202        {
203                fprintf( stderr, "%d", k1 );
204                for( k2=0; k2<lgth2; k2++ )
205                        fprintf( stderr, "%2.1f ", impmtx[k1][k2] );
206                fprintf( stderr, "\n" );
207        }
208        exit( 1 );
209#endif
210}
211
212
213void part_imp_rnaQ( int nseq1, int nseq2, char **seq1, char **seq2, double *eff1, double *eff2, RNApair ***grouprna1, RNApair ***grouprna2, int *gapmap1, int *gapmap2, RNApair *additionalpair )
214{
215        foldrna( nseq1, nseq2, seq1, seq2, eff1, eff2, grouprna1, grouprna2, impmtx, gapmap1, gapmap2, additionalpair );
216}
217
218
219void part_imp_match_initQ( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom )
220{
221        int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2;
222        static int impalloclen = 0;
223        char *pt;
224        static char *nocount1 = NULL;
225        static char *nocount2 = NULL;
226
227        if( impalloclen < lgth1 || impalloclen < lgth2 )
228        {
229                if( impmtx ) FreeFloatMtx( impmtx );
230                if( nocount1 ) free( nocount1 );
231                if( nocount2 ) free( nocount2 );
232                impalloclen = MAX( lgth1, lgth2 ) + 2;
233                impmtx = AllocateFloatMtx( impalloclen, impalloclen );
234                nocount1 = AllocateCharVec( impalloclen );
235                nocount2 = AllocateCharVec( impalloclen );
236                impalloclen -= 2;
237        }
238
239        for( i=0; i<lgth1; i++ )
240        {
241                for( j=0; j<clus1; j++ )
242                        if( seq1[j][i] == '-' ) break;
243                if( j != clus1 ) nocount1[i] = 1; 
244                else                     nocount1[i] = 0;
245        }
246        for( i=0; i<lgth2; i++ )
247        {
248                for( j=0; j<clus2; j++ )
249                        if( seq2[j][i] == '-' ) break;
250                if( j != clus2 ) nocount2[i] = 1;
251                else                     nocount2[i] = 0;
252        }
253
254#if 0
255fprintf( stderr, "nocount2 =\n" );
256for( i = 0; i<impalloclen; i++ )
257{
258        fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
259}
260#endif
261
262        for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
263                impmtx[i][j] = 0.0;
264        for( i=0; i<clus1; i++ )
265        {
266                fprintf( stderr, "i = %d, seq1 = %s\n", i, seq1[i] );
267                for( j=0; j<clus2; j++ )
268                {
269                        fprintf( stderr, "start1 = %d\n", localhom[i][j]->start1 );
270                        fprintf( stderr, "end1   = %d\n", localhom[i][j]->end1   );
271                        fprintf( stderr, "j = %d, seq2 = %s\n", j, seq2[j] );
272                        pt = seq1[i];
273                        tmpint = -1;
274                        while( *pt != 0 )
275                        {
276                                if( *pt++ != '-' ) tmpint++;
277                                if( tmpint == localhom[i][j]->start1 ) break;
278                        }
279                        start1 = pt - seq1[i] - 1;
280
281                        while( *pt != 0 )
282                        {
283//                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, localhom[i][j].end1, pt-seq1[i] );
284                                if( *pt++ != '-' ) tmpint++;
285                                if( tmpint == localhom[i][j]->end1 ) break;
286                        }
287                        end1 = pt - seq1[i] - 1;
288
289                        pt = seq2[j];
290                        tmpint = -1;
291                        while( *pt != 0 )
292                        {
293                                if( *pt++ != '-' ) tmpint++;
294                                if( tmpint == localhom[i][j]->start2 ) break;
295                        }
296                        start2 = pt - seq2[j] - 1;
297                        while( *pt != 0 )
298                        {
299                                if( *pt++ != '-' ) tmpint++;
300                                if( tmpint == localhom[i][j]->end2 ) break;
301                        }
302                        end2 = pt - seq2[j] - 1;
303//                      fprintf( stderr, "start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
304                        k1 = start1;
305                        k2 = start2;
306                        fprintf( stderr, "step 0\n" );
307                        while( k1 <= end1 && k2 <= end2 )
308                        {
309#if 0
310                                if( !nocount1[k1] && !nocount2[k2] )
311                                        impmtx[k1][k2] += localhom[i][j].wimportance * eff1[i] * eff2[j];
312                                k1++; k2++;
313#else
314                                if( !nocount1[k1] && !nocount2[k2] )
315                                        impmtx[k1][k2] += localhom[i][j]->wimportance * eff1[i] * eff2[j];
316                                k1++; k2++;
317#endif
318                        }
319
320                        dif = ( end1 - start1 ) - ( end2 - start2 );
321                        fprintf( stderr, "dif = %d\n", dif );
322                        if( dif > 0 )
323                        {
324                                do
325                                {
326                                        fprintf( stderr, "dif = %d\n", dif );
327                                        k1 = start1;
328                                        k2 = start2 - dif;
329                                        while( k1 <= end1 && k2 <= end2 )
330                                        {
331                                                if( 0 <= k2 && start2 <= k2 && !nocount1[k1] && !nocount2[k2] )
332                                                        impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
333                                                k1++; k2++;
334                                        }
335                                }
336                                while( dif-- );
337                        }
338                        else
339                        {
340                                do
341                                {
342                                        k1 = start1 + dif;
343                                        k2 = start2;
344                                        while( k1 <= end1 )
345                                        {
346                                                if( k1 >= 0 && k1 >= start1 && !nocount1[k1] && !nocount2[k2] )
347                                                        impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
348                                                k1++; k2++;
349                                        }
350                                }
351                                while( dif++ );
352                        }
353                }
354        }
355#if 0
356        fprintf( stderr, "impmtx = \n" );
357        for( k2=0; k2<lgth2; k2++ )
358                fprintf( stderr, "%6.3f ", (double)k2 );
359        fprintf( stderr, "\n" );
360        for( k1=0; k1<lgth1; k1++ )
361        {
362                fprintf( stderr, "%d", k1 );
363                for( k2=0; k2<lgth2; k2++ )
364                        fprintf( stderr, "%6.3f ", impmtx[k1][k2] );
365                fprintf( stderr, "\n" );
366        }
367        exit( 1 );
368#endif
369}
370
371static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
372{
373#if FASTMATCHCALC
374        int j, l;
375        float scarr[26];
376        float **cpmxpd = floatwork;
377        int **cpmxpdn = intwork;
378        float *matchpt, *cpmxpdpt, **cpmxpdptpt;
379        int *cpmxpdnpt, **cpmxpdnptpt;
380        if( initialize )
381        {
382                int count = 0;
383                for( j=0; j<lgth2; j++ )
384                {
385                        count = 0;
386                        for( l=0; l<26; l++ )
387                        {
388                                if( cpmx2[l][j] )
389                                {
390                                        cpmxpd[j][count] = cpmx2[l][j];
391                                        cpmxpdn[j][count] = l;
392                                        count++;
393                                }
394                        }
395                        cpmxpdn[j][count] = -1;
396                }
397        }
398
399        {
400                for( l=0; l<26; l++ )
401                {
402                        scarr[l] = 0.0;
403                        for( j=0; j<26; j++ )
404//                              scarr[l] += n_dis[j][l] * cpmx1[j][i1];
405                                scarr[l] += n_dis_consweight_multi[j][l] * cpmx1[j][i1];
406                }
407                matchpt = match;
408                cpmxpdnptpt = cpmxpdn;
409                cpmxpdptpt = cpmxpd;
410                while( lgth2-- )
411                {
412                        *matchpt = 0.0;
413                        cpmxpdnpt = *cpmxpdnptpt++;
414                        cpmxpdpt = *cpmxpdptpt++;
415                        while( *cpmxpdnpt>-1 )
416                                *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;
417                        matchpt++;
418                } 
419        }
420#else
421        int j, k, l;
422        float scarr[26];
423        float **cpmxpd = floatwork;
424        int **cpmxpdn = intwork;
425        // simple
426        if( initialize )
427        {
428                int count = 0;
429                for( j=0; j<lgth2; j++ )
430                {
431                        count = 0;
432                        for( l=0; l<26; l++ )
433                        {
434                                if( cpmx2[l][j] )
435                                {
436                                        cpmxpd[count][j] = cpmx2[l][j];
437                                        cpmxpdn[count][j] = l;
438                                        count++;
439                                }
440                        }
441                        cpmxpdn[count][j] = -1;
442                }
443        }
444        for( l=0; l<26; l++ )
445        {
446                scarr[l] = 0.0;
447                for( k=0; k<26; k++ )
448//                      scarr[l] += n_dis[k][l] * cpmx1[k][i1];
449                        scarr[l] += n_dis_consweight_multi[k][l] * cpmx1[k][i1];
450        }
451        for( j=0; j<lgth2; j++ )
452        {
453                match[j] = 0.0;
454                for( k=0; cpmxpdn[k][j]>-1; k++ )
455                        match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
456        } 
457#endif
458}
459
460static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw, 
461                                                char **seq1, char **seq2, 
462                        char **mseq1, char **mseq2, 
463                        float **cpmx1, float **cpmx2, 
464                        int **ijp, int icyc, int jcyc,
465                                                int start1, int end1, int start2, int end2,
466                                                int *gapmap1, int *gapmap2 )
467{
468        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
469//      char gap[] = "-";
470        char *gap;
471        float wm;
472        gap = newgapstr;
473        lgth1 = strlen( seq1[0] );
474        lgth2 = strlen( seq2[0] );
475
476#if 0
477        for( i=0; i<lgth1; i++ )
478        {
479                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
480        }
481#endif
482 
483        if( outgap == 1 )
484                ;
485        else
486        {
487                wm = lastverticalw[0];
488                for( i=0; i<lgth1; i++ )
489                {
490                        if( lastverticalw[i] >= wm )
491                        {
492                                wm = lastverticalw[i];
493                                iin = i; jin = lgth2-1;
494                                ijp[lgth1][lgth2] = +( lgth1 - i );
495                        }
496                }
497                for( j=0; j<lgth2; j++ )
498                {
499                        if( lasthorizontalw[j] >= wm )
500                        {
501                                wm = lasthorizontalw[j];
502                                iin = lgth1-1; jin = j;
503                                ijp[lgth1][lgth2] = -( lgth2 - j );
504                        }
505                }
506        }
507
508    for( i=0; i<lgth1+1; i++ ) 
509    {
510        ijp[i][0] = i + 1;
511    }
512    for( j=0; j<lgth2+1; j++ ) 
513    {
514        ijp[0][j] = -( j + 1 );
515    }
516
517        for( i=0; i<icyc; i++ )
518        {
519                mseq1[i] += lgth1+lgth2;
520                *mseq1[i] = 0;
521        }
522        for( j=0; j<jcyc; j++ )
523        {
524                mseq2[j] += lgth1+lgth2;
525                *mseq2[j] = 0;
526        }
527        iin = lgth1; jin = lgth2;
528        *impwmpt = 0.0;
529        for( k=0; k<=lgth1+lgth2; k++ ) 
530        {
531                if( ijp[iin][jin] < 0 ) 
532                {
533                        ifi = iin-1; jfi = jin+ijp[iin][jin];
534                }
535                else if( ijp[iin][jin] > 0 )
536                {
537                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
538                }
539                else
540                {
541                        ifi = iin-1; jfi = jin-1;
542                }
543                l = iin - ifi;
544                while( --l ) 
545                {
546                        for( i=0; i<icyc; i++ )
547                                *--mseq1[i] = seq1[i][ifi+l];
548                        for( j=0; j<jcyc; j++ ) 
549                                *--mseq2[j] = *gap;
550                        k++;
551                }
552                l= jin - jfi;
553                while( --l )
554                {
555                        for( i=0; i<icyc; i++ ) 
556                                *--mseq1[i] = *gap;
557                        for( j=0; j<jcyc; j++ ) 
558                                *--mseq2[j] = seq2[j][jfi+l];
559                        k++;
560                }
561                if( iin != lgth1 && jin != lgth2 ) // ??
562                {
563                        *impwmpt += part_imp_match_out_scQ( gapmap1[iin]+start1, gapmap2[jin]+start2 );
564//                      fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
565                }
566                if( iin <= 0 || jin <= 0 ) break;
567                for( i=0; i<icyc; i++ ) 
568                        *--mseq1[i] = seq1[i][ifi];
569                for( j=0; j<jcyc; j++ ) 
570                        *--mseq2[j] = seq2[j][jfi];
571                k++;
572                iin = ifi; jin = jfi;
573        }
574}
575static float Atracking( float *lasthorizontalw, float *lastverticalw, 
576                                                char **seq1, char **seq2, 
577                        char **mseq1, char **mseq2, 
578                        float **cpmx1, float **cpmx2, 
579                        int **ijp, int icyc, int jcyc )
580{
581        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, lastk;
582//      char gap[] = "-";
583        char *gap;
584        float wm = 0.0;
585        gap = newgapstr;
586        lgth1 = strlen( seq1[0] );
587        lgth2 = strlen( seq2[0] );
588
589#if 0
590        for( i=0; i<lgth1; i++ )
591        {
592                fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
593        }
594#endif
595 
596        if( outgap == 1 )
597                ;
598        else
599        {
600                wm = lastverticalw[0];
601                for( i=0; i<lgth1; i++ )
602                {
603                        if( lastverticalw[i] >= wm )
604                        {
605                                wm = lastverticalw[i];
606                                iin = i; jin = lgth2-1;
607                                ijp[lgth1][lgth2] = +( lgth1 - i );
608                        }
609                }
610                for( j=0; j<lgth2; j++ )
611                {
612                        if( lasthorizontalw[j] >= wm )
613                        {
614                                wm = lasthorizontalw[j];
615                                iin = lgth1-1; jin = j;
616                                ijp[lgth1][lgth2] = -( lgth2 - j );
617                        }
618                }
619        }
620
621    for( i=0; i<lgth1+1; i++ ) 
622    {
623        ijp[i][0] = i + 1;
624    }
625    for( j=0; j<lgth2+1; j++ ) 
626    {
627        ijp[0][j] = -( j + 1 );
628    }
629
630        for( i=0; i<icyc; i++ )
631        {
632                mseq1[i] += lgth1+lgth2;
633                *mseq1[i] = 0;
634        }
635        for( j=0; j<jcyc; j++ )
636        {
637                mseq2[j] += lgth1+lgth2;
638                *mseq2[j] = 0;
639        }
640        iin = lgth1; jin = lgth2;
641        lastk = lgth1+lgth2;
642        for( k=0; k<=lastk; k++ ) 
643        {
644                if( ijp[iin][jin] < 0 ) 
645                {
646                        ifi = iin-1; jfi = jin+ijp[iin][jin];
647                }
648                else if( ijp[iin][jin] > 0 )
649                {
650                        ifi = iin-ijp[iin][jin]; jfi = jin-1;
651                }
652                else
653                {
654                        ifi = iin-1; jfi = jin-1;
655                }
656                l = iin - ifi;
657                while( --l ) 
658                {
659                        for( i=0; i<icyc; i++ )
660                                *--mseq1[i] = seq1[i][ifi+l];
661                        for( j=0; j<jcyc; j++ ) 
662                                *--mseq2[j] = *gap;
663                        k++;
664                }
665                l= jin - jfi;
666                while( --l )
667                {
668                        for( i=0; i<icyc; i++ ) 
669                                *--mseq1[i] = *gap;
670                        for( j=0; j<jcyc; j++ ) 
671                                *--mseq2[j] = seq2[j][jfi+l];
672                        k++;
673                }
674                if( iin <= 0 || jin <= 0 ) break;
675                for( i=0; i<icyc; i++ ) 
676                        *--mseq1[i] = seq1[i][ifi];
677                for( j=0; j<jcyc; j++ ) 
678                        *--mseq2[j] = seq2[j][jfi];
679                k++;
680                iin = ifi; jin = jfi;
681        }
682        return( 0.0 );
683}
684
685float partQ__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, int start1, int end1, int start2, int end2, int *gapmap1, int *gapmap2, char *sgap1, char *sgap2, char *egap1, char *egap2 )
686/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
687{
688//      int k;
689        register int i, j;
690        int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
691        int lgth1, lgth2;
692        int resultlen;
693        float wm = 0.0;   /* int ?????? */
694        float g;
695        float *currentw, *previousw;
696#if 1
697        float *wtmp;
698        int *ijppt;
699        float *mjpt, *prept, *curpt;
700        int *mpjpt;
701#endif
702        static float mi, *m;
703        static int **ijp;
704        static int mpi, *mp;
705        static float *w1, *w2;
706        static float *match;
707        static float *initverticalw;    /* kufuu sureba iranai */
708        static float *lastverticalw;    /* kufuu sureba iranai */
709        static char **mseq1;
710        static char **mseq2;
711        static char **mseq;
712        static float *digf1; 
713        static float *digf2; 
714        static float *diaf1; 
715        static float *diaf2; 
716        static float *gapz1; 
717        static float *gapz2; 
718        static float *gapf1; 
719        static float *gapf2; 
720        static float *ogcp1g;
721        static float *ogcp2g;
722        static float *fgcp1g;
723        static float *fgcp2g;
724        static float *og_h_dg_n1_p;
725        static float *og_h_dg_n2_p;
726        static float *fg_h_dg_n1_p;
727        static float *fg_h_dg_n2_p;
728        static float *og_t_fg_h_dg_n1_p;
729        static float *og_t_fg_h_dg_n2_p;
730        static float *fg_t_og_h_dg_n1_p;
731        static float *fg_t_og_h_dg_n2_p;
732        static float *gapz_n1;
733        static float *gapz_n2;
734        static float **cpmx1;
735        static float **cpmx2;
736        static int **intwork;
737        static float **floatwork;
738        static int orlgth1 = 0, orlgth2 = 0;
739        float fpenalty = (float)penalty;
740#if USE_PENALTY_EX
741        float fpenalty_ex = (float)penalty_ex;
742#endif
743        float tmppenal;
744        float *fg_t_og_h_dg_n2_p_pt;
745        float *og_t_fg_h_dg_n2_p_pt;
746        float *og_h_dg_n2_p_pt;
747        float *fg_h_dg_n2_p_pt;
748        float *gapz_n2_pt0;
749        float *gapz_n2_pt1;
750        float *fgcp2pt;
751        float *ogcp2pt;
752        float fg_t_og_h_dg_n1_p_va;
753        float og_t_fg_h_dg_n1_p_va;
754        float og_h_dg_n1_p_va;
755        float fg_h_dg_n1_p_va;
756        float gapz_n1_va0;
757        float gapz_n1_va1;
758        float fgcp1va;
759        float ogcp1va;
760
761
762
763#if 0
764        fprintf( stderr, "eff in SA+++align\n" );
765        for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
766#endif
767        if( orlgth1 == 0 )
768        {
769                mseq1 = AllocateCharMtx( njob, 0 );
770                mseq2 = AllocateCharMtx( njob, 0 );
771        }
772
773
774        lgth1 = strlen( seq1[0] );
775        lgth2 = strlen( seq2[0] );
776
777        if( lgth1 > orlgth1 || lgth2 > orlgth2 )
778        {
779                int ll1, ll2;
780
781                if( orlgth1 > 0 && orlgth2 > 0 )
782                {
783                        FreeFloatVec( w1 );
784                        FreeFloatVec( w2 );
785                        FreeFloatVec( match );
786                        FreeFloatVec( initverticalw );
787                        FreeFloatVec( lastverticalw );
788
789                        FreeFloatVec( m );
790                        FreeIntVec( mp );
791
792                        FreeCharMtx( mseq );
793
794                        FreeFloatVec( digf1 );
795                        FreeFloatVec( digf2 );
796                        FreeFloatVec( diaf1 );
797                        FreeFloatVec( diaf2 );
798                        FreeFloatVec( gapz1 );
799                        FreeFloatVec( gapz2 );
800                        FreeFloatVec( gapf1 );
801                        FreeFloatVec( gapf2 );
802                        FreeFloatVec( ogcp1g );
803                        FreeFloatVec( ogcp2g );
804                        FreeFloatVec( fgcp1g );
805                        FreeFloatVec( fgcp2g );
806                        FreeFloatVec( og_h_dg_n1_p );
807                        FreeFloatVec( og_h_dg_n2_p );
808                        FreeFloatVec( fg_h_dg_n1_p );
809                        FreeFloatVec( fg_h_dg_n2_p );
810                        FreeFloatVec( og_t_fg_h_dg_n1_p );
811                        FreeFloatVec( og_t_fg_h_dg_n2_p );
812                        FreeFloatVec( fg_t_og_h_dg_n1_p );
813                        FreeFloatVec( fg_t_og_h_dg_n2_p );
814                        FreeFloatVec( gapz_n1 );
815                        FreeFloatVec( gapz_n2 );
816
817                        FreeFloatMtx( cpmx1 );
818                        FreeFloatMtx( cpmx2 );
819
820                        FreeFloatMtx( floatwork );
821                        FreeIntMtx( intwork );
822                }
823
824                ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
825                ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
826
827#if DEBUG
828                fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
829#endif
830
831                w1 = AllocateFloatVec( ll2+2 );
832                w2 = AllocateFloatVec( ll2+2 );
833                match = AllocateFloatVec( ll2+2 );
834
835                initverticalw = AllocateFloatVec( ll1+2 );
836                lastverticalw = AllocateFloatVec( ll1+2 );
837
838                m = AllocateFloatVec( ll2+2 );
839                mp = AllocateIntVec( ll2+2 );
840
841                mseq = AllocateCharMtx( njob, ll1+ll2 );
842
843                digf1 = AllocateFloatVec( ll1+2 );
844                digf2 = AllocateFloatVec( ll2+2 );
845                diaf1 = AllocateFloatVec( ll1+2 );
846                diaf2 = AllocateFloatVec( ll2+2 );
847                gapz1 = AllocateFloatVec( ll1+2 );
848                gapz2 = AllocateFloatVec( ll2+2 );
849                gapf1 = AllocateFloatVec( ll1+2 );
850                gapf2 = AllocateFloatVec( ll2+2 );
851                ogcp1g = AllocateFloatVec( ll1+2 );
852                ogcp2g = AllocateFloatVec( ll2+2 );
853                fgcp1g = AllocateFloatVec( ll1+2 );
854                fgcp2g = AllocateFloatVec( ll2+2 );
855                og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
856                og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
857                fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
858                fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
859                og_t_fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
860                og_t_fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
861                fg_t_og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );
862                fg_t_og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );
863                gapz_n1 = AllocateFloatVec( ll1+2 );
864                gapz_n2 = AllocateFloatVec( ll2+2 );
865
866                cpmx1 = AllocateFloatMtx( 26, ll1+2 );
867                cpmx2 = AllocateFloatMtx( 26, ll2+2 );
868
869#if FASTMATCHCALC
870                floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 ); 
871                intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 ); 
872#else
873                floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 ); 
874                intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 ); 
875#endif
876
877#if DEBUG
878                fprintf( stderr, "succeeded\n" );
879#endif
880
881                orlgth1 = ll1 - 100;
882                orlgth2 = ll2 - 100;
883        }
884
885
886        for( i=0; i<icyc; i++ ) mseq1[i] = mseq[i];
887        for( j=0; j<jcyc; j++ ) mseq2[j] = mseq[icyc+j];
888
889
890        if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
891        {
892                int ll1, ll2;
893
894                if( commonAlloc1 && commonAlloc2 )
895                {
896                        FreeIntMtx( commonIP );
897                }
898
899                ll1 = MAX( orlgth1, commonAlloc1 );
900                ll2 = MAX( orlgth2, commonAlloc2 );
901
902#if DEBUG
903                fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
904#endif
905
906                commonIP = AllocateIntMtx( ll1+10, ll2+10 );
907
908#if DEBUG
909                fprintf( stderr, "succeeded\n\n" );
910#endif
911
912                commonAlloc1 = ll1;
913                commonAlloc2 = ll2;
914        }
915        ijp = commonIP;
916
917        cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
918        cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
919
920        if( sgap1 )
921        {
922                new_OpeningGapCount_zure( ogcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
923                new_OpeningGapCount_zure( ogcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
924                new_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );
925                new_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
926                getdigapfreq_part( digf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
927                getdigapfreq_part( digf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
928                getdiaminofreq_part( diaf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );
929                getdiaminofreq_part( diaf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );
930                getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
931                getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
932                getgapfreq_zure_part( gapz1, icyc, seq1, eff1, lgth1, sgap1 );
933                getgapfreq_zure_part( gapz2, jcyc, seq2, eff2, lgth2, sgap1 );
934        }
935        else
936        {
937                st_OpeningGapCount( ogcp1g, icyc, seq1, eff1, lgth1 );
938                st_OpeningGapCount( ogcp2g, jcyc, seq2, eff2, lgth2 );
939                st_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1 );
940                st_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2 );
941                getdigapfreq_st( digf1, icyc, seq1, eff1, lgth1 );
942                getdigapfreq_st( digf2, jcyc, seq2, eff2, lgth2 );
943                getdiaminofreq_x( diaf1, icyc, seq1, eff1, lgth1 );
944                getdiaminofreq_x( diaf2, jcyc, seq2, eff2, lgth2 );
945                getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );
946                getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );
947                getgapfreq_zure( gapz1, icyc, seq1, eff1, lgth1 );
948                getgapfreq_zure( gapz2, jcyc, seq2, eff2, lgth2 );
949        }
950
951#if 1
952        lastj = lgth2+2;
953        for( i=0; i<lastj; i++ )
954        {
955                og_h_dg_n2_p[i] = ( 1.0-ogcp2g[i]-digf2[i] ) * fpenalty * 0.5;
956                fg_h_dg_n2_p[i] = ( 1.0-fgcp2g[i]-digf2[i] ) * fpenalty * 0.5;
957                og_t_fg_h_dg_n2_p[i] = (1.0-ogcp2g[i]+fgcp2g[i]-digf2[i]) * 0.5 * fpenalty;
958                fg_t_og_h_dg_n2_p[i] = (1.0-fgcp2g[i]+ogcp2g[i]-digf2[i]) * 0.5 * fpenalty;
959                gapz_n2[i] = (1.0-gapz2[i]);
960        }
961        lastj = lgth1+2;
962        for( i=0; i<lastj; i++ )
963        {
964                og_h_dg_n1_p[i] = ( 1.0-ogcp1g[i]-digf1[i] ) * fpenalty * 0.5;
965                fg_h_dg_n1_p[i] = ( 1.0-fgcp1g[i]-digf1[i] ) * fpenalty * 0.5;
966                og_t_fg_h_dg_n1_p[i] = (1.0-ogcp1g[i]+fgcp1g[i]-digf1[i]) * 0.5 * fpenalty;
967                fg_t_og_h_dg_n1_p[i] = (1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) * 0.5 * fpenalty;
968                gapz_n1[i] = (1.0-gapz1[i]);
969        }
970#endif
971
972        currentw = w1;
973        previousw = w2;
974
975
976        match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
977        if( localhom )
978                part_imp_match_out_vead_tate_gapmapQ( initverticalw, gapmap2[0]+start2, lgth1, start1, gapmap1 );
979
980
981        match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
982        if( localhom )
983                part_imp_match_out_vead_gapmapQ( currentw, gapmap1[0]+start1, lgth2, start2, gapmap2 );
984#if 0 // -> tbfast.c
985        if( localhom )
986                imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
987
988#endif
989
990        if( outgap == 1 )
991        {
992                g = 0.0;
993
994                g += ogcp1g[0] * og_h_dg_n2_p[0];
995
996                g += ogcp2g[0] * og_h_dg_n1_p[0];
997
998                g += fgcp1g[0] * fg_h_dg_n2_p[0];
999
1000                g += fgcp2g[0] * fg_h_dg_n1_p[0];
1001
1002                initverticalw[0] += g;
1003                currentw[0] += g;
1004
1005                for( i=1; i<lgth1+1; i++ )
1006                {
1007                        tmppenal = gapz_n2[0]*og_t_fg_h_dg_n1_p[0];
1008                        initverticalw[i] += tmppenal;
1009
1010                        tmppenal = gapz_n2[1]*fg_t_og_h_dg_n1_p[i];
1011                        initverticalw[i] += tmppenal;
1012
1013                }
1014                for( j=1; j<lgth2+1; j++ )
1015                {
1016                        tmppenal = gapz_n1[0]*og_t_fg_h_dg_n2_p[0];
1017                        currentw[j] += tmppenal;
1018
1019                        tmppenal = gapz_n1[1]*fg_t_og_h_dg_n2_p[j];
1020                        currentw[j] += tmppenal;
1021                }
1022        }
1023#if OUTGAP0TRY
1024        else
1025        {
1026                for( j=1; j<lgth2+1; j++ )
1027                                currentw[j] -= offset * j / 2.0;
1028                for( i=1; i<lgth1+1; i++ )
1029                        initverticalw[i] -= offset * i / 2.0;
1030        }
1031#endif
1032
1033        m[0] = 0.0; // iranai
1034        for( j=1; j<lgth2+1; ++j ) 
1035        {
1036                mp[j] = 0;
1037                m[j] = currentw[j-1] + 10000 * fpenalty; //iinoka?
1038        }
1039        if( lgth2 == 0 )
1040                lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari
1041        else
1042                lastverticalw[0] = currentw[lgth2-1];
1043
1044        if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1045        lastj = lgth2+1;
1046
1047#if XXXXXXX
1048fprintf( stderr, "currentw = \n" );
1049for( i=0; i<lgth1+1; i++ )
1050{
1051        fprintf( stderr, "%5.2f ", currentw[i] );
1052}
1053fprintf( stderr, "\n" );
1054fprintf( stderr, "initverticalw = \n" );
1055for( i=0; i<lgth2+1; i++ )
1056{
1057        fprintf( stderr, "%5.2f ", initverticalw[i] );
1058}
1059fprintf( stderr, "\n" );
1060fprintf( stderr, "fcgp\n" );
1061for( i=0; i<lgth1; i++ ) 
1062        fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1063for( i=0; i<lgth2; i++ ) 
1064        fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1065#endif
1066
1067        for( i=1; i<lasti; i++ )
1068        {
1069                wtmp = previousw; 
1070                previousw = currentw;
1071                currentw = wtmp;
1072
1073                previousw[0] = initverticalw[i-1];
1074
1075                match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1076#if XXXXXXX
1077fprintf( stderr, "\n" );
1078fprintf( stderr, "i=%d\n", i );
1079fprintf( stderr, "currentw = \n" );
1080for( j=0; j<lgth2; j++ )
1081{
1082        fprintf( stderr, "%5.2f ", currentw[j] );
1083}
1084fprintf( stderr, "\n" );
1085#endif
1086                if( localhom )
1087                {
1088//                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1089//                      imp_match_out_vead( currentw, i, lgth2 );
1090                        part_imp_match_out_vead_gapmapQ( currentw, gapmap1[i]+start1, lgth2, start2, gapmap2 );
1091                }
1092#if XXXXXXX
1093fprintf( stderr, "\n" );
1094fprintf( stderr, "i=%d\n", i );
1095fprintf( stderr, "currentw = \n" );
1096for( j=0; j<lgth2; j++ )
1097{
1098        fprintf( stderr, "%5.2f ", currentw[j] );
1099}
1100fprintf( stderr, "\n" );
1101#endif
1102                currentw[0] = initverticalw[i];
1103
1104                mpi = 0;
1105                mi = previousw[0] + 10000 * fpenalty;
1106
1107                ijppt = ijp[i] + 1;
1108                mjpt = m + 1;
1109                prept = previousw;
1110                curpt = currentw + 1;
1111                mpjpt = mp + 1;
1112                fg_t_og_h_dg_n2_p_pt = fg_t_og_h_dg_n2_p + 1;
1113                og_t_fg_h_dg_n2_p_pt = og_t_fg_h_dg_n2_p + 1;
1114                og_h_dg_n2_p_pt = og_h_dg_n2_p + 1;
1115                fg_h_dg_n2_p_pt = fg_h_dg_n2_p + 1;
1116                gapz_n2_pt0 = gapz_n2 + 1;
1117                gapz_n2_pt1 = gapz_n2 + 2;
1118                fgcp2pt = fgcp2g + 1;
1119                ogcp2pt = ogcp2g + 1;
1120
1121                fg_t_og_h_dg_n1_p_va = fg_t_og_h_dg_n1_p[i];
1122                og_t_fg_h_dg_n1_p_va = og_t_fg_h_dg_n1_p[i];
1123                og_h_dg_n1_p_va = og_h_dg_n1_p[i];
1124                fg_h_dg_n1_p_va = fg_h_dg_n1_p[i];
1125                gapz_n1_va0 = gapz_n1[i];
1126                gapz_n1_va1 = gapz_n1[i+1];
1127                fgcp1va = fgcp1g[i];
1128                ogcp1va = ogcp1g[i];
1129
1130                for( j=1; j<lastj; j++ )
1131                {
1132                        wm = *prept;
1133
1134                        g = ogcp1va * *og_h_dg_n2_p_pt;
1135                        wm += g;
1136
1137                        g = *ogcp2pt * og_h_dg_n1_p_va;
1138                        wm += g;
1139
1140                        g = fgcp1va * *fg_h_dg_n2_p_pt;
1141                        wm += g;
1142
1143                        g = *fgcp2pt * fg_h_dg_n1_p_va;
1144                        wm += g;
1145
1146                        *ijppt = 0;
1147
1148
1149#if 0
1150                        fprintf( stderr, "%5.0f->", wm );
1151#endif
1152                        tmppenal = gapz_n1_va1 * *fg_t_og_h_dg_n2_p_pt;
1153#if 0
1154                        fprintf( stderr, "%5.0f?", g );
1155#endif
1156                        if( (g=mi+tmppenal) > wm )
1157                        {
1158                                wm = g;
1159                                *ijppt = -( j - mpi );
1160                        }
1161
1162                        tmppenal = gapz_n1_va0 * *og_t_fg_h_dg_n2_p_pt;
1163                        if( (g=*prept+tmppenal) >= mi )
1164                        {
1165                                mi = g;
1166                                mpi = j-1;
1167                        }
1168#if USE_PENALTY_EX
1169                        mi += fpenalty_ex;
1170#endif
1171
1172                        tmppenal = *gapz_n2_pt1 * fg_t_og_h_dg_n1_p_va;
1173#if 0 
1174                        fprintf( stderr, "%5.0f?", g );
1175#endif
1176                        if( (g=*mjpt+tmppenal) > wm )
1177                        {
1178                                wm = g;
1179                                *ijppt = +( i - *mpjpt );
1180                        }
1181
1182                        tmppenal = *gapz_n2_pt0 * og_t_fg_h_dg_n1_p_va;
1183                        if( (g=*prept+tmppenal) >= *mjpt )
1184                        {
1185                                *mjpt = g;
1186                                *mpjpt = i-1;
1187                        }
1188#if USE_PENALTY_EX
1189                        m[j] += fpenalty_ex;
1190#endif
1191
1192#if 0
1193                        fprintf( stderr, "%5.0f ", wm );
1194#endif
1195                        *curpt++ += wm;
1196                        ijppt++;
1197                        mjpt++;
1198                        prept++;
1199                        mpjpt++;
1200                        fg_t_og_h_dg_n2_p_pt++;
1201                        og_t_fg_h_dg_n2_p_pt++;
1202                        og_h_dg_n2_p_pt++;
1203                        fg_h_dg_n2_p_pt++;
1204                        gapz_n2_pt0++;
1205                        gapz_n2_pt1++;
1206                        fgcp2pt++;
1207                        ogcp2pt++;
1208                }
1209                lastverticalw[i] = currentw[lgth2-1];
1210        }
1211
1212#if OUTGAP0TRY
1213        if( !outgap )
1214        {
1215                for( j=1; j<lgth2+1; j++ )
1216                        currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1217                for( i=1; i<lgth1+1; i++ )
1218                        lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1219        }
1220#endif
1221               
1222        /*
1223        fprintf( stderr, "\n" );
1224        for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1225        fprintf( stderr, "#####\n" );
1226        for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1227        fprintf( stderr, "====>" );
1228        for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1229        for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1230        */
1231        if( localhom )
1232        {
1233                Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, start1, end1, start2, end2, gapmap1, gapmap2 );
1234        }
1235        else
1236                Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1237
1238//      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1239
1240        resultlen = strlen( mseq1[0] );
1241        if( alloclen < resultlen || resultlen > N )
1242        {
1243                fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1244                ErrorExit( "LENGTH OVER!\n" );
1245        }
1246
1247
1248        for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1249        for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1250        /*
1251        fprintf( stderr, "\n" );
1252        for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1253        fprintf( stderr, "#####\n" );
1254        for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1255        */
1256
1257//      fprintf( stderr, "wm = %f\n", wm );
1258
1259
1260        return( wm );
1261}
1262
Note: See TracBrowser for help on using the repository browser.