src/dist_nodes.c

   1 /* dist_nodes.c       2012-08-14 */
   2
   3 /* Copyright 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 #define DINDEX2(i, j) i + NM * j
  11
  12 /* The algorithm is pretty simple: the tree must be in cladewise order
  13    because the edges are visited contiguously. Each edge gives trivially
  14    one distance, then by moving up along the edge matrix, one finds nodes
  15    that have already been visited and the distance matrix can be updated. */
  16
  17 void dist_nodes(int *n, int *m, int *e1, int *e2, double *el, int *N, double *D)
  18 /* n: nb of tips, m: nb of nodes, N: nb of edges */
  19 {
  20         int i, j, k, a, d, NM = *n + *m, ROOT;
  21         double x;
  22
  23         ROOT = e1[0]; d = e2[0]; /* the 2 nodes of the 1st edge */
  24         D[DINDEX2(ROOT, d)] = D[DINDEX2(d, ROOT)] = el[0]; /* the 1st edge gives the 1st distance */
  25
  26 /* go down along the edge matrix
  27    starting at the 2nd edge: */
  28         for (i = 1; i < *N; i++) {
  29                 a = e1[i]; d = e2[i]; x = el[i]; /* get the i-th nodes and branch length */
  30                 D[DINDEX2(a, d)] = D[DINDEX2(d, a)] = x;
  31                 /* then go up along the edge matrix from the i-th edge
  32                    to visit the nodes already visited and update the distances: */
  33                 for (j = i - 1; j >= 0; j--) {
  34                         k = e2[j];
  35                         if (k == a) continue;
  36                         D[DINDEX2(k, d)] = D[DINDEX2(d, k)] = D[DINDEX2(a, k)] + x;
  37                 }
  38                 if (k != ROOT)
  39                         D[DINDEX2(ROOT, d)] = D[DINDEX2(d, ROOT)] = D[DINDEX2(ROOT, a)] + x;
  40         }
  41 }