3 /*boolean leaf(node *v);
4 edge *siblingEdge(edge *e);
5 edge *depthFirstTraverse(tree *T, edge *e);
6 edge *findBottomLeft(edge *e);
7 edge *topFirstTraverse(tree *T, edge *e);
8 edge *moveUpRight(edge *e);*/
10 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger);
11 void assignBMEWeights(tree *T, double **A);
12 //void updateAveragesMatrix(tree *T, double **A, node *v,int direction);
13 void bNNItopSwitch(tree *T, edge *e, int direction, double **A);
14 int bNNIEdgeTest(edge *e, tree *T, double **A, double *weight);
15 void updatePair(double **A, edge *nearEdge, edge *farEdge, node *closer, node *further, double dcoeff, int direction);
17 int *initPerm(int size);
19 void reHeapElement(int *p, int *q, double *v, int length, int i);
20 void pushHeap(int *p, int *q, double *v, int length, int i);
21 void popHeap(int *p, int *q, double *v, int length, int i);
24 void bNNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray,
25 double *weights, int *location, int *possibleSwaps)
28 tloc = location[e->head->index+1];
29 location[e->head->index+1] =
30 bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
31 if (NONE == location[e->head->index+1])
34 popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
39 pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
41 reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
45 int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
47 void permInverse(int *p, int *q, int length);
49 void weighTree(tree *T)
53 for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
54 T->weight += e->distance;
57 //void bNNI(tree *T, double **avgDistArray, int *count)
58 void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
60 edge *e; /* , *centerEdge; */
62 int *p, *location, *q;
66 p = initPerm(T->size+1);
67 q = initPerm(T->size+1);
68 edgeArray = (edge **) malloc((T->size+1)*sizeof(double));
69 weights = (double *) malloc((T->size+1)*sizeof(double));
70 location = (int *) malloc((T->size+1)*sizeof(int));
73 for (i=0; i<numSpecies; i++)
74 for (j=0; j<numSpecies; j++)
76 epsilon = (epsilon / (numSpecies * numSpecies)) * EPSILON;
78 for (i=0;i<T->size+1;i++)
85 assignBMEWeights(T,avgDistArray);
88 e = findBottomLeft(T->root->leftEdge);
91 edgeArray[e->head->index+1] = e;
92 location[e->head->index+1] =
93 bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
94 e = depthFirstTraverse(T,e);
96 possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
97 permInverse(p,q,T->size+1);
98 /*we put the negative values of weights into a heap, indexed by p
99 with the minimum value pointed to by p[1]*/
100 /*p[i] is index (in edgeArray) of edge with i-th position
101 in the heap, q[j] is the position of edge j in the heap */
102 while (weights[p[1]] + epsilon < 0)
104 centerEdge = edgeArray[p[1]];
108 T->weight = T->weight + weights[p[1]];
109 printf("New tree weight is %lf.\n",T->weight);
111 bNNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
112 location[p[1]] = NONE;
113 weights[p[1]] = 0.0; /*after the bNNI, this edge is in optimal
115 popHeap(p,q,weights,possibleSwaps--,1);
116 /*but we must retest the other edges of T*/
117 /*CHANGE 2/28/2003 expanding retesting to _all_ edges of T*/
118 e = depthFirstTraverse(T,NULL);
121 bNNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
122 e = depthFirstTraverse(T,e);
130 assignBMEWeights(T,avgDistArray);
134 /*This function is the meat of the average distance matrix recalculation*/
135 /*Idea is: we are looking at the subtree rooted at rootEdge. The subtree
136 rooted at closer is closer to rootEdge after the NNI, while the subtree
137 rooted at further is further to rootEdge after the NNI. direction tells
138 the direction of the NNI with respect to rootEdge*/
139 void updateSubTreeAfterNNI(double **A, node *v, edge *rootEdge, node *closer, node *further,
140 double dcoeff, int direction)
145 case UP: /*rootEdge is below the center edge of the NNI*/
146 /*recursive calls to subtrees, if necessary*/
147 if (NULL != rootEdge->head->leftEdge)
148 updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,UP);
149 if (NULL != rootEdge->head->rightEdge)
150 updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,UP);
151 updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
152 sib = siblingEdge(v->parentEdge);
153 A[rootEdge->head->index][v->index] =
154 A[v->index][rootEdge->head->index] =
155 0.5*A[rootEdge->head->index][sib->head->index] +
156 0.5*A[rootEdge->head->index][v->parentEdge->tail->index];
158 case DOWN: /*rootEdge is above the center edge of the NNI*/
159 sib = siblingEdge(rootEdge);
161 updateSubTreeAfterNNI(A, v, sib, closer, further, 0.5*dcoeff, SKEW);
162 if (NULL != rootEdge->tail->parentEdge)
163 updateSubTreeAfterNNI(A, v, rootEdge->tail->parentEdge, closer, further, 0.5*dcoeff, DOWN);
164 updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, DOWN);
165 A[rootEdge->head->index][v->index] =
166 A[v->index][rootEdge->head->index] =
167 0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
168 0.5*A[rootEdge->head->index][v->rightEdge->head->index];
170 case SKEW: /*rootEdge is in subtree skew to v*/
171 if (NULL != rootEdge->head->leftEdge)
172 updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,SKEW);
173 if (NULL != rootEdge->head->rightEdge)
174 updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,SKEW);
175 updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
176 A[rootEdge->head->index][v->index] =
177 A[v->index][rootEdge->head->index] =
178 0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
179 0.5*A[rootEdge->head->index][v->rightEdge->head->index];
184 /*swapping across edge whose head is v*/
185 void bNNIupdateAverages(double **A, node *v, edge *par, edge *skew,
186 edge *swap, edge *fixed)
188 A[v->index][v->index] = 0.25*(A[fixed->head->index][par->head->index] +
189 A[fixed->head->index][swap->head->index] +
190 A[skew->head->index][par->head->index] +
191 A[skew->head->index][swap->head->index]);
192 updateSubTreeAfterNNI(A, v, fixed, skew->head, swap->head, 0.25, UP);
193 updateSubTreeAfterNNI(A, v, par, swap->head, skew->head, 0.25, DOWN);
194 updateSubTreeAfterNNI(A, v, skew, fixed->head, par->head, 0.25, UP);
195 updateSubTreeAfterNNI(A, v, swap, par->head, fixed->head, 0.25, SKEW);
199 void bNNItopSwitch(tree *T, edge *e, int direction, double **A)
201 edge *down, *swap, *fixed;
205 printf("Performing branch swap across edge %s ",e->label);
207 if (LEFT == direction)
209 else printf("right ");
210 printf("subtree.\n");
212 down = siblingEdge(e);
215 if (LEFT == direction)
217 swap = e->head->leftEdge;
218 fixed = e->head->rightEdge;
223 swap = e->head->rightEdge;
224 fixed = e->head->leftEdge;
229 if(e->tail->leftEdge == e)
233 bNNIupdateAverages(A, v, e->tail->parentEdge, down, swap, fixed);
236 double wf5(double D_AD, double D_BC, double D_AC, double D_BD,
237 double D_AB, double D_CD)
240 weight = 0.25*(D_AC + D_BD + D_AD + D_BC) + 0.5*(D_AB + D_CD);
244 int bNNIEdgeTest(edge *e, tree *T, double **A, double *weight)
247 double D_LR, D_LU, D_LD, D_RD, D_RU, D_DU;
250 printf("Branch swap: testing edge %s.\n",e->label);*/
251 if ((leaf(e->tail)) || (leaf(e->head)))
256 D_LR = A[e->head->leftEdge->head->index][e->head->rightEdge->head->index];
257 D_LU = A[e->head->leftEdge->head->index][e->tail->index];
258 D_LD = A[e->head->leftEdge->head->index][f->head->index];
259 D_RU = A[e->head->rightEdge->head->index][e->tail->index];
260 D_RD = A[e->head->rightEdge->head->index][f->head->index];
261 D_DU = A[e->tail->index][f->head->index];
263 w0 = wf5(D_RU,D_LD,D_LU,D_RD,D_DU,D_LR); /*weight of current config*/
264 w1 = wf5(D_RU,D_LD,D_DU,D_LR,D_LU,D_RD); /*weight with L<->D switch*/
265 w2 = wf5(D_DU,D_LR,D_LU,D_RD,D_RU,D_LD); /*weight with R<->D switch*/
268 if (w0 <= w2) /*w0 <= w1,w2*/
273 else /*w2 < w0 <= w1 */
278 printf("Possible swap across %s. ",e->label);
279 printf("Weight dropping by %lf.\n",w0 - w2);
280 printf("New weight would be %lf.\n",T->weight + w2 - w0);
285 else if (w2 <= w1) /*w2 <= w1 < w0*/
290 printf("Possible swap across %s. ",e->label);
291 printf("Weight dropping by %lf.\n",w0 - w2);
292 printf("New weight should be %lf.\n",T->weight + w2 - w0);
301 printf("Possible swap across %s. ",e->label);
302 printf("Weight dropping by %lf.\n",w0 - w1);
303 printf("New weight should be %lf.\n",T->weight + w1 - w0);
309 /*limitedFillTableUp fills all the entries in D associated with
310 e->head,f->head and those edges g->head above e->head, working
311 recursively and stopping when trigger is reached*/
312 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger)
315 g = f->tail->parentEdge;
317 limitedFillTableUp(e,g,A,trigger);
319 A[e->head->index][f->head->index] =
320 A[f->head->index][e->head->index] =
321 0.5*(A[e->head->index][g->head->index] + A[e->head->index][h->head->index]);