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);
23 //functions from SPR.c
24 void SPR(tree *T, double **D, double **A, int *count);
25 //functions from TBR.c
26 void TBR(tree *T, double **D, double **A);
29 void me_b(double *X, int *N, char **labels,
30 int *nni, int *spr, int *tbr,
31 int *edge1, int *edge2, double *el, char **tl)
34 set *species, *slooper;
42 species = (set *) malloc(sizeof(set));
43 species->firstNode = NULL;
44 species->secondNode = NULL;
45 D = loadMatrix(X, labels, n, species);
46 A = initDoubleMatrix(2*n - 2);
48 for(slooper = species; NULL != slooper; slooper = slooper->secondNode)
50 addNode = copyNode(slooper->firstNode);
51 T = BMEaddSpecies(T, addNode, D, A);
54 if (*nni) bNNI(T, A, &nniCount, D, n);
55 assignBMEWeights(T,A);
57 if (*spr) SPR(T, D, A, &nniCount);
58 if (*tbr) TBR(T, D, A);
60 tree2phylo(T, edge1, edge2, el, tl, n);
63 freeMatrix(A,2*n - 2);
69 void me_o(double *X, int *N, char **labels, int *nni,
70 int *edge1, int *edge2, double *el, char **tl)
73 set *species, *slooper;
81 species = (set *) malloc(sizeof(set));
82 species->firstNode = NULL;
83 species->secondNode = NULL;
85 D = loadMatrix (X, labels, n, species);
86 A = initDoubleMatrix(2 * n - 2);
88 for(slooper = species; NULL != slooper; slooper = slooper->secondNode)
90 addNode = copyNode(slooper->firstNode);
91 T = GMEaddSpecies(T,addNode,D,A);
93 makeOLSAveragesTable(T,D,A);
96 NNI(T,A,&nniCount,D,n);
97 assignOLSWeights(T,A);
99 tree2phylo(T, edge1, edge2, el, tl, n);
102 freeMatrix(A,2*n - 2);
108 /*************************************************************************
112 *************************************************************************/
114 double **initDoubleMatrix(int d)
118 A = (double **) malloc(d*sizeof(double *));
121 A[i] = (double *) malloc(d*sizeof(double));
128 double **loadMatrix (double *X, char **labels, int n, set *S)
130 char nextString[MAX_LABEL_LENGTH];
135 table = (double **) calloc(n,sizeof(double *));
137 table[i] = (double *) calloc(n,sizeof(double));
141 strncpy (nextString, labels[i], MAX_LABEL_LENGTH);
142 // ReplaceForbiddenChars (nextString, '_');
143 v = makeNewNode(nextString,-1);
146 for (j=i; j<n; j++) {
149 table[j][i] = X[XINDEX(a,b)];
150 table[i][j] = X[XINDEX(a,b)];
158 /*************************************************************************
162 *************************************************************************/
164 set *addToSet(node *v, set *X)
168 X = (set *) malloc(sizeof(set));
170 X->secondNode = NULL;
172 else if (NULL == X->firstNode)
175 X->secondNode = addToSet(v,X->secondNode);
179 node *makeNewNode(char *label, int i)
181 return(makeNode(label,NULL,i));
184 node *makeNode(char *label, edge *parentEdge, int index)
186 node *newNode; /*points to new node added to the graph*/
187 newNode = (node *) malloc(sizeof(node));
188 strncpy(newNode->label,label,NODE_LABEL_LENGTH);
189 newNode->index = index;
190 newNode->index2 = -1;
191 newNode->parentEdge = parentEdge;
192 newNode->leftEdge = NULL;
193 newNode->middleEdge = NULL;
194 newNode->rightEdge = NULL;
195 /*all fields have been initialized*/
199 /*copyNode returns a copy of v which has all of the fields identical to those
200 of v, except the node pointer fields*/
201 node *copyNode(node *v)
204 w = makeNode(v->label,NULL,v->index);
205 w->index2 = v->index2;
209 edge *siblingEdge(edge *e)
211 if(e == e->tail->leftEdge)
212 return(e->tail->rightEdge);
214 return(e->tail->leftEdge);
217 edge *makeEdge(char *label, node *tail, node *head, double weight)
220 newEdge = (edge *) malloc(sizeof(edge));
221 strncpy(newEdge->label,label,EDGE_LABEL_LENGTH);
222 newEdge->tail = tail;
223 newEdge->head = head;
224 newEdge->distance = weight;
225 newEdge->totalweight = 0.0;
232 T = (tree *) malloc(sizeof(tree));
239 void updateSizes(edge *e, int direction)
245 f = e->head->leftEdge;
248 f = e->head->rightEdge;
257 f = e->tail->parentEdge;
265 /*detrifurcate takes the (possibly trifurcated) input tree
266 and reroots the tree to a leaf*/
267 /*assumes tree is only trifurcated at root*/
268 tree *detrifurcate(tree *T)
275 if (NULL != v->parentEdge)
277 Rprintf ("Error: root %s is poorly rooted.\n",v->label);
280 for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f )
289 w->parentEdge = NULL;
295 void compareSets(tree *T, set *S)
300 e = depthFirstTraverse(T,NULL);
304 for(X = S; NULL != X; X = X->secondNode)
307 if (0 == strcmp(v->label,w->label))
309 v->index2 = w->index2;
314 e = depthFirstTraverse(T,e);
317 for(X = S; NULL != X; X = X->secondNode)
320 if (0 == strcmp(v->label,w->label))
322 v->index2 = w->index2;
329 Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
332 e = depthFirstTraverse(T,NULL);
336 if ((leaf(v)) && (-1 == v->index2))
338 Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
341 e = depthFirstTraverse(T,e);
343 for(X = S; NULL != X; X = X->secondNode)
344 if (X->firstNode->index2 > -1)
346 Rprintf("Error node %s in matrix but not a leaf in tree.\n",X->firstNode->label);
352 void partitionSizes(tree *T)
355 e = depthFirstTraverse(T,NULL);
361 e->bottomsize = e->head->leftEdge->bottomsize
362 + e->head->rightEdge->bottomsize;
363 e->topsize = (T->size + 2)/2 - e->bottomsize;
364 e = depthFirstTraverse(T,e);
368 /*************************************************************************
372 *************************************************************************/
374 edge *depthFirstTraverse(tree *T, edge *e)
375 /*depthFirstTraverse returns the edge f which is least in T according
376 to the depth-first order, but which is later than e in the search
377 pattern. If e is null, f is the least edge of T*/
382 f = T->root->leftEdge;
384 f = findBottomLeft(f);
385 return(f); /*this is the first edge of this search pattern*/
387 else /*e is non-null*/
389 if (e->tail->leftEdge == e)
390 /*if e is a left-oriented edge, we skip the entire
391 tree cut below e, and find least edge*/
393 else /*if e is a right-oriented edge, we have already looked at its
394 sibling and everything below e, so we move up*/
395 f = e->tail->parentEdge;
400 edge *findBottomLeft(edge *e)
401 /*findBottomLeft searches by gottom down in the tree and to the left.*/
405 while (NULL != f->head->leftEdge)
406 f = f->head->leftEdge;
410 edge *moveRight(edge *e)
413 f = e->tail->rightEdge; /*this step moves from a left-oriented edge
414 to a right-oriented edge*/
416 f = findBottomLeft(f);
420 edge *topFirstTraverse(tree *T, edge *e)
421 /*topFirstTraverse starts from the top of T, and from there moves stepwise
422 down, left before right*/
423 /*assumes tree has been detrifurcated*/
427 return(T->root->leftEdge); /*first Edge searched*/
428 else if (!(leaf(e->head)))
429 return(e->head->leftEdge); /*down and to the left is preferred*/
430 else /*e->head is a leaf*/
437 edge *moveUpRight(edge *e)
441 while ((NULL != f) && ( f->tail->leftEdge != f))
442 f = f->tail->parentEdge;
443 /*go up the tree until f is a leftEdge*/
445 return(f); /*triggered at end of search*/
447 return(f->tail->rightEdge);
448 /*and then go right*/
451 /*************************************************************************
455 *************************************************************************/
457 void freeMatrix(double **D, int size)
468 freeSet(S->secondNode);
472 void freeTree(tree *T)
476 if (NULL != v->leftEdge)
477 freeSubTree(v->leftEdge);
482 void freeSubTree(edge *e)