1 /* bNNI.c 2007-09-05 */
3 /* Copyright 2007 Vincent Lefort */
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
10 /*boolean leaf(node *v);
11 edge *siblingEdge(edge *e);
12 edge *depthFirstTraverse(tree *T, edge *e);
13 edge *findBottomLeft(edge *e);
14 edge *topFirstTraverse(tree *T, edge *e);
15 edge *moveUpRight(edge *e);*/
17 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger);
18 void assignBMEWeights(tree *T, double **A);
19 //void updateAveragesMatrix(tree *T, double **A, node *v,int direction);
20 void bNNItopSwitch(tree *T, edge *e, int direction, double **A);
21 int bNNIEdgeTest(edge *e, tree *T, double **A, double *weight);
22 void updatePair(double **A, edge *nearEdge, edge *farEdge, node *closer, node *further, double dcoeff, int direction);
24 int *initPerm(int size);
26 void reHeapElement(int *p, int *q, double *v, int length, int i);
27 void pushHeap(int *p, int *q, double *v, int length, int i);
28 void popHeap(int *p, int *q, double *v, int length, int i);
31 void bNNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray,
32 double *weights, int *location, int *possibleSwaps)
35 tloc = location[e->head->index+1];
36 location[e->head->index+1] =
37 bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
38 if (NONE == location[e->head->index+1])
41 popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
46 pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
48 reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
52 int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
54 void permInverse(int *p, int *q, int length);
56 void weighTree(tree *T)
60 for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
61 T->weight += e->distance;
64 //void bNNI(tree *T, double **avgDistArray, int *count)
65 void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
69 int *p, *location, *q;
73 p = initPerm(T->size+1);
74 q = initPerm(T->size+1);
75 edgeArray = (edge **) malloc((T->size+1)*sizeof(double));
76 weights = (double *) malloc((T->size+1)*sizeof(double));
77 location = (int *) malloc((T->size+1)*sizeof(int));
80 for (i=0; i<numSpecies; i++)
81 for (j=0; j<numSpecies; j++)
83 epsilon = (epsilon / (numSpecies * numSpecies)) * EPSILON;
85 for (i=0;i<T->size+1;i++)
92 assignBMEWeights(T,avgDistArray);
95 e = findBottomLeft(T->root->leftEdge);
98 edgeArray[e->head->index+1] = e;
99 location[e->head->index+1] =
100 bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
101 e = depthFirstTraverse(T,e);
103 possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
104 permInverse(p,q,T->size+1);
105 /*we put the negative values of weights into a heap, indexed by p
106 with the minimum value pointed to by p[1]*/
107 /*p[i] is index (in edgeArray) of edge with i-th position
108 in the heap, q[j] is the position of edge j in the heap */
109 while (weights[p[1]] + epsilon < 0)
111 centerEdge = edgeArray[p[1]];
115 T->weight = T->weight + weights[p[1]];
116 printf("New tree weight is %lf.\n",T->weight);
118 bNNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
119 location[p[1]] = NONE;
120 weights[p[1]] = 0.0; /*after the bNNI, this edge is in optimal
122 popHeap(p,q,weights,possibleSwaps--,1);
123 /*but we must retest the other edges of T*/
124 /*CHANGE 2/28/2003 expanding retesting to _all_ edges of T*/
125 e = depthFirstTraverse(T,NULL);
128 bNNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
129 e = depthFirstTraverse(T,e);
137 assignBMEWeights(T,avgDistArray);
141 /*This function is the meat of the average distance matrix recalculation*/
142 /*Idea is: we are looking at the subtree rooted at rootEdge. The subtree
143 rooted at closer is closer to rootEdge after the NNI, while the subtree
144 rooted at further is further to rootEdge after the NNI. direction tells
145 the direction of the NNI with respect to rootEdge*/
146 void updateSubTreeAfterNNI(double **A, node *v, edge *rootEdge, node *closer, node *further,
147 double dcoeff, int direction)
152 case UP: /*rootEdge is below the center edge of the NNI*/
153 /*recursive calls to subtrees, if necessary*/
154 if (NULL != rootEdge->head->leftEdge)
155 updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,UP);
156 if (NULL != rootEdge->head->rightEdge)
157 updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,UP);
158 updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
159 sib = siblingEdge(v->parentEdge);
160 A[rootEdge->head->index][v->index] =
161 A[v->index][rootEdge->head->index] =
162 0.5*A[rootEdge->head->index][sib->head->index] +
163 0.5*A[rootEdge->head->index][v->parentEdge->tail->index];
165 case DOWN: /*rootEdge is above the center edge of the NNI*/
166 sib = siblingEdge(rootEdge);
168 updateSubTreeAfterNNI(A, v, sib, closer, further, 0.5*dcoeff, SKEW);
169 if (NULL != rootEdge->tail->parentEdge)
170 updateSubTreeAfterNNI(A, v, rootEdge->tail->parentEdge, closer, further, 0.5*dcoeff, DOWN);
171 updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, DOWN);
172 A[rootEdge->head->index][v->index] =
173 A[v->index][rootEdge->head->index] =
174 0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
175 0.5*A[rootEdge->head->index][v->rightEdge->head->index];
177 case SKEW: /*rootEdge is in subtree skew to v*/
178 if (NULL != rootEdge->head->leftEdge)
179 updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,SKEW);
180 if (NULL != rootEdge->head->rightEdge)
181 updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,SKEW);
182 updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
183 A[rootEdge->head->index][v->index] =
184 A[v->index][rootEdge->head->index] =
185 0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
186 0.5*A[rootEdge->head->index][v->rightEdge->head->index];
191 /*swapping across edge whose head is v*/
192 void bNNIupdateAverages(double **A, node *v, edge *par, edge *skew,
193 edge *swap, edge *fixed)
195 A[v->index][v->index] = 0.25*(A[fixed->head->index][par->head->index] +
196 A[fixed->head->index][swap->head->index] +
197 A[skew->head->index][par->head->index] +
198 A[skew->head->index][swap->head->index]);
199 updateSubTreeAfterNNI(A, v, fixed, skew->head, swap->head, 0.25, UP);
200 updateSubTreeAfterNNI(A, v, par, swap->head, skew->head, 0.25, DOWN);
201 updateSubTreeAfterNNI(A, v, skew, fixed->head, par->head, 0.25, UP);
202 updateSubTreeAfterNNI(A, v, swap, par->head, fixed->head, 0.25, SKEW);
206 void bNNItopSwitch(tree *T, edge *e, int direction, double **A)
208 edge *down, *swap, *fixed;
212 printf("Performing branch swap across edge %s ",e->label);
214 if (LEFT == direction)
216 else printf("right ");
217 printf("subtree.\n");
219 down = siblingEdge(e);
222 if (LEFT == direction)
224 swap = e->head->leftEdge;
225 fixed = e->head->rightEdge;
230 swap = e->head->rightEdge;
231 fixed = e->head->leftEdge;
236 if(e->tail->leftEdge == e)
240 bNNIupdateAverages(A, v, e->tail->parentEdge, down, swap, fixed);
243 double wf5(double D_AD, double D_BC, double D_AC, double D_BD,
244 double D_AB, double D_CD)
247 weight = 0.25*(D_AC + D_BD + D_AD + D_BC) + 0.5*(D_AB + D_CD);
251 int bNNIEdgeTest(edge *e, tree *T, double **A, double *weight)
254 double D_LR, D_LU, D_LD, D_RD, D_RU, D_DU;
257 printf("Branch swap: testing edge %s.\n",e->label);*/
258 if ((leaf(e->tail)) || (leaf(e->head)))
263 D_LR = A[e->head->leftEdge->head->index][e->head->rightEdge->head->index];
264 D_LU = A[e->head->leftEdge->head->index][e->tail->index];
265 D_LD = A[e->head->leftEdge->head->index][f->head->index];
266 D_RU = A[e->head->rightEdge->head->index][e->tail->index];
267 D_RD = A[e->head->rightEdge->head->index][f->head->index];
268 D_DU = A[e->tail->index][f->head->index];
270 w0 = wf5(D_RU,D_LD,D_LU,D_RD,D_DU,D_LR); /*weight of current config*/
271 w1 = wf5(D_RU,D_LD,D_DU,D_LR,D_LU,D_RD); /*weight with L<->D switch*/
272 w2 = wf5(D_DU,D_LR,D_LU,D_RD,D_RU,D_LD); /*weight with R<->D switch*/
275 if (w0 <= w2) /*w0 <= w1,w2*/
280 else /*w2 < w0 <= w1 */
285 printf("Possible swap across %s. ",e->label);
286 printf("Weight dropping by %lf.\n",w0 - w2);
287 printf("New weight would be %lf.\n",T->weight + w2 - w0);
292 else if (w2 <= w1) /*w2 <= w1 < w0*/
297 printf("Possible swap across %s. ",e->label);
298 printf("Weight dropping by %lf.\n",w0 - w2);
299 printf("New weight should be %lf.\n",T->weight + w2 - w0);
308 printf("Possible swap across %s. ",e->label);
309 printf("Weight dropping by %lf.\n",w0 - w1);
310 printf("New weight should be %lf.\n",T->weight + w1 - w0);
316 /*limitedFillTableUp fills all the entries in D associated with
317 e->head,f->head and those edges g->head above e->head, working
318 recursively and stopping when trigger is reached*/
319 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger)
322 g = f->tail->parentEdge;
324 limitedFillTableUp(e,g,A,trigger);
326 A[e->head->index][f->head->index] =
327 A[f->head->index][e->head->index] =
328 0.5*(A[e->head->index][g->head->index] + A[e->head->index][h->head->index]);