src/reorder_phylo.c

   1 /* reorder_phylo.c       2012-09-03 */
   2
   3 /* Copyright 2008-2012 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
  10 static int iii;
  11
  12 void foo_reorder(int node, int n, int m, int *e1, int *e2, int *neworder, int *L, int *pos)
  13 {
  14         int i = node - n - 1, j, k;
  15
  16 /* 'i' is the C index corresponding to 'node' */
  17
  18         for (j = 0; j < pos[i]; j++) {
  19                 k = L[i + m * j];
  20                 neworder[iii++] = k + 1;
  21                 if (e2[k] > n) /* is it an internal edge? */
  22                         foo_reorder(e2[k], n, m, e1, e2, neworder, L, pos);
  23         }
  24 }
  25
  26 void bar_reorder(int node, int n, int m, int *e1, int *e2, int *neworder, int *L, int *pos)
  27 {
  28         int i = node - n - 1, j, k;
  29
  30         for (j = pos[i] - 1; j >= 0; j--)
  31                 neworder[iii--] = L[i + m * j] + 1;
  32
  33         for (j = 0; j < pos[i]; j++) {
  34                 k = e2[L[i + m * j]];
  35                 if (k > n)
  36                         bar_reorder(k, n, m, e1, e2, neworder, L, pos);
  37         }
  38 }
  39
  40 void neworder_phylo(int *n, int *e1, int *e2, int *N, int *neworder, int *order)
  41 /* n: nb of tips
  42    m: nb of nodes
  43    N: nb of edges */
  44 {
  45         int i, j, k, *L, *pos, m = *N - *n + 1, degrmax = *n - m + 1;
  46
  47 /* degrmax is the largest value that a node degree can be */
  48
  49 /* L is a 1-d array storing, for each node, the C indices of the rows of
  50    the edge matrix where the node is ancestor (i.e., present in the 1st
  51    column). It is used in the same way than a matrix (which is actually
  52    a vector) is used in R as a 2-d structure. */
  53
  54         L = (int*)R_alloc(m * degrmax, sizeof(int));
  55
  56 /* pos gives the position for each 'row' of L, that is the number of elements
  57    which have already been stored for that 'row'. */
  58
  59         pos = (int*)R_alloc(m, sizeof(int));
  60         memset(pos, 0, m * sizeof(int));
  61
  62 /* we now go down along the edge matrix */
  63
  64         for (i = 0; i < *N; i++) {
  65                 k = e1[i] - *n - 1; /* k is the 'row' index in L corresponding to node e1[i] */
  66                 j = pos[k]; /* the current 'column' position corresping to k */
  67                 pos[k]++; /* increment in case the same node is found in another row of the edge matrix */
  68                 L[k + m * j] = i;
  69         }
  70
  71 /* L is now ready: we can start the recursive calls. */
  72 /* We start with the root 'n + 1': its index will be changed into
  73    the corresponding C index inside the recursive function. */
  74
  75         switch(*order) {
  76         case 1 : iii = 0;
  77                 foo_reorder(*n + 1, *n, m, e1, e2, neworder, L, pos);
  78                 break;
  79         case 2 : iii = *N - 1;
  80                 bar_reorder(*n + 1, *n, m, e1, e2, neworder, L, pos);
  81                 break;
  82         }
  83 }
  84
  85 #define DO_NODE_PRUNING\
  86     /* go back down in `edge' to set `neworder' */\
  87     for (j = 0; j <= i; j++) {\
  88         /* if find the edge where `node' is */\
  89         /* the descendant, make as ready */\
  90         if (edge2[j] == node) ready[j] = 1;\
  91         if (edge1[j] != node) continue;\
  92         neworder[nextI] = j + 1;\
  93         ready[j] = 0; /* mark the edge as done */\
  94         nextI++;\
  95     }
  96
  97 void neworder_pruningwise(int *ntip, int *nnode, int *edge1,
  98                           int *edge2, int *nedge, int *neworder)
  99 {
 100     int *ready, *Ndegr, i, j, node, nextI, n;
 101
 102     nextI = *ntip +  *nnode;
 103     Ndegr = (int*)R_alloc(nextI, sizeof(int));
 104     memset(Ndegr, 0, nextI*sizeof(int));
 105     for (i = 0; i < *nedge; i++) (Ndegr[edge1[i] - 1])++;
 106
 107     ready = (int*)R_alloc(*nedge, sizeof(int));
 108
 109     /* `ready' indicates whether an edge is ready to be */
 110     /* collected; only the terminal edges are initially ready */
 111     for (i = 0; i < *nedge; i++)
 112             ready[i] = (edge2[i] <= *ntip) ? 1 : 0;
 113
 114     /* `n' counts the number of times a node has been seen. */
 115     /* This algo will work if the tree is in cladewise order, */
 116     /* so that the nodes of "cherries" will be contiguous in `edge'. */
 117     n = 0;
 118     nextI = 0;
 119     while (nextI < *nedge - Ndegr[*ntip]) {
 120         for (i = 0; i < *nedge; i++) {
 121             if (!ready[i]) continue;
 122             if (!n) {
 123                 /* if found an edge ready, initialize `node' and start counting */
 124                 node = edge1[i];
 125                 n = 1;
 126             } else { /* else counting has already started */
 127                 if (edge1[i] == node) n++;
 128                 else {
 129                     /* if the node has changed we checked that all edges */
 130                     /* from `node' have been found */
 131                     if (n == Ndegr[node - 1]) {
 132                         DO_NODE_PRUNING
 133                     }
 134                     /* in all cases reset `n' and `node' and carry on */
 135                     node = edge1[i];
 136                     n = 1;
 137                 }
 138             } /* go to the next edge */
 139             /* if at the end of `edge', check that we can't do a node */
 140             if (n == Ndegr[node - 1]) {
 141                 DO_NODE_PRUNING
 142                 n = 0;
 143             }
 144         }
 145     }
 146     for (i = 0; i < *nedge; i++) {
 147         if (!ready[i]) continue;
 148         neworder[nextI] = i + 1;
 149         nextI++;
 150     }
 151 }
 152