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

Go to the source code of this file.

Functions

static R_INLINE SEXP paste0 (SEXP a, SEXP b, SEXP c)
 
static SEXP pomp_default_rinit (SEXP params, SEXP Pnames, int npar, int nrep, int nsim)
 
static R_INLINE SEXP add_args (SEXP args, SEXP Pnames, SEXP Cnames)
 
static R_INLINE SEXP eval_call (SEXP fn, SEXP args, double *t0, double *p, int npar, double *c, int ncov)
 
static R_INLINE SEXP ret_array (int m, int n, SEXP names)
 
SEXP do_rinit (SEXP object, SEXP params, SEXP t0, SEXP nsim, SEXP gnsi)
 

Function Documentation

◆ add_args()

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

Definition at line 21 of file rinit.c.

22 {
23  SEXP var;
24  int v;
25 
26  PROTECT(args = VectorToPairList(args));
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  // Time
47  var = NEW_NUMERIC(1);
48  args = LCONS(var,args);
49  UNPROTECT(1);
50  PROTECT(args);
51  SET_TAG(args,install("t0"));
52 
53  UNPROTECT(1);
54  return args;
55 
56 }
SEXP args
Definition: trajectory.c:139
Here is the caller graph for this function:

◆ do_rinit()

SEXP do_rinit ( SEXP  object,
SEXP  params,
SEXP  t0,
SEXP  nsim,
SEXP  gnsi 
)

Definition at line 91 of file rinit.c.

92 {
93 
94  SEXP Pnames, Cnames, Snames, pcnames;
95  SEXP x = R_NilValue;
96  SEXP pompfun, fn, args;
98  lookup_table_t covariate_table;
99  SEXP cvec;
100  double *cov;
101  int *dim;
102  int npar, nrep, nvar, ncovars, nsims, ns;
103 
104  nsims = *(INTEGER(AS_INTEGER(nsim)));
105  PROTECT(params = as_matrix(params));
106  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
107  PROTECT(pcnames = GET_COLNAMES(GET_DIMNAMES(params)));
108 
109  dim = INTEGER(GET_DIM(params));
110  npar = dim[0]; nrep = dim[1];
111  ns = nsims*nrep;
112 
113  // set up the covariate table
114  covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
115  PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
116  PROTECT(cvec = NEW_NUMERIC(ncovars));
117  cov = REAL(cvec);
118 
119  table_lookup(&covariate_table,*(REAL(t0)),cov);
120 
121  // extract userdata
122  PROTECT(args = GET_SLOT(object,install("userdata")));
123 
124  PROTECT(pompfun = GET_SLOT(object,install("rinit")));
125  PROTECT(Snames = GET_SLOT(pompfun,install("statenames")));
126 
127  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
128 
129  int nprotect = 9;
130 
131  switch (mode) {
132  case Rfun: {
133 
134  SEXP ans;
135  double *time = REAL(AS_NUMERIC(t0));
136  double *ps = REAL(params);
137  double *xs, *xt = NULL;
138  int *midx;
139  int j;
140 
141  PROTECT(args = add_args(args,Pnames,Cnames));
142  PROTECT(ans = eval_call(fn,args,time,ps,npar,cov,ncovars));
143  PROTECT(ans = AS_NUMERIC(ans));
144  PROTECT(Snames = GET_NAMES(ans));
145 
146  if (invalid_names(Snames))
147  err("user 'rinit' must return a named numeric vector.");
148 
149  nvar = LENGTH(ans);
150  xs = REAL(ans);
151  midx = INTEGER(PROTECT(match(Pnames,Snames,0)));
152 
153  for (j = 0; j < nvar; j++) {
154  if (midx[j] != 0)
155  err("a state variable and a parameter share the name: '%s'.",CHAR(STRING_ELT(Snames,j)));
156  }
157 
158  PROTECT(x = ret_array(nvar,ns,Snames));
159  xt = REAL(x);
160 
161  memcpy(xt,xs,nvar*sizeof(double));
162 
163  nprotect += 6;
164 
165  for (j = 1, xt += nvar; j < ns; j++, xt += nvar) {
166  PROTECT(ans = eval_call(fn,args,time,ps+npar*(j%nrep),npar,cov,ncovars));
167  xs = REAL(ans);
168  if (LENGTH(ans) != nvar)
169  err("user 'rinit' returns vectors of variable length.");
170  memcpy(xt,xs,nvar*sizeof(double));
171  UNPROTECT(1);
172  }
173 
174  }
175 
176  break;
177 
178  case native: case regNative: {
179 
180  int *sidx, *pidx, *cidx;
181  double *xt, *ps, time;
182  pomp_rinit *ff = NULL;
183  int j;
184 
185  nvar = *INTEGER(GET_SLOT(object,install("nstatevars")));
186  PROTECT(x = ret_array(nvar,ns,Snames)); nprotect++;
187 
188  sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
189  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
190  cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
191 
192  // address of native routine
193  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
194 
195  GetRNGstate();
196 
197  time = *(REAL(t0));
198 
199  // loop over replicates
200  for (j = 0, xt = REAL(x), ps = REAL(params); j < ns; j++, xt += nvar)
201  (*ff)(xt,ps+npar*(j%nrep),time,sidx,pidx,cidx,cov);
202 
203  PutRNGstate();
204 
205  }
206 
207  break;
208 
209  default: {
210 
211  PROTECT(x = pomp_default_rinit(params,Pnames,npar,nrep,ns)); nprotect++;
212 
213  }
214 
215  break;
216 
217  }
218 
219  // now add column names
220  if (nrep > 1) {
221  SEXP dn, xn;
222  int k, *p;
223 
224  if (isNull(pcnames)) {
225  PROTECT(pcnames = NEW_INTEGER(nrep)); nprotect++;
226  for (k = 0, p = INTEGER(pcnames); k < nrep; k++, p++) *p = k+1;
227  }
228 
229  if (nsims > 1) {
230  int k, *sp;
231  SEXP us;
232  PROTECT(us = mkString("_"));
233  PROTECT(xn = NEW_INTEGER(ns));
234  for (k = 0, sp = INTEGER(xn); k < ns; k++, sp++) *sp = (k/nrep)+1;
235  PROTECT(xn = paste0(pcnames,us,xn));
236  PROTECT(dn = GET_DIMNAMES(x));
237  nprotect += 4;
238  SET_ELEMENT(dn,1,xn);
239  SET_DIMNAMES(x,dn);
240 
241  } else {
242 
243  PROTECT(dn = GET_DIMNAMES(x)); nprotect++;
244  SET_ELEMENT(dn,1,pcnames);
245  SET_DIMNAMES(x,dn);
246 
247  }
248  }
249 
250  UNPROTECT(nprotect);
251  return x;
252 }
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
void pomp_rinit(double *x, const double *p, double t0, const int *stateindex, const int *parindex, const int *covindex, const double *covars)
Definition: pomp.h:40
#define err(...)
Definition: pomp.h:21
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
static R_INLINE SEXP as_matrix(SEXP x)
Definition: pomp_defines.h:145
static R_INLINE SEXP ret_array(int m, int n, SEXP names)
Definition: rinit.c:79
static R_INLINE SEXP paste0(SEXP a, SEXP b, SEXP c)
Definition: rinit.c:10
static SEXP pomp_default_rinit(SEXP params, SEXP Pnames, int npar, int nrep, int nsim)
Definition: rinit.c:254
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t0, double *p, int npar, double *c, int ncov)
Definition: rinit.c:59
static R_INLINE SEXP add_args(SEXP args, SEXP Pnames, SEXP Cnames)
Definition: rinit.c:21
SEXP Snames
Definition: trajectory.c:140
int ncovars
Definition: trajectory.c:133
SEXP params
Definition: trajectory.c:128
pompfunmode mode
Definition: trajectory.c:126
SEXP cov
Definition: trajectory.c:129
SEXP fn
Definition: trajectory.c:138
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 *  t0,
double *  p,
int  npar,
double *  c,
int  ncov 
)
static

Definition at line 58 of file rinit.c.

62  {
63 
64  SEXP var = args, ans, ob;
65  int v;
66 
67  *(REAL(CAR(var))) = *t0; var = CDR(var);
68  for (v = 0; v < npar; v++, p++, var=CDR(var)) *(REAL(CAR(var))) = *p;
69  for (v = 0; v < ncov; v++, c++, var=CDR(var)) *(REAL(CAR(var))) = *c;
70 
71  PROTECT(ob = LCONS(fn,args));
72  PROTECT(ans = eval(ob,CLOENV(fn)));
73 
74  UNPROTECT(2);
75  return ans;
76 
77 }
Here is the caller graph for this function:

◆ paste0()

static R_INLINE SEXP paste0 ( SEXP  a,
SEXP  b,
SEXP  c 
)
static

Definition at line 10 of file rinit.c.

10  {
11  SEXP p, v;
12  PROTECT(p = lang4(install("paste0"),a,b,c));
13  PROTECT(v = eval(p,R_BaseEnv));
14  UNPROTECT(2);
15  return v;
16 }
Here is the caller graph for this function:

◆ pomp_default_rinit()

static SEXP pomp_default_rinit ( SEXP  params,
SEXP  Pnames,
int  npar,
int  nrep,
int  nsim 
)
static

Definition at line 254 of file rinit.c.

256 {
257 
258  SEXP fcall, pat, ivpnames, statenames, x;
259  int *pidx;
260  int nvar, j, k;
261  double *xp, *pp;
262 
263  // extract names of IVPs using 'grep'
264  PROTECT(pat = mkString("[\\_\\.]0$"));
265  PROTECT(fcall = LCONS(ScalarLogical(1),R_NilValue));
266  SET_TAG(fcall,install("value"));
267  PROTECT(fcall = LCONS(Pnames,fcall));
268  SET_TAG(fcall,install("x"));
269  PROTECT(fcall = LCONS(pat,fcall));
270  SET_TAG(fcall,install("pattern"));
271  PROTECT(fcall = LCONS(install("grep"),fcall));
272  PROTECT(ivpnames = eval(fcall,R_BaseEnv));
273 
274  nvar = LENGTH(ivpnames);
275  if (nvar < 1)
276  warn("in default 'rinit': there are no parameters with suffix '.0' or '_0'. See '?rinit_spec'.");
277 
278  pidx = INTEGER(PROTECT(match(Pnames,ivpnames,0)));
279  for (k = 0; k < nvar; k++) pidx[k]--;
280 
281  // construct names of state variables using 'sub'
282  PROTECT(fcall = LCONS(ivpnames,R_NilValue));
283  SET_TAG(fcall,install("x"));
284  PROTECT(fcall = LCONS(mkString(""),fcall));
285  SET_TAG(fcall,install("replacement"));
286  PROTECT(fcall = LCONS(pat,fcall));
287  SET_TAG(fcall,install("pattern"));
288  PROTECT(fcall = LCONS(install("sub"),fcall));
289  PROTECT(statenames = eval(fcall,R_BaseEnv));
290 
291  PROTECT(x = ret_array(nvar,nsim,statenames));
292 
293  for (j = 0, xp = REAL(x); j < nsim; j++) {
294  pp = REAL(params) + npar*(j%nrep);
295  for (k = 0; k < nvar; k++, xp++) *xp = pp[pidx[k]];
296  }
297 
298  UNPROTECT(13);
299  return x;
300 }
#define warn(...)
Definition: pomp.h:22
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ret_array()

static R_INLINE SEXP ret_array ( int  m,
int  n,
SEXP  names 
)
static

Definition at line 79 of file rinit.c.

80 {
81  int dim[2] = {m, n};
82  const char *dimnm[2] = {"name",".id"};
83  SEXP X;
84  PROTECT(X = makearray(2,dim));
85  fillrownames(X,names);
86  fixdimnames(X,dimnm,2);
87  UNPROTECT(1);
88  return X;
89 }
#define X
Definition: gompertz.c:14
static R_INLINE void fixdimnames(SEXP x, const char **names, int n)
Definition: pomp_defines.h:129
static R_INLINE void fillrownames(SEXP x, SEXP names)
Definition: pomp_defines.h:87
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: