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