1 /* #include <stdio.h> */
2 /* #include <stdlib.h> */
3 /* #include <math.h> */
4 /* #include "graph.h" */
5 /* #include "fastme.h" */
8 /*functions from bNNI.c*/
9 void makeBMEAveragesTable(tree *T, double **D, double **A);
10 void assignBMEWeights(tree *T, double **A);
13 edge *siblingEdge(edge *e);
14 void weighTree(tree *T);
15 void freeMatrix(double **D, int size);
16 edge *depthFirstTraverse(tree *T, edge *e);
17 double **initDoubleMatrix(int d);
20 node *indexedNode(tree *T, int i);
21 edge *indexedEdge(tree *T, int i);
22 void assignSPRWeights(node *v, double **A, double ***swapWeights);
23 void SPRTopShift(tree *T, node *vmove, edge *esplit, int UpOrDown);
24 void assignDownWeightsUp(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights);
25 void assignDownWeightsSkew(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights);
26 void assignDownWeightsDown(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights);
27 void assignUpWeights(edge *etest, node *vtest, node *va, edge *back, node *cprve, double oldD_AB, double coeff, double **A, double ***swapWeights);
29 void zero3DMatrix(double ***X, int h, int l, int w)
39 void findTableMin(int *imin, int *jmin, int *kmin, int n, double ***X, double *min)
46 if (X[i][j][k] < *min)
57 void SPR(tree *T, double **D, double **A, int *count)
63 /* char filename[MAX_LABEL_LENGTH];*/
64 double ***swapWeights;
65 double swapValue = 0.0;
66 swapWeights = (double ***)malloc(2*sizeof(double **));
67 makeBMEAveragesTable(T,D,A);
68 assignBMEWeights(T,A);
71 printf("Before SPRs: tree length is %lf.\n",T->weight);*/
73 swapWeights[i] = initDoubleMatrix(T->size);
77 zero3DMatrix(swapWeights,2,T->size,T->size);
79 for(e=depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
80 assignSPRWeights(e->head,A,swapWeights);
81 findTableMin(&i,&j,&k,T->size,swapWeights,&swapValue);
82 swapValue = swapWeights[i][j][k];
83 if (swapValue < -EPSILON)
86 // printf("New tree weight should be %lf.\n",T->weight + 0.25*swapValue);
90 // printf("Swapping tree below %s to split edge %s with head %s and tail %s\n",
91 // v->parentEdge->label,f->label,f->head->label,f->tail->label);
92 SPRTopShift(T,v,f,2-i);
93 makeBMEAveragesTable(T,D,A);
94 assignBMEWeights(T,A);
97 /*sprintf(filename,"tree%d.new",*count);*/
99 // printf("After %d SPRs, tree weight is %lf.\n\n",*count,T->weight);
100 /*treefile = fopen(filename,"w");
101 NewickPrintTree(T,treefile);
104 } while (swapValue < -EPSILON);
106 freeMatrix(swapWeights[i],T->size);
112 /*assigns values to array swapWeights*/
113 /*swapWeights[0][j][k] will be the value of removing the tree below the edge whose head node has index j
114 and reattaching it to split the edge whose head has the index k*/
115 /*swapWeights[1][j][k] will be the value of removing the tree above the edge whose head node has index j
116 and reattaching it to split the edge whose head has the index k*/
117 void assignSPRWeights(node *vtest, double **A, double ***swapWeights)
119 edge *etest, *left, *right, *sib, *par;
120 etest = vtest->parentEdge;
121 left = vtest->leftEdge;
122 right = vtest->rightEdge;
123 par = etest->tail->parentEdge;
124 sib = siblingEdge(etest);
126 assignDownWeightsUp(par,vtest,sib->head,NULL,NULL,0.0,1.0,A,swapWeights);
128 assignDownWeightsSkew(sib,vtest,sib->tail,NULL,NULL,0.0,1.0,A,swapWeights);
129 /*assigns values for moving subtree rooted at vtest, starting with edge
130 parental to tail of edge parental to vtest*/
133 assignUpWeights(left,vtest,right->head,NULL,NULL,0.0,1.0,A,swapWeights);
134 assignUpWeights(right,vtest,left->head,NULL,NULL,0.0,1.0,A,swapWeights);
139 /*recall NNI formula: change in tree length from AB|CD split to AC|BD split is
140 proportional to D_AC + D_BD - D_AB - D_CD*/
141 /*in our case B is the tree being moved (below vtest), A is the tree backwards below back, but
142 with the vtest subtree removed, C is the sibling tree of back and D is the tree above etest*/
143 /*use va to denote the root of the sibling tree to B in the original tree*/
144 /*please excuse the multiple uses of the same letters: A,D, etc.*/
145 void assignDownWeightsUp(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights)
147 edge *par, *sib, *skew;
148 double D_AC, D_BD, D_AB, D_CD;
149 par = etest->tail->parentEdge;
150 skew = siblingEdge(etest);
151 if (NULL == back) /*first recursive call*/
155 else /*start the process of assigning weights recursively*/
157 assignDownWeightsUp(par,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
158 assignDownWeightsSkew(skew,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
161 else /*second or later recursive call*/
163 sib = siblingEdge(back);
164 D_BD = A[vtest->index][etest->head->index]; /*straightforward*/
165 D_CD = A[sib->head->index][etest->head->index]; /*this one too*/
166 D_AC = A[sib->head->index][back->head->index] + coeff*(A[sib->head->index][va->index]
167 - A[sib->head->index][vtest->index]);
168 D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
169 swapWeights[0][vtest->index][etest->head->index] = swapWeights[0][vtest->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
172 assignDownWeightsUp(par,vtest,va,etest,sib->head,D_AB,0.5*coeff,A,swapWeights);
173 assignDownWeightsSkew(skew,vtest,va,etest,sib->head,D_AB,0.5*coeff,A,swapWeights);
178 void assignDownWeightsSkew(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights)
180 /*same general idea as assignDownWeights, except needing to keep track of things a bit differently*/
181 edge *par, *left, *right;
182 /*par here = sib before
183 left, right here = par, skew before*/
184 double D_AB, D_CD, D_AC, D_BD;
185 /*B is subtree being moved - below vtest
186 A is subtree remaining fixed - below va, unioned with all trees already passed by B*/
187 /*C is subtree being passed by B, in this case above par
188 D is subtree below etest, fixed on other side*/
189 par = etest->tail->parentEdge;
190 left = etest->head->leftEdge;
191 right = etest->head->rightEdge;
198 assignDownWeightsDown(left,vtest,va,etest,etest->tail,A[vtest->index][etest->tail->index],0.5,A,swapWeights);
199 assignDownWeightsDown(right,vtest,va,etest,etest->tail,A[vtest->index][etest->tail->index],0.5,A,swapWeights);
204 D_BD = A[vtest->index][etest->head->index];
205 D_CD = A[par->head->index][etest->head->index];
206 D_AC = A[back->head->index][par->head->index] + coeff*(A[va->index][par->head->index] - A[vtest->index][par->head->index]);
207 D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
208 swapWeights[0][vtest->index][etest->head->index] = swapWeights[0][vtest->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
211 assignDownWeightsDown(left,vtest, va, etest, etest->tail, D_AB, 0.5*coeff, A, swapWeights);
212 assignDownWeightsDown(right,vtest, va, etest, etest->tail, D_AB, 0.5*coeff, A, swapWeights);
217 void assignDownWeightsDown(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights)
219 /*again the same general idea*/
220 edge *sib, *left, *right;
221 /*sib here = par in assignDownWeightsSkew
222 rest is the same as assignDownWeightsSkew*/
223 double D_AB, D_CD, D_AC, D_BD;
224 /*B is below vtest, A is below va unioned with all trees already passed by B*/
225 /*C is subtree being passed - below sib*/
226 /*D is tree below etest*/
227 sib = siblingEdge(etest);
228 left = etest->head->leftEdge;
229 right = etest->head->rightEdge;
230 D_BD = A[vtest->index][etest->head->index];
231 D_CD = A[sib->head->index][etest->head->index];
232 D_AC = A[sib->head->index][back->head->index] + coeff*(A[sib->head->index][va->index] - A[sib->head->index][vtest->index]);
233 D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
234 swapWeights[0][vtest->index][etest->head->index] = swapWeights[0][vtest->index][back->head->index] + ( D_AC + D_BD - D_AB - D_CD);
237 assignDownWeightsDown(left,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
238 assignDownWeightsDown(right,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
242 void assignUpWeights(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A,
243 double ***swapWeights)
245 /*SPR performed on tree above vtest...*/
246 /*same idea as above, with appropriate selections of edges and nodes*/
247 edge *sib, *left, *right;
248 /*B is above vtest, A is other tree below vtest unioned with trees in path to vtest*/
249 /*sib is tree C being passed by B*/
250 /*D is tree below etest*/
251 double D_AB, D_CD, D_AC, D_BD;
253 sib = siblingEdge(etest);
254 left = etest->head->leftEdge;
255 right = etest->head->rightEdge;
256 if (NULL == back) /*first recursive call*/
260 else /*start the process of assigning weights recursively*/
262 assignUpWeights(left,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
263 assignUpWeights(right,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
266 else /*second or later recursive call*/
268 D_BD = A[vtest->index][etest->head->index];
269 D_CD = A[sib->head->index][etest->head->index];
270 D_AC = A[back->head->index][sib->head->index] + coeff*(A[va->index][sib->head->index] - A[vtest->index][sib->head->index]);
271 D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
272 thisWeight = swapWeights[1][vtest->index][etest->head->index] = swapWeights[1][vtest->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
275 assignUpWeights(left,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
276 assignUpWeights(right,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
281 void pruneSubtree(edge *p, edge *u, edge *d)
282 /*starting with edge u above edges p, d*/
283 /*removes p, d from tree, u connects to d->head to compensate*/
285 p->tail->parentEdge = NULL;/*remove p subtree*/
287 d->head->parentEdge = u; /*u connected to d->head*/
288 d->head = NULL; /*d removed from tree*/
291 void SPRsplitEdge(edge *e, edge *p, edge *d)
292 /*splits edge e to make it parental to p,d. d is parental to what
293 previously was below e*/
297 p->tail->parentEdge = e;
298 d->head->parentEdge = d;
302 /*topological shift function*/
303 /*removes subtree rooted at v and re-inserts to spilt e*/
304 void SPRDownShift(tree *T, node *v, edge *e)
306 edge *vup, *vdown, *vpar;
307 vpar = v->parentEdge;
308 vdown = siblingEdge(vpar);
309 vup = vpar->tail->parentEdge;
310 /*topological shift*/
311 pruneSubtree(vpar,vup,vdown);
312 /*removes v subtree and vdown, extends vup*/
313 SPRsplitEdge(e,vpar,vdown); /*splits e to make e sibling edge to vpar,
317 void SPRUpShift(tree *T, node *vmove, edge *esplit)
318 /*an inelegant iterative version*/
326 for(f=esplit->tail->parentEdge,pathLength=1;f->tail != vmove;f=f->tail->parentEdge)
328 /*count number of edges to vmove*/
329 /*note that pathLength > 0 will hold*/
331 EPath = (edge **)malloc(pathLength*sizeof(edge *));
332 v = (node **)malloc(pathLength*sizeof(edge *));
333 sib = (edge **)malloc((pathLength+1)*sizeof(edge *));
334 /*there are pathLength + 1 side trees, one at the head and tail of each edge in the path*/
336 sib[pathLength] = siblingEdge(esplit);
338 f = esplit->tail->parentEdge;
343 sib[i] = siblingEdge(f);
345 f = f->tail->parentEdge;
347 /*indexed so head of Epath is v[i], tail is v[i-1] and sibling edge is sib[i]*/
348 /*need to assign head, tail of each edge in path
349 as well as have proper values for the left and right fields*/
351 if (esplit == esplit->tail->leftEdge)
353 vmove->leftEdge = esplit;
354 vmove->rightEdge = EPath[pathLength-1];
358 vmove->rightEdge = esplit;
359 vmove->leftEdge = EPath[pathLength-1];
361 esplit->tail = vmove;
362 /*espilt->head remains unchanged*/
363 /*vmove has the proper fields for left, right, and parentEdge*/
365 for(i=0;i<(pathLength-1);i++)
366 EPath[i]->tail = v[i+1];
368 /*this bit flips the orientation along the path
369 the tail of Epath[i] is now v[i+1] instead of v[i-1]*/
371 EPath[pathLength-1]->tail = vmove;
373 for(i=1;i<pathLength;i++)
375 if (sib[i+1] == v[i]->leftEdge)
376 v[i]->rightEdge = EPath[i-1];
378 v[i]->leftEdge = EPath[i-1];
380 if (sib[1] == v[0]->leftEdge)
381 v[0]->rightEdge = sib[0];
383 v[0]->leftEdge = sib[0];
391 void SPRTopShift(tree *T, node *vmove, edge *esplit, int UpOrDown)
393 if (DOWN == UpOrDown)
394 SPRDownShift(T,vmove,esplit);
396 SPRUpShift(T,vmove,esplit);
399 node *indexedNode(tree *T, int i)
402 for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
403 if (i == e->head->index)
405 if (i == T->root->index)
410 edge *indexedEdge(tree *T, int i)
413 for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
414 if (i == e->head->index)