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