X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=src%2Freorder_phylo.c;h=7a2f80eb0070d1f0739035d43199bfa6d5587852;hb=b0d1251527d8dd48ca1703b1cfaf217f413eda0e;hp=c54ea52ee9103d28eef2ca301903503dca327e04;hpb=5f07c2ca4b34dafc7ecc22e028c10d2d7eadffef;p=ape.git

diff --git a/src/reorder_phylo.c b/src/reorder_phylo.c
index c54ea52..7a2f80e 100644
--- a/src/reorder_phylo.c
+++ b/src/reorder_phylo.c
@@ -1,49 +1,85 @@
-/* reorder_phylo.c       2008-03-17 */
+/* reorder_phylo.c       2012-09-03 */
 
-/* Copyright 2008 Emmanuel Paradis */
+/* Copyright 2008-2012 Emmanuel Paradis */
 
 /* This file is part of the R-package `ape'. */
 /* See the file ../COPYING for licensing issues. */
 
 #include <R.h>
-#include <R_ext/Applic.h>
 
-void neworder_cladewise(int *n, int *edge1, int *edge2,
-			int *N, int *neworder)
-/* n: nb of tips, N: nb of edges */
+static int iii;
+
+void foo_reorder(int node, int n, int m, int *e1, int *e2, int *neworder, int *L, int *pos)
 {
-    int i, j, k, node, *done, dn, *node_back, eb;
-    done = &dn;
-    node_back = &eb;
-
-    /* done: indicates whether an edge has been collected
-       node_back: the series of node from the root to `node'
-       node: the current node */
-
-    done = (int*)R_alloc(*N, sizeof(int));
-    node_back = (int*)R_alloc(*N + 2 - *n, sizeof(int));
-    memset(done, 0, *N * sizeof(int));
-
-    j = k = 0;
-    node = *n + 1;
-    while (j < *N) {
-        for (i = 0; i < *N; i++) {
-  	    if (done[i] || edge1[i] != node) continue;
-	    neworder[j] = i + 1;
-	    j++;
-	    done[i] = 1;
-	    if (edge2[i] > *n) {
-	        node_back[k] = node;
-		k++;
-		node = edge2[i];
-		/* if found a new node, reset the loop */
-		i = -1;
-	    }
+	int i = node - n - 1, j, k;
+
+/* 'i' is the C index corresponding to 'node' */
+
+	for (j = 0; j < pos[i]; j++) {
+		k = L[i + m * j];
+		neworder[iii++] = k + 1;
+		if (e2[k] > n) /* is it an internal edge? */
+			foo_reorder(e2[k], n, m, e1, e2, neworder, L, pos);
+	}
+}
+
+void bar_reorder(int node, int n, int m, int *e1, int *e2, int *neworder, int *L, int *pos)
+{
+	int i = node - n - 1, j, k;
+
+	for (j = pos[i] - 1; j >= 0; j--)
+		neworder[iii--] = L[i + m * j] + 1;
+
+	for (j = 0; j < pos[i]; j++) {
+		k = e2[L[i + m * j]];
+		if (k > n)
+			bar_reorder(k, n, m, e1, e2, neworder, L, pos);
+	}
+}
+
+void neworder_phylo(int *n, int *e1, int *e2, int *N, int *neworder, int *order)
+/* n: nb of tips
+   m: nb of nodes
+   N: nb of edges */
+{
+	int i, j, k, *L, *pos, m = *N - *n + 1, degrmax = *n - m + 1;
+
+/* degrmax is the largest value that a node degree can be */
+
+/* L is a 1-d array storing, for each node, the C indices of the rows of
+   the edge matrix where the node is ancestor (i.e., present in the 1st
+   column). It is used in the same way than a matrix (which is actually
+   a vector) is used in R as a 2-d structure. */
+
+	L = (int*)R_alloc(m * degrmax, sizeof(int));
+
+/* pos gives the position for each 'row' of L, that is the number of elements
+   which have already been stored for that 'row'. */
+
+	pos = (int*)R_alloc(m, sizeof(int));
+	memset(pos, 0, m * sizeof(int));
+
+/* we now go down along the edge matrix */
+
+	for (i = 0; i < *N; i++) {
+		k = e1[i] - *n - 1; /* k is the 'row' index in L corresponding to node e1[i] */
+		j = pos[k]; /* the current 'column' position corresping to k */
+		pos[k]++; /* increment in case the same node is found in another row of the edge matrix */
+		L[k + m * j] = i;
+	}
+
+/* L is now ready: we can start the recursive calls. */
+/* We start with the root 'n + 1': its index will be changed into
+   the corresponding C index inside the recursive function. */
+
+	switch(*order) {
+	case 1 : iii = 0;
+		foo_reorder(*n + 1, *n, m, e1, e2, neworder, L, pos);
+		break;
+	case 2 : iii = *N - 1;
+		bar_reorder(*n + 1, *n, m, e1, e2, neworder, L, pos);
+		break;
 	}
-	/* if arrived at the end of `edge', go down one node */
-	k--;
-	node = node_back[k];
-    }
 }
 
 #define DO_NODE_PRUNING\
@@ -61,23 +97,19 @@ void neworder_cladewise(int *n, int *edge1, int *edge2,
 void neworder_pruningwise(int *ntip, int *nnode, int *edge1,
 			  int *edge2, int *nedge, int *neworder)
 {
-    int *Ndegr, degree, *ready, rdy, i, j, node, nextI, n;
-    Ndegr = &degree;
-    ready = &rdy;
-
-    ready = (int*)R_alloc(*nedge, sizeof(int));
+    int *ready, *Ndegr, i, j, node, nextI, n;
 
-    /* use `nextI' temporarily because need an address for R_tabulate */
     nextI = *ntip +  *nnode;
     Ndegr = (int*)R_alloc(nextI, sizeof(int));
     memset(Ndegr, 0, nextI*sizeof(int));
-    R_tabulate(edge1, nedge, &nextI, Ndegr);
+    for (i = 0; i < *nedge; i++) (Ndegr[edge1[i] - 1])++;
+
+    ready = (int*)R_alloc(*nedge, sizeof(int));
 
     /* `ready' indicates whether an edge is ready to be */
     /* collected; only the terminal edges are initially ready */
     for (i = 0; i < *nedge; i++)
-      if (edge2[i] <= *ntip) ready[i] = 1;
-      else ready[i] = 0;
+	    ready[i] = (edge2[i] <= *ntip) ? 1 : 0;
 
     /* `n' counts the number of times a node has been seen. */
     /* This algo will work if the tree is in cladewise order, */
@@ -117,3 +149,4 @@ void neworder_pruningwise(int *ntip, int *nnode, int *edge1,
 	nextI++;
     }
 }
+