3 void BalWFext(edge *e, double **A) /*works except when e is the one edge
4 inserted to new vertex v by firstInsert*/
7 if ((leaf(e->head)) && (leaf(e->tail)))
8 e->distance = A[e->head->index][e->head->index];
9 else if (leaf(e->head))
11 f = e->tail->parentEdge;
13 e->distance = 0.5*(A[e->head->index][g->head->index]
14 + A[e->head->index][f->head->index]
15 - A[g->head->index][f->head->index]);
19 f = e->head->leftEdge;
20 g = e->head->rightEdge;
21 e->distance = 0.5*(A[g->head->index][e->head->index]
22 + A[f->head->index][e->head->index]
23 - A[f->head->index][g->head->index]);
27 void BalWFint(edge *e, double **A)
29 int up, down, left, right;
31 down = (siblingEdge(e))->head->index;
32 left = e->head->leftEdge->head->index;
33 right = e->head->rightEdge->head->index;
34 e->distance = 0.25*(A[up][left] + A[up][right] + A[left][down] + A[right][down]) - 0.5*(A[down][up] + A[left][right]);
37 void assignBMEWeights(tree *T, double **A)
40 e = depthFirstTraverse(T,NULL);
42 if ((leaf(e->head)) || (leaf(e->tail)))
46 e = depthFirstTraverse(T,e);
50 void BMEcalcDownAverage(tree *T, node *v, edge *e, double **D, double **A)
54 A[e->head->index][v->index] = D[v->index2][e->head->index2];
57 left = e->head->leftEdge;
58 right = e->head->rightEdge;
59 A[e->head->index][v->index] = 0.5 * A[left->head->index][v->index]
60 + 0.5 * A[right->head->index][v->index];
64 void BMEcalcUpAverage(tree *T, node *v, edge *e, double **D, double **A)
67 if (T->root == e->tail)
68 A[v->index][e->head->index] = D[v->index2][e->tail->index2];
69 /*for now, use convention
70 v->index first => looking up
71 v->index second => looking down */
74 up = e->tail->parentEdge;
75 down = siblingEdge(e);
76 A[v->index][e->head->index] = 0.5 * A[v->index][up->head->index]
77 +0.5 * A[down->head->index][v->index];
82 void BMEcalcNewvAverages(tree *T, node *v, double **D, double **A)
85 /*depth-first search*/
88 e = depthFirstTraverse(T,e); /*the downward averages need to be
89 calculated from bottom to top */
92 BMEcalcDownAverage(T,v,e,D,A);
93 e = depthFirstTraverse(T,e);
96 e = topFirstTraverse(T,e); /*the upward averages need to be calculated
100 BMEcalcUpAverage(T,v,e,D,A);
101 e = topFirstTraverse(T,e);
106 /*update Pair updates A[nearEdge][farEdge] and makes recursive call to subtree
108 /*root is head or tail of edge being split, depending on direction toward
110 void updatePair(double **A, edge *nearEdge, edge *farEdge, node *v,
111 node *root, double dcoeff, int direction)
114 switch(direction) /*the various cases refer to where the new vertex has
115 been inserted, in relation to the edge nearEdge*/
117 case UP: /*this case is called when v has been inserted above
119 /*do recursive calls first!*/
120 if (NULL != farEdge->head->leftEdge)
121 updatePair(A,nearEdge,farEdge->head->leftEdge,v,root,dcoeff,UP);
122 if (NULL != farEdge->head->rightEdge)
123 updatePair(A,nearEdge,farEdge->head->rightEdge,v,root,dcoeff,UP);
124 A[farEdge->head->index][nearEdge->head->index] =
125 A[nearEdge->head->index][farEdge->head->index]
126 = A[farEdge->head->index][nearEdge->head->index]
127 + dcoeff*A[farEdge->head->index][v->index]
128 - dcoeff*A[farEdge->head->index][root->index];
130 case DOWN: /*called when v has been inserted below farEdge*/
131 if (NULL != farEdge->tail->parentEdge)
132 updatePair(A,nearEdge,farEdge->tail->parentEdge,v,root,dcoeff,DOWN);
133 sib = siblingEdge(farEdge);
135 updatePair(A,nearEdge,sib,v,root,dcoeff,UP);
136 A[farEdge->head->index][nearEdge->head->index] =
137 A[nearEdge->head->index][farEdge->head->index]
138 = A[farEdge->head->index][nearEdge->head->index]
139 + dcoeff*A[v->index][farEdge->head->index]
140 - dcoeff*A[farEdge->head->index][root->index];
144 void updateSubTree(double **A, edge *nearEdge, node *v, node *root,
145 node *newNode, double dcoeff, int direction)
150 case UP: /*newNode is above the edge nearEdge*/
151 A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
152 A[newNode->index][nearEdge->head->index] =
153 A[nearEdge->head->index][newNode->index] =
154 A[nearEdge->head->index][root->index];
155 if (NULL != nearEdge->head->leftEdge)
156 updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff, UP);
157 if (NULL != nearEdge->head->rightEdge)
158 updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff, UP);
159 updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
161 case DOWN: /*newNode is below the edge nearEdge*/
162 A[nearEdge->head->index][v->index] = A[v->index][nearEdge->head->index];
163 A[newNode->index][nearEdge->head->index] =
164 A[nearEdge->head->index][newNode->index] =
165 0.5*(A[nearEdge->head->index][root->index]
166 + A[v->index][nearEdge->head->index]);
167 sib = siblingEdge(nearEdge);
169 updateSubTree(A, sib, v, root, newNode, 0.5*dcoeff, SKEW);
170 if (NULL != nearEdge->tail->parentEdge)
171 updateSubTree(A, nearEdge->tail->parentEdge, v, root, newNode, 0.5*dcoeff, DOWN);
172 updatePair(A, nearEdge, nearEdge, v, root, dcoeff, DOWN);
174 case SKEW: /*newNode is neither above nor below nearEdge*/
175 A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
176 A[newNode->index][nearEdge->head->index] =
177 A[nearEdge->head->index][newNode->index] =
178 0.5*(A[nearEdge->head->index][root->index] +
179 A[nearEdge->head->index][v->index]);
180 if (NULL != nearEdge->head->leftEdge)
181 updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff,SKEW);
182 if (NULL != nearEdge->head->rightEdge)
183 updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff,SKEW);
184 updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
189 /*we update all the averages for nodes (u1,u2), where the insertion point of
190 v is in "direction" from both u1 and u2 */
191 /*The general idea is to proceed in a direction from those edges already corrected
194 /*r is the root of the tree relative to the inserted node*/
196 void BMEupdateAveragesMatrix(double **A, edge *e, node *v,node *newNode)
198 edge *sib, *par, *left, *right;
199 /*first, update the v,newNode entries*/
200 A[newNode->index][newNode->index] = 0.5*(A[e->head->index][e->head->index]
201 + A[v->index][e->head->index]);
202 A[v->index][newNode->index] = A[newNode->index][v->index] =
203 A[v->index][e->head->index];
204 A[v->index][v->index] =
205 0.5*(A[e->head->index][v->index] + A[v->index][e->head->index]);
206 left = e->head->leftEdge;
207 right = e->head->rightEdge;
209 updateSubTree(A,left,v,e->head,newNode,0.25,UP); /*updates left and below*/
211 updateSubTree(A,right,v,e->head,newNode,0.25,UP); /*updates right and below*/
212 sib = siblingEdge(e);
214 updateSubTree(A,sib,v,e->head,newNode,0.25,SKEW); /*updates sib and below*/
215 par = e->tail->parentEdge;
217 updateSubTree(A,par,v,e->head,newNode,0.25,DOWN); /*updates par and above*/
219 /*must change values A[e->head][*] last, as they are used to update
220 the rest of the matrix*/
221 A[newNode->index][e->head->index] = A[e->head->index][newNode->index]
222 = A[e->head->index][e->head->index];
223 A[v->index][e->head->index] = A[e->head->index][v->index];
225 updatePair(A,e,e,v,e->head,0.5,UP); /*updates e->head fields only*/
228 /*A is tree below sibling, B is tree below edge, C is tree above edge*/
229 double wf3(double D_AB, double D_AC, double D_kB, double D_kC)
231 return(D_AC + D_kB - D_AB - D_kC);
234 void BMEtestEdge(edge *e, node *v, double **A)
237 down = siblingEdge(e);
238 up = e->tail->parentEdge;
239 e->totalweight = wf3(A[e->head->index][down->head->index],
240 A[down->head->index][e->tail->index],
241 A[e->head->index][v->index],
242 A[v->index][e->tail->index])
246 void BMEsplitEdge(tree *T, node *v, edge *e, double **A)
248 edge *newPendantEdge;
249 edge *newInternalEdge;
251 char nodeLabel[NODE_LABEL_LENGTH];
252 char edgeLabel1[EDGE_LABEL_LENGTH];
253 char edgeLabel2[EDGE_LABEL_LENGTH];
254 snprintf(nodeLabel,1,"");
255 //sprintf(edgeLabel1,"E%d",T->size);
256 //sprintf(edgeLabel2,"E%d",T->size+1);
257 snprintf(edgeLabel1,EDGE_LABEL_LENGTH,"E%d",T->size);
258 snprintf(edgeLabel2,EDGE_LABEL_LENGTH,"E%d",T->size+1);
260 /*make the new node and edges*/
261 newNode = makeNewNode(nodeLabel,T->size+1);
262 newPendantEdge = makeEdge(edgeLabel1,newNode,v,0.0);
263 newInternalEdge = makeEdge(edgeLabel2,newNode,e->head,0.0);
265 /*update the matrix of average distances*/
266 BMEupdateAveragesMatrix(A,e,v,newNode);
268 /*put them in the correct topology*/
269 newNode->parentEdge = e;
270 e->head->parentEdge = newInternalEdge;
271 v->parentEdge = newPendantEdge;
274 T->size = T->size + 2;
276 if (e->tail->leftEdge == e)
277 /*actually this is totally arbitrary and probably unnecessary*/
279 newNode->leftEdge = newInternalEdge;
280 newNode->rightEdge = newPendantEdge;
284 newNode->leftEdge = newInternalEdge;
285 newNode->rightEdge = newPendantEdge;
289 tree *BMEaddSpecies(tree *T,node *v, double **D, double **A)
290 /*the key function of the program addSpeices inserts
291 the node v to the tree T. It uses testEdge to see what the relative
292 weight would be if v split a particular edge. Once insertion point
293 is found, v is added to T, and A is updated. Edge weights
294 are not assigned until entire tree is build*/
298 edge *e; /*loop variable*/
299 edge *e_min; /*points to best edge seen thus far*/
300 double w_min = 0.0; /*used to keep track of tree weights*/
302 /*initialize variables as necessary*/
304 /*CASE 1: T is empty, v is the first node*/
305 if (NULL == T) /*create a tree with v as only vertex, no edges*/
309 /*note that we are rooting T arbitrarily at a leaf.
310 T->root is not the phylogenetic root*/
315 /*CASE 2: T is a single-vertex tree*/
319 e = makeEdge("",T->root,v,0.0);
320 //sprintf(e->label,"E1");
321 snprintf(e->label,EDGE_LABEL_LENGTH,"E1");
322 A[v->index][v->index] = D[v->index2][T->root->index2];
323 T->root->leftEdge = v->parentEdge = e;
327 /*CASE 3: T has at least two nodes and an edge. Insert new node
328 by breaking one of the edges*/
331 BMEcalcNewvAverages(T,v,D,A);
332 /*calcNewvAverages will update A for the row and column
333 include the node v. Will do so using pre-existing averages in T and
334 information from A,D*/
335 e_min = T->root->leftEdge;
336 e = e_min->head->leftEdge;
340 /*testEdge tests weight of tree if loop variable
341 e is the edge split, places this value in the e->totalweight field */
342 if (e->totalweight < w_min)
345 w_min = e->totalweight;
347 e = topFirstTraverse(T,e);
349 /*e_min now points at the edge we want to split*/
351 printf("Inserting %s between %s and %s on %s\n",v->label,e_min->tail->label,
352 e_min->head->label,e_min->label);*/
353 BMEsplitEdge(T,v,e_min,A);
357 /*calcUpAverages will ensure that A[e->head->index][f->head->index] is
358 filled for any f >= g. Works recursively*/
359 void calcUpAverages(double **D, double **A, edge *e, edge *g)
363 if (!(leaf(g->tail)))
365 calcUpAverages(D,A,e,g->tail->parentEdge);
369 A[e->head->index][g->head->index] = A[g->head->index][e->head->index]
370 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
374 void makeBMEAveragesTable(tree *T, double **D, double **A)
376 edge *e, *f, *exclude;
378 /*first, let's deal with the averages involving the root of T*/
379 e = T->root->leftEdge;
380 f = depthFirstTraverse(T,NULL);
383 A[e->head->index][f->head->index] = A[f->head->index][e->head->index]
384 = D[e->tail->index2][f->head->index2];
388 u = f->head->leftEdge->head;
389 v = f->head->rightEdge->head;
390 A[e->head->index][f->head->index] = A[f->head->index][e->head->index]
391 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
393 f = depthFirstTraverse(T,f);
395 e = depthFirstTraverse(T,NULL);
396 while (T->root->leftEdge != e) {
398 while (T->root->leftEdge != f) {
400 exclude = exclude->tail->parentEdge;
401 else if (leaf(e->head))
404 A[e->head->index][f->head->index] =
405 A[f->head->index][e->head->index]
406 = D[e->head->index2][f->head->index2];
409 u = f->head->leftEdge->head; /*since f is chosen using a
410 depth-first search, other values
411 have been calculated*/
412 v = f->head->rightEdge->head;
413 A[e->head->index][f->head->index]
414 = A[f->head->index][e->head->index]
415 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
420 u = e->head->leftEdge->head;
421 v = e->head->rightEdge->head;
422 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]);
424 f = depthFirstTraverse(T,f);
426 e = depthFirstTraverse(T,e);
428 e = depthFirstTraverse(T,NULL);
429 while (T->root->leftEdge != e)
431 calcUpAverages(D,A,e,e); /*calculates averages for
432 A[e->head->index][g->head->index] for
433 any edge g in path from e to root of tree*/
434 e = depthFirstTraverse(T,e);
436 } /*makeAveragesMatrix*/