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 bNNI.c*/
11 void makeBMEAveragesTable(tree *T, double **D, double **A);
12 void assignBMEWeights(tree *T, double **A);
15 edge *siblingEdge(edge *e);
16 void weighTree(tree *T);
17 void freeMatrix(double **D, int size);
18 edge *depthFirstTraverse(tree *T, edge *e);
19 double **initDoubleMatrix(int d);
22 node *indexedNode(tree *T, int i);
23 edge *indexedEdge(tree *T, int i);
24 void assignSPRWeights(node *v, double **A, double ***swapWeights);
25 void SPRTopShift(tree *T, node *vmove, edge *esplit, int UpOrDown);
26 void assignDownWeightsUp(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights);
27 void assignDownWeightsSkew(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights);
28 void assignDownWeightsDown(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights);
29 void assignUpWeights(edge *etest, node *vtest, node *va, edge *back, node *cprve, double oldD_AB, double coeff, double **A, double ***swapWeights);
31 void zero3DMatrix(double ***X, int h, int l, int w)
41 void findTableMin(int *imin, int *jmin, int *kmin, int n, double ***X, double *min)
48 if (X[i][j][k] < *min)
59 void SPR(tree *T, double **D, double **A, int *count)
65 /* char filename[MAX_LABEL_LENGTH];*/
66 double ***swapWeights;
67 double swapValue = 0.0;
68 swapWeights = (double ***)malloc(2*sizeof(double **));
69 makeBMEAveragesTable(T,D,A);
70 assignBMEWeights(T,A);
73 printf("Before SPRs: tree length is %lf.\n",T->weight);*/
75 swapWeights[i] = initDoubleMatrix(T->size);
79 zero3DMatrix(swapWeights,2,T->size,T->size);
81 for(e=depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
82 assignSPRWeights(e->head,A,swapWeights);
83 findTableMin(&i,&j,&k,T->size,swapWeights,&swapValue);
84 swapValue = swapWeights[i][j][k];
85 if (swapValue < -EPSILON)
88 // printf("New tree weight should be %lf.\n",T->weight + 0.25*swapValue);
92 // printf("Swapping tree below %s to split edge %s with head %s and tail %s\n",
93 // v->parentEdge->label,f->label,f->head->label,f->tail->label);
94 SPRTopShift(T,v,f,2-i);
95 makeBMEAveragesTable(T,D,A);
96 assignBMEWeights(T,A);
99 /*sprintf(filename,"tree%d.new",*count);*/
101 // printf("After %d SPRs, tree weight is %lf.\n\n",*count,T->weight);
102 /*treefile = fopen(filename,"w");
103 NewickPrintTree(T,treefile);
106 } while (swapValue < -EPSILON);
108 freeMatrix(swapWeights[i],T->size);
114 /*assigns values to array swapWeights*/
115 /*swapWeights[0][j][k] will be the value of removing the tree below the edge whose head node has index j
116 and reattaching it to split the edge whose head has the index k*/
117 /*swapWeights[1][j][k] will be the value of removing the tree above the edge whose head node has index j
118 and reattaching it to split the edge whose head has the index k*/
119 void assignSPRWeights(node *vtest, double **A, double ***swapWeights)
121 edge *etest, *left, *right, *sib, *par;
122 etest = vtest->parentEdge;
123 left = vtest->leftEdge;
124 right = vtest->rightEdge;
125 par = etest->tail->parentEdge;
126 sib = siblingEdge(etest);
128 assignDownWeightsUp(par,vtest,sib->head,NULL,NULL,0.0,1.0,A,swapWeights);
130 assignDownWeightsSkew(sib,vtest,sib->tail,NULL,NULL,0.0,1.0,A,swapWeights);
131 /*assigns values for moving subtree rooted at vtest, starting with edge
132 parental to tail of edge parental to vtest*/
135 assignUpWeights(left,vtest,right->head,NULL,NULL,0.0,1.0,A,swapWeights);
136 assignUpWeights(right,vtest,left->head,NULL,NULL,0.0,1.0,A,swapWeights);
141 /*recall NNI formula: change in tree length from AB|CD split to AC|BD split is
142 proportional to D_AC + D_BD - D_AB - D_CD*/
143 /*in our case B is the tree being moved (below vtest), A is the tree backwards below back, but
144 with the vtest subtree removed, C is the sibling tree of back and D is the tree above etest*/
145 /*use va to denote the root of the sibling tree to B in the original tree*/
146 /*please excuse the multiple uses of the same letters: A,D, etc.*/
147 void assignDownWeightsUp(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights)
149 edge *par, *sib, *skew;
150 double D_AC, D_BD, D_AB, D_CD;
151 par = etest->tail->parentEdge;
152 skew = siblingEdge(etest);
153 if (NULL == back) /*first recursive call*/
157 else /*start the process of assigning weights recursively*/
159 assignDownWeightsUp(par,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
160 assignDownWeightsSkew(skew,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
163 else /*second or later recursive call*/
165 sib = siblingEdge(back);
166 D_BD = A[vtest->index][etest->head->index]; /*straightforward*/
167 D_CD = A[sib->head->index][etest->head->index]; /*this one too*/
168 D_AC = A[sib->head->index][back->head->index] + coeff*(A[sib->head->index][va->index]
169 - A[sib->head->index][vtest->index]);
170 D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
171 swapWeights[0][vtest->index][etest->head->index] = swapWeights[0][vtest->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
174 assignDownWeightsUp(par,vtest,va,etest,sib->head,D_AB,0.5*coeff,A,swapWeights);
175 assignDownWeightsSkew(skew,vtest,va,etest,sib->head,D_AB,0.5*coeff,A,swapWeights);
180 void assignDownWeightsSkew(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights)
182 /*same general idea as assignDownWeights, except needing to keep track of things a bit differently*/
183 edge *par, *left, *right;
184 /*par here = sib before
185 left, right here = par, skew before*/
186 double D_AB, D_CD, D_AC, D_BD;
187 /*B is subtree being moved - below vtest
188 A is subtree remaining fixed - below va, unioned with all trees already passed by B*/
189 /*C is subtree being passed by B, in this case above par
190 D is subtree below etest, fixed on other side*/
191 par = etest->tail->parentEdge;
192 left = etest->head->leftEdge;
193 right = etest->head->rightEdge;
200 assignDownWeightsDown(left,vtest,va,etest,etest->tail,A[vtest->index][etest->tail->index],0.5,A,swapWeights);
201 assignDownWeightsDown(right,vtest,va,etest,etest->tail,A[vtest->index][etest->tail->index],0.5,A,swapWeights);
206 D_BD = A[vtest->index][etest->head->index];
207 D_CD = A[par->head->index][etest->head->index];
208 D_AC = A[back->head->index][par->head->index] + coeff*(A[va->index][par->head->index] - A[vtest->index][par->head->index]);
209 D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
210 swapWeights[0][vtest->index][etest->head->index] = swapWeights[0][vtest->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
213 assignDownWeightsDown(left,vtest, va, etest, etest->tail, D_AB, 0.5*coeff, A, swapWeights);
214 assignDownWeightsDown(right,vtest, va, etest, etest->tail, D_AB, 0.5*coeff, A, swapWeights);
219 void assignDownWeightsDown(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights)
221 /*again the same general idea*/
222 edge *sib, *left, *right;
223 /*sib here = par in assignDownWeightsSkew
224 rest is the same as assignDownWeightsSkew*/
225 double D_AB, D_CD, D_AC, D_BD;
226 /*B is below vtest, A is below va unioned with all trees already passed by B*/
227 /*C is subtree being passed - below sib*/
228 /*D is tree below etest*/
229 sib = siblingEdge(etest);
230 left = etest->head->leftEdge;
231 right = etest->head->rightEdge;
232 D_BD = A[vtest->index][etest->head->index];
233 D_CD = A[sib->head->index][etest->head->index];
234 D_AC = A[sib->head->index][back->head->index] + coeff*(A[sib->head->index][va->index] - A[sib->head->index][vtest->index]);
235 D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
236 swapWeights[0][vtest->index][etest->head->index] = swapWeights[0][vtest->index][back->head->index] + ( D_AC + D_BD - D_AB - D_CD);
239 assignDownWeightsDown(left,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
240 assignDownWeightsDown(right,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
244 void assignUpWeights(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A,
245 double ***swapWeights)
247 /*SPR performed on tree above vtest...*/
248 /*same idea as above, with appropriate selections of edges and nodes*/
249 edge *sib, *left, *right;
250 /*B is above vtest, A is other tree below vtest unioned with trees in path to vtest*/
251 /*sib is tree C being passed by B*/
252 /*D is tree below etest*/
253 double D_AB, D_CD, D_AC, D_BD;
255 sib = siblingEdge(etest);
256 left = etest->head->leftEdge;
257 right = etest->head->rightEdge;
258 if (NULL == back) /*first recursive call*/
262 else /*start the process of assigning weights recursively*/
264 assignUpWeights(left,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
265 assignUpWeights(right,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
268 else /*second or later recursive call*/
270 D_BD = A[vtest->index][etest->head->index];
271 D_CD = A[sib->head->index][etest->head->index];
272 D_AC = A[back->head->index][sib->head->index] + coeff*(A[va->index][sib->head->index] - A[vtest->index][sib->head->index]);
273 D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
274 thisWeight = swapWeights[1][vtest->index][etest->head->index] = swapWeights[1][vtest->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
277 assignUpWeights(left,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
278 assignUpWeights(right,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
283 void pruneSubtree(edge *p, edge *u, edge *d)
284 /*starting with edge u above edges p, d*/
285 /*removes p, d from tree, u connects to d->head to compensate*/
287 p->tail->parentEdge = NULL;/*remove p subtree*/
289 d->head->parentEdge = u; /*u connected to d->head*/
290 d->head = NULL; /*d removed from tree*/
293 void SPRsplitEdge(edge *e, edge *p, edge *d)
294 /*splits edge e to make it parental to p,d. d is parental to what
295 previously was below e*/
299 p->tail->parentEdge = e;
300 d->head->parentEdge = d;
304 /*topological shift function*/
305 /*removes subtree rooted at v and re-inserts to spilt e*/
306 void SPRDownShift(tree *T, node *v, edge *e)
308 edge *vup, *vdown, *vpar;
309 vpar = v->parentEdge;
310 vdown = siblingEdge(vpar);
311 vup = vpar->tail->parentEdge;
312 /*topological shift*/
313 pruneSubtree(vpar,vup,vdown);
314 /*removes v subtree and vdown, extends vup*/
315 SPRsplitEdge(e,vpar,vdown); /*splits e to make e sibling edge to vpar,
319 void SPRUpShift(tree *T, node *vmove, edge *esplit)
320 /*an inelegant iterative version*/
328 for(f=esplit->tail->parentEdge,pathLength=1;f->tail != vmove;f=f->tail->parentEdge)
330 /*count number of edges to vmove*/
331 /*note that pathLength > 0 will hold*/
333 EPath = (edge **)malloc(pathLength*sizeof(edge *));
334 v = (node **)malloc(pathLength*sizeof(edge *));
335 sib = (edge **)malloc((pathLength+1)*sizeof(edge *));
336 /*there are pathLength + 1 side trees, one at the head and tail of each edge in the path*/
338 sib[pathLength] = siblingEdge(esplit);
340 f = esplit->tail->parentEdge;
345 sib[i] = siblingEdge(f);
347 f = f->tail->parentEdge;
349 /*indexed so head of Epath is v[i], tail is v[i-1] and sibling edge is sib[i]*/
350 /*need to assign head, tail of each edge in path
351 as well as have proper values for the left and right fields*/
353 if (esplit == esplit->tail->leftEdge)
355 vmove->leftEdge = esplit;
356 vmove->rightEdge = EPath[pathLength-1];
360 vmove->rightEdge = esplit;
361 vmove->leftEdge = EPath[pathLength-1];
363 esplit->tail = vmove;
364 /*espilt->head remains unchanged*/
365 /*vmove has the proper fields for left, right, and parentEdge*/
367 for(i=0;i<(pathLength-1);i++)
368 EPath[i]->tail = v[i+1];
370 /*this bit flips the orientation along the path
371 the tail of Epath[i] is now v[i+1] instead of v[i-1]*/
373 EPath[pathLength-1]->tail = vmove;
375 for(i=1;i<pathLength;i++)
377 if (sib[i+1] == v[i]->leftEdge)
378 v[i]->rightEdge = EPath[i-1];
380 v[i]->leftEdge = EPath[i-1];
382 if (sib[1] == v[0]->leftEdge)
383 v[0]->rightEdge = sib[0];
385 v[0]->leftEdge = sib[0];
393 void SPRTopShift(tree *T, node *vmove, edge *esplit, int UpOrDown)
395 if (DOWN == UpOrDown)
396 SPRDownShift(T,vmove,esplit);
398 SPRUpShift(T,vmove,esplit);
401 node *indexedNode(tree *T, int i)
404 for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
405 if (i == e->head->index)
407 if (i == T->root->index)
412 edge *indexedEdge(tree *T, int i)
415 for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
416 if (i == e->head->index)