pomp
Inference for partially observed Markov processes
•All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
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}
args
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;
97 pompfunmode mode = undef;
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}
make_covariate_table
lookup_table_t make_covariate_table(SEXP, int *)
Definition lookup_table.c:11
pomp_fun_handler
SEXP pomp_fun_handler(SEXP, SEXP, pompfunmode *, SEXP, SEXP, SEXP, SEXP)
Definition pomp_fun.c:30
table_lookup
void table_lookup(lookup_table_t *, double, double *)
Definition lookup_table.c:53
get_covariate_names
SEXP get_covariate_names(SEXP)
Definition lookup_table.c:7
pomp_rinit
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
err
#define err(...)
Definition pomp.h:21
pompfunmode
pompfunmode
Definition pomp_defines.h:16
Rfun
@ Rfun
Definition pomp_defines.h:16
native
@ native
Definition pomp_defines.h:16
undef
@ undef
Definition pomp_defines.h:16
regNative
@ regNative
Definition pomp_defines.h:16
invalid_names
static R_INLINE int invalid_names(SEXP names)
Definition pomp_defines.h:55
as_matrix
static R_INLINE SEXP as_matrix(SEXP x)
Definition pomp_defines.h:145
ret_array
static R_INLINE SEXP ret_array(int m, int n, SEXP names)
Definition rinit.c:79
paste0
static R_INLINE SEXP paste0(SEXP a, SEXP b, SEXP c)
Definition rinit.c:10
pomp_default_rinit
static SEXP pomp_default_rinit(SEXP params, SEXP Pnames, int npar, int nrep, int nsim)
Definition rinit.c:254
eval_call
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t0, double *p, int npar, double *c, int ncov)
Definition rinit.c:59
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Pnames, SEXP Cnames)
Definition rinit.c:21
lookup_table_t
Definition pomp_defines.h:20
Snames
SEXP Snames
Definition trajectory.c:140
ncovars
int ncovars
Definition trajectory.c:133
params
SEXP params
Definition trajectory.c:128
mode
pompfunmode mode
Definition trajectory.c:126
cov
SEXP cov
Definition trajectory.c:129
fn
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}
warn
#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}
X
#define X
Definition gompertz.c:14
fixdimnames
static R_INLINE void fixdimnames(SEXP x, const char **names, int n)
Definition pomp_defines.h:129
fillrownames
static R_INLINE void fillrownames(SEXP x, SEXP names)
Definition pomp_defines.h:87
makearray
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: