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