5 #define EF_THREEWAY 1.0
16 void checkMinusLength( int nseq, double **len )
19 for( i=0; i<nseq-1; i++ ) for( j=0; j<2; j++ )
20 if( len[i][j] < MINLEN ) len[i][j] = MINLEN;
23 void negativeMember2( int *mem, int *query, int locnseq )
30 tmp = AllocateCharVec( locnseq );
32 for( i=0; i<locnseq; i++ ) tmp[i] = 0;
33 while( (n=*query++) != -1 ) tmp[n] = 1;
36 for( i=0; i<locnseq; i++ )
47 int *negativeMember( int *query, int locnseq )
49 int *bk, *value = NULL;
54 tmp = AllocateCharVec( locnseq );
55 bk = value = AllocateIntVec( locnseq );
56 if( !value ) ErrorExit( "Cannot allocate value" );
58 for( i=0; i<locnseq; i++ ) tmp[i] = 0;
59 while( (n=*query++) != -1 ) tmp[n] = 1;
61 for( i=0; i<locnseq; i++ )
65 fprintf( stderr, "%3d ", i );
69 fprintf( stderr, "\n" );
75 int IntExistsInVec( int query, int *vector )
77 while( *vector != -1 )
78 if( query == *vector++ ) return( 1 );
82 NodeInCub searchParent( int top, int ***topol, int Start, int End )
86 for( i=Start; i<End; i++ )
90 if( IntExistsInVec( top, topol[i][j] ) )
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
105 void stopolInit( int n, Node *stopol )
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;
121 while( --numintvec >= 0 )
123 free( tmpintvec[numintvec] );
130 void treeCnv( Node *stopol, int locnseq, int ***topol, double **len, double **bw )
138 static int **tmpintvec = NULL;
139 static int numintvec = 0;
141 count = AllocateIntVec( 2 * locnseq ); /* oome */
142 if( !count ) ErrorExit( "Cannot allocate count.\n" );
144 checkMinusLength( locnseq, len ); /* uwagaki */
146 stopolInit( locnseq * 2, stopol );
147 for( i=0; i<locnseq * 2; i++ ) count[i] = 0;
149 for( i=locnseq; i<locnseq*2; i++ )
152 parent = searchParent( rep, topol, 0, locnseq-1 );
154 fprintf( stderr, "Parent of node No.%d ( Seq No.%d ) = %d - %d\n", i, i-locnseq, parent.step, parent.LorR );
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]++;
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];
173 for( i=0; i<locnseq-2; i++ )
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]++;
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;
191 stopol[i].members[ccount] = negativeMember( topol[parent.step][parent.LorR], locnseq );
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];
202 fprintf( stderr, "Parent of node No.%d = %d - %d\n", i, parent.step, parent.LorR );
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
214 fprintf( stderr, "BEFORE MODIFY\n" );
215 for( i=0; i<locnseq*2; i++ )
219 fprintf( stderr, "stopol[%d].tmpChildren[%d] = %d, children[%d] = %d \n", i, j, stopol[i].tmpChildren[j], j, stopol[i].children[j] - stopol );
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" );
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];
233 tmpint = (int)( stopol[locnseq-2].children[i] - stopol );
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;
242 for( i=0; i<locnseq*2; i++ )
246 fprintf( stderr, "stopol[%d].tmpChildren[%d] = %d, children[%d] = %d \n", i, j, stopol[i].tmpChildren[j], j, stopol[i].children[j] - stopol );
250 for( i=0; i<locnseq*2; i++ )
252 fprintf( stderr, "-- stopol[%d]\n", i );
255 if( !stopol[i].members[j] )
257 fprintf( stderr, "LEAF\n" );
260 fprintf( stderr, " group %d are \n", j );
261 for( k=0; (n=stopol[i].members[j][k]) != -1; k++ )
263 fprintf( stderr, "%#5d", n );
265 fprintf( stderr, "\n" );
267 fprintf( stderr, "\n" );
277 int isLeaf( Node node )
279 if( node.children[1] ) return( 0 );
283 double syntheticLength( Node *ob, Node *oppositeNode )
287 int dir_pa = -10; // by katoh
288 double value, tmpvalue0, tmpvalue1;
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 );
300 fprintf( stderr, "LEAF\n\n" );
302 return( ob->length[0] );
305 for( i=0, count=0; i<3; i++ )
308 fprintf( stderr, "ob->tmpChildren[%d] = %d\n", i, ob->tmpChildren[i] );
310 if( oppositeNode != ob->children[i] ) dir_ch[count++] = i;
314 fprintf( stderr, "\n" );
319 fprintf( stderr, "Node No.%d has no child like No.%d \n", ob-stopol_g, oppositeNode-stopol_g );
321 ErrorExit( "Invalid call\n" );
324 tmpvalue0 = syntheticLength( ob->children[dir_ch[0]], ob );
325 tmpvalue1 = syntheticLength( ob->children[dir_ch[1]], ob );
328 fprintf( stderr, "tmpvalue0 = %f\n", tmpvalue0 );
329 fprintf( stderr, "tmpvalue1 = %f\n", tmpvalue1 );
331 if( tmpvalue0 ) tmpvalue0 = 1.0 / tmpvalue0;
333 if( tmpvalue1 ) tmpvalue1 = 1.0 / tmpvalue1;
336 if( nanflag ) value = 0.0;
339 value = tmpvalue0 + tmpvalue1;
342 value += ob->length[dir_pa];
344 fprintf( stderr, "value = %f\n", value );
350 double calcW( Node *ob, Node *op )
354 int dir_pa = -10; // by katoh
355 double a, b, c, f, s;
361 for( i=0, count=0; i<3; i++ )
363 if( op != ob->children[i] ) dir_ch[count++] = i;
366 if( count != 2 ) ErrorExit( "Invalid call of calcW\n" );
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" );
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 );
383 fprintf( stderr, "a = %f\n", a );
384 fprintf( stderr, "b = %f\n", b );
385 fprintf( stderr, "c = %f\n", c );
388 if( !c ) return( MAXBW );
389 if ( !a || !b ) return( MINBW ); /* ? */
392 s = ( b*c + c*a + a*b );
394 value = a*b*(c+a)*(c+b) / ( c*(a+b) * f * s );
396 value = sqrt( value );
401 void calcBranchWeight( double **bw, int locnseq, Node *stopol, int ***topol, double **len )
406 Node *topNode, *btmNode;
409 for( i=locnseq; i<locnseq*2; i++ )
412 parent = searchParent( rep, topol, 0, locnseq-1 );
413 if( parent.step == locnseq - 2 ) continue;
415 topNode = stopol+parent.step; btmNode = stopol+i;
417 fprintf( stderr, "In calcBranchWeight, topNode=%d, btmNode=%d\n", topNode-stopol_g, btmNode-stopol_g );
419 topW = calcW( topNode, btmNode );
420 btmW = calcW( btmNode, topNode );
421 bw[parent.step][parent.LorR] = topW * btmW;
423 for( i=0; i<locnseq-3; i++ )
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;
431 fprintf( stderr, "In calcBranchWeight, topNode=%d, btmNode=%d\n", topNode-stopol_g, btmNode-stopol_g );
433 topW = calcW( topNode, btmNode );
434 btmW = calcW( btmNode, topNode );
435 bw[parent.step][parent.LorR] = topW * btmW;
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;
446 void branchWeightToPairWeight( int locnseq, int ***topol, double **pw, double **bw )
450 double wFromLeaf[locnseq];
452 static double *wFromLeaf = NULL;
453 if( wFromLeaf == NULL )
454 wFromLeaf = AllocateDoubleVec( locnseq );
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] );
462 for( i=0; i<locnseq; i++ ) wFromLeaf[i] = 1.0;
463 for( i=0; i<locnseq; i++ ) for( j=0; j<locnseq; j++ )
465 for( i=0; i<locnseq-1; i++ )
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];
478 void weightFromABranch_rec( double *result, Node *ob, Node *op )
481 int dir_ch[3], dir_pa;
484 fprintf( stderr, "In weightFromABranch_rec, ob = %d\n", ob - stopol_g );
486 if( isLeaf( *ob ) ) return;
487 for( i=0, count=0; i<3; i++ )
489 if( ob->children[i] != op ) dir_ch[count++] = i;
495 fprintf( stderr, "Node No.%d has no child like No.%d \n", ob-stopol_g, op-stopol_g );
497 ErrorExit( "Incorrect call of weightFromABranch_rec" );
499 for( i=0; (n=ob->members[dir_ch[0]][i])!=-1; i++ )
500 result[n] *= *ob->weightptr[dir_ch[0]];
501 weightFromABranch_rec( result, ob->children[dir_ch[0]], ob );
503 for( i=0; (n=ob->members[dir_ch[1]][i])!=-1; i++ )
504 result[n] *= *ob->weightptr[dir_ch[1]];
505 weightFromABranch_rec( result, ob->children[dir_ch[1]], ob );
508 void weightFromABranch( int nseq, double *result, Node *stopol, int ***topol, int step, int LorR )
510 Node *topNode, *btmNode;
513 if( step == nseq - 2 )
515 topNode = stopol[nseq-2].children[0];
516 btmNode = stopol + nseq-3;
518 fprintf( stderr, "Now step == nseq-3, topNode = %d, btmNode = %d\n", topNode - stopol_g, btmNode-stopol_g );
526 if( stopol[step].members[i][0] == topol[step][LorR][0] )
529 if( i== 3 ) ErrorExit( "Incorrect call of weightFromABranch." );
530 btmNode = stopol[step].children[i];
531 topNode = stopol+step;
534 for( i=0; i<nseq; i++ ) result[i] = 1.0;
535 weightFromABranch_rec( result, btmNode, topNode );
536 weightFromABranch_rec( result, topNode, btmNode );
538 for( i=0; i<nseq; i++ )
539 fprintf( stdout, "w[%d] = %f\n", i, result[i] );
541 // fprintf( stderr, "new weight!\n" );
542 // for( i=0; i<nseq; i++ )
543 // result[i] *= result[i];
547 void assignstrweight_rec( double *strweight, Node *ob, Node *op, char *kozoari, double *seqweight )
549 int i, n, count, lastkozo;
550 int dir_ch[3], dir_pa;
554 fprintf( stderr, "In weightFromABranch_rec, ob = %d\n", ob - stopol_g );
558 // fprintf( stderr, "Leaf!\n" );
561 for( i=0, count=0; i<3; i++ )
563 if( ob->children[i] != op ) dir_ch[count++] = i;
569 fprintf( stderr, "Node No.%d has no child like No.%d \n", ob-stopol_g, op-stopol_g );
571 ErrorExit( "Incorrect call of weightFromABranch_rec" );
575 // fprintf( stderr, "\n" );
579 for( i=0; (n=ob->members[dir_ch[0]][i])!=-1; i++ )
581 // fprintf( stderr, "member1! n=%d\n", n );
582 sumweight += seqweight[n];
589 for( i=0; (n=ob->members[dir_ch[1]][i])!=-1; i++ )
591 // fprintf( stderr, "member2! n=%d\n", n );
592 sumweight += seqweight[n];
600 // fprintf( stderr, "count = %d\n", count );
603 strweight[lastkozo] = sumweight;
606 assignstrweight_rec( strweight, ob->children[dir_ch[0]], ob, kozoari, seqweight );
607 assignstrweight_rec( strweight, ob->children[dir_ch[1]], ob, kozoari, seqweight );
611 void assignstrweight( int nseq, double *strweight, Node *stopol, int ***topol, int step, int LorR, char *kozoari, double *seqweight )
613 Node *topNode, *btmNode;
616 if( step == nseq - 2 )
618 topNode = stopol[nseq-2].children[0];
619 btmNode = stopol + nseq-3;
621 fprintf( stderr, "Now step == nseq-3, topNode = %d, btmNode = %d\n", topNode - stopol_g, btmNode-stopol_g );
629 if( stopol[step].members[i][0] == topol[step][LorR][0] )
632 if( i== 3 ) ErrorExit( "Incorrect call of weightFromABranch." );
633 btmNode = stopol[step].children[i];
634 topNode = stopol+step;
637 for( i=0; i<nseq; i++ ) strweight[i] = 0.0;
638 for( i=0; i<nseq; i++ ) if( kozoari[i] ) strweight[i] = seqweight[i];
639 // fprintf( stderr, "calling _rec (1)\n" );
640 assignstrweight_rec( strweight, btmNode, topNode, kozoari, seqweight );
641 // fprintf( stderr, "calling _rec (2)\n" );
642 assignstrweight_rec( strweight, topNode, btmNode, kozoari, seqweight );
644 #if 1 // nazeka kokowo tobasuto seido ga sagaru ?????
645 fprintf( stderr, "STEP %d\n", step );
646 for( i=0; topol[step][0][i]>-1; i++ )
647 fprintf( stderr, "%3d ", topol[step][0][i] );
648 fprintf( stderr, "\n" );
649 for( i=0; topol[step][1][i]>-1; i++ )
650 fprintf( stderr, "%3d ", topol[step][1][i] );
651 fprintf( stderr, "\n" );
652 for( i=0; i<nseq; i++ )
653 fprintf( stderr, "seqweight[%d] = %f\n", i, seqweight[i] );
654 for( i=0; i<nseq; i++ )
655 fprintf( stderr, "strweight[%d] = %f\n", i, strweight[i] );
656 fprintf( stderr, "\n" );