3 /* Copyright 2007-2008 Olivier Gascuel, Rick Desper,
4 R port by Vincent Lefort, me_*() below modified by Emmanuel Paradis */
6 /* This file is part of the R-package `ape'. */
7 /* See the file ../COPYING for licensing issues. */
11 //functions from me_balanced.c
12 tree *BMEaddSpecies(tree *T, node *v, double **D, double **A);
13 void assignBMEWeights(tree *T, double **A);
14 void makeBMEAveragesTable(tree *T, double **D, double **A);
15 //functions from me_ols.c
16 tree *GMEaddSpecies(tree *T, node *v, double **D, double **A);
17 void assignOLSWeights(tree *T, double **A);
18 void makeOLSAveragesTable(tree *T, double **D, double **A);
19 //functions from bNNI.c
20 void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies);
21 //functions from NNI.c
22 void NNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies);
25 void me_b(double *X, int *N, char **labels, int *nni,
26 int *edge1, int *edge2, double *el, char **tl)
29 set *species, *slooper;
37 species = (set *) malloc(sizeof(set));
38 species->firstNode = NULL;
39 species->secondNode = NULL;
40 D = loadMatrix(X, labels, n, species);
41 A = initDoubleMatrix(2*n - 2);
43 for(slooper = species; NULL != slooper; slooper = slooper->secondNode)
45 addNode = copyNode(slooper->firstNode);
46 T = BMEaddSpecies(T, addNode, D, A);
50 bNNI(T, A, &nniCount, D, n);
51 assignBMEWeights(T,A);
53 tree2phylo(T, edge1, edge2, el, tl, n);
56 freeMatrix(A,2*n - 2);
62 void me_o(double *X, int *N, char **labels, int *nni,
63 int *edge1, int *edge2, double *el, char **tl)
66 set *species, *slooper;
74 species = (set *) malloc(sizeof(set));
75 species->firstNode = NULL;
76 species->secondNode = NULL;
78 D = loadMatrix (X, labels, n, species);
79 A = initDoubleMatrix(2 * n - 2);
81 for(slooper = species; NULL != slooper; slooper = slooper->secondNode)
83 addNode = copyNode(slooper->firstNode);
84 T = GMEaddSpecies(T,addNode,D,A);
86 makeOLSAveragesTable(T,D,A);
89 NNI(T,A,&nniCount,D,n);
90 assignOLSWeights(T,A);
92 tree2phylo(T, edge1, edge2, el, tl, n);
95 freeMatrix(A,2*n - 2);
101 /*************************************************************************
105 *************************************************************************/
107 double **initDoubleMatrix(int d)
111 A = (double **) malloc(d*sizeof(double *));
114 A[i] = (double *) malloc(d*sizeof(double));
121 double **loadMatrix (double *X, char **labels, int n, set *S)
123 char nextString[MAX_LABEL_LENGTH];
128 table = (double **) calloc(n,sizeof(double *));
130 table[i] = (double *) calloc(n,sizeof(double));
134 strncpy (nextString, labels[i], MAX_LABEL_LENGTH);
135 // ReplaceForbiddenChars (nextString, '_');
136 v = makeNewNode(nextString,-1);
139 for (j=i; j<n; j++) {
142 table[j][i] = X[XINDEX(a,b)];
143 table[i][j] = X[XINDEX(a,b)];
151 /*************************************************************************
155 *************************************************************************/
157 set *addToSet(node *v, set *X)
161 X = (set *) malloc(sizeof(set));
163 X->secondNode = NULL;
165 else if (NULL == X->firstNode)
168 X->secondNode = addToSet(v,X->secondNode);
172 node *makeNewNode(char *label, int i)
174 return(makeNode(label,NULL,i));
177 node *makeNode(char *label, edge *parentEdge, int index)
179 node *newNode; /*points to new node added to the graph*/
180 newNode = (node *) malloc(sizeof(node));
181 strncpy(newNode->label,label,NODE_LABEL_LENGTH);
182 newNode->index = index;
183 newNode->index2 = -1;
184 newNode->parentEdge = parentEdge;
185 newNode->leftEdge = NULL;
186 newNode->middleEdge = NULL;
187 newNode->rightEdge = NULL;
188 /*all fields have been initialized*/
192 /*copyNode returns a copy of v which has all of the fields identical to those
193 of v, except the node pointer fields*/
194 node *copyNode(node *v)
197 w = makeNode(v->label,NULL,v->index);
198 w->index2 = v->index2;
202 edge *siblingEdge(edge *e)
204 if(e == e->tail->leftEdge)
205 return(e->tail->rightEdge);
207 return(e->tail->leftEdge);
210 edge *makeEdge(char *label, node *tail, node *head, double weight)
213 newEdge = (edge *) malloc(sizeof(edge));
214 strncpy(newEdge->label,label,EDGE_LABEL_LENGTH);
215 newEdge->tail = tail;
216 newEdge->head = head;
217 newEdge->distance = weight;
218 newEdge->totalweight = 0.0;
225 T = (tree *) malloc(sizeof(tree));
232 void updateSizes(edge *e, int direction)
238 f = e->head->leftEdge;
241 f = e->head->rightEdge;
250 f = e->tail->parentEdge;
258 /*detrifurcate takes the (possibly trifurcated) input tree
259 and reroots the tree to a leaf*/
260 /*assumes tree is only trifurcated at root*/
261 tree *detrifurcate(tree *T)
268 if (NULL != v->parentEdge)
270 Rprintf ("Error: root %s is poorly rooted.\n",v->label);
273 for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f )
282 w->parentEdge = NULL;
288 void compareSets(tree *T, set *S)
293 e = depthFirstTraverse(T,NULL);
297 for(X = S; NULL != X; X = X->secondNode)
300 if (0 == strcmp(v->label,w->label))
302 v->index2 = w->index2;
307 e = depthFirstTraverse(T,e);
310 for(X = S; NULL != X; X = X->secondNode)
313 if (0 == strcmp(v->label,w->label))
315 v->index2 = w->index2;
322 Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
325 e = depthFirstTraverse(T,NULL);
329 if ((leaf(v)) && (-1 == v->index2))
331 Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
334 e = depthFirstTraverse(T,e);
336 for(X = S; NULL != X; X = X->secondNode)
337 if (X->firstNode->index2 > -1)
339 Rprintf("Error node %s in matrix but not a leaf in tree.\n",X->firstNode->label);
345 void partitionSizes(tree *T)
348 e = depthFirstTraverse(T,NULL);
354 e->bottomsize = e->head->leftEdge->bottomsize
355 + e->head->rightEdge->bottomsize;
356 e->topsize = (T->size + 2)/2 - e->bottomsize;
357 e = depthFirstTraverse(T,e);
361 /*************************************************************************
365 *************************************************************************/
367 edge *depthFirstTraverse(tree *T, edge *e)
368 /*depthFirstTraverse returns the edge f which is least in T according
369 to the depth-first order, but which is later than e in the search
370 pattern. If e is null, f is the least edge of T*/
375 f = T->root->leftEdge;
377 f = findBottomLeft(f);
378 return(f); /*this is the first edge of this search pattern*/
380 else /*e is non-null*/
382 if (e->tail->leftEdge == e)
383 /*if e is a left-oriented edge, we skip the entire
384 tree cut below e, and find least edge*/
386 else /*if e is a right-oriented edge, we have already looked at its
387 sibling and everything below e, so we move up*/
388 f = e->tail->parentEdge;
393 edge *findBottomLeft(edge *e)
394 /*findBottomLeft searches by gottom down in the tree and to the left.*/
398 while (NULL != f->head->leftEdge)
399 f = f->head->leftEdge;
403 edge *moveRight(edge *e)
406 f = e->tail->rightEdge; /*this step moves from a left-oriented edge
407 to a right-oriented edge*/
409 f = findBottomLeft(f);
413 edge *topFirstTraverse(tree *T, edge *e)
414 /*topFirstTraverse starts from the top of T, and from there moves stepwise
415 down, left before right*/
416 /*assumes tree has been detrifurcated*/
420 return(T->root->leftEdge); /*first Edge searched*/
421 else if (!(leaf(e->head)))
422 return(e->head->leftEdge); /*down and to the left is preferred*/
423 else /*e->head is a leaf*/
430 edge *moveUpRight(edge *e)
434 while ((NULL != f) && ( f->tail->leftEdge != f))
435 f = f->tail->parentEdge;
436 /*go up the tree until f is a leftEdge*/
438 return(f); /*triggered at end of search*/
440 return(f->tail->rightEdge);
441 /*and then go right*/
444 /*************************************************************************
448 *************************************************************************/
450 void freeMatrix(double **D, int size)
461 freeSet(S->secondNode);
465 void freeTree(tree *T)
469 if (NULL != v->leftEdge)
470 freeSubTree(v->leftEdge);
475 void freeSubTree(edge *e)