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         double alphaT, M, S;
  16
  17         switch(*model) {
  18         case 1 : for (i = *Nedge - 1; i >= 0; i--) {
  19                         GetRNGstate();
  20                         x[edge2[i]] = x[edge1[i]] + sqrt(el[i]) * sigma[i] * norm_rand();
  21                         PutRNGstate();
  22                 }
  23                 break;
  24         case 2 : for (i = *Nedge - 1; i >= 0; i--) {
  25                         if (alpha[i]) {
  26                                 alphaT = alpha[i] * el[i];
  27                                 M = exp(-alphaT);
  28                                 S = sigma[i] * sqrt((1 - exp(-2 * alphaT))/(2 * alpha[i]));
  29                         } else { /* same than if (alpha[i] ==0) */
  30                                 M = 1;
  31                                 S = sqrt(el[i]) * sigma[i];
  32                         }
  33                         GetRNGstate();
  34                         x[edge2[i]] = x[edge1[i]] * M +  theta[i] * (1 - M) + S * norm_rand();
  35                         PutRNGstate();
  36                 }
  37                 break;
  38         }
  39 }