[Pomp-commits] r390 - in pkg: . src

[noreply at r-forge.r-project.org]

Author: kingaa
Date: 2010-10-20 13:20:08 +0200 (Wed, 20 Oct 2010)
New Revision: 390

Added:
   pkg/src/blowfly.c
Modified:
   pkg/DESCRIPTION
Log:

- add C codes for blowfly example


Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2010-10-19 13:49:03 UTC (rev 389)
+++ pkg/DESCRIPTION	2010-10-20 11:20:08 UTC (rev 390)
@@ -2,7 +2,7 @@
 Type: Package
 Title: Statistical inference for partially observed Markov processes
 Version: 0.34-3
-Date: 2010-10-19
+Date: 2010-10-20
 Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, 
 	Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
 Maintainer: Aaron A. King <kingaa at umich.edu>

Added: pkg/src/blowfly.c
===================================================================
--- pkg/src/blowfly.c	                        (rev 0)
+++ pkg/src/blowfly.c	2010-10-20 11:20:08 UTC (rev 390)
@@ -0,0 +1,53 @@
+// dear emacs, please treat this as -*- C++ -*-
+
+#include <Rmath.h>
+
+#include "pomp.h"
+
+#define LOG_P       (p[parindex[0]]) // growth rate
+#define LOG_NZERO   (p[parindex[1]]) // density-dependence parameter
+#define LOG_DELTA   (p[parindex[2]]) // survival parameter
+#define LOG_SIGMAP  (p[parindex[3]]) // recruitment noise SD
+#define LOG_SIGMAD  (p[parindex[4]]) // survivorship noise SD
+#define TAU         (p[parindex[5]]) // delay
+
+#define N          (&x[stateindex[0]]) // total population
+#define R           (x[stateindex[1]]) // recruits
+#define S           (x[stateindex[2]]) // survivors
+#define E           (x[stateindex[3]]) // recruitment noise
+#define EPS         (x[stateindex[4]]) // survival noise
+
+// Ricker model with log-normal process noise
+void _blowfly_simulator (double *x, const double *p, 
+			 const int *stateindex, const int *parindex, const int *covindex,
+			 int covdim, const double *covar, 
+			 double t, double dt)
+{
+  double var_p = exp(2*LOG_SIGMAP);
+  double var_d = exp(2*LOG_SIGMAD);
+  double e = (var_p > 0.0) ? rgamma(1.0/var_p,var_p) : 1.0;
+  double eps = (var_d > 0.0) ? rgamma(1.0/var_d,var_d) : 1.0;
+  double P = exp(LOG_P);
+  int tau = (int) TAU;
+  int k;
+  R = rpois(P*N[tau]*exp(-N[tau]/exp(LOG_NZERO))*e);
+  S = rbinom(N[0],exp(-exp(LOG_DELTA)*eps));
+  E = e;
+  EPS = eps;
+  for (k = tau; k > 0; k--) N[k] = N[k-1];
+  N[0] = R+S;
+}
+
+#undef N
+#undef R
+#undef S
+#undef E
+#undef EPS
+
+#undef LOG_P
+#undef LOG_NZERO
+#undef LOG_DELTA
+#undef LOG_SIGMAP
+#undef LOG_SIGMAD
+#undef TAU
+