src/reorder_phylo.c

   1 /* reorder_phylo.c       2012-08-01 */
   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 void neworder_cladewise(int *n, int *edge1, int *edge2,
  11                         int *N, int *neworder)
  12 /* n: nb of tips, N: nb of edges */
  13 {
  14     int i, j, k, node, *done, dn, *node_back, eb;
  15     done = &dn;
  16     node_back = &eb;
  17
  18     /* done: indicates whether an edge has been collected
  19        node_back: the series of node from the root to `node'
  20        node: the current node */
  21
  22     done = (int*)R_alloc(*N, sizeof(int));
  23     node_back = (int*)R_alloc(*N + 2 - *n, sizeof(int));
  24     memset(done, 0, *N * sizeof(int));
  25
  26     j = k = 0;
  27     node = *n + 1;
  28     while (j < *N) {
  29         for (i = 0; i < *N; i++) {
  30             if (done[i] || edge1[i] != node) continue;
  31             neworder[j] = i + 1;
  32             j++;
  33             done[i] = 1;
  34             if (edge2[i] > *n) {
  35                 node_back[k] = node;
  36                 k++;
  37                 node = edge2[i];
  38                 /* if found a new node, reset the loop */
  39                 i = -1;
  40             }
  41         }
  42         /* if arrived at the end of `edge', go down one node */
  43         k--;
  44         node = node_back[k];
  45     }
  46 }
  47
  48 #define DO_NODE_PRUNING\
  49     /* go back down in `edge' to set `neworder' */\
  50     for (j = 0; j <= i; j++) {\
  51         /* if find the edge where `node' is */\
  52         /* the descendant, make as ready */\
  53         if (edge2[j] == node) ready[j] = 1;\
  54         if (edge1[j] != node) continue;\
  55         neworder[nextI] = j + 1;\
  56         ready[j] = 0; /* mark the edge as done */\
  57         nextI++;\
  58     }
  59
  60 void neworder_pruningwise(int *ntip, int *nnode, int *edge1,
  61                           int *edge2, int *nedge, int *neworder)
  62 {
  63     int *ready, *Ndegr, i, j, node, nextI, n;
  64
  65     nextI = *ntip +  *nnode;
  66     Ndegr = (int*)R_alloc(nextI, sizeof(int));
  67     memset(Ndegr, 0, nextI*sizeof(int));
  68     for (i = 0; i < *nedge; i++) (Ndegr[edge1[i] - 1])++;
  69
  70     ready = (int*)R_alloc(*nedge, sizeof(int));
  71
  72     /* `ready' indicates whether an edge is ready to be */
  73     /* collected; only the terminal edges are initially ready */
  74     for (i = 0; i < *nedge; i++)
  75             ready[i] = (edge2[i] <= *ntip) ? 1 : 0;
  76
  77     /* `n' counts the number of times a node has been seen. */
  78     /* This algo will work if the tree is in cladewise order, */
  79     /* so that the nodes of "cherries" will be contiguous in `edge'. */
  80     n = 0;
  81     nextI = 0;
  82     while (nextI < *nedge - Ndegr[*ntip]) {
  83         for (i = 0; i < *nedge; i++) {
  84             if (!ready[i]) continue;
  85             if (!n) {
  86                 /* if found an edge ready, initialize `node' and start counting */
  87                 node = edge1[i];
  88                 n = 1;
  89             } else { /* else counting has already started */
  90                 if (edge1[i] == node) n++;
  91                 else {
  92                     /* if the node has changed we checked that all edges */
  93                     /* from `node' have been found */
  94                     if (n == Ndegr[node - 1]) {
  95                         DO_NODE_PRUNING
  96                     }
  97                     /* in all cases reset `n' and `node' and carry on */
  98                     node = edge1[i];
  99                     n = 1;
 100                 }
 101             } /* go to the next edge */
 102             /* if at the end of `edge', check that we can't do a node */
 103             if (n == Ndegr[node - 1]) {
 104                 DO_NODE_PRUNING
 105                 n = 0;
 106             }
 107         }
 108     }
 109     for (i = 0; i < *nedge; i++) {
 110         if (!ready[i]) continue;
 111         neworder[nextI] = i + 1;
 112         nextI++;
 113     }
 114 }