Mercurial > repos > nick > duplex

comparison mafft/core/tddis.c @ 18:e4d75f9efb90 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty
author nick
date Thu, 02 Feb 2017 18:44:31 -0500
parents
children
comparison
equal deleted inserted replaced
17:836fa4fe9494 18:e4d75f9efb90
1 #include "mltaln.h"
2
3 #define DEBUG 0
4
5 #if 0
6 void mdfymtx( char pair[njob][njob], int s1, double **partialmtx, double **mtx )
7 #else
8 void 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
52 float 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
131 void 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<nalphabets; 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
143 void 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<nalphabets; 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
159 void 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 = nalphabets;
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 }
186 void 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 = nalphabets;
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<nalphabets; 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
222 void 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
290 void 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
306 int 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
308 int 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
349 void 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
359 int 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
392 int 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
451 int fastconjuction_noname( int *memlist, char **seq, char **aseq, double *peff, double *eff, char *d, double mineff )
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 if( eff[m] < mineff )
475 peff[k] = mineff;
476 else
477 peff[k] = eff[m];
478
479 total += peff[k];
480 }
481 #if 1
482 for( m=0; m<k; m++ )
483 {
484 // fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
485 peff[m] /= total;
486 }
487 #endif
488 return( k );
489 }
490
491 int fastconjuction( int *memlist, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
492 {
493 int m, k, dln;
494 char b[B];
495 double total;
496
497 #if DEBUG
498 fprintf( stderr, "s = %d\n", s );
499 #endif
500
501 total = 0.0;
502 d[0] = 0;
503 dln = 0;
504 for( k=0; *memlist!=-1; memlist++, k++ )
505 {
506 m = *memlist;
507 dln += sprintf( b, " %d", m+1 );
508 #if 1
509 if( dln < 100 ) strcat( d, b );
510 #else
511 strcat( d, b );
512 #endif
513 aseq[k] = seq[m];
514 peff[k] = eff[m];
515 total += peff[k];
516 #if 0
517 strcpy( aname[k], name[m] );
518 #endif
519 }
520 #if 1
521 for( m=0; m<k; m++ )
522 peff[m] /= total;
523 #endif
524 return( k );
525 }
526
527
528
529 int conjuctionfortbfast_old( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char *d )
530 {
531 int m, k;
532 char *b;
533 double total;
534
535 b = calloc( B, sizeof( char ) );
536 #if DEBUG
537 fprintf( stderr, "s = %d\n", s );
538 #endif
539
540 total = 0.0;
541 d[0] = 0;
542 for( m=s, k=0; m<njob; m++ )
543 {
544 if( pair[s][m] != 0 )
545 {
546 sprintf( b, " %d", m+1 );
547 #if 1
548 if( strlen( d ) < 100 ) strcat( d, b );
549 #else
550 strcat( d, b );
551 #endif
552 aseq[k] = seq[m];
553 peff[k] = eff[m];
554 total += peff[k];
555 #if 0
556 strcpy( aname[k], name[m] );
557 #endif
558 k++;
559 }
560 }
561 #if 1
562 for( m=0; m<k; m++ )
563 peff[m] /= total;
564 #endif
565 free( b );
566 return( k );
567 }
568
569 int conjuction( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char **name, char **aname, char *d )
570 {
571 int m, k;
572 char b[B];
573 double total;
574
575 #if DEBUG
576 fprintf( stderr, "s = %d\n", s );
577 #endif
578
579 total = 0.0;
580 d[0] = 0;
581 for( m=s, k=0; m<njob; m++ )
582 {
583 if( pair[s][m] != 0 )
584 {
585 sprintf( b, " %d", m+1 );
586 #if 1
587 if( strlen( d ) < 100 ) strcat( d, b );
588 #else
589 strcat( d, b );
590 #endif
591 aseq[k] = seq[m];
592 peff[k] = eff[m];
593 total += peff[k];
594 #if 0
595 strcpy( aname[k], name[m] );
596 #endif
597 k++;
598 }
599 }
600 #if 0
601 for( m=0; m<k; m++ )
602 peff[m] /= total;
603 #endif
604 return( k );
605 }
606
607 void floatdelete( float **cpmx, int d, int len )
608 {
609 int i, j;
610
611 for( i=d; i<len-1; i++ )
612 {
613 for( j=0; j<nalphabets; j++ )
614 {
615 cpmx[j][i] = cpmx[j][i+1];
616 }
617 }
618 }
619
620 void chardelete( char *seq, int d )
621 {
622 char b[N];
623
624 strcpy( b, seq+d+1 );
625 strcpy( seq+d, b );
626 }
627
628 int RootBranchNode( int nseq, int ***topol, int step, int branch )
629 {
630 int i, j, s1, s2, value;
631
632 value = 1;
633 for( i=step+1; i<nseq-2; i++ )
634 {
635 for( j=0; (s1=topol[i][0][j])>-1; j++ )
636 if( s1 == topol[step][branch][0] ) value++;
637 for( j=0; (s2=topol[i][1][j])>-1; j++ )
638 if( s2 == topol[step][branch][0] ) value++;
639 }
640 return( value );
641 }
642
643 void BranchLeafNode( int nseq, int ***topol, int *node, int step, int branch )
644 {
645 int i, j, k, s;
646
647 for( i=0; i<nseq; i++ ) node[i] = 0;
648 for( i=0; i<step-1; i++ )
649 for( k=0; k<2; k++ )
650 for( j=0; (s=topol[i][k][j])>-1; j++ )
651 node[s]++;
652 for( k=0; k<branch+1; k++ )
653 for( j=0; (s=topol[step][k][j])>-1; j++ )
654 node[s]++;
655 }
656
657 void RootLeafNode( int nseq, int ***topol, int *node )
658 {
659 int i, j, k, s;
660 for( i=0; i<nseq; i++ ) node[i] = 0;
661 for( i=0; i<nseq-2; i++ )
662 for( k=0; k<2; k++ )
663 for( j=0; (s=topol[i][k][j])>-1; j++ )
664 node[s]++;
665 }
666
667 void nodeFromABranch( int nseq, int *result, int **pairwisenode, int ***topol, double **len, int step, int num )
668 {
669 int i, s, count;
670 int *innergroup;
671 int *outergroup1;
672 #if 0
673 int outergroup2[nseq];
674 int table[nseq];
675 #else
676 static int *outergroup2 = NULL;
677 static int *table = NULL;
678 if( outergroup2 == NULL )
679 {
680 outergroup2 = AllocateIntVec( nseq );
681 table = AllocateIntVec( nseq );
682 }
683 #endif
684 innergroup = topol[step][num];
685 outergroup1 = topol[step][!num];
686 for( i=0; i<nseq; i++ ) table[i] = 1;
687 for( i=0; (s=innergroup[i])>-1; i++ ) table[s] = 0;
688 for( i=0; (s=outergroup1[i])>-1; i++ ) table[s] = 0;
689 for( i=0, count=0; i<nseq; i++ )
690 {
691 if( table[i] )
692 {
693 outergroup2[count] = i;
694 count++;
695 }
696 }
697 outergroup2[count] = -1;
698
699 for( i=0; (s=innergroup[i])>-1; i++ )
700 {
701 result[s] = pairwisenode[s][outergroup1[0]]
702 + pairwisenode[s][outergroup2[0]]
703 - pairwisenode[outergroup1[0]][outergroup2[0]] - 1;
704 result[s] /= 2;
705 }
706 for( i=0; (s=outergroup1[i])>-1; i++ )
707 {
708 result[s] = pairwisenode[s][outergroup2[0]]
709 + pairwisenode[s][innergroup[0]]
710 - pairwisenode[innergroup[0]][outergroup2[0]] + 1;
711 result[s] /= 2;
712 }
713 for( i=0; (s=outergroup2[i])>-1; i++ )
714 {
715 result[s] = pairwisenode[s][outergroup1[0]]
716 + pairwisenode[s][innergroup[0]]
717 - pairwisenode[innergroup[0]][outergroup1[0]] + 1;
718 result[s] /= 2;
719 }
720
721 #if 0
722 for( i=0; i<nseq; i++ )
723 fprintf( stderr, "result[%d] = %d\n", i+1, result[i] );
724 #endif
725 }
726
727
728
729
730
731
732
733
734
735
736
737
738
739 void OneClusterAndTheOther_fast( int locnjob, int *memlist1, int *memlist2, int *s1, int *s2, char *pair, int ***topol, int step, int branch, double **smalldistmtx, double **distmtx )
740 {
741 int i, k, j;
742 int r1;
743 // char *pair;
744
745 // pair = calloc( locnjob, sizeof( char ) );
746
747 for( i=0; i<locnjob; i++ ) pair[i] = 0;
748 for( i=0, k=0; (r1=topol[step][branch][i])>-1; i++ )
749 {
750 pair[r1] = 1;
751 memlist1[k++] = r1;
752 }
753 memlist1[k] = -1;
754
755 for( i=0, k=0; i<locnjob; i++ )
756 {
757 if( !pair[i] )
758 {
759 memlist2[k++] = i;
760 }
761 }
762 memlist2[k] = -1;
763
764 *s1 = memlist1[0];
765 *s2 = memlist2[0];
766
767 if( smalldistmtx )
768 {
769 int im, jm;
770 for( i=0; (im=memlist1[i])!=-1; i++ ) for( j=0; (jm=memlist2[j])!=-1; j++ )
771 {
772 smalldistmtx[i][j] = distmtx[MIN(im,jm)][MAX(im,jm)];
773 // fprintf( stderr, "#### %d-%d, %f\n", im, jm, smalldistmtx[i][j] );
774 }
775 }
776 // free( pair );
777 }
778
779
780 void makeEffMtx( int nseq, double **mtx, double *vec )
781 {
782 int i, j;
783 for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ )
784 mtx[i][j] = vec[i] * vec[j];
785 }
786
787 void node_eff( int nseq, double *eff, int *node )
788 {
789 int i;
790 extern double ipower( double, int );
791 for( i=0; i<nseq; i++ )
792 eff[i] = ipower( 0.5, node[i] ) + geta2;
793 /*
794 eff[i] = ipower( 0.5, node[i] ) + 0;
795 */
796 #if DEBUG
797 for( i=0; i<nseq; i++ )
798 #endif
799 }
800
801
802 int shrinklocalhom( char **pair, int s1, int s2, LocalHom **localhom, LocalHom ***localhomshrink )
803 {
804 int m1, k1, m2, k2;
805
806 for( m1=s1, k1=0; m1<njob; m1++ )
807 {
808 if( pair[s1][m1] != 0 )
809 {
810 for( m2=s2, k2=0; m2<njob; m2++ )
811 {
812 if( pair[s2][m2] != 0 )
813 {
814 if( localhom[m1][m2].opt == -1 )
815 localhomshrink[k1][k2] = NULL;
816 else
817 localhomshrink[k1][k2] = localhom[m1]+m2;
818 k2++;
819 }
820 }
821 k1++;
822 }
823 }
824 return( 0 );
825 }
826
827 int msshrinklocalhom_fast( int *memlist1, int *memlist2, LocalHom **localhom, LocalHom ***localhomshrink )
828 {
829 int m1, k1, m2, k2;
830
831 for( k1=0; (m1=memlist1[k1])!=-1; k1++ )
832 {
833 for( k2=0; (m2=memlist2[k2])!=-1; k2++ )
834 {
835 if( localhom[m1][m2].opt == -1 )
836 localhomshrink[k1][k2] = NULL;
837 else
838 localhomshrink[k1][k2] = localhom[m1]+m2;
839 }
840 }
841 return( 0 );
842 }
843 int fastshrinklocalhom_one( int *mem1, int *mem2, int norg, LocalHom **localhom, LocalHom ***localhomshrink )
844 {
845 int k1, k2;
846 int *intpt1, *intpt2;
847
848
849 for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
850 {
851 for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
852 {
853 if( *intpt2 != norg )
854 {
855 fprintf( stderr, "ERROR! *intpt2 = %d\n", *intpt2 );
856 exit( 1 );
857 }
858 if( localhom[*intpt1][0].opt == -1 )
859 localhomshrink[k1][k2] = NULL;
860 else
861 localhomshrink[k1][k2] = localhom[*intpt1];
862 }
863 }
864 return( 0 );
865 }
866
867 int fastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
868 {
869 int k1, k2;
870 int *intpt1, *intpt2;
871
872
873 for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
874 {
875 for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
876 {
877 if( localhom[*intpt1][*intpt2].opt == -1 )
878 localhomshrink[k1][k2] = NULL;
879 else
880 localhomshrink[k1][k2] = localhom[*intpt1]+*intpt2;
881 }
882 }
883 return( 0 );
884 }
885
886 int msfastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
887 {
888 int k1, k2;
889 int *intpt1, *intpt2;
890 int m1, m2;
891
892 for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
893 {
894 for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
895 {
896 m1 = MIN(*intpt1,*intpt2); m2 = MAX(*intpt1,*intpt2);
897 if( localhom[m1][m2].opt == -1 )
898 localhomshrink[k1][k2] = NULL;
899 else
900 localhomshrink[k1][k2] = localhom[m1]+m2;
901 }
902 }
903 return( 0 );
904 }
905