1 /* me_balanced.c 2012-04-30 */
3 /* Copyright 2007 Vincent Lefort
4 BMEsplitEdge() modified by Emmanuel Paradis */
6 /* This file is part of the R-package `ape'. */
7 /* See the file ../COPYING for licensing issues. */
11 void BalWFext(edge *e, double **A) /*works except when e is the one edge
12 inserted to new vertex v by firstInsert*/
15 if ((leaf(e->head)) && (leaf(e->tail)))
16 e->distance = A[e->head->index][e->head->index];
17 else if (leaf(e->head))
19 f = e->tail->parentEdge;
21 e->distance = 0.5*(A[e->head->index][g->head->index]
22 + A[e->head->index][f->head->index]
23 - A[g->head->index][f->head->index]);
27 f = e->head->leftEdge;
28 g = e->head->rightEdge;
29 e->distance = 0.5*(A[g->head->index][e->head->index]
30 + A[f->head->index][e->head->index]
31 - A[f->head->index][g->head->index]);
35 void BalWFint(edge *e, double **A)
37 int up, down, left, right;
39 down = (siblingEdge(e))->head->index;
40 left = e->head->leftEdge->head->index;
41 right = e->head->rightEdge->head->index;
42 e->distance = 0.25*(A[up][left] + A[up][right] + A[left][down] + A[right][down]) - 0.5*(A[down][up] + A[left][right]);
45 void assignBMEWeights(tree *T, double **A)
48 e = depthFirstTraverse(T,NULL);
50 if ((leaf(e->head)) || (leaf(e->tail)))
54 e = depthFirstTraverse(T,e);
58 void BMEcalcDownAverage(tree *T, node *v, edge *e, double **D, double **A)
62 A[e->head->index][v->index] = D[v->index2][e->head->index2];
65 left = e->head->leftEdge;
66 right = e->head->rightEdge;
67 A[e->head->index][v->index] = 0.5 * A[left->head->index][v->index]
68 + 0.5 * A[right->head->index][v->index];
72 void BMEcalcUpAverage(tree *T, node *v, edge *e, double **D, double **A)
75 if (T->root == e->tail)
76 A[v->index][e->head->index] = D[v->index2][e->tail->index2];
77 /*for now, use convention
78 v->index first => looking up
79 v->index second => looking down */
82 up = e->tail->parentEdge;
83 down = siblingEdge(e);
84 A[v->index][e->head->index] = 0.5 * A[v->index][up->head->index]
85 +0.5 * A[down->head->index][v->index];
90 void BMEcalcNewvAverages(tree *T, node *v, double **D, double **A)
93 /*depth-first search*/
96 e = depthFirstTraverse(T,e); /*the downward averages need to be
97 calculated from bottom to top */
100 BMEcalcDownAverage(T,v,e,D,A);
101 e = depthFirstTraverse(T,e);
104 e = topFirstTraverse(T,e); /*the upward averages need to be calculated
105 from top to bottom */
108 BMEcalcUpAverage(T,v,e,D,A);
109 e = topFirstTraverse(T,e);
114 /*update Pair updates A[nearEdge][farEdge] and makes recursive call to subtree
116 /*root is head or tail of edge being split, depending on direction toward
118 void updatePair(double **A, edge *nearEdge, edge *farEdge, node *v,
119 node *root, double dcoeff, int direction)
122 switch(direction) /*the various cases refer to where the new vertex has
123 been inserted, in relation to the edge nearEdge*/
125 case UP: /*this case is called when v has been inserted above
127 /*do recursive calls first!*/
128 if (NULL != farEdge->head->leftEdge)
129 updatePair(A,nearEdge,farEdge->head->leftEdge,v,root,dcoeff,UP);
130 if (NULL != farEdge->head->rightEdge)
131 updatePair(A,nearEdge,farEdge->head->rightEdge,v,root,dcoeff,UP);
132 A[farEdge->head->index][nearEdge->head->index] =
133 A[nearEdge->head->index][farEdge->head->index]
134 = A[farEdge->head->index][nearEdge->head->index]
135 + dcoeff*A[farEdge->head->index][v->index]
136 - dcoeff*A[farEdge->head->index][root->index];
138 case DOWN: /*called when v has been inserted below farEdge*/
139 if (NULL != farEdge->tail->parentEdge)
140 updatePair(A,nearEdge,farEdge->tail->parentEdge,v,root,dcoeff,DOWN);
141 sib = siblingEdge(farEdge);
143 updatePair(A,nearEdge,sib,v,root,dcoeff,UP);
144 A[farEdge->head->index][nearEdge->head->index] =
145 A[nearEdge->head->index][farEdge->head->index]
146 = A[farEdge->head->index][nearEdge->head->index]
147 + dcoeff*A[v->index][farEdge->head->index]
148 - dcoeff*A[farEdge->head->index][root->index];
152 void updateSubTree(double **A, edge *nearEdge, node *v, node *root,
153 node *newNode, double dcoeff, int direction)
158 case UP: /*newNode is above the edge nearEdge*/
159 A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
160 A[newNode->index][nearEdge->head->index] =
161 A[nearEdge->head->index][newNode->index] =
162 A[nearEdge->head->index][root->index];
163 if (NULL != nearEdge->head->leftEdge)
164 updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff, UP);
165 if (NULL != nearEdge->head->rightEdge)
166 updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff, UP);
167 updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
169 case DOWN: /*newNode is below the edge nearEdge*/
170 A[nearEdge->head->index][v->index] = A[v->index][nearEdge->head->index];
171 A[newNode->index][nearEdge->head->index] =
172 A[nearEdge->head->index][newNode->index] =
173 0.5*(A[nearEdge->head->index][root->index]
174 + A[v->index][nearEdge->head->index]);
175 sib = siblingEdge(nearEdge);
177 updateSubTree(A, sib, v, root, newNode, 0.5*dcoeff, SKEW);
178 if (NULL != nearEdge->tail->parentEdge)
179 updateSubTree(A, nearEdge->tail->parentEdge, v, root, newNode, 0.5*dcoeff, DOWN);
180 updatePair(A, nearEdge, nearEdge, v, root, dcoeff, DOWN);
182 case SKEW: /*newNode is neither above nor below nearEdge*/
183 A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
184 A[newNode->index][nearEdge->head->index] =
185 A[nearEdge->head->index][newNode->index] =
186 0.5*(A[nearEdge->head->index][root->index] +
187 A[nearEdge->head->index][v->index]);
188 if (NULL != nearEdge->head->leftEdge)
189 updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff,SKEW);
190 if (NULL != nearEdge->head->rightEdge)
191 updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff,SKEW);
192 updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
197 /*we update all the averages for nodes (u1,u2), where the insertion point of
198 v is in "direction" from both u1 and u2 */
199 /*The general idea is to proceed in a direction from those edges already corrected
202 /*r is the root of the tree relative to the inserted node*/
204 void BMEupdateAveragesMatrix(double **A, edge *e, node *v,node *newNode)
206 edge *sib, *par, *left, *right;
207 /*first, update the v,newNode entries*/
208 A[newNode->index][newNode->index] = 0.5*(A[e->head->index][e->head->index]
209 + A[v->index][e->head->index]);
210 A[v->index][newNode->index] = A[newNode->index][v->index] =
211 A[v->index][e->head->index];
212 A[v->index][v->index] =
213 0.5*(A[e->head->index][v->index] + A[v->index][e->head->index]);
214 left = e->head->leftEdge;
215 right = e->head->rightEdge;
217 updateSubTree(A,left,v,e->head,newNode,0.25,UP); /*updates left and below*/
219 updateSubTree(A,right,v,e->head,newNode,0.25,UP); /*updates right and below*/
220 sib = siblingEdge(e);
222 updateSubTree(A,sib,v,e->head,newNode,0.25,SKEW); /*updates sib and below*/
223 par = e->tail->parentEdge;
225 updateSubTree(A,par,v,e->head,newNode,0.25,DOWN); /*updates par and above*/
227 /*must change values A[e->head][*] last, as they are used to update
228 the rest of the matrix*/
229 A[newNode->index][e->head->index] = A[e->head->index][newNode->index]
230 = A[e->head->index][e->head->index];
231 A[v->index][e->head->index] = A[e->head->index][v->index];
233 updatePair(A,e,e,v,e->head,0.5,UP); /*updates e->head fields only*/
236 /*A is tree below sibling, B is tree below edge, C is tree above edge*/
237 double wf3(double D_AB, double D_AC, double D_kB, double D_kC)
239 return(D_AC + D_kB - D_AB - D_kC);
242 void BMEtestEdge(edge *e, node *v, double **A)
245 down = siblingEdge(e);
246 up = e->tail->parentEdge;
247 e->totalweight = wf3(A[e->head->index][down->head->index],
248 A[down->head->index][e->tail->index],
249 A[e->head->index][v->index],
250 A[v->index][e->tail->index])
254 void BMEsplitEdge(tree *T, node *v, edge *e, double **A)
256 edge *newPendantEdge;
257 edge *newInternalEdge;
259 int nodeLabel = 0;//char nodeLabel[NODE_LABEL_LENGTH];
260 char edgeLabel1[EDGE_LABEL_LENGTH];
261 char edgeLabel2[EDGE_LABEL_LENGTH];
262 //snprintf(nodeLabel,1,"");
263 //sprintf(edgeLabel1,"E%d",T->size);
264 //sprintf(edgeLabel2,"E%d",T->size+1);
265 snprintf(edgeLabel1,EDGE_LABEL_LENGTH,"E%d",T->size);
266 snprintf(edgeLabel2,EDGE_LABEL_LENGTH,"E%d",T->size+1);
268 /*make the new node and edges*/
269 newNode = makeNewNode(nodeLabel,T->size+1);
270 newPendantEdge = makeEdge(edgeLabel1,newNode,v,0.0);
271 newInternalEdge = makeEdge(edgeLabel2,newNode,e->head,0.0);
273 /*update the matrix of average distances*/
274 BMEupdateAveragesMatrix(A,e,v,newNode);
276 /*put them in the correct topology*/
277 newNode->parentEdge = e;
278 e->head->parentEdge = newInternalEdge;
279 v->parentEdge = newPendantEdge;
282 T->size = T->size + 2;
284 if (e->tail->leftEdge == e)
285 /*actually this is totally arbitrary and probably unnecessary*/
287 newNode->leftEdge = newInternalEdge;
288 newNode->rightEdge = newPendantEdge;
292 newNode->leftEdge = newInternalEdge;
293 newNode->rightEdge = newPendantEdge;
297 tree *BMEaddSpecies(tree *T,node *v, double **D, double **A)
298 /*the key function of the program addSpeices inserts
299 the node v to the tree T. It uses testEdge to see what the relative
300 weight would be if v split a particular edge. Once insertion point
301 is found, v is added to T, and A is updated. Edge weights
302 are not assigned until entire tree is build*/
306 edge *e; /*loop variable*/
307 edge *e_min; /*points to best edge seen thus far*/
308 double w_min = 0.0; /*used to keep track of tree weights*/
310 /*initialize variables as necessary*/
312 /*CASE 1: T is empty, v is the first node*/
313 if (NULL == T) /*create a tree with v as only vertex, no edges*/
317 /*note that we are rooting T arbitrarily at a leaf.
318 T->root is not the phylogenetic root*/
323 /*CASE 2: T is a single-vertex tree*/
327 e = makeEdge("",T->root,v,0.0);
328 //sprintf(e->label,"E1");
329 snprintf(e->label,EDGE_LABEL_LENGTH,"E1");
330 A[v->index][v->index] = D[v->index2][T->root->index2];
331 T->root->leftEdge = v->parentEdge = e;
335 /*CASE 3: T has at least two nodes and an edge. Insert new node
336 by breaking one of the edges*/
339 BMEcalcNewvAverages(T,v,D,A);
340 /*calcNewvAverages will update A for the row and column
341 include the node v. Will do so using pre-existing averages in T and
342 information from A,D*/
343 e_min = T->root->leftEdge;
344 e = e_min->head->leftEdge;
348 /*testEdge tests weight of tree if loop variable
349 e is the edge split, places this value in the e->totalweight field */
350 if (e->totalweight < w_min)
353 w_min = e->totalweight;
355 e = topFirstTraverse(T,e);
357 /*e_min now points at the edge we want to split*/
359 printf("Inserting %s between %s and %s on %s\n",v->label,e_min->tail->label,
360 e_min->head->label,e_min->label);*/
361 BMEsplitEdge(T,v,e_min,A);
365 /*calcUpAverages will ensure that A[e->head->index][f->head->index] is
366 filled for any f >= g. Works recursively*/
367 void calcUpAverages(double **D, double **A, edge *e, edge *g)
371 if (!(leaf(g->tail)))
373 calcUpAverages(D,A,e,g->tail->parentEdge);
377 A[e->head->index][g->head->index] = A[g->head->index][e->head->index]
378 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
382 void makeBMEAveragesTable(tree *T, double **D, double **A)
384 edge *e, *f, *exclude;
386 /*first, let's deal with the averages involving the root of T*/
387 e = T->root->leftEdge;
388 f = depthFirstTraverse(T,NULL);
391 A[e->head->index][f->head->index] = A[f->head->index][e->head->index]
392 = D[e->tail->index2][f->head->index2];
396 u = f->head->leftEdge->head;
397 v = f->head->rightEdge->head;
398 A[e->head->index][f->head->index] = A[f->head->index][e->head->index]
399 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
401 f = depthFirstTraverse(T,f);
403 e = depthFirstTraverse(T,NULL);
404 while (T->root->leftEdge != e) {
406 while (T->root->leftEdge != f) {
408 exclude = exclude->tail->parentEdge;
409 else if (leaf(e->head))
412 A[e->head->index][f->head->index] =
413 A[f->head->index][e->head->index]
414 = D[e->head->index2][f->head->index2];
417 u = f->head->leftEdge->head; /*since f is chosen using a
418 depth-first search, other values
419 have been calculated*/
420 v = f->head->rightEdge->head;
421 A[e->head->index][f->head->index]
422 = A[f->head->index][e->head->index]
423 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
428 u = e->head->leftEdge->head;
429 v = e->head->rightEdge->head;
430 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]);
432 f = depthFirstTraverse(T,f);
434 e = depthFirstTraverse(T,e);
436 e = depthFirstTraverse(T,NULL);
437 while (T->root->leftEdge != e)
439 calcUpAverages(D,A,e,e); /*calculates averages for
440 A[e->head->index][g->head->index] for
441 any edge g in path from e to root of tree*/
442 e = depthFirstTraverse(T,e);
444 } /*makeAveragesMatrix*/