src/reorder_phylo.c

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