6 void BalWFext(edge *e, double **A) /*works except when e is the one edge
7 inserted to new vertex v by firstInsert*/
10 if ((leaf(e->head)) && (leaf(e->tail)))
11 e->distance = A[e->head->index][e->head->index];
12 else if (leaf(e->head))
14 f = e->tail->parentEdge;
16 e->distance = 0.5*(A[e->head->index][g->head->index]
17 + A[e->head->index][f->head->index]
18 - A[g->head->index][f->head->index]);
22 f = e->head->leftEdge;
23 g = e->head->rightEdge;
24 e->distance = 0.5*(A[g->head->index][e->head->index]
25 + A[f->head->index][e->head->index]
26 - A[f->head->index][g->head->index]);
30 void BalWFint(edge *e, double **A)
32 int up, down, left, right;
34 down = (siblingEdge(e))->head->index;
35 left = e->head->leftEdge->head->index;
36 right = e->head->rightEdge->head->index;
37 e->distance = 0.25*(A[up][left] + A[up][right] + A[left][down] + A[right][down]) - 0.5*(A[down][up] + A[left][right]);
40 void assignBMEWeights(tree *T, double **A)
43 e = depthFirstTraverse(T,NULL);
45 if ((leaf(e->head)) || (leaf(e->tail)))
49 e = depthFirstTraverse(T,e);
53 void BMEcalcDownAverage(tree *T, node *v, edge *e, double **D, double **A)
57 A[e->head->index][v->index] = D[v->index2][e->head->index2];
60 left = e->head->leftEdge;
61 right = e->head->rightEdge;
62 A[e->head->index][v->index] = 0.5 * A[left->head->index][v->index]
63 + 0.5 * A[right->head->index][v->index];
67 void BMEcalcUpAverage(tree *T, node *v, edge *e, double **D, double **A)
70 if (T->root == e->tail)
71 A[v->index][e->head->index] = D[v->index2][e->tail->index2];
72 /*for now, use convention
73 v->index first => looking up
74 v->index second => looking down */
77 up = e->tail->parentEdge;
78 down = siblingEdge(e);
79 A[v->index][e->head->index] = 0.5 * A[v->index][up->head->index]
80 +0.5 * A[down->head->index][v->index];
85 void BMEcalcNewvAverages(tree *T, node *v, double **D, double **A)
88 /*depth-first search*/
91 e = depthFirstTraverse(T,e); /*the downward averages need to be
92 calculated from bottom to top */
95 BMEcalcDownAverage(T,v,e,D,A);
96 e = depthFirstTraverse(T,e);
99 e = topFirstTraverse(T,e); /*the upward averages need to be calculated
100 from top to bottom */
103 BMEcalcUpAverage(T,v,e,D,A);
104 e = topFirstTraverse(T,e);
109 /*update Pair updates A[nearEdge][farEdge] and makes recursive call to subtree
111 /*root is head or tail of edge being split, depending on direction toward
113 void updatePair(double **A, edge *nearEdge, edge *farEdge, node *v,
114 node *root, double dcoeff, int direction)
117 switch(direction) /*the various cases refer to where the new vertex has
118 been inserted, in relation to the edge nearEdge*/
120 case UP: /*this case is called when v has been inserted above
122 /*do recursive calls first!*/
123 if (NULL != farEdge->head->leftEdge)
124 updatePair(A,nearEdge,farEdge->head->leftEdge,v,root,dcoeff,UP);
125 if (NULL != farEdge->head->rightEdge)
126 updatePair(A,nearEdge,farEdge->head->rightEdge,v,root,dcoeff,UP);
127 A[farEdge->head->index][nearEdge->head->index] =
128 A[nearEdge->head->index][farEdge->head->index]
129 = A[farEdge->head->index][nearEdge->head->index]
130 + dcoeff*A[farEdge->head->index][v->index]
131 - dcoeff*A[farEdge->head->index][root->index];
133 case DOWN: /*called when v has been inserted below farEdge*/
134 if (NULL != farEdge->tail->parentEdge)
135 updatePair(A,nearEdge,farEdge->tail->parentEdge,v,root,dcoeff,DOWN);
136 sib = siblingEdge(farEdge);
138 updatePair(A,nearEdge,sib,v,root,dcoeff,UP);
139 A[farEdge->head->index][nearEdge->head->index] =
140 A[nearEdge->head->index][farEdge->head->index]
141 = A[farEdge->head->index][nearEdge->head->index]
142 + dcoeff*A[v->index][farEdge->head->index]
143 - dcoeff*A[farEdge->head->index][root->index];
147 void updateSubTree(double **A, edge *nearEdge, node *v, node *root,
148 node *newNode, double dcoeff, int direction)
153 case UP: /*newNode is above the edge nearEdge*/
154 A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
155 A[newNode->index][nearEdge->head->index] =
156 A[nearEdge->head->index][newNode->index] =
157 A[nearEdge->head->index][root->index];
158 if (NULL != nearEdge->head->leftEdge)
159 updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff, UP);
160 if (NULL != nearEdge->head->rightEdge)
161 updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff, UP);
162 updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
164 case DOWN: /*newNode is below the edge nearEdge*/
165 A[nearEdge->head->index][v->index] = A[v->index][nearEdge->head->index];
166 A[newNode->index][nearEdge->head->index] =
167 A[nearEdge->head->index][newNode->index] =
168 0.5*(A[nearEdge->head->index][root->index]
169 + A[v->index][nearEdge->head->index]);
170 sib = siblingEdge(nearEdge);
172 updateSubTree(A, sib, v, root, newNode, 0.5*dcoeff, SKEW);
173 if (NULL != nearEdge->tail->parentEdge)
174 updateSubTree(A, nearEdge->tail->parentEdge, v, root, newNode, 0.5*dcoeff, DOWN);
175 updatePair(A, nearEdge, nearEdge, v, root, dcoeff, DOWN);
177 case SKEW: /*newNode is neither above nor below nearEdge*/
178 A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
179 A[newNode->index][nearEdge->head->index] =
180 A[nearEdge->head->index][newNode->index] =
181 0.5*(A[nearEdge->head->index][root->index] +
182 A[nearEdge->head->index][v->index]);
183 if (NULL != nearEdge->head->leftEdge)
184 updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff,SKEW);
185 if (NULL != nearEdge->head->rightEdge)
186 updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff,SKEW);
187 updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
192 /*we update all the averages for nodes (u1,u2), where the insertion point of
193 v is in "direction" from both u1 and u2 */
194 /*The general idea is to proceed in a direction from those edges already corrected
197 /*r is the root of the tree relative to the inserted node*/
199 void BMEupdateAveragesMatrix(double **A, edge *e, node *v,node *newNode)
201 edge *sib, *par, *left, *right;
202 /*first, update the v,newNode entries*/
203 A[newNode->index][newNode->index] = 0.5*(A[e->head->index][e->head->index]
204 + A[v->index][e->head->index]);
205 A[v->index][newNode->index] = A[newNode->index][v->index] =
206 A[v->index][e->head->index];
207 A[v->index][v->index] =
208 0.5*(A[e->head->index][v->index] + A[v->index][e->head->index]);
209 left = e->head->leftEdge;
210 right = e->head->rightEdge;
212 updateSubTree(A,left,v,e->head,newNode,0.25,UP); /*updates left and below*/
214 updateSubTree(A,right,v,e->head,newNode,0.25,UP); /*updates right and below*/
215 sib = siblingEdge(e);
217 updateSubTree(A,sib,v,e->head,newNode,0.25,SKEW); /*updates sib and below*/
218 par = e->tail->parentEdge;
220 updateSubTree(A,par,v,e->head,newNode,0.25,DOWN); /*updates par and above*/
222 /*must change values A[e->head][*] last, as they are used to update
223 the rest of the matrix*/
224 A[newNode->index][e->head->index] = A[e->head->index][newNode->index]
225 = A[e->head->index][e->head->index];
226 A[v->index][e->head->index] = A[e->head->index][v->index];
228 updatePair(A,e,e,v,e->head,0.5,UP); /*updates e->head fields only*/
231 /*A is tree below sibling, B is tree below edge, C is tree above edge*/
232 double wf3(double D_AB, double D_AC, double D_kB, double D_kC)
234 return(D_AC + D_kB - D_AB - D_kC);
237 void BMEtestEdge(edge *e, node *v, double **A)
240 down = siblingEdge(e);
241 up = e->tail->parentEdge;
242 e->totalweight = wf3(A[e->head->index][down->head->index],
243 A[down->head->index][e->tail->index],
244 A[e->head->index][v->index],
245 A[v->index][e->tail->index])
249 void BMEsplitEdge(tree *T, node *v, edge *e, double **A)
251 edge *newPendantEdge;
252 edge *newInternalEdge;
254 char nodeLabel[NODE_LABEL_LENGTH];
255 char edgeLabel1[EDGE_LABEL_LENGTH];
256 char edgeLabel2[EDGE_LABEL_LENGTH];
257 snprintf(nodeLabel,1,"");
258 //sprintf(edgeLabel1,"E%d",T->size);
259 //sprintf(edgeLabel2,"E%d",T->size+1);
260 snprintf(edgeLabel1,EDGE_LABEL_LENGTH,"E%d",T->size);
261 snprintf(edgeLabel2,EDGE_LABEL_LENGTH,"E%d",T->size+1);
264 /*make the new node and edges*/
265 newNode = makeNewNode(nodeLabel,T->size+1);
266 newPendantEdge = makeEdge(edgeLabel1,newNode,v,0.0);
267 newInternalEdge = makeEdge(edgeLabel2,newNode,e->head,0.0);
269 /*update the matrix of average distances*/
270 BMEupdateAveragesMatrix(A,e,v,newNode);
272 /*put them in the correct topology*/
273 newNode->parentEdge = e;
274 e->head->parentEdge = newInternalEdge;
275 v->parentEdge = newPendantEdge;
278 T->size = T->size + 2;
280 if (e->tail->leftEdge == e)
281 /*actually this is totally arbitrary and probably unnecessary*/
283 newNode->leftEdge = newInternalEdge;
284 newNode->rightEdge = newPendantEdge;
288 newNode->leftEdge = newInternalEdge;
289 newNode->rightEdge = newPendantEdge;
293 tree *BMEaddSpecies(tree *T,node *v, double **D, double **A)
294 /*the key function of the program addSpeices inserts
295 the node v to the tree T. It uses testEdge to see what the relative
296 weight would be if v split a particular edge. Once insertion point
297 is found, v is added to T, and A is updated. Edge weights
298 are not assigned until entire tree is build*/
302 edge *e; /*loop variable*/
303 edge *e_min; /*points to best edge seen thus far*/
304 double w_min = 0.0; /*used to keep track of tree weights*/
306 /*initialize variables as necessary*/
308 /*CASE 1: T is empty, v is the first node*/
309 if (NULL == T) /*create a tree with v as only vertex, no edges*/
313 /*note that we are rooting T arbitrarily at a leaf.
314 T->root is not the phylogenetic root*/
319 /*CASE 2: T is a single-vertex tree*/
323 e = makeEdge("",T->root,v,0.0);
324 //sprintf(e->label,"E1");
325 snprintf(e->label,EDGE_LABEL_LENGTH,"E1");
326 A[v->index][v->index] = D[v->index2][T->root->index2];
327 T->root->leftEdge = v->parentEdge = e;
331 /*CASE 3: T has at least two nodes and an edge. Insert new node
332 by breaking one of the edges*/
335 BMEcalcNewvAverages(T,v,D,A);
336 /*calcNewvAverages will update A for the row and column
337 include the node v. Will do so using pre-existing averages in T and
338 information from A,D*/
339 e_min = T->root->leftEdge;
340 e = e_min->head->leftEdge;
344 /*testEdge tests weight of tree if loop variable
345 e is the edge split, places this value in the e->totalweight field */
346 if (e->totalweight < w_min)
349 w_min = e->totalweight;
351 e = topFirstTraverse(T,e);
353 /*e_min now points at the edge we want to split*/
355 printf("Inserting %s between %s and %s on %s\n",v->label,e_min->tail->label,
356 e_min->head->label,e_min->label);*/
357 BMEsplitEdge(T,v,e_min,A);
361 /*calcUpAverages will ensure that A[e->head->index][f->head->index] is
362 filled for any f >= g. Works recursively*/
363 void calcUpAverages(double **D, double **A, edge *e, edge *g)
367 if (!(leaf(g->tail)))
369 calcUpAverages(D,A,e,g->tail->parentEdge);
373 A[e->head->index][g->head->index] = A[g->head->index][e->head->index]
374 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
378 void makeBMEAveragesTable(tree *T, double **D, double **A)
380 edge *e, *f, *exclude;
382 /*first, let's deal with the averages involving the root of T*/
383 e = T->root->leftEdge;
384 f = depthFirstTraverse(T,NULL);
387 A[e->head->index][f->head->index] = A[f->head->index][e->head->index]
388 = D[e->tail->index2][f->head->index2];
392 u = f->head->leftEdge->head;
393 v = f->head->rightEdge->head;
394 A[e->head->index][f->head->index] = A[f->head->index][e->head->index]
395 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
397 f = depthFirstTraverse(T,f);
399 e = depthFirstTraverse(T,NULL);
400 while (T->root->leftEdge != e) {
402 while (T->root->leftEdge != f) {
404 exclude = exclude->tail->parentEdge;
405 else if (leaf(e->head))
408 A[e->head->index][f->head->index] =
409 A[f->head->index][e->head->index]
410 = D[e->head->index2][f->head->index2];
413 u = f->head->leftEdge->head; /*since f is chosen using a
414 depth-first search, other values
415 have been calculated*/
416 v = f->head->rightEdge->head;
417 A[e->head->index][f->head->index]
418 = A[f->head->index][e->head->index]
419 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
424 u = e->head->leftEdge->head;
425 v = e->head->rightEdge->head;
426 A[e->head->index][f->head->index] = A[f->head->index][e->head->index] = 0.5*(A[f->head->index][u->index] + A[f->head->index][v->index]);
428 f = depthFirstTraverse(T,f);
430 e = depthFirstTraverse(T,e);
432 e = depthFirstTraverse(T,NULL);
433 while (T->root->leftEdge != e)
435 calcUpAverages(D,A,e,e); /*calculates averages for
436 A[e->head->index][g->head->index] for
437 any edge g in path from e to root of tree*/
438 e = depthFirstTraverse(T,e);
440 } /*makeAveragesMatrix*/