source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/core/tddis.c

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

updated mafft version. Added extensions (no svn ignore, yet)

File size: 17.8 KB
Line 
1#include "mltaln.h"
2
3#define DEBUG 0
4
5#if 0
6void mdfymtx( char pair[njob][njob], int s1, double **partialmtx, double **mtx )
7#else
8void mdfymtx( char **pair, int s1, double **partialmtx, double **mtx )
9#endif
10{
11        int i, j;
12        int icount, jcount;
13#if DEBUG
14        FILE *fp;
15        static char name[M][B];
16
17        for( i=0; i<M; i++ ) name[i][0] = 0;
18        fprintf( stdout, "s1 = %d\n", s1 );
19        for( i=0; i<njob; i++ ) 
20        {
21                for( j=0; j<njob; j++ ) 
22                {
23                        printf( "%#2d", pair[i][j] );
24                }
25                printf( "\n" );
26        }
27#endif
28
29        for( i=0, icount=0; i<njob-1; i++ )
30        {
31                if( !pair[s1][i] ) continue;
32                for( j=i+1, jcount=icount+1; j<njob; j++ ) 
33                {
34                        if( !pair[s1][j] ) continue;
35                        partialmtx[icount][jcount] = mtx[i][j];
36                        jcount++;
37                }
38                icount++;
39        }
40#if DEBUG
41        fp = fopen( "hat2.org", "w" );
42        WriteHat2( fp, njob, name, mtx );
43        fclose( fp );
44        fp = fopen( "hat2.mdf", "w" );
45        WriteHat2( fp, icount, name, partialmtx );
46        fclose( fp );
47#endif
48               
49}
50
51               
52float score_calc( char **seq, int s )  /* method 3  */
53{
54    int i, j, k, c;
55    int len = strlen( seq[0] );
56    float score;
57    int tmpscore;
58    char *mseq1, *mseq2;
59
60    score = 0.0;
61    for( i=0; i<s-1; i++ )
62    {
63        for( j=i+1; j<s; j++ )
64        {
65            mseq1 = seq[i];
66            mseq2 = seq[j];
67            tmpscore = 0;
68            c = 0;
69            for( k=0; k<len; k++ )
70            {
71                if( mseq1[k] == '-' && mseq2[k] == '-' ) continue;
72                c++;
73                tmpscore += amino_dis[(int)mseq1[k]][(int)mseq2[k]];
74                if( mseq1[k] == '-' )
75                {
76                    tmpscore += penalty;
77                    while( mseq1[++k] == '-' )
78                        ;
79                    k--;
80                    if( k > len-2 ) break;
81                    continue;
82                }
83                if( mseq2[k] == '-' )
84                {
85                    tmpscore += penalty;
86                    while( mseq2[++k] == '-' )
87                        ;
88                    k--;
89                    if( k > len-2 ) break;
90                    continue;
91                }
92            }
93            /*
94            if( mseq1[0] == '-' || mseq2[0] == '-' )
95            {
96                for( k=0; k<len; k++ )
97                {
98                    if( mseq1[k] == '-' && mseq2[k] == '-' ) continue;
99                    if( !( mseq1[k] != '-' && mseq2[k] != '-' ) )
100                    {
101                        c--;
102                        tmpscore -= penalty;
103                        break;
104                    }
105                    else break;
106                }
107            }
108            if( mseq1[len-1] == '-' || mseq2[len-1] == '-' )
109            {
110                for( k=0; k<len; k++ )
111                {
112                    if( mseq1[k] == '-' && mseq2[k] == '-' ) continue;
113                    if( !( mseq1[k] != '-' && mseq2[k] != '-' ) )
114                    {
115                        c--;
116                        tmpscore -= penalty;
117                        break;
118                    }
119                    else break;
120                }
121            }
122            */
123            score += (double)tmpscore / (double)c;
124        }
125    }
126    score = (float)score / ( ( (double)s * ((double)s-1.0) ) / 2.0 );
127        fprintf( stderr, "score in score_calc = %f\n", score );
128    return( score );
129}
130
131void cpmx_calc( char **seq, float **cpmx, double *eff, int lgth, int clus )
132{
133        int  i, j, k;
134        double totaleff = 0.0;
135
136        for( i=0; i<clus; i++ ) totaleff += eff[i]; 
137        for( i=0; i<26; i++ ) for( j=0; j<lgth; j++ ) cpmx[i][j] = 0.0;
138        for( j=0; j<lgth; j++ ) for( k=0; k<clus; k++ )
139                        cpmx[(int)amino_n[(int)seq[k][j]]][j] += (float)eff[k] / totaleff;
140}
141
142
143void cpmx_calc_new_bk( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
144{
145    int  i, j, k;
146    float feff;
147
148    for( i=0; i<26; i++ ) for( j=0; j<lgth; j++ ) cpmx[i][j] = 0.0;
149    for( k=0; k<clus; k++ )
150    {
151        feff = (float)eff[k];
152        for( j=0; j<lgth; j++ ) 
153        {
154            cpmx[(int)amino_n[(int)seq[k][j]]][j] += feff;
155        }
156    }
157}
158
159void cpmx_calc_new( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
160{
161        int  i, j, k;
162        float feff;
163        float *cpmxpt, **cpmxptpt;
164        char *seqpt;
165
166        j = 26;
167        cpmxptpt = cpmx;
168        while( j-- )
169        {
170                cpmxpt = *cpmxptpt++;
171                i = lgth;
172                while( i-- )
173                        *cpmxpt++ = 0.0;
174        }
175        for( k=0; k<clus; k++ )
176        {
177                feff = (float)eff[k];
178                seqpt = seq[k];
179//              fprintf( stderr, "seqpt = %s, lgth=%d\n", seqpt, lgth );
180                for( j=0; j<lgth; j++ )
181                {
182                        cpmx[(int)amino_n[(int)*seqpt++]][j] += feff;
183                }
184        }
185}
186void MScpmx_calc_new( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
187{
188        int  i, j, k;
189        float feff;
190        float **cpmxptpt, *cpmxpt;
191        char *seqpt;
192
193        j = lgth;
194        cpmxptpt = cpmx;
195        while( j-- )
196        {
197                cpmxpt = *cpmxptpt++;
198                i = 26;
199                while( i-- )
200                        *cpmxpt++ = 0.0;
201        }
202        for( k=0; k<clus; k++ )
203        {
204                feff = (float)eff[k];
205                seqpt = seq[k];
206                cpmxptpt = cpmx;
207                j = lgth;
208                while( j-- )
209                        (*cpmxptpt++)[(int)amino_n[(int)*seqpt++]] += feff;
210        }
211#if 0
212        for( j=0; j<lgth; j++ ) for( i=0; i<26; i++ ) cpmx[j][i] = 0.0;
213        for( k=0; k<clus; k++ )
214        {
215                feff = (float)eff[k];
216                for( j=0; j<lgth; j++ )
217                        cpmx[j][(int)amino_n[(int)seq[k][j]]] += feff;
218        }
219#endif
220}
221
222void cpmx_ribosum( char **seq, char **seqr, char *dir, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
223{
224        int  i, j, k;
225        float feff;
226        float **cpmxptpt, *cpmxpt;
227        char *seqpt, *seqrpt, *dirpt;
228        int code, code1, code2;
229
230        j = lgth;
231        cpmxptpt = cpmx;
232        while( j-- )
233        {
234                cpmxpt = *cpmxptpt++;
235                i = 37;
236                while( i-- )
237                        *cpmxpt++ = 0.0;
238        }
239        for( k=0; k<clus; k++ )
240        {
241                feff = (float)eff[k];
242                seqpt = seq[k];
243                seqrpt = seqr[k];
244                dirpt = dir;
245                cpmxptpt = cpmx;
246                j = lgth;
247                while( j-- )
248                {
249#if 0
250                        code1 = amino_n[(int)*seqpt];
251                        if( code1 > 3 ) code = 36;
252                        else
253                                code = code1;
254#else
255                        code1 = amino_n[(int)*seqpt];
256                        code2 = amino_n[(int)*seqrpt];
257                        if( code1 > 3 ) 
258                        {
259                                code = 36;
260                        }
261                        else if( code2 > 3 )
262                        {
263                                code = code1;
264                        }
265                        else if( *dirpt == '5' ) 
266                        {
267                                code = 4 + code2 * 4  + code1;
268                        }
269                        else if( *dirpt == '3' ) 
270                        {
271                                code = 20 + code2 * 4  + code1;
272                        }
273                        else // if( *dirpt == 'o' ) // nai
274                        {
275                                code = code1;
276                        }
277#endif
278
279//                      fprintf( stderr, "%c -> code=%d toa=%d, tog=%d, toc=%d, tot=%d, ton=%d, efee=%f\n", *seqpt, code%4, ribosumdis[code][4+0], ribosumdis[code][4+1], ribosumdis[code][4+2], ribosumdis[code][20+3], ribosumdis[code][36], feff );
280
281                        seqpt++;
282                        seqrpt++;
283                        dirpt++;
284                       
285                        (*cpmxptpt++)[code] += feff;
286                }
287        }
288}
289
290void mseqcat( char **seq1, char **seq2, double **eff, double *effarr1, double *effarr2, char name1[M][B], char name2[M][B], int clus1, int clus2 )
291{
292        int i, j;
293    for( i=0; i<clus2; i++ )
294        seq1[clus1+i] = seq2[i];
295    for( i=0; i<clus2; i++ ) strcpy( name1[clus1+i], name2[i] );
296
297        for( i=0; i<clus1; i++ ) for( j=0; j<clus1; j++ ) eff[i][j] = effarr1[i]* effarr1[j]; 
298        for( i=0; i<clus1; i++ ) for( j=clus1; j<clus1+clus2; j++ ) eff[i][j] = effarr1[i]* effarr2[j-clus1]; 
299        for( i=clus1; i<clus1+clus2; i++ ) for( j=0; j<clus1; j++ ) eff[i][j] = effarr2[i-clus1] * effarr1[j]; 
300        for( i=clus1; i<clus1+clus2; i++ ) for( j=clus1; j<clus1+clus2; j++ ) eff[i][j] = effarr2[i-clus1] * effarr2[j-clus1]; 
301}
302
303       
304
305#if 0
306int conjuction( char pair[njob][njob], int s, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
307#else
308int conjuctionforgaln( int s0, int s1, char **seq, char **aseq, double *peff, double *eff, char **name, char **aname, char *d )
309#endif
310{
311        int m, k;
312        char b[B];
313        double total;
314
315#if 0
316        fprintf( stderr, "s0 = %d, s1 = %d\n", s0, s1 );
317#endif
318
319        total = 0.0;
320        d[0] = 0;
321        for( m=s0, k=0; m<s1; m++ )
322        {
323                {
324                        sprintf( b, " %d", m+1 ); 
325#if 1
326                        if( strlen( d ) < 100 ) strcat( d, b );
327#else
328                        strcat( d, b );
329#endif
330                        aseq[k] = seq[m];
331                        peff[k] = eff[m];
332                        total += peff[k];
333#if 0
334                        strcpy( aname[k], name[m] );
335#endif
336                        k++;
337                }
338        }
339#if 1
340        for( m=0; m<k; m++ )
341        {
342                peff[m] /= total;
343//              fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
344        }
345#endif
346        return( k );
347}
348
349void makegrouprna( RNApair ***group, RNApair ***all, int *memlist )
350{
351        int k, m;
352        for( k=0; (m=*memlist)!=-1; memlist++, k++ )
353        {
354                group[k] = all[m];
355        }
356}
357
358
359int fastconjuction_noweight( int *memlist, char **seq, char **aseq, double *peff, char *d )
360{
361        int m, k, dln;
362        char b[B];
363        double total;
364
365#if DEBUG
366        fprintf( stderr, "s = %d\n", s );
367#endif
368
369        total = 0.0;
370        d[0] = 0;
371        dln = 0;
372        for( k=0; *memlist!=-1; memlist++, k++ )
373        {
374                m = *memlist;
375                dln += sprintf( b, " %d", m+1 ); 
376#if 1
377                if( dln < 100 ) strcat( d, b );
378#else
379                strcat( d, b );
380#endif
381                aseq[k] = seq[m];
382                peff[k] = 1.0;
383                total += peff[k];
384        }
385#if 1
386        for( m=0; m<k; m++ )
387                peff[m] /= total;
388#endif
389        return( k );
390}
391
392int fastconjuction_noname_kozo( int *memlist, char **seq, char **aseq, double *peff, double *eff, double *peff_kozo, double *eff_kozo, char *d )
393{
394        int m, k, dln;
395        char b[B];
396        double total;
397        double total_kozo;
398
399#if DEBUG
400        fprintf( stderr, "s = %d\n", s );
401#endif
402
403        total = 0.0;
404        total_kozo = 0.0;
405        d[0] = 0;
406        dln = 0;
407        for( k=0; *memlist!=-1; memlist++, k++ )
408        {
409                m = *memlist;
410                dln += sprintf( b, " %d", m+1 ); 
411#if 1
412                if( dln < 100 ) strcat( d, b );
413#else
414                strcat( d, b );
415#endif
416                aseq[k] = seq[m];
417                peff[k] = eff[m];
418                peff_kozo[k] = eff_kozo[m];
419                total += peff[k];
420                total_kozo += peff_kozo[k];
421        }
422#if 1
423        for( m=0; m<k; m++ )
424        {
425//              fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
426                peff[m] /= total;
427        }
428        if( total_kozo ) 
429        {
430                for( m=0; m<k; m++ )
431                {
432//                      fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
433                        peff_kozo[m] /= total_kozo;
434                        if( peff_kozo[m] > 0.0 ) peff_kozo[m] += peff[m];
435                }
436        }
437        else  //iranai
438        {
439                for( m=0; m<k; m++ )
440                {
441                        peff_kozo[m] = 0.0;
442                }
443        }
444#endif
445
446//      fprintf( stderr, "\n\ntbfast_total_kozo = %f\n\n", total_kozo );
447        return( k );
448}
449
450
451int fastconjuction_noname( int *memlist, char **seq, char **aseq, double *peff, double *eff, char *d )
452{
453        int m, k, dln;
454        char b[B];
455        double total;
456
457#if DEBUG
458        fprintf( stderr, "s = %d\n", s );
459#endif
460
461        total = 0.0;
462        d[0] = 0;
463        dln = 0;
464        for( k=0; *memlist!=-1; memlist++, k++ )
465        {
466                m = *memlist;
467                dln += sprintf( b, " %d", m+1 ); 
468#if 1
469                if( dln < 100 ) strcat( d, b );
470#else
471                strcat( d, b );
472#endif
473                aseq[k] = seq[m];
474                peff[k] = eff[m];
475                total += peff[k];
476        }
477#if 1
478        for( m=0; m<k; m++ )
479        {
480//              fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
481                peff[m] /= total;
482        }
483#endif
484        return( k );
485}
486
487int fastconjuction( int *memlist, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
488{
489        int m, k, dln;
490        char b[B];
491        double total;
492
493#if DEBUG
494        fprintf( stderr, "s = %d\n", s );
495#endif
496
497        total = 0.0;
498        d[0] = 0;
499        dln = 0;
500        for( k=0; *memlist!=-1; memlist++, k++ )
501        {
502                m = *memlist;
503                dln += sprintf( b, " %d", m+1 ); 
504#if 1
505                if( dln < 100 ) strcat( d, b );
506#else
507                strcat( d, b );
508#endif
509                aseq[k] = seq[m];
510                peff[k] = eff[m];
511                total += peff[k];
512#if 0
513                        strcpy( aname[k], name[m] );
514#endif
515        }
516#if 1
517        for( m=0; m<k; m++ )
518                peff[m] /= total;
519#endif
520        return( k );
521}
522
523
524
525int conjuctionfortbfast_old( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char *d )
526{
527        int m, k;
528        char *b;
529        double total;
530
531        b = calloc( B, sizeof( char ) );
532#if DEBUG
533        fprintf( stderr, "s = %d\n", s );
534#endif
535
536        total = 0.0;
537        d[0] = 0;
538        for( m=s, k=0; m<njob; m++ )
539        {
540                if( pair[s][m] != 0 )
541                {
542                        sprintf( b, " %d", m+1 ); 
543#if 1
544                        if( strlen( d ) < 100 ) strcat( d, b );
545#else
546                        strcat( d, b );
547#endif
548                        aseq[k] = seq[m];
549                        peff[k] = eff[m];
550                        total += peff[k];
551#if 0
552                        strcpy( aname[k], name[m] );
553#endif
554                        k++;
555                }
556        }
557#if 1
558        for( m=0; m<k; m++ )
559                peff[m] /= total;
560#endif
561        free( b );
562        return( k );
563}
564
565int conjuction( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char **name, char **aname, char *d )
566{
567        int m, k;
568        char b[B];
569        double total;
570
571#if DEBUG
572        fprintf( stderr, "s = %d\n", s );
573#endif
574
575        total = 0.0;
576        d[0] = 0;
577        for( m=s, k=0; m<njob; m++ )
578        {
579                if( pair[s][m] != 0 )
580                {
581                        sprintf( b, " %d", m+1 ); 
582#if 1
583                        if( strlen( d ) < 100 ) strcat( d, b );
584#else
585                        strcat( d, b );
586#endif
587                        aseq[k] = seq[m];
588                        peff[k] = eff[m];
589                        total += peff[k];
590#if 0
591                        strcpy( aname[k], name[m] );
592#endif
593                        k++;
594                }
595        }
596#if 0
597        for( m=0; m<k; m++ )
598                peff[m] /= total;
599#endif
600        return( k );
601}
602                       
603void floatdelete( float **cpmx, int d, int len )
604{
605        int i, j;
606
607        for( i=d; i<len-1; i++ )
608        {
609                for( j=0; j<26; j++ )
610                {
611                        cpmx[j][i] = cpmx[j][i+1]; 
612                }
613        }
614}
615               
616void chardelete( char *seq, int d )
617{
618        char b[N]; 
619
620                strcpy( b, seq+d+1 );
621                strcpy( seq+d, b );
622}
623
624int RootBranchNode( int nseq, int ***topol, int step, int branch )
625{
626        int i, j, s1, s2, value;
627
628        value = 1;
629        for( i=step+1; i<nseq-2; i++ ) 
630        {
631                for( j=0; (s1=topol[i][0][j])>-1; j++ ) 
632                        if( s1 == topol[step][branch][0] ) value++;
633                for( j=0; (s2=topol[i][1][j])>-1; j++ ) 
634                        if( s2 == topol[step][branch][0] ) value++;
635        }
636        return( value );
637}
638
639void BranchLeafNode( int nseq, int ***topol, int *node, int step, int branch )
640{
641        int i, j, k, s;
642
643        for( i=0; i<nseq; i++ ) node[i] = 0;
644        for( i=0; i<step-1; i++ )
645                for( k=0; k<2; k++ ) 
646                        for( j=0; (s=topol[i][k][j])>-1; j++ ) 
647                                node[s]++;
648        for( k=0; k<branch+1; k++ ) 
649                for( j=0; (s=topol[step][k][j])>-1; j++ )
650                        node[s]++;
651}
652
653void RootLeafNode( int nseq, int ***topol, int *node )
654{
655        int i, j, k, s;
656        for( i=0; i<nseq; i++ ) node[i] = 0;
657        for( i=0; i<nseq-2; i++ )
658                for( k=0; k<2; k++ ) 
659                        for( j=0; (s=topol[i][k][j])>-1; j++ ) 
660                                node[s]++;
661}
662               
663void nodeFromABranch( int nseq, int *result, int **pairwisenode, int ***topol, double **len, int step, int num )
664{
665        int i, s, count;
666        int *innergroup;
667        int *outergroup1;
668#if 0
669        int outergroup2[nseq];
670        int table[nseq];
671#else
672        static int *outergroup2 = NULL;
673        static int *table = NULL;
674        if( outergroup2 == NULL )
675        {
676                outergroup2 = AllocateIntVec( nseq );
677                table = AllocateIntVec( nseq );
678        }
679#endif
680        innergroup = topol[step][num];
681        outergroup1 = topol[step][!num];
682        for( i=0; i<nseq; i++ ) table[i] = 1;
683        for( i=0; (s=innergroup[i])>-1; i++ ) table[s] = 0;
684        for( i=0; (s=outergroup1[i])>-1; i++ ) table[s] = 0;
685        for( i=0, count=0; i<nseq; i++ ) 
686        {
687                if( table[i] )
688                {
689                        outergroup2[count] = i;
690                        count++;
691                }
692        }
693        outergroup2[count] = -1;
694
695        for( i=0; (s=innergroup[i])>-1; i++ )
696        {
697                result[s] = pairwisenode[s][outergroup1[0]]
698                                  + pairwisenode[s][outergroup2[0]]
699                                  - pairwisenode[outergroup1[0]][outergroup2[0]] - 1;
700                result[s] /= 2;
701        }
702        for( i=0; (s=outergroup1[i])>-1; i++ ) 
703        {
704                result[s] = pairwisenode[s][outergroup2[0]]
705                                  + pairwisenode[s][innergroup[0]]
706                                  - pairwisenode[innergroup[0]][outergroup2[0]] + 1;
707                result[s] /= 2;
708        }
709        for( i=0; (s=outergroup2[i])>-1; i++ ) 
710        {
711                result[s] = pairwisenode[s][outergroup1[0]]
712                                  + pairwisenode[s][innergroup[0]]
713                                  - pairwisenode[innergroup[0]][outergroup1[0]] + 1;
714                result[s] /= 2;
715        }
716
717#if 0
718        for( i=0; i<nseq; i++ )
719                fprintf( stderr, "result[%d] = %d\n", i+1, result[i] );
720#endif
721}
722               
723
724
725
726
727       
728
729               
730               
731                                       
732                                       
733
734                               
735void OneClusterAndTheOther_fast( int locnjob, int *memlist1, int *memlist2, int *s1, int *s2, char *pair, int ***topol, int step, int branch )
736{
737        int i, k;
738        int r1;
739//      char *pair;
740
741//      pair = calloc( locnjob, sizeof( char ) );
742
743        for( i=0; i<locnjob; i++ ) pair[i] = 0;
744    for( i=0, k=0; (r1=topol[step][branch][i])>-1; i++ ) 
745        {
746        pair[r1] = 1;
747                memlist1[k++] = r1;
748        }
749        memlist1[k] = -1;
750
751    for( i=0, k=0; i<locnjob; i++ ) 
752    {
753        if( !pair[i] ) 
754        {
755                        memlist2[k++] = i;
756        }
757    }
758        memlist2[k] = -1;
759
760        *s1 = memlist1[0];
761        *s2 = memlist2[0];
762//      free( pair );
763}
764               
765
766void makeEffMtx( int nseq, double **mtx, double *vec )
767{
768        int i, j;
769        for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ ) 
770                mtx[i][j] = vec[i] * vec[j];
771}
772       
773void node_eff( int nseq, double *eff, int *node )
774{
775        int i;
776        extern double ipower( double, int );
777        for( i=0; i<nseq; i++ ) 
778                eff[i] = ipower( 0.5, node[i] ) + geta2;
779        /*
780                eff[i] = ipower( 0.5, node[i] ) + 0;
781        */
782#if DEBUG
783        for( i=0; i<nseq; i++ ) 
784#endif
785}
786
787
788int shrinklocalhom( char **pair, int s1, int s2, LocalHom **localhom, LocalHom ***localhomshrink )
789{
790        int m1, k1, m2, k2;
791
792        for( m1=s1, k1=0; m1<njob; m1++ )
793        {
794                if( pair[s1][m1] != 0 )
795                {
796                        for( m2=s2, k2=0; m2<njob; m2++ )
797                        {
798                                if( pair[s2][m2] != 0 )
799                                {
800                                        if( localhom[m1][m2].opt == -1 )
801                                                localhomshrink[k1][k2] = NULL;
802                                        else
803                                                localhomshrink[k1][k2] = localhom[m1]+m2;
804                                        k2++;
805                                }
806                        }
807                        k1++;
808                }
809        }
810        return( 0 );
811}
812
813int msshrinklocalhom_fast( int *memlist1, int *memlist2, LocalHom **localhom, LocalHom ***localhomshrink )
814{
815        int m1, k1, m2, k2;
816
817        for( k1=0; (m1=memlist1[k1])!=-1; k1++ )
818        {
819                for( k2=0; (m2=memlist2[k2])!=-1; k2++ )
820                {
821                        if( localhom[m1][m2].opt == -1 )
822                                localhomshrink[k1][k2] = NULL;
823                        else
824                                localhomshrink[k1][k2] = localhom[m1]+m2;
825                }
826        }
827        return( 0 );
828}
829int fastshrinklocalhom_one( int *mem1, int *mem2, int norg, LocalHom **localhom, LocalHom ***localhomshrink )
830{
831        int k1, k2;
832        int *intpt1, *intpt2;
833
834       
835        for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
836        {
837                for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
838                {
839                        if( *intpt2 != norg ) 
840                        {
841                                fprintf( stderr, "ERROR! *intpt2 = %d\n", *intpt2 );
842                                exit( 1 );
843                        }
844                        if( localhom[*intpt1][0].opt == -1 )
845                                localhomshrink[k1][k2] = NULL;
846                        else
847                                localhomshrink[k1][k2] = localhom[*intpt1];
848                }
849        }
850        return( 0 );
851}
852
853int fastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
854{
855        int k1, k2;
856        int *intpt1, *intpt2;
857
858       
859        for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
860        {
861                for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
862                {
863                        if( localhom[*intpt1][*intpt2].opt == -1 )
864                                localhomshrink[k1][k2] = NULL;
865                        else
866                                localhomshrink[k1][k2] = localhom[*intpt1]+*intpt2;
867                }
868        }
869        return( 0 );
870}
871
872int msfastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
873{
874        int k1, k2;
875        int *intpt1, *intpt2;
876        int m1, m2;
877       
878        for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
879        {
880                for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
881                {
882                        m1 = MIN(*intpt1,*intpt2); m2 = MAX(*intpt1,*intpt2);
883                        if( localhom[m1][m2].opt == -1 )
884                                localhomshrink[k1][k2] = NULL;
885                        else
886                                localhomshrink[k1][k2] = localhom[m1]+m2;
887                }
888        }
889        return( 0 );
890}
891
Note: See TracBrowser for help on using the repository browser.