3 /* Copyright 2007-2008 Olivier Gascuel, Rick Desper,
4 R port by Vincent Lefort and 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, int *labels,
30 int *nni, int *spr, int *tbr,
31 int *edge1, int *edge2, double *el)
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, labels, n);
63 freeMatrix(A,2*n - 2);
69 void me_o(double *X, int *N, int *labels, int *nni,
70 int *edge1, int *edge2, double *el)
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, labels, n);
102 freeMatrix(A,2*n - 2);
110 -- MATRIX FUNCTIONS --
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)
129 double **loadMatrix (double *X, int *labels, int n, set *S)
131 // char nextString[MAX_LABEL_LENGTH];
136 table = (double **) calloc(n,sizeof(double *));
138 table[i] = (double *) calloc(n,sizeof(double));
142 // strncpy (nextString, labels[i], MAX_LABEL_LENGTH);
143 // ReplaceForbiddenChars (nextString, '_');
144 // v = makeNewNode(nextString,-1);
145 v = makeNewNode(labels[i], -1);
148 for (j=i; j<n; j++) {
151 table[j][i] = X[XINDEX(a,b)];
152 table[i][j] = X[XINDEX(a,b)];
162 -- GRAPH FUNCTIONS --
166 set *addToSet(node *v, set *X)
170 X = (set *) malloc(sizeof(set));
172 X->secondNode = NULL;
174 else if (NULL == X->firstNode)
177 X->secondNode = addToSet(v,X->secondNode);
181 //node *makeNewNode(char *label, int i)
182 node *makeNewNode(int label, int i)
184 return(makeNode(label,NULL,i));
187 //node *makeNode(char *label, edge *parentEdge, int index)
188 node *makeNode(int label, edge *parentEdge, int index)
190 node *newNode; /*points to new node added to the graph*/
191 newNode = (node *) malloc(sizeof(node));
192 // strncpy(newNode->label,label,NODE_LABEL_LENGTH);
193 newNode->label = label;
194 newNode->index = index;
195 newNode->index2 = -1;
196 newNode->parentEdge = parentEdge;
197 newNode->leftEdge = NULL;
198 newNode->middleEdge = NULL;
199 newNode->rightEdge = NULL;
200 /*all fields have been initialized*/
204 /*copyNode returns a copy of v which has all of the fields identical to those
205 of v, except the node pointer fields*/
206 node *copyNode(node *v)
209 w = makeNode(v->label,NULL,v->index);
210 w->index2 = v->index2;
214 edge *siblingEdge(edge *e)
216 if(e == e->tail->leftEdge)
217 return(e->tail->rightEdge);
219 return(e->tail->leftEdge);
222 edge *makeEdge(char *label, node *tail, node *head, double weight)
225 newEdge = (edge *) malloc(sizeof(edge));
226 strncpy(newEdge->label,label,EDGE_LABEL_LENGTH);
227 newEdge->tail = tail;
228 newEdge->head = head;
229 newEdge->distance = weight;
230 newEdge->totalweight = 0.0;
237 T = (tree *) malloc(sizeof(tree));
244 void updateSizes(edge *e, int direction)
250 f = e->head->leftEdge;
253 f = e->head->rightEdge;
262 f = e->tail->parentEdge;
270 /*detrifurcate takes the (possibly trifurcated) input tree
271 and reroots the tree to a leaf*/
272 /*assumes tree is only trifurcated at root*/
273 tree *detrifurcate(tree *T)
280 if (NULL != v->parentEdge)
282 error("root %d is poorly rooted.", v->label);
284 for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f )
293 w->parentEdge = NULL;
299 void compareSets(tree *T, set *S)
304 e = depthFirstTraverse(T,NULL);
308 for(X = S; NULL != X; X = X->secondNode)
311 // if (0 == strcmp(v->label,w->label))
312 if (v->label == w->label)
314 v->index2 = w->index2;
319 e = depthFirstTraverse(T,e);
322 for(X = S; NULL != X; X = X->secondNode)
325 // if (0 == strcmp(v->label,w->label))
326 if (v->label == w->label)
328 v->index2 = w->index2;
335 error("leaf %d in tree not in distance matrix.", v->label);
337 e = depthFirstTraverse(T,NULL);
341 if ((leaf(v)) && (-1 == v->index2))
343 error("leaf %d in tree not in distance matrix.", v->label);
345 e = depthFirstTraverse(T,e);
347 for(X = S; NULL != X; X = X->secondNode)
348 if (X->firstNode->index2 > -1)
350 error("node %d in matrix but not a leaf in tree.", X->firstNode->label);
355 void partitionSizes(tree *T)
358 e = depthFirstTraverse(T,NULL);
364 e->bottomsize = e->head->leftEdge->bottomsize
365 + e->head->rightEdge->bottomsize;
366 e->topsize = (T->size + 2)/2 - e->bottomsize;
367 e = depthFirstTraverse(T,e);
371 /*************************************************************************
375 *************************************************************************/
377 edge *depthFirstTraverse(tree *T, edge *e)
378 /*depthFirstTraverse returns the edge f which is least in T according
379 to the depth-first order, but which is later than e in the search
380 pattern. If e is null, f is the least edge of T*/
385 f = T->root->leftEdge;
387 f = findBottomLeft(f);
388 return(f); /*this is the first edge of this search pattern*/
390 else /*e is non-null*/
392 if (e->tail->leftEdge == e)
393 /*if e is a left-oriented edge, we skip the entire
394 tree cut below e, and find least edge*/
396 else /*if e is a right-oriented edge, we have already looked at its
397 sibling and everything below e, so we move up*/
398 f = e->tail->parentEdge;
403 edge *findBottomLeft(edge *e)
404 /*findBottomLeft searches by gottom down in the tree and to the left.*/
408 while (NULL != f->head->leftEdge)
409 f = f->head->leftEdge;
413 edge *moveRight(edge *e)
416 f = e->tail->rightEdge; /*this step moves from a left-oriented edge
417 to a right-oriented edge*/
419 f = findBottomLeft(f);
423 edge *topFirstTraverse(tree *T, edge *e)
424 /*topFirstTraverse starts from the top of T, and from there moves stepwise
425 down, left before right*/
426 /*assumes tree has been detrifurcated*/
430 return(T->root->leftEdge); /*first Edge searched*/
431 else if (!(leaf(e->head)))
432 return(e->head->leftEdge); /*down and to the left is preferred*/
433 else /*e->head is a leaf*/
440 edge *moveUpRight(edge *e)
444 while ((NULL != f) && ( f->tail->leftEdge != f))
445 f = f->tail->parentEdge;
446 /*go up the tree until f is a leftEdge*/
448 return(f); /*triggered at end of search*/
450 return(f->tail->rightEdge);
451 /*and then go right*/
454 /*************************************************************************
458 *************************************************************************/
460 void freeMatrix(double **D, int size)
471 freeSet(S->secondNode);
475 void freeTree(tree *T)
479 if (NULL != v->leftEdge)
480 freeSubTree(v->leftEdge);
485 void freeSubTree(edge *e)