src/rTrait.c

   1 /* rTrait.c       2011-06-25 */
   2
   3 /* Copyright 2010-2011 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 rTraitCont(int *model, int *Nedge, int *edge1, int *edge2, double *el,
  11                 double *sigma, double *alpha, double *theta, double *x)
  12 {
  13 /* The tree must be in pruningwise order */
  14         int i;
  15
  16         switch(*model) {
  17         case 1 : for (i = *Nedge - 1; i >= 0; i--) {
  18                         GetRNGstate();
  19                         x[edge2[i]] = x[edge1[i]] + sqrt(el[i]) * sigma[i] * norm_rand();
  20                         PutRNGstate();
  21                 }
  22                 break;
  23         case 2 : for (i = *Nedge - 1; i >= 0; i--) {
  24                         if (alpha[i]) {
  25                                 alphaT = alpha[i] * el[i];
  26                                 M = exp(-alphaT);
  27                                 S = sigma[i] * sqrt((1 - exp(-2 * alphaT))/(2 * alpha[i]));
  28                         } else {
  29                                 M = 1;
  30                                 S = sqrt(el[i]) * sigma[i];
  31                         }
  32                         GetRNGstate();
  33                         x[edge2[i]] = x[edge1[i]] * M +  theta[i] * (1 - M) + S * norm_rand();
  34                         PutRNGstate();
  35                 }
  36                 break;
  37         }
  38 }