src/TBR.c

   1 /* TBR.c    2009-01-12 */
   2
   3 /* Copyright 2009 Richard Desper */
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 #include "me.h"
   9
  10 /*functions from me_balanced.c*/
  11 void makeBMEAveragesTable(tree *T, double **D, double **A);
  12 void assignBMEWeights(tree *T, double **A);
  13
  14 /*from me.c*/
  15 edge *siblingEdge(edge *e);
  16 double **initDoubleMatrix(int d);
  17 void freeMatrix(double **D, int size);
  18 edge *depthFirstTraverse(tree *T, edge *e);
  19 edge *findBottomLeft(edge *e);
  20
  21 /*from bnni.c*/
  22 void weighTree(tree *T);
  23
  24 void freeTree(tree *T);
  25
  26 /*from SPR.c*/
  27 void zero3DMatrix(double ***X, int h, int l, int w);
  28
  29 void assignTBRDownWeightsUp(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights,
  30                             double *bestWeight, edge **eBestSplit, edge **eBestTop, edge **eBestBottom);
  31 void assignTBRDownWeightsSkew(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights,
  32                               double *bestWeight, edge **eBestSplit, edge **eBestTop, edge **eBestBototm);
  33 void assignTBRDownWeightsDown(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights,
  34                               double *bestWeight, edge **eBestSplit, edge **eBestTop, edge **eBestBottom);
  35
  36 void assignTBRUpWeights(edge *ebottom, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A,
  37                                                 double **dXTop, double ***swapWeights, edge *etop, double *bestWeight,
  38                                                 edge **bestSplit, edge **bestTop, edge **bestBottom)
  39                                                 /*function assigns the value for etop if the tree below vtest is moved to be below etop*/
  40                                                 /*and SPR for tree bottom tree splits ebottom*/
  41                                                 /*recursive function searching over values of ebottom*/
  42                                                 /*minor variant of SPR.c's assignUpWeights
  43                                                 difference is the index assignment in the array swapWeights, which has a different meaning
  44                                                 for the TBR routines*/
  45 /*also using dXTop to assign value of average distance to tree above vtest*/
  46 { /*SPR performed on tree above vtest...*/
  47   edge *sib, *left, *right;
  48   /*B is above vtest, A is other tree below vtest unioned with trees in path to vtest*/
  49   /*sib is tree C being passed by B*/
  50   /*D is tree below etest*/
  51   double D_AB, D_CD, D_AC, D_BD;
  52   sib = siblingEdge(ebottom);
  53   left = ebottom->head->leftEdge;
  54   right = ebottom->head->rightEdge;
  55   if (NULL == etop)
  56     {
  57       if (NULL == back) /*first recursive call*/
  58         {
  59           if (NULL == left)
  60             return; /*no subtree below for SPR*/
  61           else       /*NULL == back and NULL == etop*/
  62             {
  63               assignTBRUpWeights(left,vtest,va,ebottom,va,A[va->index][vtest->index],0.5,A,dXTop,swapWeights,NULL,bestWeight,bestSplit,bestTop,bestBottom);
  64               assignTBRUpWeights(right,vtest,va,ebottom,va,A[va->index][vtest->index],0.5,A,dXTop,swapWeights,NULL,bestWeight,bestSplit,bestTop,bestBottom);
  65             }
  66         }
  67       else /*NULL != back*/
  68         {
  69           D_BD = A[ebottom->head->index][vtest->index];
  70           /*B is tree above vtest, D is tree below ebottom*/
  71           D_CD = A[sib->head->index][ebottom->head->index]; /*C is tree below sib*/
  72           D_AC = A[back->head->index][sib->head->index] +
  73             coeff*(A[va->index][sib->head->index] - A[vtest->index][sib->head->index]);
  74           /*va is root of subtree skew back at vtest*/
  75           /*A is union of tree below va and all subtrees already passed in path from vtest to ebottom*/
  76           D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
  77           swapWeights[vtest->index][ebottom->head->index][ebottom->head->index] = swapWeights[vtest->index][back->head->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
  78           if (swapWeights[vtest->index][ebottom->head->index][ebottom->head->index] < *bestWeight)
  79             {
  80               *bestSplit = vtest->parentEdge;
  81               *bestTop = NULL;
  82               *bestBottom = ebottom;
  83               *bestWeight = swapWeights[vtest->index][ebottom->head->index][ebottom->head->index];
  84             }
  85           if (NULL != left)
  86             {
  87               assignTBRUpWeights(left,vtest,va,ebottom,sib->head,D_AB,0.5*coeff,A,dXTop,swapWeights,NULL,bestWeight,bestSplit,bestTop,bestBottom);
  88               assignTBRUpWeights(right,vtest,va,ebottom,sib->head,D_AB,0.5*coeff,A,dXTop,swapWeights,NULL,bestWeight,bestSplit,bestTop,bestBottom);
  89             }
  90         }
  91     }
  92   else /*NULL != etop*/
  93     {
  94       if (NULL == back) /*first recursive call*/
  95         {
  96           if (swapWeights[vtest->index][etop->head->index][etop->head->index]< *bestWeight)
  97             /*this represents value of SPR from esplit to etop, with no SPR in bottom tree*/
  98             {
  99               *bestSplit = vtest->parentEdge;
 100               *bestTop = etop;
 101               *bestBottom = NULL;
 102               *bestWeight = swapWeights[vtest->index][etop->head->index][etop->head->index];
 103             }
 104           if (NULL == left)
 105             return; /*no subtree below for SPR*/
 106           else if (NULL != etop)/*start the process of assigning weights recursively*/
 107             {
 108               assignTBRUpWeights(left,vtest,va,ebottom,va,dXTop[va->index][etop->head->index],0.5,A,dXTop,swapWeights,etop,bestWeight,bestSplit,bestTop,bestBottom);
 109               assignTBRUpWeights(right,vtest,va,ebottom,va,dXTop[va->index][etop->head->index],0.5,A,dXTop,swapWeights,etop,bestWeight,bestSplit,bestTop,bestBottom);
 110             }
 111         } /*NULL == back*/
 112       /*in following bit, any average distance of form A[vtest->index][x->index] is
 113         replaced by dXTop[x->index][etop->head->index]*/
 114       else /*second or later recursive call, NULL != etop*/
 115         {
 116           D_BD = dXTop[ebottom->head->index][etop->head->index]; /*B is tree above vtest - it is in configuration
 117                                                                    indexed by etop*/
 118           /*D is tree below ebottom*/
 119           D_CD = A[sib->head->index][ebottom->head->index]; /*C is tree below sib*/
 120           D_AC = A[back->head->index][sib->head->index] +
 121             coeff*(A[va->index][sib->head->index] - A[sib->head->index][vtest->index]);
 122           /*it is correct to use A[][] here because the bad average distances involving B from the first term will
 123             be cancelled by the bad average distances involving B in the subtracted part*/
 124           /*va is root of subtree skew back at vtest*/
 125           /*A is union of tree below va and all subtrees already passed in path from vtest to ebottom*/
 126           D_AB = 0.5*(oldD_AB + dXTop[cprev->index][etop->head->index]);
 127           swapWeights[vtest->index][ebottom->head->index][etop->head->index] = swapWeights[vtest->index][back->head->index][etop->head->index]  + (D_AC + D_BD - D_AB - D_CD);
 128           if (swapWeights[vtest->index][ebottom->head->index][etop->head->index] + swapWeights[vtest->index][etop->head->index][etop->head->index]< *bestWeight)
 129             /*first term is contribution of second SPR, second term is contribution of first SPR*/
 130             {
 131               *bestSplit = vtest->parentEdge;
 132               *bestTop = etop;
 133               *bestBottom = ebottom;
 134               *bestWeight = swapWeights[vtest->index][ebottom->head->index][etop->head->index] + swapWeights[vtest->index][etop->head->index][etop->head->index];
 135             }
 136           if (NULL != left)
 137             {
 138               assignTBRUpWeights(left,vtest, va, ebottom, sib->head, D_AB, 0.5*coeff, A, dXTop, swapWeights,etop,bestWeight,bestSplit,bestTop,bestBottom);
 139               assignTBRUpWeights(right,vtest, va, ebottom, sib->head, D_AB, 0.5*coeff, A, dXTop, swapWeights,etop,bestWeight,bestSplit,bestTop,bestBottom);
 140             }
 141         } /*else NULL != back, etop*/
 142     }
 143 }
 144
 145 /*recall NNI formula: change in tree length from AB|CD split to AC|BD split is
 146 proportional to D_AC + D_BD - D_AB - D_CD*/
 147 /*in our case B is the tree being moved (below vtest), A is the tree backwards below back, but
 148 with the vtest subtree removed, C is the sibling tree of back and D is the tree above vtest*/
 149 /*use va to denote the root of the sibling tree to B in the original tree*/
 150 /*please excuse the multiple uses of the same letters: A,D, etc.*/
 151 void assignTBRDownWeightsUp(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights,
 152                             double *bestWeight, edge **bestSplitEdge, edge **bestTop, edge **bestBottom)
 153 {
 154         edge *par, *sib, *skew;
 155         double D_AC, D_BD, D_AB, D_CD;
 156         par = etest->tail->parentEdge;
 157         skew = siblingEdge(etest);
 158         if (NULL == back) /*first recursive call*/
 159           {
 160           if (NULL == par)
 161             return;
 162           else /*start the process of assigning weights recursively*/
 163             {
 164               assignTBRDownWeightsUp(par,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 165               assignTBRDownWeightsSkew(skew,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 166             }
 167         }
 168         else /*second or later recursive call*/
 169           {
 170             sib = siblingEdge(back);
 171             D_BD = A[vtest->index][etest->head->index]; /*straightforward*/
 172             D_CD = A[sib->head->index][etest->head->index]; /*this one too*/
 173             D_AC = A[sib->head->index][back->head->index] + coeff*(A[sib->head->index][va->index]
 174                                                                    - A[sib->head->index][vtest->index]);
 175             D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
 176             swapWeights[vtest->index][etest->head->index][etest->head->index] = swapWeights[vtest->index][back->head->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
 177             /*using diagonal to store values for SPR swaps above the split edge*/
 178             /*this is value of swapping tree below vtest to break etest*/
 179             if (swapWeights[vtest->index][etest->head->index][etest->head->index] < *bestWeight)
 180               {
 181                 *bestWeight = swapWeights[vtest->index][etest->head->index][etest->head->index];
 182                 *bestSplitEdge = vtest->parentEdge;
 183                 *bestTop = etest;
 184                 *bestBottom = NULL;
 185               }
 186             if (NULL != par)
 187               {
 188                 assignTBRDownWeightsUp(par,vtest,va,etest,sib->head,D_AB,0.5*coeff,A,swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 189                 assignTBRDownWeightsSkew(skew,vtest,va,etest,sib->head,D_AB,0.5*coeff,A,swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 190               }
 191           }
 192 }
 193
 194
 195 void assignTBRDownWeightsSkew(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights,
 196                               double *bestWeight, edge **bestSplitEdge, edge **bestTop, edge **bestBottom)
 197 {
 198   /*same general idea as assignDownWeights, except needing to keep track of things a bit differently*/
 199   edge *par, *left, *right;
 200   /*par here = sib before
 201     left, right here = par, skew before*/
 202   double D_AB, D_CD, D_AC, D_BD;
 203   /*B is subtree being moved - below vtest
 204     A is subtree remaining fixed - below va, unioned with all trees already passed by B*/
 205   /*C is subtree being passed by B, in this case above par
 206     D is subtree below etest, fixed on other side*/
 207   par = etest->tail->parentEdge;
 208   left = etest->head->leftEdge;
 209   right = etest->head->rightEdge;
 210   if (NULL == back)
 211     {
 212       if (NULL == left)
 213         return;
 214       else
 215         {
 216           assignTBRDownWeightsDown(left,vtest,va,etest,etest->tail,A[vtest->index][etest->tail->index],0.5,A,swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 217           assignTBRDownWeightsDown(right,vtest,va,etest,etest->tail,A[vtest->index][etest->tail->index],0.5,A,swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 218         }
 219     }
 220   else
 221     {
 222       D_BD = A[vtest->index][etest->head->index];
 223       D_CD = A[par->head->index][etest->head->index];
 224       D_AC = A[back->head->index][par->head->index] + coeff*(A[va->index][par->head->index] - A[vtest->index][par->head->index]);
 225       D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
 226       swapWeights[vtest->index][etest->head->index][etest->head->index] = swapWeights[vtest->index][back->head->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
 227       if (swapWeights[vtest->index][etest->head->index][etest->head->index] < *bestWeight)
 228         {
 229           *bestWeight = swapWeights[vtest->index][etest->head->index][etest->head->index];
 230           *bestSplitEdge = vtest->parentEdge;
 231           *bestTop = etest;
 232           *bestBottom = NULL;
 233         }
 234       if (NULL != left)
 235         {
 236           assignTBRDownWeightsDown(left,vtest, va, etest, etest->tail, D_AB, 0.5*coeff, A, swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 237           assignTBRDownWeightsDown(right,vtest, va, etest, etest->tail, D_AB, 0.5*coeff, A, swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 238         }
 239     }
 240 }
 241
 242 void assignTBRDownWeightsDown(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights,
 243                               double *bestWeight, edge **bestSplitEdge, edge **bestTop, edge **bestBottom)
 244 {
 245   /*again the same general idea*/
 246   edge *sib, *left, *right;
 247   /*sib here = par in assignDownWeightsSkew
 248     rest is the same as assignDownWeightsSkew*/
 249   double D_AB, D_CD, D_AC, D_BD;
 250   /*B is below vtest, A is below va unioned with all trees already passed by B*/
 251   /*C is subtree being passed - below sib*/
 252   /*D is tree below etest*/
 253   sib = siblingEdge(etest);
 254   left = etest->head->leftEdge;
 255   right = etest->head->rightEdge;
 256   D_BD = A[vtest->index][etest->head->index];
 257   D_CD = A[sib->head->index][etest->head->index];
 258   D_AC = A[sib->head->index][back->head->index] + coeff*(A[sib->head->index][va->index] - A[sib->head->index][vtest->index]);
 259   D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
 260   swapWeights[vtest->index][etest->head->index][etest->head->index] = swapWeights[vtest->index][back->head->index][back->head->index] + ( D_AC + D_BD - D_AB - D_CD);
 261   if (swapWeights[vtest->index][etest->head->index][etest->head->index] < *bestWeight)
 262     {
 263       *bestWeight = swapWeights[vtest->index][etest->head->index][etest->head->index];
 264       *bestSplitEdge = vtest->parentEdge;
 265       *bestTop = etest;
 266       *bestBottom = NULL;
 267     }
 268   if (NULL != left)
 269     {
 270       assignTBRDownWeightsDown(left,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 271       assignTBRDownWeightsDown(right,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights,bestWeight,bestSplitEdge,bestTop,bestBottom);
 272     }
 273 }
 274
 275 /*general idea is to have a double loop for a given edge, testing all SPRs for the subtrees above and below a given edge.
 276   Then that function loops over all the edges of a tree*/
 277
 278 void TBRswitch(tree *T, edge *e1, edge *e2, edge *e3);
 279
 280 /*vbottom is node below esplit for average calculations in matrix dXTop, A is matrix of average
 281   distances from original tree, dXout is average distance from vbottom to tree rooted at far edge
 282   of eback, if SPR breaking eback, UpOrDown indicates whether etop is in path above split edge
 283   (Up) or not (Down)*/
 284 void calcTBRTopBottomAverage(node *vbottom, double **A, double **dXTop, double dXOut,
 285                              edge *esplit, edge *etop, edge *eback, int UpOrDown)
 286 {
 287   edge *enew1, *enew2, *emove;
 288   double newdXOut;
 289   if (esplit == eback) /*top level call - means trivial SPR*/
 290     dXTop[vbottom->index][etop->head->index] = A[vbottom->index][esplit->head->index];
 291   else
 292     dXTop[vbottom->index][etop->head->index] = dXTop[vbottom->index][eback->head->index] +
 293       0.25*(A[vbottom->index][etop->head->index] - dXOut);
 294   /*by moving etop past the vbottom tree, everything in the etop tree is closer by coefficient of
 295     0.25, while everything in the old back tree is further by a coefficient of 0.25*/
 296   /*everything in the tree that is being moved (emove) keeps the same relative weight in the average
 297     distance calculation*/
 298   if (UP == UpOrDown)
 299     {
 300       enew1 = etop->tail->parentEdge;
 301       if (NULL != enew1) /*need to do recursive calls*/
 302         {
 303           enew2 = siblingEdge(etop);
 304           emove = siblingEdge(eback);    /*emove is third edge meeting at vertex with eback, etest*/
 305           if (esplit == eback)
 306             newdXOut = A[vbottom->index][emove->head->index];
 307           else
 308             newdXOut = 0.5*(dXOut + A[vbottom->index][emove->head->index]);
 309           calcTBRTopBottomAverage(vbottom,A,dXTop,newdXOut,esplit, enew1,etop,UP); /*old etop is new value for back*/
 310           calcTBRTopBottomAverage(vbottom,A,dXTop,newdXOut,esplit, enew2,etop,DOWN);
 311         }
 312     }
 313   else /*moving down*/
 314     {
 315       enew1 = etop->head->leftEdge;
 316       if (NULL != enew1)
 317         {
 318           enew2 = etop->head->rightEdge;
 319           if (eback == siblingEdge(etop))
 320             emove = etop->tail->parentEdge;
 321           else
 322             emove = siblingEdge(etop);
 323           if (esplit == eback)
 324             newdXOut = A[vbottom->index][emove->head->index];
 325           else
 326             newdXOut = 0.5*(dXOut + A[vbottom->index][emove->head->index]);
 327           calcTBRTopBottomAverage(vbottom,A,dXTop,newdXOut,esplit,enew1,etop,DOWN);
 328           calcTBRTopBottomAverage(vbottom,A,dXTop,newdXOut,esplit,enew2,etop,DOWN);
 329         }
 330     }
 331 }
 332
 333 void calcTBRaverages(tree *T, edge *esplit, double **A, double **dXTop)
 334 {
 335   edge *ebottom, *par, *sib;
 336   for (ebottom = findBottomLeft(esplit); ebottom != esplit; ebottom = depthFirstTraverse(T,ebottom))
 337     {
 338       par = esplit->tail->parentEdge;
 339       sib = siblingEdge(esplit);
 340       calcTBRTopBottomAverage(ebottom->head,A, dXTop, 0.0, esplit, par,esplit,UP);
 341       calcTBRTopBottomAverage(ebottom->head,A, dXTop, 0.0, esplit, sib,esplit,DOWN);
 342     }
 343 }
 344
 345 void TBR(tree *T, double **D, double **A)
 346 {
 347   int i;
 348   edge *esplit, *etop, *eBestTop, *eBestBottom, *eBestSplit;
 349   edge *eout, *block;
 350   edge *left, *right, *par, *sib;
 351   double **dXTop; /*dXTop[i][j] is average distance from subtree rooted at i to tree above split edge, if
 352                     SPR above split edge cuts edge whose head has index j*/
 353   double bestWeight;
 354   double ***TBRWeights;
 355   dXTop = initDoubleMatrix(T->size);
 356   weighTree(T);
 357   TBRWeights = (double ***)calloc(T->size,sizeof(double **));
 358   for(i=0;i<T->size;i++)
 359     TBRWeights[i] = initDoubleMatrix(T->size);
 360   do
 361     {
 362       zero3DMatrix(TBRWeights,T->size,T->size,T->size);
 363       bestWeight = 0.0;
 364       eBestSplit = eBestTop = eBestBottom = NULL;
 365       for(esplit=depthFirstTraverse(T,NULL);NULL!=esplit;esplit=depthFirstTraverse(T,esplit))
 366         {
 367           par = esplit->tail->parentEdge;
 368           if (NULL != par)
 369             {
 370               sib = siblingEdge(esplit);
 371               /*next two lines calculate the possible improvements for any SPR above esplit*/
 372               assignTBRDownWeightsUp(par,esplit->head,sib->head,NULL,NULL,0.0,1.0,A,TBRWeights,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 373               assignTBRDownWeightsSkew(sib,esplit->head,sib->tail,NULL,NULL,0.0,1.0,A,TBRWeights,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 374               calcTBRaverages(T,esplit,A,dXTop); /*calculates the average distance from any subtree
 375                                                    below esplit to the entire subtree above esplit,
 376                                                    after any possible SPR above*/
 377               /*for etop above esplit, we loop using information from above to calculate values
 378                 for all possible SPRs below esplit*/
 379             }
 380
 381           right = esplit->head->rightEdge;
 382           if (NULL != right)
 383             {
 384               left = esplit->head->leftEdge;
 385               /*test case: etop = null means only do bottom SPR*/
 386               assignTBRUpWeights(left,esplit->head,right->head,NULL,NULL,0.0,1.0,A,dXTop,TBRWeights,NULL,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 387               assignTBRUpWeights(right,esplit->head,left->head,NULL,NULL,0.0,1.0,A,dXTop,TBRWeights,NULL,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 388
 389               block = esplit;
 390               while (NULL != block)
 391                 {
 392                   if (block != esplit)
 393                     {
 394                       etop = block;
 395                       assignTBRUpWeights(left,esplit->head,right->head,NULL,NULL,0.0,1.0,A,dXTop,TBRWeights,etop,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 396                       assignTBRUpWeights(right,esplit->head,left->head,NULL,NULL,0.0,1.0,A,dXTop,TBRWeights,etop,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 397                     }
 398                   eout = siblingEdge(block);
 399                   if (NULL != eout)
 400                     {
 401                       etop = findBottomLeft(eout);
 402                       while (etop->tail != eout->tail)
 403                         {
 404                           /*for ebottom below esplit*/
 405
 406                           assignTBRUpWeights(left,esplit->head,right->head,NULL,NULL,0.0,1.0,A,dXTop,TBRWeights,etop,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 407                           assignTBRUpWeights(right,esplit->head,left->head,NULL,NULL,0.0,1.0,A,dXTop,TBRWeights,etop,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 408                           etop = depthFirstTraverse(T,etop);
 409                         }
 410
 411                       /*etop == eout*/
 412
 413                       assignTBRUpWeights(left,esplit->head,right->head,NULL,NULL,0.0,1.0,A,dXTop,TBRWeights,etop,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 414                       assignTBRUpWeights(right,esplit->head,left->head,NULL,NULL,0.0,1.0,A,dXTop,TBRWeights,etop,&bestWeight,&eBestSplit,&eBestTop,&eBestBottom);
 415                     }
 416                   block = block->tail->parentEdge;
 417                 }
 418             } /*if NULL != right*/
 419         } /*for esplit*/
 420       /*find bestWeight, best split edge, etc.*/
 421       if (bestWeight < -EPSILON)
 422         {
 423 //        if (verbose)
 424 //          {
 425 //            printf("TBR #%d: Splitting edge %s: top edge %s, bottom edge %s\n",*count+1,
 426 //                   eBestSplit->label, eBestTop->label,eBestBottom->label);
 427 //            printf("Old tree weight is %lf, new tree weight should be %lf\n",T->weight, T->weight + 0.25*bestWeight);
 428 //          }
 429           TBRswitch(T,eBestSplit,eBestTop,eBestBottom);
 430           makeBMEAveragesTable(T,D,A);
 431           assignBMEWeights(T,A);
 432           weighTree(T);
 433 //        if (verbose)
 434 //          printf("TBR: new tree weight is %lf\n\n",T->weight);<
 435 //        (*count)++;
 436         }
 437       else
 438         bestWeight = 1.0;
 439     } while (bestWeight < -EPSILON);
 440   for(i=0;i<T->size;i++)
 441     freeMatrix(TBRWeights[i],T->size);
 442   freeMatrix(dXTop,T->size);
 443 }
 444
 445 void SPRTopShift(tree *T, node *v, edge *e, int UpOrDown);
 446
 447 void TBRswitch(tree *T, edge *es, edge *et, edge *eb)
 448 {
 449         if (NULL != et)
 450                 SPRTopShift(T,es->head,et,DOWN); /*DOWN because tree being moved is below split edge*/
 451         if (NULL != eb)
 452                 SPRTopShift(T,es->head,eb,UP);   /*UP because tree being moved is above split edge*/
 453 }