pomp
Inference for partially observed Markov processes
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
euler.c File Reference
#include <R.h>
#include <Rmath.h>
#include <Rdefines.h>
#include "internal.h"
Include dependency graph for euler.c:

Go to the source code of this file.

Functions

static R_INLINE SEXP add_args (SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
 
static R_INLINE SEXP eval_call (SEXP fn, SEXP args, double *t, double *dt, double *x, int nvar, double *p, int npar, double *c, int ncov)
 
static R_INLINE SEXP ret_array (int n, int nreps, int ntimes, SEXP names)
 
SEXP euler_simulator (SEXP func, SEXP xstart, SEXP tstart, SEXP times, SEXP params, double deltat, rprocmode method, SEXP accumvars, SEXP covar, SEXP args, SEXP gnsi)
 
int num_euler_steps (double t1, double t2, double *deltat)
 
int num_map_steps (double t1, double t2, double deltat)
 

Function Documentation

◆ add_args()

static R_INLINE SEXP add_args ( SEXP  args,
SEXP  Snames,
SEXP  Pnames,
SEXP  Cnames 
)
static

Definition at line 10 of file euler.c.

11 {
12 
13  SEXP var;
14  int v;
15 
16  PROTECT(args = VectorToPairList(args));
17 
18  // we construct the call from end to beginning
19  // delta.t, covariates, parameter, states, then time
20 
21  // 'delta.t'
22  var = NEW_NUMERIC(1);
23  args = LCONS(var,args);
24  UNPROTECT(1);
25  PROTECT(args);
26  SET_TAG(args,install("delta.t"));
27 
28  // Covariates
29  for (v = LENGTH(Cnames)-1; v >= 0; v--) {
30  var = NEW_NUMERIC(1);
31  args = LCONS(var,args);
32  UNPROTECT(1);
33  PROTECT(args);
34  SET_TAG(args,installChar(STRING_ELT(Cnames,v)));
35  }
36 
37  // Parameters
38  for (v = LENGTH(Pnames)-1; v >= 0; v--) {
39  var = NEW_NUMERIC(1);
40  args = LCONS(var,args);
41  UNPROTECT(1);
42  PROTECT(args);
43  SET_TAG(args,installChar(STRING_ELT(Pnames,v)));
44  }
45 
46  // Latent state variables
47  for (v = LENGTH(Snames)-1; v >= 0; v--) {
48  var = NEW_NUMERIC(1);
49  args = LCONS(var,args);
50  UNPROTECT(1);
51  PROTECT(args);
52  SET_TAG(args,installChar(STRING_ELT(Snames,v)));
53  }
54 
55  // Time
56  var = NEW_NUMERIC(1);
57  args = LCONS(var,args);
58  UNPROTECT(1);
59  PROTECT(args);
60  SET_TAG(args,install("t"));
61 
62  UNPROTECT(1);
63  return args;
64 
65 }
SEXP Snames
Definition: trajectory.c:140
SEXP args
Definition: trajectory.c:139
Here is the caller graph for this function:

◆ euler_simulator()

SEXP euler_simulator ( SEXP  func,
SEXP  xstart,
SEXP  tstart,
SEXP  times,
SEXP  params,
double  deltat,
rprocmode  method,
SEXP  accumvars,
SEXP  covar,
SEXP  args,
SEXP  gnsi 
)

Definition at line 107 of file euler.c.

112  {
113 
115  int nvars, npars, nreps, ntimes, nzeros, ncovars;
116  double *cov, t0;
117  SEXP cvec, X, fn;
118 
119  if (deltat <= 0) err("'delta.t' should be a positive number."); // #nocov
120 
121  int *dim;
122  dim = INTEGER(GET_DIM(xstart)); nvars = dim[0]; nreps = dim[1];
123  dim = INTEGER(GET_DIM(params)); npars = dim[0];
124  ntimes = LENGTH(times);
125 
126  PROTECT(tstart = AS_NUMERIC(tstart));
127  PROTECT(times = AS_NUMERIC(times));
128  t0 = *(REAL(tstart));
129  if (t0 > *(REAL(times))) err("'t0' must be no later than 'times[1]'.");
130 
131  SEXP Snames, Pnames, Cnames;
132  PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(xstart)));
133  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
134  PROTECT(Cnames = get_covariate_names(covar));
135 
136  // set up the covariate table
137  lookup_table_t covariate_table = make_covariate_table(covar,&ncovars);
138  PROTECT(cvec = NEW_NUMERIC(ncovars));
139  cov = REAL(cvec);
140 
141  // indices of accumulator variables
142  nzeros = LENGTH(accumvars);
143  int *zidx = INTEGER(PROTECT(matchnames(Snames,accumvars,"state variables")));
144 
145  // extract user function
146  PROTECT(fn = pomp_fun_handler(func,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
147 
148  // array to hold results
149  PROTECT(X = ret_array(nvars,nreps,ntimes,Snames));
150 
151  // copy the start values into the result array
152  memcpy(REAL(X),REAL(xstart),nvars*nreps*sizeof(double));
153 
154  // set up
155 
156  int nprotect = 9;
157  int *pidx = 0, *sidx = 0, *cidx = 0;
158  pomp_onestep_sim *ff = NULL;
159 
160  switch (mode) {
161 
162  case Rfun:
163  // construct list of all arguments
164  PROTECT(args = add_args(args,Snames,Pnames,Cnames)); nprotect++;
165  break;
166 
167  case native: case regNative:
168  // construct state, parameter, covariate indices
169  sidx = INTEGER(GET_SLOT(func,install("stateindex")));
170  pidx = INTEGER(GET_SLOT(func,install("paramindex")));
171  cidx = INTEGER(GET_SLOT(func,install("covarindex")));
172 
173  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
174 
175  GetRNGstate();
176  break;
177 
178  default: // #nocov
179  err("unrecognized 'mode' %d",mode); // #nocov
180  break; // #nocov
181  }
182 
183  // main computation loop
184  int step;
185  double *xt, *time, t;
186  int first = 1;
187 
188  for (
189  step = 0, xt = REAL(X), time = REAL(times), t = t0;
190  step < ntimes;
191  step++, xt += nvars*nreps
192  ) {
193 
194  double dt;
195  int nstep = 0;
196  int i, j, k;
197 
198  R_CheckUserInterrupt();
199 
200  if (t > time[step]) err("'times' must be an increasing sequence."); // #nocov
201 
202  // set accumulator variables to zero
203  for (j = 0; j < nreps; j++)
204  for (i = 0; i < nzeros; i++)
205  xt[zidx[i]+nvars*j] = 0.0;
206 
207  // determine size and number of time-steps
208  switch (method) {
209  case onestep: default: // one step
210  dt = time[step]-t;
211  nstep = 1;
212  break;
213  case discrete: // fixed step
214  dt = deltat;
215  nstep = num_map_steps(t,time[step],dt);
216  break;
217  case euler: // Euler method
218  dt = deltat;
219  nstep = num_euler_steps(t,time[step],&dt);
220  break;
221  }
222 
223  // loop over individual steps
224  for (k = 0; k < nstep; k++) {
225 
226  // interpolate the covar functions for the covariates
227  table_lookup(&covariate_table,t,cov);
228 
229  // loop over replicates
230  double *ap, *pm, *xm, *ps = REAL(params);
231 
232  for (j = 0, pm = ps, xm = xt; j < nreps; j++, pm += npars, xm += nvars) {
233 
234  switch (mode) {
235 
236  case Rfun:
237  {
238  SEXP ans, nm;
239 
240  if (first) {
241 
242  PROTECT(ans = eval_call(fn,args,&t,&dt,xm,nvars,pm,npars,cov,ncovars));
243 
244  PROTECT(nm = GET_NAMES(ans));
245  if (invalid_names(nm))
246  err("'rprocess' must return a named numeric vector.");
247  pidx = INTEGER(PROTECT(matchnames(nm,Snames,"state variables")));
248 
249  nprotect += 3;
250 
251  ap = REAL(AS_NUMERIC(ans));
252  for (i = 0; i < nvars; i++) xm[pidx[i]] = ap[i];
253 
254  first = 0;
255 
256  } else {
257 
258  PROTECT(ans = eval_call(fn,args,&t,&dt,xm,nvars,pm,npars,cov,ncovars));
259 
260  ap = REAL(AS_NUMERIC(ans));
261  for (i = 0; i < nvars; i++) xm[pidx[i]] = ap[i];
262 
263  UNPROTECT(1);
264 
265  }
266  }
267  break;
268 
269  case native: case regNative:
270  (*ff)(xm,pm,sidx,pidx,cidx,cov,t,dt);
271  break;
272 
273  default: // #nocov
274  err("unrecognized 'mode' %d",mode); // #nocov
275  break; // #nocov
276 
277  }
278 
279  }
280 
281  t += dt;
282 
283  if ((method == euler) && (k == nstep-2)) { // penultimate step
284  dt = time[step]-t;
285  t = time[step]-dt;
286  }
287 
288  }
289 
290  if (step < ntimes-1)
291  memcpy(xt+nvars*nreps,xt,nreps*nvars*sizeof(double));
292 
293  }
294 
295  // clean up
296  switch (mode) {
297 
298  case native: case regNative: {
299 
300  PutRNGstate();
301 
302  }
303 
304  break;
305 
306  case Rfun: default:
307 
308  break;
309 
310  }
311 
312  UNPROTECT(nprotect);
313  return X;
314 
315 }
lookup_table_t make_covariate_table(SEXP, int *)
Definition: lookup_table.c:11
SEXP pomp_fun_handler(SEXP, SEXP, pompfunmode *, SEXP, SEXP, SEXP, SEXP)
Definition: pomp_fun.c:30
void table_lookup(lookup_table_t *, double, double *)
Definition: lookup_table.c:53
SEXP get_covariate_names(SEXP)
Definition: lookup_table.c:7
int num_map_steps(double t1, double t2, double deltat)
Definition: euler.c:342
int num_euler_steps(double t1, double t2, double *deltat)
Definition: euler.c:317
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *dt, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition: euler.c:67
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition: euler.c:10
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition: euler.c:92
#define X
Definition: gompertz.c:14
void pomp_onestep_sim(double *x, const double *p, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t, double dt)
Definition: pomp.h:59
#define err(...)
Definition: pomp.h:21
@ discrete
Definition: pomp_defines.h:17
@ onestep
Definition: pomp_defines.h:17
@ euler
Definition: pomp_defines.h:17
static R_INLINE SEXP matchnames(SEXP provided, SEXP needed, const char *where)
Definition: pomp_defines.h:60
pompfunmode
Definition: pomp_defines.h:16
@ Rfun
Definition: pomp_defines.h:16
@ native
Definition: pomp_defines.h:16
@ undef
Definition: pomp_defines.h:16
@ regNative
Definition: pomp_defines.h:16
static R_INLINE int invalid_names(SEXP names)
Definition: pomp_defines.h:55
int npars
Definition: trajectory.c:132
int ncovars
Definition: trajectory.c:133
SEXP params
Definition: trajectory.c:128
pompfunmode mode
Definition: trajectory.c:126
int nreps
Definition: trajectory.c:134
SEXP cov
Definition: trajectory.c:129
SEXP fn
Definition: trajectory.c:138
int nvars
Definition: trajectory.c:131
Here is the call graph for this function:
Here is the caller graph for this function:

◆ eval_call()

static R_INLINE SEXP eval_call ( SEXP  fn,
SEXP  args,
double *  t,
double *  dt,
double *  x,
int  nvar,
double *  p,
int  npar,
double *  c,
int  ncov 
)
static

Definition at line 67 of file euler.c.

73 {
74 
75  SEXP var = args, ans, ob;
76  int v;
77 
78  *(REAL(CAR(var))) = *t; var = CDR(var);
79  for (v = 0; v < nvar; v++, x++, var=CDR(var)) *(REAL(CAR(var))) = *x;
80  for (v = 0; v < npar; v++, p++, var=CDR(var)) *(REAL(CAR(var))) = *p;
81  for (v = 0; v < ncov; v++, c++, var=CDR(var)) *(REAL(CAR(var))) = *c;
82  *(REAL(CAR(var))) = *dt; var = CDR(var);
83 
84  PROTECT(ob = LCONS(fn,args));
85  PROTECT(ans = eval(ob,CLOENV(fn)));
86 
87  UNPROTECT(2);
88  return ans;
89 
90 }
Here is the caller graph for this function:

◆ num_euler_steps()

int num_euler_steps ( double  t1,
double  t2,
double *  deltat 
)

Definition at line 317 of file euler.c.

317  {
318  double tol = sqrt(DBL_EPSILON);
319  int nstep;
320  // nstep will be the number of Euler steps to take in going from t1 to t2.
321  // note also that the stepsize changes.
322  // this choice is meant to be conservative
323  // (i.e., so that the actual deltat does not exceed the specified deltat
324  // by more than the relative tolerance 'tol')
325  // and to counteract roundoff error.
326  // It seems to work well, but is not guaranteed:
327  // suggestions would be appreciated.
328 
329  if (t1 >= t2) {
330  *deltat = 0.0;
331  nstep = 0;
332  } else if (t1+*deltat >= t2) {
333  *deltat = t2-t1;
334  nstep = 1;
335  } else {
336  nstep = (int) ceil((t2-t1)/(*deltat)/(1+tol));
337  *deltat = (t2-t1)/((double) nstep);
338  }
339  return nstep;
340 }
Here is the caller graph for this function:

◆ num_map_steps()

int num_map_steps ( double  t1,
double  t2,
double  deltat 
)

Definition at line 342 of file euler.c.

342  {
343  double tol = sqrt(DBL_EPSILON);
344  int nstep;
345  // nstep will be the number of discrete-time steps to take in going from t1 to t2.
346  nstep = (int) floor((t2-t1)/deltat/(1-tol));
347  return (nstep > 0) ? nstep : 0;
348 }
Here is the caller graph for this function:

◆ ret_array()

static R_INLINE SEXP ret_array ( int  n,
int  nreps,
int  ntimes,
SEXP  names 
)
static

Definition at line 92 of file euler.c.

93 {
94  int dim[3] = {n, nreps, ntimes};
95  const char *dimnm[3] = {"name", ".id", "time"};
96  SEXP Y;
97 
98  PROTECT(Y = makearray(3,dim));
99  setrownames(Y,names,3);
100  fixdimnames(Y,dimnm,3);
101 
102  UNPROTECT(1);
103  return Y;
104 
105 }
#define Y
Definition: gompertz.c:13
static R_INLINE void fixdimnames(SEXP x, const char **names, int n)
Definition: pomp_defines.h:129
static R_INLINE void setrownames(SEXP x, SEXP names, int rank)
Definition: pomp_defines.h:110
static R_INLINE SEXP makearray(int rank, const int *dim)
Definition: pomp_defines.h:40
Here is the call graph for this function:
Here is the caller graph for this function: