3 /* Copyright 2009 Richard Desper */
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
10 /*functions from me_balanced.c*/
11 void makeBMEAveragesTable(tree *T, double **D, double **A);
12 void assignBMEWeights(tree *T, double **A);
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);
22 void weighTree(tree *T);
24 void freeTree(tree *T);
27 void zero3DMatrix(double ***X, int h, int l, int w);
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);
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;
57 if (NULL == back) /*first recursive call*/
60 return; /*no subtree below for SPR*/
61 else /*NULL == back and NULL == etop*/
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);
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)
80 *bestSplit = vtest->parentEdge;
82 *bestBottom = ebottom;
83 *bestWeight = swapWeights[vtest->index][ebottom->head->index][ebottom->head->index];
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);
94 if (NULL == back) /*first recursive call*/
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*/
99 *bestSplit = vtest->parentEdge;
102 *bestWeight = swapWeights[vtest->index][etop->head->index][etop->head->index];
105 return; /*no subtree below for SPR*/
106 else if (NULL != etop)/*start the process of assigning weights recursively*/
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);
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*/
116 D_BD = dXTop[ebottom->head->index][etop->head->index]; /*B is tree above vtest - it is in configuration
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*/
131 *bestSplit = vtest->parentEdge;
133 *bestBottom = ebottom;
134 *bestWeight = swapWeights[vtest->index][ebottom->head->index][etop->head->index] + swapWeights[vtest->index][etop->head->index][etop->head->index];
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);
141 } /*else NULL != back, etop*/
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)
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*/
162 else /*start the process of assigning weights recursively*/
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);
168 else /*second or later recursive call*/
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)
181 *bestWeight = swapWeights[vtest->index][etest->head->index][etest->head->index];
182 *bestSplitEdge = vtest->parentEdge;
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);
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)
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;
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);
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)
229 *bestWeight = swapWeights[vtest->index][etest->head->index][etest->head->index];
230 *bestSplitEdge = vtest->parentEdge;
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);
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)
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)
263 *bestWeight = swapWeights[vtest->index][etest->head->index][etest->head->index];
264 *bestSplitEdge = vtest->parentEdge;
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);
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*/
278 void TBRswitch(tree *T, edge *e1, edge *e2, edge *e3);
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
284 void calcTBRTopBottomAverage(node *vbottom, double **A, double **dXTop, double dXOut,
285 edge *esplit, edge *etop, edge *eback, int UpOrDown)
287 edge *enew1, *enew2, *emove;
289 if (esplit == eback) /*top level call - means trivial SPR*/
290 dXTop[vbottom->index][etop->head->index] = A[vbottom->index][esplit->head->index];
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*/
300 enew1 = etop->tail->parentEdge;
301 if (NULL != enew1) /*need to do recursive calls*/
303 enew2 = siblingEdge(etop);
304 emove = siblingEdge(eback); /*emove is third edge meeting at vertex with eback, etest*/
306 newdXOut = A[vbottom->index][emove->head->index];
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);
315 enew1 = etop->head->leftEdge;
318 enew2 = etop->head->rightEdge;
319 if (eback == siblingEdge(etop))
320 emove = etop->tail->parentEdge;
322 emove = siblingEdge(etop);
324 newdXOut = A[vbottom->index][emove->head->index];
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);
333 void calcTBRaverages(tree *T, edge *esplit, double **A, double **dXTop)
335 edge *ebottom, *par, *sib;
336 for (ebottom = findBottomLeft(esplit); ebottom != esplit; ebottom = depthFirstTraverse(T,ebottom))
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);
345 void TBR(tree *T, double **D, double **A)
348 edge *esplit, *etop, *eBestTop, *eBestBottom, *eBestSplit;
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*/
354 double ***TBRWeights;
355 dXTop = initDoubleMatrix(T->size);
357 TBRWeights = (double ***)calloc(T->size,sizeof(double **));
358 for(i=0;i<T->size;i++)
359 TBRWeights[i] = initDoubleMatrix(T->size);
362 zero3DMatrix(TBRWeights,T->size,T->size,T->size);
364 eBestSplit = eBestTop = eBestBottom = NULL;
365 for(esplit=depthFirstTraverse(T,NULL);NULL!=esplit;esplit=depthFirstTraverse(T,esplit))
367 par = esplit->tail->parentEdge;
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*/
381 right = esplit->head->rightEdge;
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);
390 while (NULL != 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);
398 eout = siblingEdge(block);
401 etop = findBottomLeft(eout);
402 while (etop->tail != eout->tail)
404 /*for ebottom below esplit*/
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);
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);
416 block = block->tail->parentEdge;
418 } /*if NULL != right*/
420 /*find bestWeight, best split edge, etc.*/
421 if (bestWeight < -EPSILON)
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);
429 TBRswitch(T,eBestSplit,eBestTop,eBestBottom);
430 makeBMEAveragesTable(T,D,A);
431 assignBMEWeights(T,A);
434 // printf("TBR: new tree weight is %lf\n\n",T->weight);<
439 } while (bestWeight < -EPSILON);
440 for(i=0;i<T->size;i++)
441 freeMatrix(TBRWeights[i],T->size);
442 freeMatrix(dXTop,T->size);
445 void SPRTopShift(tree *T, node *v, edge *e, int UpOrDown);
447 void TBRswitch(tree *T, edge *es, edge *et, edge *eb)
450 SPRTopShift(T,es->head,et,DOWN); /*DOWN because tree being moved is below split edge*/
452 SPRTopShift(T,es->head,eb,UP); /*UP because tree being moved is above split edge*/