src/tree_build.c

   1 /* tree_build.c    2008-03-09 */
   2
   3 /* Copyright 2008 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 - 2; i >= 0; i--, k *= 10)
  16                 ans += ((int)x[i] - 48) * k;
  17
  18         return ans;
  19 }
  20
  21 void decode_edge(const char *x, int a, int b, int *node, double *w)
  22 {
  23         int i, k, co = a;
  24         char *endstr, str[100];
  25
  26         while (x[co] != ':') co++;
  27         if (a == co) *node = 0;
  28         else {
  29                 for (i = a, k = 0; i < co; i++, k++) str[k] = x[i];
  30                 str[k] = '\0';
  31                 *node = str2int(str, k + 1);
  32         }
  33         for (i = co + 1, k = 0; i <= b; i++, k++) str[k] = x[i];
  34         str[k] = '\0';
  35         *w = R_strtod(str, &endstr);
  36 }
  37
  38 #define ADD_INTERNAL_EDGE\
  39     e[j] = curnode;\
  40     e[j + nedge] = curnode = ++node;\
  41     j++
  42
  43 #define ADD_TERMINAL_EDGE\
  44     e[j] = curnode;\
  45     decode_edge(x, pr + 1, ps - 1, &tmpi, &tmpd);\
  46     e[j + nedge] = tmpi;\
  47     el[j] = tmpd;\
  48     j++
  49
  50 #define GO_DOWN\
  51     l = j - 1;\
  52     while (e[l + nedge] != curnode) l--;\
  53     decode_edge(x, ps + 1, pt - 1, &tmpi, &tmpd);\
  54     nl[curnode - ntip - 1] = tmpi;\
  55     el[l] = tmpd;\
  56     curnode = e[l]
  57
  58 SEXP treeBuildWithTokens(SEXP nwk)
  59 {
  60         const char *x;
  61         int n, i, ntip = 1, nnode = 0, nedge, *e, *nl, curnode, node, j, *skeleton, nsk = 0, ps, pr, pt, tmpi, l;
  62         double *el, tmpd;
  63         SEXP node_label, edge, edge_length, Nnode, phy;
  64
  65         PROTECT(nwk = coerceVector(nwk, STRSXP));
  66         x = CHAR(STRING_ELT(nwk, 0));
  67         n = strlen(x);
  68         skeleton = (int *)R_alloc(n, sizeof(int *));
  69         for (i = 0; i < n; i++) {
  70                 if (x[i] == '(' || x[i] == ',' || x[i] == ')') {
  71                         skeleton[nsk] = i;
  72                         nsk++;
  73                         switch(x[i]) {
  74                         case '(': break;
  75                         case ',': ntip++; break;
  76                         case ')': nnode++; break;
  77                         }
  78                 }
  79         }
  80         nedge = ntip + nnode - 1;
  81
  82         PROTECT(node_label = allocVector(INTSXP, nnode));
  83         PROTECT(Nnode = allocVector(INTSXP, 1));
  84         PROTECT(edge = allocVector(INTSXP, nedge*2));
  85         PROTECT(edge_length = allocVector(REALSXP, nedge));
  86         INTEGER(Nnode)[0] = nnode;
  87
  88         nl = INTEGER(node_label);
  89         memset(nl, 0, nnode*sizeof(int));
  90         e = INTEGER(edge);
  91         el = REAL(edge_length);
  92
  93         curnode = node = ntip + 1;
  94         j = 0;
  95
  96         for (i = 1; i < nsk - 1; i++) {
  97                 ps = skeleton[i];
  98                 Rprintf(""); /* <- again !!!! */
  99                 if (x[ps] == '(') {
 100                         ADD_INTERNAL_EDGE;
 101                         continue;
 102                 }
 103                 pr = skeleton[i - 1];
 104                 if (x[ps] == ',') {
 105                         if (x[pr] != ')') {
 106                                 /* !!! accolades indispensables !!! */
 107                                 ADD_TERMINAL_EDGE;
 108                         }
 109                         continue;
 110                 }
 111                 if (x[ps] == ')') {
 112                         pt = skeleton[i + 1];
 113                         if (x[pr] == ',') {
 114                                 ADD_TERMINAL_EDGE;
 115                                 GO_DOWN;
 116                                 continue;
 117                         }
 118                         if (x[pr] == ')') {
 119                                 GO_DOWN;
 120                         }
 121                 }
 122         }
 123
 124         pr = skeleton[nsk - 2];
 125         ps = skeleton[nsk - 1];
 126 /* is the last edge terminal? */
 127         if (x[pr] == ',' && x[ps] == ')') {
 128                 ADD_TERMINAL_EDGE;
 129         }
 130
 131 /* is there a root edge? */
 132         if (ps < n - 2) {
 133                 PROTECT(phy = allocVector(VECSXP, 5));
 134                 SEXP root_edge;
 135                 decode_edge(x, ps + 1, n - 2, &tmpi, &tmpd);
 136                 PROTECT(root_edge = allocVector(REALSXP, 1));
 137                 nl[0] = tmpi;
 138                 REAL(root_edge)[0] = tmpd;
 139                 SET_VECTOR_ELT(phy, 4, root_edge);
 140                 UNPROTECT(1);
 141         } else PROTECT(phy = allocVector(VECSXP, 4));
 142
 143         SET_VECTOR_ELT(phy, 0, edge);
 144         SET_VECTOR_ELT(phy, 1, edge_length);
 145         SET_VECTOR_ELT(phy, 2, Nnode);
 146         SET_VECTOR_ELT(phy, 3, node_label);
 147
 148         UNPROTECT(6);
 149         return phy;
 150 }
 151
 152 #undef ADD_INTERNAL_EDGE
 153 #undef ADD_TERMINAL_EDGE
 154 #undef GO_DOWN