src/pic.c

   1 /* pic.c       2006-11-13 */
   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
  10 void pic(int *ntip, int *nnode, int *edge1, int *edge2,
  11          double *edge_len, double *phe, double *contr,
  12          double *var_contr, int *var, int *scaled)
  13 {
  14 /* The tree must be in pruningwise order */
  15     int anc, d1, d2, ic, i, j, k;
  16     double sumbl;
  17
  18     for (i = 0; i < *ntip * 2 - 3; i += 2) {
  19         j = i + 1;
  20         anc = edge1[i];
  21         d1 = edge2[i] - 1;
  22         d2 = edge2[j] - 1;
  23         sumbl = edge_len[i] + edge_len[j];
  24         ic = anc - *ntip - 1;
  25         contr[ic] = phe[d1] - phe[d2];
  26         if (*scaled) contr[ic] = contr[ic]/sqrt(sumbl);
  27         if (*var) var_contr[ic] = sumbl;
  28         phe[anc - 1] = (phe[d1]*edge_len[j] + phe[d2]*edge_len[i])/sumbl;
  29         /* find the edge where `anc' is a descendant (except if at the root):
  30            it is obviously below the j'th edge */
  31         if (j != *ntip * 2 - 3) {
  32             k = j + 1;
  33             while (edge2[k] != anc) k++;
  34             edge_len[k] = edge_len[k] + edge_len[i]*edge_len[j]/sumbl;
  35         }
  36     }
  37 }