src/tree_build.c

   1 /* tree_build.c    2011-06-23 */
   2
   3 /* Copyright 2008-2011 Emmanuel Paradis */
   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 static int str2int(char *x, int n)
  12 {
  13         int i, k = 1, ans = 0;
  14
  15         for (i = n - 1; i >= 0; i--, k *= 10)
  16                 ans += ((int)x[i] - 48) * k;
  17
  18         return ans;
  19 }
  20
  21 void extract_portion_Newick(const char *x, int a, int b, char *y)
  22 {
  23         int i, j;
  24
  25         for (i = a, j = 0; i <= b; i++, j++) y[j] = x[i];
  26
  27         y[j] = '\0';
  28 }
  29
  30 void decode_terminal_edge_token(const char *x, int a, int b, int *node, double *w)
  31 {
  32         int co = a;
  33         char *endstr, str[100];
  34
  35         while (x[co] != ':') co++;
  36
  37         extract_portion_Newick(x, a, co - 1, str);
  38         *node = str2int(str, co - a);
  39         extract_portion_Newick(x, co + 1, b, str);
  40         *w = R_strtod(str, &endstr);
  41 }
  42
  43 void decode_internal_edge(const char *x, int a, int b, char *lab, double *w)
  44 {
  45         int co = a;
  46         char *endstr, str[100];
  47
  48         while (x[co] != ':') co++;
  49
  50         if (a == co) lab[0] = '\0'; /* if no node label */
  51         else extract_portion_Newick(x, a, co - 1, lab);
  52
  53         extract_portion_Newick(x, co + 1, b, str);
  54         *w = R_strtod(str, &endstr);
  55 }
  56
  57 #define ADD_INTERNAL_EDGE            \
  58     e[j] = curnode;                  \
  59     e[j + nedge] = curnode = ++node; \
  60     stack_internal[k++] = j;         \
  61     j++
  62
  63 #define ADD_TERMINAL_EDGE                                        \
  64     e[j] = curnode;                                              \
  65     decode_terminal_edge_token(x, pr + 1, ps - 1, &tmpi, &tmpd); \
  66     e[j + nedge] = tmpi;                                         \
  67     el[j] = tmpd;                                                \
  68     j++
  69
  70 #define GO_DOWN                                                  \
  71     decode_internal_edge(x, ps + 1, pt - 1, lab, &tmpd);         \
  72     SET_STRING_ELT(node_label, curnode - 1 - ntip, mkChar(lab)); \
  73     l = stack_internal[--k];                                     \
  74     el[l] = tmpd;                                                \
  75     curnode = e[l]
  76
  77 SEXP treeBuildWithTokens(SEXP nwk)
  78 {
  79         const char *x;
  80         int n, i, ntip = 1, nnode = 0, nedge, *e, curnode, node, j, *skeleton, nsk = 0, ps, pr, pt, tmpi, l, k, stack_internal[10000];
  81         double *el, tmpd;
  82         char lab[512];
  83         SEXP edge, edge_length, Nnode, node_label, phy;
  84
  85         PROTECT(nwk = coerceVector(nwk, STRSXP));
  86         x = CHAR(STRING_ELT(nwk, 0));
  87         n = strlen(x);
  88         skeleton = (int *)R_alloc(n, sizeof(int *));
  89
  90 /* first pass on the Newick string to localize parentheses and commas */
  91         for (i = 0; i < n; i++) {
  92                 if (x[i] == '(') {
  93                         skeleton[nsk] = i;
  94                         nsk++;
  95                         continue;
  96                 }
  97                 if (x[i] == ',') {
  98                         skeleton[nsk] = i;
  99                         nsk++;
 100                         ntip++;
 101                         continue;
 102                 }
 103                 if (x[i] == ')') {
 104                         skeleton[nsk] = i;
 105                         nsk++;
 106                         nnode++;
 107                 }
 108         }
 109
 110         nedge = ntip + nnode - 1;
 111
 112         PROTECT(Nnode = allocVector(INTSXP, 1));
 113         PROTECT(edge = allocVector(INTSXP, nedge*2));
 114         PROTECT(edge_length = allocVector(REALSXP, nedge));
 115         PROTECT(node_label = allocVector(STRSXP, nnode));
 116         INTEGER(Nnode)[0] = nnode;
 117
 118         e = INTEGER(edge);
 119         el = REAL(edge_length);
 120
 121         curnode = node = ntip + 1;
 122         k = j = 0;
 123 /* j: index of the current position in the edge matrix */
 124 /* k: index of the current position in stack_internal */
 125 /* stack_internal is a simple array storing the indices of the
 126    successive internal edges from the root; it's a stack so it is
 127    incremented every time an internal edge is added, and decremented
 128    every GO_DOWN step. This makes easy to the index of the
 129    subtending edge. */
 130
 131 /* second pass on the Newick string to build the "phylo" object elements */
 132         for (i = 1; i < nsk - 1; i++) {
 133                 ps = skeleton[i];
 134 //              Rprintf(""); /* <- again !!!! */
 135                 if (x[ps] == '(') {
 136                         ADD_INTERNAL_EDGE;
 137                         continue;
 138                 }
 139                 pr = skeleton[i - 1];
 140                 if (x[ps] == ',') {
 141                         if (x[pr] != ')') {
 142                                 /* !!! accolades indispensables !!! */
 143                                 ADD_TERMINAL_EDGE;
 144                         }
 145                         continue;
 146                 }
 147                 if (x[ps] == ')') {
 148                         pt = skeleton[i + 1]; // <- utile ???
 149                         if (x[pr] == ',') {
 150                                 ADD_TERMINAL_EDGE;
 151                                 GO_DOWN;
 152                                 continue;
 153                         }
 154                         if (x[pr] == ')') {
 155                                 GO_DOWN;
 156                         }
 157                 }
 158         }
 159
 160         pr = skeleton[nsk - 2];
 161         ps = skeleton[nsk - 1];
 162 /* is the last edge terminal? */
 163         if (x[pr] == ',' && x[ps] == ')') {
 164                 ADD_TERMINAL_EDGE;
 165         }
 166
 167 /* is there a root edge and/or root label? */
 168         if (ps < n - 2) {
 169                 i = ps + 1;
 170                 while (i < n - 2 && x[i] != ':') i++;
 171                 if (i < n - 2) {
 172                         PROTECT(phy = allocVector(VECSXP, 5));
 173                         SEXP root_edge;
 174                         decode_internal_edge(x, ps + 1, n - 2, lab, &tmpd);
 175                         PROTECT(root_edge = allocVector(REALSXP, 1));
 176                         REAL(root_edge)[0] = tmpd;
 177                         SET_VECTOR_ELT(phy, 4, root_edge);
 178                         UNPROTECT(1);
 179                         SET_STRING_ELT(node_label, 0, mkChar(lab));
 180                 } else {
 181                         extract_portion_Newick(x, ps + 1, n - 2, lab);
 182                         SET_STRING_ELT(node_label, 0, mkChar(lab));
 183                         PROTECT(phy = allocVector(VECSXP, 4));
 184                 }
 185         } else PROTECT(phy = allocVector(VECSXP, 4));
 186
 187         SET_VECTOR_ELT(phy, 0, edge);
 188         SET_VECTOR_ELT(phy, 1, edge_length);
 189         SET_VECTOR_ELT(phy, 2, Nnode);
 190         SET_VECTOR_ELT(phy, 3, node_label);
 191
 192         UNPROTECT(6);
 193         return phy;
 194 }
 195
 196 #undef ADD_INTERNAL_EDGE
 197 #undef ADD_TERMINAL_EDGE
 198 #undef GO_DOWN