src/bipartition.c

   1 /* bipartition.c    2011-06-23 */
   2
   3 /* Copyright 2005-2011 Emmanuel Paradis, and 2007 R Development Core Team */
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 #include <R.h>
   9 #include <Rinternals.h>
  10
  11 SEXP seq_root2tip(SEXP edge, SEXP nbtip, SEXP nbnode)
  12 {
  13     int i, j, k, Nedge, *x, *done, dn, sumdone, lt, ROOT, Ntip, Nnode;
  14     SEXP ans, seqnod, tmp_vec;
  15
  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;
  25     ROOT = Ntip + 1;
  26
  27     PROTECT(ans = allocVector(VECSXP, Ntip));
  28     PROTECT(seqnod = allocVector(VECSXP, Nnode));
  29
  30     done = &dn;
  31     done = (int*)R_alloc(Nnode, sizeof(int));
  32     for (i = 0; i < Nnode; i++) done[i] = 0;
  33
  34     tmp_vec = allocVector(INTSXP, 1);
  35     INTEGER(tmp_vec)[0] = ROOT; /* sure ? */
  36     SET_VECTOR_ELT(seqnod, 0, tmp_vec);
  37     sumdone = 0;
  38
  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);
  56             }
  57             done[i] = 1;
  58             sumdone++;
  59         }
  60     }
  61
  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);
  73     }
  74
  75     UNPROTECT(5);
  76     return ans;
  77 } /* EOF seq_root2tip */
  78
  79 SEXP bipartition(SEXP edge, SEXP nbtip, SEXP nbnode)
  80 {
  81     int i, j, k, lt, lt2, inod, Ntip, Nnode;
  82     SEXP ans, seqnod, tmp_vec;
  83
  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);
  89
  90     PROTECT(ans = allocVector(VECSXP, Nnode));
  91     PROTECT(seqnod = seq_root2tip(edge, nbtip, nbnode));
  92
  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;
 100             } else {
 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;
 106             }
 107             SET_VECTOR_ELT(ans, inod, tmp_vec);
 108         }
 109     }
 110
 111     UNPROTECT(5);
 112     return ans;
 113 } /* bipartition */
 114
 115 /* From R-ext manual
 116    (not the same than in library/stats/src/nls.c) */
 117 SEXP getListElement(SEXP list, char *str)
 118 {
 119     SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);
 120     int i;
 121
 122     for (i = 0; i < length(list); i++)
 123       if(strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {
 124           elmt = VECTOR_ELT(list, i);
 125           break;
 126       }
 127     return elmt;
 128 }
 129
 130 int SameClade(SEXP clade1, SEXP clade2)
 131 {
 132     int i, n = LENGTH(clade1), *c1, *c2;
 133
 134     if (n != LENGTH(clade2)) return 0;
 135
 136     c1 = INTEGER(clade1);
 137     c2 = INTEGER(clade2);
 138     for (i = 0; i < n; i++)
 139       if (c1[i] != c2[i]) return 0;
 140
 141     return 1;
 142 }
 143
 144 SEXP prop_part(SEXP TREES, SEXP nbtree, SEXP keep_partitions)
 145 {
 146     int i, j, k, KeepPartition, Ntree, Ntip, Nnode, Npart, NpartCurrent, *no;
 147     SEXP bp, ans, nbtip, nbnode, number;
 148
 149     PROTECT(nbtree = coerceVector(nbtree, INTSXP));
 150     PROTECT(keep_partitions = coerceVector(keep_partitions, INTSXP));
 151     Ntree = *INTEGER(nbtree);
 152     KeepPartition = *INTEGER(keep_partitions);
 153
 154
 155     Ntip = LENGTH(getListElement(VECTOR_ELT(TREES, 0), "tip.label"));
 156     Nnode = *INTEGER(getListElement(VECTOR_ELT(TREES, 0), "Nnode"));
 157
 158     PROTECT(nbtip = allocVector(INTSXP, 1));
 159     PROTECT(nbnode = allocVector(INTSXP, 1));
 160     INTEGER(nbtip)[0] = Ntip;
 161     INTEGER(nbnode)[0] = Nnode;
 162
 163     if (KeepPartition) Npart = Ntree*(Nnode - 1) + 1;
 164     else Npart = Nnode;
 165
 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: */
 170     no[0] = Ntree;
 171     /* The partitions in the first tree are obviously observed once: */
 172     for (i = 1; i < Nnode; i++) no[i] = 1;
 173
 174     if (KeepPartition) {
 175         for (i = Nnode; i < Npart; i++) no[i] = 0;
 176
 177         PROTECT(ans = allocVector(VECSXP, Npart));
 178         PROTECT(bp = bipartition(getListElement(VECTOR_ELT(TREES, 0), "edge"),
 179                                  nbtip, nbnode));
 180         for (i = 0; i < Nnode; i++)
 181           SET_VECTOR_ELT(ans, i, VECTOR_ELT(bp, i));
 182         UNPROTECT(1);
 183     } else {
 184         PROTECT(ans = bipartition(getListElement(VECTOR_ELT(TREES, 0), "edge"),
 185                                   nbtip, nbnode));
 186     }
 187
 188     NpartCurrent = Nnode;
 189
 190     /* We start on the 2nd tree: */
 191     for (k = 1; k < Ntree; k++) {
 192         PROTECT(bp = bipartition(getListElement(VECTOR_ELT(TREES, k), "edge"),
 193                                  nbtip, nbnode));
 194         for (i = 1; i < Nnode; i++) {
 195             j = 1;
 196 next_j:
 197             if (SameClade(VECTOR_ELT(bp, i), VECTOR_ELT(ans, j))) {
 198                 no[j]++;
 199                 continue;
 200             }
 201             j++;
 202             if (j < NpartCurrent) goto next_j;
 203             if (KeepPartition) {
 204                 no[NpartCurrent]++;
 205                 SET_VECTOR_ELT(ans, NpartCurrent, VECTOR_ELT(bp, i));
 206                 NpartCurrent++;
 207             }
 208         }
 209         UNPROTECT(1);
 210     }
 211
 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);
 217         UNPROTECT(7);
 218         return bp;
 219     } else {
 220         setAttrib(ans, install("number"), number);
 221         UNPROTECT(6);
 222         return ans;
 223     }
 224 } /* prop_part */