9 /*boolean leaf(node *v);
10 edge *siblingEdge(edge *e);
11 edge *depthFirstTraverse(tree *T, edge *e);
12 edge *findBottomLeft(edge *e);
13 edge *topFirstTraverse(tree *T, edge *e);
14 edge *moveUpRight(edge *e);*/
16 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger);
17 void assignBMEWeights(tree *T, double **A);
18 //void updateAveragesMatrix(tree *T, double **A, node *v,int direction);
19 void bNNItopSwitch(tree *T, edge *e, int direction, double **A);
20 int bNNIEdgeTest(edge *e, tree *T, double **A, double *weight);
21 void updatePair(double **A, edge *nearEdge, edge *farEdge, node *closer, node *further, double dcoeff, int direction);
23 int *initPerm(int size);
25 void reHeapElement(int *p, int *q, double *v, int length, int i);
26 void pushHeap(int *p, int *q, double *v, int length, int i);
27 void popHeap(int *p, int *q, double *v, int length, int i);
30 void bNNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray,
31 double *weights, int *location, int *possibleSwaps)
34 tloc = location[e->head->index+1];
35 location[e->head->index+1] =
36 bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
37 if (NONE == location[e->head->index+1])
40 popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
45 pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
47 reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
51 int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
53 void permInverse(int *p, int *q, int length);
55 void weighTree(tree *T)
59 for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
60 T->weight += e->distance;
63 //void bNNI(tree *T, double **avgDistArray, int *count)
64 void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
68 int *p, *location, *q;
72 p = initPerm(T->size+1);
73 q = initPerm(T->size+1);
74 edgeArray = (edge **) malloc((T->size+1)*sizeof(double));
75 weights = (double *) malloc((T->size+1)*sizeof(double));
76 location = (int *) malloc((T->size+1)*sizeof(int));
79 for (i=0; i<numSpecies; i++)
80 for (j=0; j<numSpecies; j++)
82 epsilon = (epsilon / (numSpecies * numSpecies)) * EPSILON;
84 for (i=0;i<T->size+1;i++)
91 assignBMEWeights(T,avgDistArray);
94 e = findBottomLeft(T->root->leftEdge);
97 edgeArray[e->head->index+1] = e;
98 location[e->head->index+1] =
99 bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
100 e = depthFirstTraverse(T,e);
102 possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
103 permInverse(p,q,T->size+1);
104 /*we put the negative values of weights into a heap, indexed by p
105 with the minimum value pointed to by p[1]*/
106 /*p[i] is index (in edgeArray) of edge with i-th position
107 in the heap, q[j] is the position of edge j in the heap */
108 while (weights[p[1]] + epsilon < 0)
110 centerEdge = edgeArray[p[1]];
114 T->weight = T->weight + weights[p[1]];
115 printf("New tree weight is %lf.\n",T->weight);
117 bNNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
118 location[p[1]] = NONE;
119 weights[p[1]] = 0.0; /*after the bNNI, this edge is in optimal
121 popHeap(p,q,weights,possibleSwaps--,1);
122 /*but we must retest the other edges of T*/
123 /*CHANGE 2/28/2003 expanding retesting to _all_ edges of T*/
124 e = depthFirstTraverse(T,NULL);
127 bNNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
128 e = depthFirstTraverse(T,e);
136 assignBMEWeights(T,avgDistArray);
140 /*This function is the meat of the average distance matrix recalculation*/
141 /*Idea is: we are looking at the subtree rooted at rootEdge. The subtree
142 rooted at closer is closer to rootEdge after the NNI, while the subtree
143 rooted at further is further to rootEdge after the NNI. direction tells
144 the direction of the NNI with respect to rootEdge*/
145 void updateSubTreeAfterNNI(double **A, node *v, edge *rootEdge, node *closer, node *further,
146 double dcoeff, int direction)
151 case UP: /*rootEdge is below the center edge of the NNI*/
152 /*recursive calls to subtrees, if necessary*/
153 if (NULL != rootEdge->head->leftEdge)
154 updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,UP);
155 if (NULL != rootEdge->head->rightEdge)
156 updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,UP);
157 updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
158 sib = siblingEdge(v->parentEdge);
159 A[rootEdge->head->index][v->index] =
160 A[v->index][rootEdge->head->index] =
161 0.5*A[rootEdge->head->index][sib->head->index] +
162 0.5*A[rootEdge->head->index][v->parentEdge->tail->index];
164 case DOWN: /*rootEdge is above the center edge of the NNI*/
165 sib = siblingEdge(rootEdge);
167 updateSubTreeAfterNNI(A, v, sib, closer, further, 0.5*dcoeff, SKEW);
168 if (NULL != rootEdge->tail->parentEdge)
169 updateSubTreeAfterNNI(A, v, rootEdge->tail->parentEdge, closer, further, 0.5*dcoeff, DOWN);
170 updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, DOWN);
171 A[rootEdge->head->index][v->index] =
172 A[v->index][rootEdge->head->index] =
173 0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
174 0.5*A[rootEdge->head->index][v->rightEdge->head->index];
176 case SKEW: /*rootEdge is in subtree skew to v*/
177 if (NULL != rootEdge->head->leftEdge)
178 updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,SKEW);
179 if (NULL != rootEdge->head->rightEdge)
180 updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,SKEW);
181 updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
182 A[rootEdge->head->index][v->index] =
183 A[v->index][rootEdge->head->index] =
184 0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
185 0.5*A[rootEdge->head->index][v->rightEdge->head->index];
190 /*swapping across edge whose head is v*/
191 void bNNIupdateAverages(double **A, node *v, edge *par, edge *skew,
192 edge *swap, edge *fixed)
194 A[v->index][v->index] = 0.25*(A[fixed->head->index][par->head->index] +
195 A[fixed->head->index][swap->head->index] +
196 A[skew->head->index][par->head->index] +
197 A[skew->head->index][swap->head->index]);
198 updateSubTreeAfterNNI(A, v, fixed, skew->head, swap->head, 0.25, UP);
199 updateSubTreeAfterNNI(A, v, par, swap->head, skew->head, 0.25, DOWN);
200 updateSubTreeAfterNNI(A, v, skew, fixed->head, par->head, 0.25, UP);
201 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);
317 /*limitedFillTableUp fills all the entries in D associated with
318 e->head,f->head and those edges g->head above e->head, working
319 recursively and stopping when trigger is reached*/
320 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger)
323 g = f->tail->parentEdge;
325 limitedFillTableUp(e,g,A,trigger);
327 A[e->head->index][f->head->index] =
328 A[f->head->index][e->head->index] =
329 0.5*(A[e->head->index][g->head->index] + A[e->head->index][h->head->index]);