binaries/src/mafft/core/treeOperation.c

   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
  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
  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
  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
  75 int IntExistsInVec( int query, int *vector )
  76 {
  77         while( *vector != -1 )
  78                 if( query == *vector++ ) return( 1 );
  79         return( 0 );
  80 }
  81
  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
 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
 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
 277 int isLeaf( Node node )
 278 {
 279         if( node.children[1] ) return( 0 );
 280         else                   return( 1 );
 281 }
 282
 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
 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
 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
 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
 478 void weightFromABranch_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 weightFromABranch_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 weightFromABranch_rec" );
 498         }
 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 );
 502
 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 );
 506 }
 507
 508 void weightFromABranch( int nseq, double *result, Node *stopol, int ***topol, int step, int LorR )
 509 {
 510         Node *topNode, *btmNode;
 511         int i;
 512
 513         if( step == nseq - 2 )
 514         {
 515                 topNode = stopol[nseq-2].children[0];
 516                 btmNode = stopol + nseq-3;
 517 #if DEBUG
 518                 fprintf( stderr, "Now step == nseq-3, topNode = %d, btmNode = %d\n", topNode - stopol_g, btmNode-stopol_g );
 519 #endif
 520         }
 521
 522         else
 523         {
 524                 for( i=0; i<3; i++ )
 525                 {
 526                         if( stopol[step].members[i][0] == topol[step][LorR][0] )
 527                         break;
 528                 }
 529                 if( i== 3 ) ErrorExit( "Incorrect call of weightFromABranch." );
 530                 btmNode = stopol[step].children[i];
 531                 topNode = stopol+step;
 532         }
 533
 534         for( i=0; i<nseq; i++ ) result[i] = 1.0;
 535         weightFromABranch_rec( result, btmNode, topNode );
 536         weightFromABranch_rec( result, topNode, btmNode );
 537 #if 0
 538         for( i=0; i<nseq; i++ )
 539                 fprintf( stdout, "w[%d] = %f\n", i, result[i] );
 540 #endif
 541 //      fprintf( stderr, "new weight!\n" );
 542 //      for( i=0; i<nseq; i++ )
 543 //              result[i] *= result[i];
 544
 545
 546 }
 547 void assignstrweight_rec( double *strweight, Node *ob, Node *op, char *kozoari, double *seqweight )
 548 {
 549         int i, n, count, lastkozo;
 550         int dir_ch[3], dir_pa;
 551         double sumweight;
 552
 553 #if DEBUG
 554         fprintf( stderr, "In weightFromABranch_rec, ob = %d\n", ob - stopol_g );
 555 #endif
 556         if( isLeaf( *ob ) )
 557         {
 558 //              fprintf( stderr, "Leaf!\n" );
 559                 return;
 560         }
 561         for( i=0, count=0; i<3; i++ )
 562         {
 563                 if( ob->children[i] != op ) dir_ch[count++] = i;
 564                 else dir_pa = i;
 565         }
 566         if( count != 2 )
 567         {
 568 #if DEBUG
 569                 fprintf( stderr, "Node No.%d has no child like No.%d \n", ob-stopol_g, op-stopol_g );
 570 #endif
 571                 ErrorExit( "Incorrect call of weightFromABranch_rec" );
 572         }
 573
 574
 575 //      fprintf( stderr, "\n" );
 576         sumweight = 0.0;
 577         count = 0;
 578         lastkozo = -1;
 579         for( i=0; (n=ob->members[dir_ch[0]][i])!=-1; i++ )
 580         {
 581 //              fprintf( stderr, "member1! n=%d\n", n );
 582                 sumweight += seqweight[n];
 583                 if( kozoari[n] )
 584                 {
 585                         count++;
 586                         lastkozo = n;
 587                 }
 588         }
 589         for( i=0; (n=ob->members[dir_ch[1]][i])!=-1; i++ )
 590         {
 591 //              fprintf( stderr, "member2! n=%d\n", n );
 592                 sumweight += seqweight[n];
 593                 if( kozoari[n] )
 594                 {
 595                         count++;
 596                         lastkozo = n;
 597                 }
 598         }
 599
 600 //      fprintf( stderr, "count = %d\n", count );
 601
 602         if( count == 1 )
 603                 strweight[lastkozo] = sumweight;
 604         else if( count > 1 )
 605         {
 606                 assignstrweight_rec( strweight, ob->children[dir_ch[0]], ob, kozoari, seqweight );
 607                 assignstrweight_rec( strweight, ob->children[dir_ch[1]], ob, kozoari, seqweight );
 608         }
 609 }
 610
 611 void assignstrweight( int nseq, double *strweight, Node *stopol, int ***topol, int step, int LorR, char *kozoari, double *seqweight )
 612 {
 613         Node *topNode, *btmNode;
 614         int i;
 615
 616         if( step == nseq - 2 )
 617         {
 618                 topNode = stopol[nseq-2].children[0];
 619                 btmNode = stopol + nseq-3;
 620 #if DEBUG
 621                 fprintf( stderr, "Now step == nseq-3, topNode = %d, btmNode = %d\n", topNode - stopol_g, btmNode-stopol_g );
 622 #endif
 623         }
 624
 625         else
 626         {
 627                 for( i=0; i<3; i++ )
 628                 {
 629                         if( stopol[step].members[i][0] == topol[step][LorR][0] )
 630                         break;
 631                 }
 632                 if( i== 3 ) ErrorExit( "Incorrect call of weightFromABranch." );
 633                 btmNode = stopol[step].children[i];
 634                 topNode = stopol+step;
 635         }
 636
 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 );
 643
 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" );
 657 #endif
 658 }