1 /* bipartition.c 2007-06-29 */
3 /* Copyright 2005-2007 Emmanuel Paradis, and 2007 R Development Core Team */
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
9 #include <Rinternals.h>
11 SEXP seq_root2tip(SEXP edge, SEXP nbtip, SEXP nbnode)
13 int i, j, k, Nedge, *x, *done, dn, sumdone, lt, ROOT, Ntip, Nnode;
14 SEXP ans, seqnod, tmp_vec;
16 /* The following is needed only if we are not sure
17 that the storage mode of `edge' is "integer". */
18 PROTECT(edge = coerceVector(edge, INTSXP));
19 PROTECT(nbtip = coerceVector(nbtip, INTSXP));
20 PROTECT(nbnode = coerceVector(nbnode, INTSXP));
21 x = INTEGER(edge); /* copy the pointer */
22 Ntip = *INTEGER(nbtip);
23 Nnode = *INTEGER(nbnode);
24 Nedge = LENGTH(edge)/2;
27 PROTECT(ans = allocVector(VECSXP, Ntip));
28 PROTECT(seqnod = allocVector(VECSXP, Nnode));
31 done = (int*)R_alloc(Nnode, sizeof(int));
32 for (i = 0; i < Nnode; i++) done[i] = 0;
34 tmp_vec = allocVector(INTSXP, 1);
35 INTEGER(tmp_vec)[0] = ROOT; /* sure ? */
36 SET_VECTOR_ELT(seqnod, 0, tmp_vec);
39 while (sumdone < Nnode) {
40 for (i = 0; i < Nnode; i++) { /* loop through all nodes */
41 /* if the vector is not empty and its */
42 /* descendants are not yet found */
43 if (VECTOR_ELT(seqnod, i) == R_NilValue || done[i]) continue;
44 /* look for the descendants in 'edge': */
45 for (j = 0; j < Nedge; j++) {
46 /* skip the terminal edges, we look only for nodes now */
47 if (x[j] - Ntip != i + 1 || x[j + Nedge] <= Ntip) continue;
48 /* can now make the sequence from */
49 /* the root to the current node */
50 lt = LENGTH(VECTOR_ELT(seqnod, i));
51 tmp_vec = allocVector(INTSXP, lt + 1);
52 for (k = 0; k < lt; k++)
53 INTEGER(tmp_vec)[k] = INTEGER(VECTOR_ELT(seqnod, i))[k];
54 INTEGER(tmp_vec)[lt] = x[j + Nedge];
55 SET_VECTOR_ELT(seqnod, x[j + Nedge] - Ntip - 1, tmp_vec);
62 /* build the sequence from root to tip */
63 /* by simply looping through 'edge' */
64 for (i = 0; i < Nedge; i++) {
65 /* skip the internal edges */
66 if (x[i + Nedge] > Ntip) continue;
67 lt = LENGTH(VECTOR_ELT(seqnod, x[i] - Ntip - 1));
68 tmp_vec = allocVector(INTSXP, lt + 1);
69 for (j = 0; j < lt; j++)
70 INTEGER(tmp_vec)[j] = INTEGER(VECTOR_ELT(seqnod, x[i] - Ntip - 1))[j];
71 INTEGER(tmp_vec)[lt] = x[i + Nedge];
72 SET_VECTOR_ELT(ans, x[i + Nedge] - 1, tmp_vec);
77 } /* EOF seq_root2tip */
79 SEXP bipartition(SEXP edge, SEXP nbtip, SEXP nbnode)
81 int i, j, k, lt, lt2, inod, Ntip, Nnode;
82 SEXP ans, seqnod, tmp_vec;
84 PROTECT(edge = coerceVector(edge, INTSXP));
85 PROTECT(nbtip = coerceVector(nbtip, INTSXP));
86 PROTECT(nbnode = coerceVector(nbnode, INTSXP));
87 Ntip = *INTEGER(nbtip);
88 Nnode = *INTEGER(nbnode);
90 PROTECT(ans = allocVector(VECSXP, Nnode));
91 PROTECT(seqnod = seq_root2tip(edge, nbtip, nbnode));
93 for (i = 0; i < LENGTH(seqnod); i++) { /* for each tip */
94 lt = LENGTH(VECTOR_ELT(seqnod, i));
95 for (j = 0; j < lt - 1; j++) {
96 inod = INTEGER(VECTOR_ELT(seqnod, i))[j] - Ntip - 1;
97 if (VECTOR_ELT(ans, inod) == R_NilValue) {
98 tmp_vec = allocVector(INTSXP, 1);
99 INTEGER(tmp_vec)[0] = i + 1;
101 lt2 = LENGTH(VECTOR_ELT(ans, inod));
102 tmp_vec = allocVector(INTSXP, lt2 + 1);
103 for (k = 0; k < lt2; k++)
104 INTEGER(tmp_vec)[k] = INTEGER(VECTOR_ELT(ans, inod))[k];
105 INTEGER(tmp_vec)[lt2] = i + 1;
107 SET_VECTOR_ELT(ans, inod, tmp_vec);
116 (not the same than in library/stats/src/nls.c) */
117 SEXP getListElement(SEXP list, char *str)
119 SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);
122 for (i = 0; i < length(list); i++)
123 if(strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {
124 elmt = VECTOR_ELT(list, i);
130 int SameClade(SEXP clade1, SEXP clade2)
132 int i, n = LENGTH(clade1), *c1, *c2;
134 if (n != LENGTH(clade2)) return 0;
136 c1 = INTEGER(clade1);
137 c2 = INTEGER(clade2);
138 for (i = 0; i < n; i++)
139 if (c1[i] != c2[i]) return 0;
144 SEXP prop_part(SEXP TREES, SEXP nbtree, SEXP keep_partitions)
146 int i, j, k, l, KeepPartition, Ntree, Ntip, Nnode, Npart, NpartCurrent, *no;
147 SEXP bp, ans, nbtip, nbnode, number;
149 PROTECT(nbtree = coerceVector(nbtree, INTSXP));
150 PROTECT(keep_partitions = coerceVector(keep_partitions, INTSXP));
151 Ntree = *INTEGER(nbtree);
152 KeepPartition = *INTEGER(keep_partitions);
155 Ntip = LENGTH(getListElement(VECTOR_ELT(TREES, 0), "tip.label"));
156 Nnode = *INTEGER(getListElement(VECTOR_ELT(TREES, 0), "Nnode"));
158 PROTECT(nbtip = allocVector(INTSXP, 1));
159 PROTECT(nbnode = allocVector(INTSXP, 1));
160 INTEGER(nbtip)[0] = Ntip;
161 INTEGER(nbnode)[0] = Nnode;
163 if (KeepPartition) Npart = Ntree*(Nnode - 1) + 1;
166 PROTECT(number = allocVector(INTSXP, Npart));
167 no = INTEGER(number); /* copy the pointer */
168 /* The first partition in the returned list has all tips,
169 so it is observed in all trees: */
171 /* The partitions in the first tree are obviously observed once: */
172 for (i = 1; i < Nnode; i++) no[i] = 1;
175 for (i = Nnode; i < Npart; i++) no[i] = 0;
177 PROTECT(ans = allocVector(VECSXP, Npart));
178 PROTECT(bp = bipartition(getListElement(VECTOR_ELT(TREES, 0), "edge"),
180 for (i = 0; i < Nnode; i++)
181 SET_VECTOR_ELT(ans, i, VECTOR_ELT(bp, i));
184 PROTECT(ans = bipartition(getListElement(VECTOR_ELT(TREES, 0), "edge"),
188 NpartCurrent = Nnode;
190 /* We start on the 2nd tree: */
191 for (k = 1; k < Ntree; k++) {
192 PROTECT(bp = bipartition(getListElement(VECTOR_ELT(TREES, k), "edge"),
194 for (i = 1; i < Nnode; i++) {
197 if (SameClade(VECTOR_ELT(bp, i), VECTOR_ELT(ans, j))) {
202 if (j < NpartCurrent) goto next_j;
205 SET_VECTOR_ELT(ans, NpartCurrent, VECTOR_ELT(bp, i));
212 if (KeepPartition && NpartCurrent < Npart) {
213 PROTECT(bp = allocVector(VECSXP, NpartCurrent));
214 for (i = 0; i < NpartCurrent; i++)
215 SET_VECTOR_ELT(bp, i, VECTOR_ELT(ans, i));
216 setAttrib(bp, install("number"), number);
220 setAttrib(ans, install("number"), number);