src/tree_build.c

   1 /* tree_build.c    2009-11-21 */
   2
   3 /* Copyright 2008-2009 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     el[l] = tmpd;\
  55     curnode = e[l]
  56
  57 SEXP treeBuildWithTokens(SEXP nwk)
  58 {
  59         const char *x;
  60         int n, i, ntip = 1, nnode = 0, nedge, *e, curnode, node, j, *skeleton, nsk = 0, ps, pr, pt, tmpi, l;
  61         double *el, tmpd;
  62         SEXP edge, edge_length, Nnode, phy;
  63
  64         PROTECT(nwk = coerceVector(nwk, STRSXP));
  65         x = CHAR(STRING_ELT(nwk, 0));
  66         n = strlen(x);
  67         skeleton = (int *)R_alloc(n, sizeof(int *));
  68         for (i = 0; i < n; i++) {
  69                 if (x[i] == '(' || x[i] == ',' || x[i] == ')') {
  70                         skeleton[nsk] = i;
  71                         nsk++;
  72                         switch(x[i]) {
  73                         case '(': break;
  74                         case ',': ntip++; break;
  75                         case ')': nnode++; break;
  76                         }
  77                 }
  78         }
  79         nedge = ntip + nnode - 1;
  80
  81         PROTECT(Nnode = allocVector(INTSXP, 1));
  82         PROTECT(edge = allocVector(INTSXP, nedge*2));
  83         PROTECT(edge_length = allocVector(REALSXP, nedge));
  84         INTEGER(Nnode)[0] = nnode;
  85
  86         e = INTEGER(edge);
  87         el = REAL(edge_length);
  88
  89         curnode = node = ntip + 1;
  90         j = 0;
  91
  92         for (i = 1; i < nsk - 1; i++) {
  93                 ps = skeleton[i];
  94                 Rprintf(""); /* <- again !!!! */
  95                 if (x[ps] == '(') {
  96                         ADD_INTERNAL_EDGE;
  97                         continue;
  98                 }
  99                 pr = skeleton[i - 1];
 100                 if (x[ps] == ',') {
 101                         if (x[pr] != ')') {
 102                                 /* !!! accolades indispensables !!! */
 103                                 ADD_TERMINAL_EDGE;
 104                         }
 105                         continue;
 106                 }
 107                 if (x[ps] == ')') {
 108                         pt = skeleton[i + 1];
 109                         if (x[pr] == ',') {
 110                                 ADD_TERMINAL_EDGE;
 111                                 GO_DOWN;
 112                                 continue;
 113                         }
 114                         if (x[pr] == ')') {
 115                                 GO_DOWN;
 116                         }
 117                 }
 118         }
 119
 120         pr = skeleton[nsk - 2];
 121         ps = skeleton[nsk - 1];
 122 /* is the last edge terminal? */
 123         if (x[pr] == ',' && x[ps] == ')') {
 124                 ADD_TERMINAL_EDGE;
 125         }
 126
 127 /* is there a root edge? */
 128         if (ps < n - 2) {
 129                 PROTECT(phy = allocVector(VECSXP, 4));
 130                 SEXP root_edge;
 131                 decode_edge(x, ps + 1, n - 2, &tmpi, &tmpd);
 132                 PROTECT(root_edge = allocVector(REALSXP, 1));
 133                 REAL(root_edge)[0] = tmpd;
 134                 SET_VECTOR_ELT(phy, 3, root_edge);
 135                 UNPROTECT(1);
 136         } else PROTECT(phy = allocVector(VECSXP, 3));
 137
 138         SET_VECTOR_ELT(phy, 0, edge);
 139         SET_VECTOR_ELT(phy, 1, edge_length);
 140         SET_VECTOR_ELT(phy, 2, Nnode);
 141
 142         UNPROTECT(5);
 143         return phy;
 144 }
 145
 146 #undef ADD_INTERNAL_EDGE
 147 #undef ADD_TERMINAL_EDGE
 148 #undef GO_DOWN