1 #include "mltaln.h"
2
3 #define DEBUG 0
4
5 #define EF_THREEWAY 1.0
6 #define MAXBW 1.0
7 #define MINBW 0.01
8
9 #define MINLEN 0.001
10
11 #if DEBUG
12 Node *stopol_g;
13 #endif
14
15
checkMinusLength(int nseq,double ** len)16 void checkMinusLength( int nseq, double **len )
17 {
18 int i, j;
19 for( i=0; i<nseq-1; i++ ) for( j=0; j<2; j++ )
20 if( len[i][j] < MINLEN ) len[i][j] = MINLEN;
21 }
22
negativeMember2(int * mem,int * query,int locnseq)23 void negativeMember2( int *mem, int *query, int locnseq )
24 {
25 int *ptr;
26 char *tmp;
27 int i;
28 int n;
29
30 tmp = AllocateCharVec( locnseq );
31
32 for( i=0; i<locnseq; i++ ) tmp[i] = 0;
33 while( (n=*query++) != -1 ) tmp[n] = 1;
34
35 ptr = mem;
36 for( i=0; i<locnseq; i++ )
37 {
38 if( !tmp[i] )
39 {
40 *ptr++ = i;
41 }
42 }
43 *ptr = -1;
44 free( tmp );
45 }
46
negativeMember(int * query,int locnseq)47 int *negativeMember( int *query, int locnseq )
48 {
49 int *bk, *value = NULL;
50 char *tmp;
51 int i;
52 int n;
53
54 tmp = AllocateCharVec( locnseq );
55 bk = value = AllocateIntVec( locnseq );
56 if( !value ) ErrorExit( "Cannot allocate value" );
57
58 for( i=0; i<locnseq; i++ ) tmp[i] = 0;
59 while( (n=*query++) != -1 ) tmp[n] = 1;
60
61 for( i=0; i<locnseq; i++ )
62 {
63 if( !tmp[i] )
64 {
65 fprintf( stderr, "%3d ", i );
66 *value++ = i;
67 }
68 }
69 fprintf( stderr, "\n" );
70 *value = -1;
71 free( tmp );
72 return( bk );
73 }
74
IntExistsInVec(int query,int * vector)75 int IntExistsInVec( int query, int *vector )
76 {
77 while( *vector != -1 )
78 if( query == *vector++ ) return( 1 );
79 return( 0 );
80 }
81
searchParent(int top,int *** topol,int Start,int End)82 NodeInCub searchParent( int top, int ***topol, int Start, int End )
83 {
84 int i, j;
85 NodeInCub value;
86 for( i=Start; i<End; i++ )
87 {
88 for( j=0; j<2; j++ )
89 {
90 if( IntExistsInVec( top, topol[i][j] ) )
91 {
92 value.step = i;
93 value.LorR = j;
94 return( value );
95 }
96 }
97 }
98 fprintf( stderr, "ERROR!!!\n" );
99 ErrorExit( "Error in searchParent" );
100 value.step=0; // by D.Mathog, katoh
101 value.LorR=0; // by D.Mathog, katoh
102 return( value );
103 }
104
stopolInit(int n,Node * stopol)105 void stopolInit( int n, Node *stopol )
106 {
107 int i, j;
108 for( i=0; i<n; i++ )
109 {
110 for( j=0; j<3; j++ )
111 {
112 stopol[i].length[j] = 0.0;
113 stopol[i].children[j] = NULL;
114 stopol[i].tmpChildren[j] = -1;
115 stopol[i].top[j] = -1;
116 stopol[i].members[j] = NULL;
117 stopol[i].weightptr[j] = NULL;
118 }
119 }
120 #if 0
121 while( --numintvec >= 0 )
122 {
123 free( tmpintvec[numintvec] );
124 }
125 free( tmpintvec );
126 numintvec = 0;
127 #endif
128 }
129
treeCnv(Node * stopol,int locnseq,int *** topol,double ** len,double ** bw)130 void treeCnv( Node *stopol, int locnseq, int ***topol, double **len, double **bw )
131 {
132 int i;
133 NodeInCub parent;
134 int *count;
135 int ccount;
136 int rep;
137 int tmpint;
138 static int **tmpintvec = NULL;
139 static int numintvec = 0;
140
141 count = AllocateIntVec( 2 * locnseq ); /* oome */
142 if( !count ) ErrorExit( "Cannot allocate count.\n" );
143
144 checkMinusLength( locnseq, len ); /* uwagaki */
145
146 stopolInit( locnseq * 2, stopol );
147 for( i=0; i<locnseq * 2; i++ ) count[i] = 0;
148
149 for( i=locnseq; i<locnseq*2; i++ )
150 {
151 rep = i - locnseq;
152 parent = searchParent( rep, topol, 0, locnseq-1 );
153 #if DEBUG
154 fprintf( stderr, "Parent of node No.%d ( Seq No.%d ) = %d - %d\n", i, i-locnseq, parent.step, parent.LorR );
155 #endif
156
157 ccount = count[parent.step];
158 stopol[parent.step].length[ccount] = len[parent.step][parent.LorR];
159 stopol[parent.step].weightptr[ccount] = &(bw[parent.step][parent.LorR]);
160 stopol[parent.step].children[ccount] = &stopol[i];
161 stopol[parent.step].tmpChildren[ccount] = i;
162 stopol[parent.step].members[ccount] = topol[parent.step][parent.LorR];
163 count[parent.step]++;
164
165 ccount = count[i];
166 stopol[i].length[ccount] = len[parent.step][parent.LorR];
167 stopol[i].weightptr[ccount] = &(bw[parent.step][parent.LorR]);
168 stopol[i].children[ccount] = &stopol[parent.step];
169 stopol[i].tmpChildren[ccount] = parent.step;
170 stopol[i].members[ccount] = topol[parent.step][parent.LorR];
171 count[i]++;
172 }
173 for( i=0; i<locnseq-2; i++ )
174 {
175 rep = MIN( topol[i][0][0], topol[i][1][0] );
176 parent = searchParent( rep, topol, i+1, locnseq-1 );
177 ccount = count[parent.step];
178 stopol[parent.step].length[ccount] = len[parent.step][parent.LorR];
179 stopol[parent.step].weightptr[ccount] = &(bw[parent.step][parent.LorR]);
180 stopol[parent.step].children[ccount] = &stopol[i];
181 stopol[parent.step].tmpChildren[ccount] = i;
182 stopol[parent.step].members[ccount] = topol[parent.step][parent.LorR];
183 count[parent.step]++;
184
185 ccount = count[i];
186 stopol[i].length[ccount] = len[parent.step][parent.LorR];
187 stopol[i].weightptr[ccount] = &(bw[parent.step][parent.LorR]);
188 stopol[i].children[ccount] = &stopol[parent.step];
189 stopol[i].tmpChildren[ccount] = parent.step;
190 #if 0
191 stopol[i].members[ccount] = negativeMember( topol[parent.step][parent.LorR], locnseq );
192 #else
193 // fprintf( stderr, "allocating numintvec = %d\n", numintvec );
194 tmpintvec = (int **)realloc( (void *)tmpintvec, (numintvec+1) * sizeof( int * ) );
195 tmpintvec[numintvec] = (int *)calloc( locnseq, sizeof( int ) );
196 negativeMember2( tmpintvec[numintvec], topol[parent.step][parent.LorR], locnseq );
197 stopol[i].members[ccount] = tmpintvec[numintvec];
198 numintvec++;
199 #endif
200 count[i]++;
201 #if DEBUG
202 fprintf( stderr, "Parent of node No.%d = %d - %d\n", i, parent.step, parent.LorR );
203 #endif
204 }
205 /*
206 Unrooted tree.
207 locnseq-2 no children no nakade,
208 locnseq-3 wo sashiteinai mono wo sagashite,
209 locnseq-3 no children ni kuwae,
210 locnseq-2 wo sashiteita node no chilren wo
211 locnseq-3 ni kaeru.
212 */
213 #if DEBUG
214 fprintf( stderr, "BEFORE MODIFY\n" );
215 for( i=0; i<locnseq*2; i++ )
216 {
217 for( j=0; j<3; j++ )
218 {
219 fprintf( stderr, "stopol[%d].tmpChildren[%d] = %d, children[%d] = %d \n", i, j, stopol[i].tmpChildren[j], j, stopol[i].children[j] - stopol );
220 }
221 }
222 #endif
223
224 if ( stopol[locnseq-2].children[0] == &stopol[locnseq-3] ) i = 1;
225 else if( stopol[locnseq-2].children[1] == &stopol[locnseq-3] ) i = 0;
226 else ErrorExit( "?\n" );
227
228 stopol[locnseq-3].length[2] = len[locnseq-2][0] + len[locnseq-2][1];
229 stopol[locnseq-3].weightptr[2] = &bw[locnseq-2][0];
230 stopol[locnseq-3].children[2] = stopol[locnseq-2].children[i];
231 stopol[locnseq-3].tmpChildren[2] = stopol[locnseq-2].tmpChildren[i];
232
233 tmpint = (int)( stopol[locnseq-2].children[i] - stopol );
234
235 stopol[tmpint].children[2] = &stopol[locnseq-3];
236 stopol[tmpint].length[2] = len[locnseq-2][0] + len[locnseq-2][1];
237 stopol[tmpint].weightptr[2] = &bw[locnseq-2][0];
238 stopol[tmpint].tmpChildren[2] = locnseq-3;
239
240
241 #if DEBUG
242 for( i=0; i<locnseq*2; i++ )
243 {
244 for( j=0; j<3; j++ )
245 {
246 fprintf( stderr, "stopol[%d].tmpChildren[%d] = %d, children[%d] = %d \n", i, j, stopol[i].tmpChildren[j], j, stopol[i].children[j] - stopol );
247 }
248 }
249
250 for( i=0; i<locnseq*2; i++ )
251 {
252 fprintf( stderr, "-- stopol[%d]\n", i );
253 for( j=0; j<3; j++ )
254 {
255 if( !stopol[i].members[j] )
256 {
257 fprintf( stderr, "LEAF\n" );
258 break;
259 }
260 fprintf( stderr, " group %d are \n", j );
261 for( k=0; (n=stopol[i].members[j][k]) != -1; k++ )
262 {
263 fprintf( stderr, "%#5d", n );
264 }
265 fprintf( stderr, "\n" );
266 }
267 fprintf( stderr, "\n" );
268 }
269 #endif
270
271 #if DEBUG
272 stopol_g = stopol;
273 #endif
274 free( count );
275 }
276
isLeaf(Node node)277 int isLeaf( Node node )
278 {
279 if( node.children[1] ) return( 0 );
280 else return( 1 );
281 }
282
syntheticLength(Node * ob,Node * oppositeNode)283 double syntheticLength( Node *ob, Node *oppositeNode )
284 {
285 int i, count;
286 int dir_ch[3];
287 int dir_pa = -10; // by katoh
288 double value, tmpvalue0, tmpvalue1;
289 int nanflag = 0;
290
291 #if DEBUG
292 fprintf( stderr, "In syntheticLength\n" );
293 fprintf( stderr, "ob - stopol_g = %d\n", ob - stopol_g );
294 fprintf( stderr, "op - stopol_g = %d\n", oppositeNode - stopol_g );
295 #endif
296
297 if( isLeaf( *ob ) )
298 {
299 #if DEBUG
300 fprintf( stderr, "LEAF\n\n" );
301 #endif
302 return( ob->length[0] );
303 }
304
305 for( i=0, count=0; i<3; i++ )
306 {
307 #if DEBUG
308 fprintf( stderr, "ob->tmpChildren[%d] = %d\n", i, ob->tmpChildren[i] );
309 #endif
310 if( oppositeNode != ob->children[i] ) dir_ch[count++] = i;
311 else dir_pa = i;
312 }
313 #if DEBUG
314 fprintf( stderr, "\n" );
315 #endif
316 if( count != 2 )
317 {
318 #if DEBUG
319 fprintf( stderr, "Node No.%d has no child like No.%d \n", ob-stopol_g, oppositeNode-stopol_g );
320 #endif
321 ErrorExit( "Invalid call\n" );
322 }
323
324 tmpvalue0 = syntheticLength( ob->children[dir_ch[0]], ob );
325 tmpvalue1 = syntheticLength( ob->children[dir_ch[1]], ob );
326
327 #if DEBUG
328 fprintf( stderr, "tmpvalue0 = %f\n", tmpvalue0 );
329 fprintf( stderr, "tmpvalue1 = %f\n", tmpvalue1 );
330 #endif
331 if( tmpvalue0 ) tmpvalue0 = 1.0 / tmpvalue0;
332 else nanflag = 1;
333 if( tmpvalue1 ) tmpvalue1 = 1.0 / tmpvalue1;
334 else nanflag = 1;
335
336 if( nanflag ) value = 0.0;
337 else
338 {
339 value = tmpvalue0 + tmpvalue1;
340 value = 1.0 / value;
341 }
342 value += ob->length[dir_pa];
343 #if DEBUG
344 fprintf( stderr, "value = %f\n", value );
345 #endif
346
347 return( value );
348 }
349
calcW(Node * ob,Node * op)350 double calcW( Node *ob, Node *op )
351 {
352 int i, count;
353 int dir_ch[3];
354 int dir_pa = -10; // by katoh
355 double a, b, c, f, s;
356 double value;
357
358 if( isLeaf( *ob ) )
359 return( 1.0 );
360
361 for( i=0, count=0; i<3; i++ )
362 {
363 if( op != ob->children[i] ) dir_ch[count++] = i;
364 else dir_pa = i;
365 }
366 if( count != 2 ) ErrorExit( "Invalid call of calcW\n" );
367
368 #if DEBUG
369 fprintf( stderr, "In calcW\n" );
370 fprintf( stderr, "ob = %d\n", ob - stopol_g );
371 fprintf( stderr, "op = %d\n", op - stopol_g );
372 fprintf( stderr, "ob->children[c1] = %d\n", ob->children[dir_ch[0]] - stopol_g );
373 fprintf( stderr, "ob->children[c2] = %d\n", ob->children[dir_ch[1]] - stopol_g );
374 fprintf( stderr, "ob->children[pa] = %d\n", ob->children[dir_pa] - stopol_g );
375 fprintf( stderr, "\n" );
376 #endif
377
378 a = syntheticLength( ob->children[dir_ch[0]], ob );
379 b = syntheticLength( ob->children[dir_ch[1]], ob );
380 c = syntheticLength( ob->children[dir_pa], ob );
381
382 #if DEBUG
383 fprintf( stderr, "a = %f\n", a );
384 fprintf( stderr, "b = %f\n", b );
385 fprintf( stderr, "c = %f\n", c );
386 #endif
387
388 if( !c ) return( MAXBW );
389 if ( !a || !b ) return( MINBW ); /* ? */
390
391 f = EF_THREEWAY;
392 s = ( b*c + c*a + a*b );
393
394 value = a*b*(c+a)*(c+b) / ( c*(a+b) * f * s );
395
396 value = sqrt( value );
397
398 return( value );
399 }
400
calcBranchWeight(double ** bw,int locnseq,Node * stopol,int *** topol,double ** len)401 void calcBranchWeight( double **bw, int locnseq, Node *stopol, int ***topol, double **len )
402 {
403 NodeInCub parent;
404 int i;
405 int rep;
406 Node *topNode, *btmNode;
407 double topW, btmW;
408
409 for( i=locnseq; i<locnseq*2; i++ )
410 {
411 rep = i - locnseq;
412 parent = searchParent( rep, topol, 0, locnseq-1 );
413 if( parent.step == locnseq - 2 ) continue;
414
415 topNode = stopol+parent.step; btmNode = stopol+i;
416 #if DEBUG
417 fprintf( stderr, "In calcBranchWeight, topNode=%d, btmNode=%d\n", topNode-stopol_g, btmNode-stopol_g );
418 #endif
419 topW = calcW( topNode, btmNode );
420 btmW = calcW( btmNode, topNode );
421 bw[parent.step][parent.LorR] = topW * btmW;
422 }
423 for( i=0; i<locnseq-3; i++ )
424 {
425 rep = MIN( topol[i][0][0], topol[i][1][0] );
426 parent = searchParent( rep, topol, i+1, locnseq-1 );
427 if( parent.step == locnseq - 2 ) continue;
428 topNode = stopol+parent.step;
429 btmNode = stopol+i;
430 #if DEBUG
431 fprintf( stderr, "In calcBranchWeight, topNode=%d, btmNode=%d\n", topNode-stopol_g, btmNode-stopol_g );
432 #endif
433 topW = calcW( topNode, btmNode );
434 btmW = calcW( btmNode, topNode );
435 bw[parent.step][parent.LorR] = topW * btmW;
436 }
437
438 topNode = stopol[locnseq-3].children[2];
439 btmNode = stopol + i;
440 topW = calcW( topNode, btmNode );
441 btmW = calcW( btmNode, topNode );
442 bw[locnseq-2][0] = topW * btmW;
443 bw[locnseq-2][1] = 1.0;
444 }
445
branchWeightToPairWeight(int locnseq,int *** topol,double ** pw,double ** bw)446 void branchWeightToPairWeight( int locnseq, int ***topol, double **pw, double **bw )
447 {
448 int i, j, k, n0, n1;
449 #if 0
450 double wFromLeaf[locnseq];
451 #else
452 static double *wFromLeaf = NULL;
453 if( wFromLeaf == NULL )
454 wFromLeaf = AllocateDoubleVec( locnseq );
455 #endif
456
457 #if DEBUG
458 for( i=0; i<locnseq-1; i++ ) for( j=0; j<2; j++ )
459 fprintf( stderr, "pw[%d][%d] = %f\n", i, j, bw[i][j] );
460 #endif
461
462 for( i=0; i<locnseq; i++ ) wFromLeaf[i] = 1.0;
463 for( i=0; i<locnseq; i++ ) for( j=0; j<locnseq; j++ )
464 pw[i][j] = 0.0;
465 for( i=0; i<locnseq-1; i++ )
466 {
467 for( j=0; (n0=topol[i][0][j])!=-1; j++ )
468 for( k=0; (n1=topol[i][1][k])!=-1; k++ )
469 pw[MIN( n0, n1 )][MAX( n0, n1 )]
470 = wFromLeaf[n0] * wFromLeaf[n1] * bw[i][0] * bw[i][1];
471 for( j=0; (n0=topol[i][0][j])!=-1; j++ )
472 wFromLeaf[n0] *= bw[i][0];
473 for( j=0; (n1=topol[i][1][j])!=-1; j++ )
474 wFromLeaf[n1] *= bw[i][1];
475 }
476 }
477
distFromABranch_rec(double * result,Node * ob,Node * op)478 static void distFromABranch_rec( double *result, Node *ob, Node *op )
479 {
480 int i, n, count;
481 int dir_ch[3], dir_pa;
482
483 #if DEBUG
484 fprintf( stderr, "In distFromABranch_rec, ob = %d\n", ob - stopol_g );
485 #endif
486 if( isLeaf( *ob ) ) return;
487 for( i=0, count=0; i<3; i++ )
488 {
489 if( ob->children[i] != op ) dir_ch[count++] = i;
490 else dir_pa = i;
491 }
492 if( count != 2 )
493 {
494 #if DEBUG
495 fprintf( stderr, "Node No.%d has no child like No.%d \n", ob-stopol_g, op-stopol_g );
496 #endif
497 ErrorExit( "Incorrect call of distFromABranch_rec" );
498 }
499 for( i=0; (n=ob->members[dir_ch[0]][i])!=-1; i++ )
500 {
501 result[n] += ob->length[dir_ch[0]];
502 }
503 distFromABranch_rec( result, ob->children[dir_ch[0]], ob );
504
505 for( i=0; (n=ob->members[dir_ch[1]][i])!=-1; i++ )
506 {
507 result[n] += ob->length[dir_ch[1]];
508 }
509 distFromABranch_rec( result, ob->children[dir_ch[1]], ob );
510 }
511
distFromABranch(int nseq,double * result,Node * stopol,int *** topol,double ** len,int step,int LorR)512 void distFromABranch( int nseq, double *result, Node *stopol, int ***topol, double **len, int step, int LorR )
513 {
514 Node *topNode, *btmNode;
515 int i;
516
517 if( nseq == 2 )
518 {
519 result[0] = len[0][0];
520 result[1] = len[0][1];
521 // reporterr( "result[0] = %f\n", result[0] );
522 // reporterr( "result[1] = %f\n", result[1] );
523 return;
524 }
525
526 if( step == nseq - 2 )
527 {
528 topNode = stopol[nseq-2].children[0];
529 btmNode = stopol + nseq-3;
530 #if DEBUG
531 fprintf( stderr, "Now step == nseq-3, topNode = %d, btmNode = %d\n", topNode - stopol_g, btmNode-stopol_g );
532 #endif
533 }
534
535 else
536 {
537 for( i=0; i<3; i++ )
538 {
539 if( stopol[step].members[i][0] == topol[step][LorR][0] )
540 break;
541 }
542 if( i== 3 ) ErrorExit( "Incorrect call of distFromABranch." );
543 btmNode = stopol[step].children[i];
544 topNode = stopol+step;
545 }
546
547 for( i=0; i<nseq; i++ ) result[i] = 0.0;
548 distFromABranch_rec( result, btmNode, topNode );
549 distFromABranch_rec( result, topNode, btmNode );
550 #if 0
551 for( i=0; i<nseq; i++ )
552 fprintf( stdout, "w[%d] = %30.20f\n", i, result[i] );
553 #endif
554 // fprintf( stderr, "new weight!\n" );
555 // for( i=0; i<nseq; i++ )
556 // result[i] *= result[i];
557
558
559 }
560
weightFromABranch_rec(double * result,Node * ob,Node * op)561 static void weightFromABranch_rec( double *result, Node *ob, Node *op )
562 {
563 int i, n, count;
564 int dir_ch[3], dir_pa;
565
566
567 #if DEBUG
568 fprintf( stderr, "In weightFromABranch_rec, ob = %d\n", ob - stopol_g );
569 #endif
570 if( isLeaf( *ob ) ) return;
571 for( i=0, count=0; i<3; i++ )
572 {
573 if( ob->children[i] != op ) dir_ch[count++] = i;
574 else dir_pa = i;
575 }
576 if( count != 2 )
577 {
578 #if DEBUG
579 fprintf( stderr, "Node No.%d has no child like No.%d \n", ob-stopol_g, op-stopol_g );
580 #endif
581 ErrorExit( "Incorrect call of weightFromABranch_rec" );
582 }
583 for( i=0; (n=ob->members[dir_ch[0]][i])!=-1; i++ )
584 result[n] *= *ob->weightptr[dir_ch[0]];
585 weightFromABranch_rec( result, ob->children[dir_ch[0]], ob );
586
587 for( i=0; (n=ob->members[dir_ch[1]][i])!=-1; i++ )
588 result[n] *= *ob->weightptr[dir_ch[1]];
589 weightFromABranch_rec( result, ob->children[dir_ch[1]], ob );
590 }
591
weightFromABranch(int nseq,double * result,Node * stopol,int *** topol,int step,int LorR)592 void weightFromABranch( int nseq, double *result, Node *stopol, int ***topol, int step, int LorR )
593 {
594 Node *topNode, *btmNode;
595 int i;
596
597 if( step == nseq - 2 )
598 {
599 topNode = stopol[nseq-2].children[0];
600 btmNode = stopol + nseq-3;
601 #if DEBUG
602 fprintf( stderr, "Now step == nseq-3, topNode = %d, btmNode = %d\n", topNode - stopol_g, btmNode-stopol_g );
603 #endif
604 }
605
606 else
607 {
608 for( i=0; i<3; i++ )
609 {
610 if( stopol[step].members[i][0] == topol[step][LorR][0] )
611 break;
612 }
613 if( i== 3 ) ErrorExit( "Incorrect call of weightFromABranch." );
614 btmNode = stopol[step].children[i];
615 topNode = stopol+step;
616 }
617
618 for( i=0; i<nseq; i++ ) result[i] = 1.0;
619 weightFromABranch_rec( result, btmNode, topNode );
620 weightFromABranch_rec( result, topNode, btmNode );
621 #if 0
622 for( i=0; i<nseq; i++ )
623 fprintf( stdout, "w[%d] = %30.20f\n", i, result[i] );
624 #endif
625 // fprintf( stderr, "new weight!\n" );
626 // for( i=0; i<nseq; i++ )
627 // result[i] *= result[i];
628
629
630 }
assignstrweight_rec(double * strweight,Node * ob,Node * op,char * kozoari,double * seqweight)631 void assignstrweight_rec( double *strweight, Node *ob, Node *op, char *kozoari, double *seqweight )
632 {
633 int i, n, count, lastkozo;
634 int dir_ch[3], dir_pa;
635 double sumweight;
636
637 #if DEBUG
638 fprintf( stderr, "In weightFromABranch_rec, ob = %d\n", ob - stopol_g );
639 #endif
640 if( isLeaf( *ob ) )
641 {
642 // fprintf( stderr, "Leaf!\n" );
643 return;
644 }
645 for( i=0, count=0; i<3; i++ )
646 {
647 if( ob->children[i] != op ) dir_ch[count++] = i;
648 else dir_pa = i;
649 }
650 if( count != 2 )
651 {
652 #if DEBUG
653 fprintf( stderr, "Node No.%d has no child like No.%d \n", ob-stopol_g, op-stopol_g );
654 #endif
655 ErrorExit( "Incorrect call of weightFromABranch_rec" );
656 }
657
658
659 // fprintf( stderr, "\n" );
660 sumweight = 0.0;
661 count = 0;
662 lastkozo = -1;
663 for( i=0; (n=ob->members[dir_ch[0]][i])!=-1; i++ )
664 {
665 // fprintf( stderr, "member1! n=%d\n", n );
666 sumweight += seqweight[n];
667 if( kozoari[n] )
668 {
669 count++;
670 lastkozo = n;
671 }
672 }
673 for( i=0; (n=ob->members[dir_ch[1]][i])!=-1; i++ )
674 {
675 // fprintf( stderr, "member2! n=%d\n", n );
676 sumweight += seqweight[n];
677 if( kozoari[n] )
678 {
679 count++;
680 lastkozo = n;
681 }
682 }
683
684 // fprintf( stderr, "count = %d\n", count );
685
686 if( count == 1 )
687 strweight[lastkozo] = sumweight;
688 else if( count > 1 )
689 {
690 assignstrweight_rec( strweight, ob->children[dir_ch[0]], ob, kozoari, seqweight );
691 assignstrweight_rec( strweight, ob->children[dir_ch[1]], ob, kozoari, seqweight );
692 }
693 }
694
assignstrweight(int nseq,double * strweight,Node * stopol,int *** topol,int step,int LorR,char * kozoari,double * seqweight)695 void assignstrweight( int nseq, double *strweight, Node *stopol, int ***topol, int step, int LorR, char *kozoari, double *seqweight )
696 {
697 Node *topNode, *btmNode;
698 int i;
699
700 if( step == nseq - 2 )
701 {
702 topNode = stopol[nseq-2].children[0];
703 btmNode = stopol + nseq-3;
704 #if DEBUG
705 fprintf( stderr, "Now step == nseq-3, topNode = %d, btmNode = %d\n", topNode - stopol_g, btmNode-stopol_g );
706 #endif
707 }
708
709 else
710 {
711 for( i=0; i<3; i++ )
712 {
713 if( stopol[step].members[i][0] == topol[step][LorR][0] )
714 break;
715 }
716 if( i== 3 ) ErrorExit( "Incorrect call of weightFromABranch." );
717 btmNode = stopol[step].children[i];
718 topNode = stopol+step;
719 }
720
721 for( i=0; i<nseq; i++ ) strweight[i] = 0.0;
722 for( i=0; i<nseq; i++ ) if( kozoari[i] ) strweight[i] = seqweight[i];
723 // fprintf( stderr, "calling _rec (1)\n" );
724 assignstrweight_rec( strweight, btmNode, topNode, kozoari, seqweight );
725 // fprintf( stderr, "calling _rec (2)\n" );
726 assignstrweight_rec( strweight, topNode, btmNode, kozoari, seqweight );
727
728 #if 1 // nazeka kokowo tobasuto seido ga sagaru ?????
729 fprintf( stderr, "STEP %d\n", step );
730 for( i=0; topol[step][0][i]>-1; i++ )
731 fprintf( stderr, "%3d ", topol[step][0][i] );
732 fprintf( stderr, "\n" );
733 for( i=0; topol[step][1][i]>-1; i++ )
734 fprintf( stderr, "%3d ", topol[step][1][i] );
735 fprintf( stderr, "\n" );
736 for( i=0; i<nseq; i++ )
737 fprintf( stderr, "seqweight[%d] = %f\n", i, seqweight[i] );
738 for( i=0; i<nseq; i++ )
739 fprintf( stderr, "strweight[%d] = %f\n", i, strweight[i] );
740 fprintf( stderr, "\n" );
741 #endif
742 }
743