src/TBR.c

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