pomp
Inference for partially observed Markov processes
•All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
emeasure.c File Reference
#include <R.h>
#include <Rmath.h>
#include <Rdefines.h>
#include "internal.h"
Include dependency graph for emeasure.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 *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 do_emeasure (SEXP object, SEXP x, SEXP times, SEXP params, SEXP gnsi)
 

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 emeasure.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 // covariates, parameter, states, then time
20
21 // Covariates
22 for (v = LENGTH(Cnames)-1; v >= 0; v--) {
23 var = NEW_NUMERIC(1);
24 args = LCONS(var,args);
25 UNPROTECT(1);
26 PROTECT(args);
27 SET_TAG(args,installChar(STRING_ELT(Cnames,v)));
28 }
29
30 // Parameters
31 for (v = LENGTH(Pnames)-1; v >= 0; v--) {
32 var = NEW_NUMERIC(1);
33 args = LCONS(var,args);
34 UNPROTECT(1);
35 PROTECT(args);
36 SET_TAG(args,installChar(STRING_ELT(Pnames,v)));
37 }
38
39 // Latent state variables
40 for (v = LENGTH(Snames)-1; v >= 0; v--) {
41 var = NEW_NUMERIC(1);
42 args = LCONS(var,args);
43 UNPROTECT(1);
44 PROTECT(args);
45 SET_TAG(args,installChar(STRING_ELT(Snames,v)));
46 }
47
48 // Time
49 var = NEW_NUMERIC(1);
50 args = LCONS(var,args);
51 UNPROTECT(1);
52 PROTECT(args);
53 SET_TAG(args,install("t"));
54
55 UNPROTECT(1);
56 return args;
57
58}
Snames
SEXP Snames
Definition trajectory.c:140
args
SEXP args
Definition trajectory.c:139
Here is the caller graph for this function:

◆ do_emeasure()

SEXP do_emeasure ( SEXP  object,
SEXP  x,
SEXP  times,
SEXP  params,
SEXP  gnsi 
)

Definition at line 98 of file emeasure.c.

99{
100 pompfunmode mode = undef;
101 int ntimes, nvars, npars, ncovars, nreps, nrepsx, nrepsp;
102 int nobs = 0;
103 SEXP Snames, Pnames, Cnames, Onames = R_NilValue;
104 SEXP fn, args;
105 SEXP pompfun;
106 SEXP Y = R_NilValue;
107 int *dim;
108 lookup_table_t covariate_table;
109 SEXP cvec;
110 double *cov;
111
112 PROTECT(times = AS_NUMERIC(times));
113 ntimes = length(times);
114 if (ntimes < 1)
115 err("length('times') = 0, no work to do.");
116
117 PROTECT(x = as_state_array(x));
118 dim = INTEGER(GET_DIM(x));
119 nvars = dim[0]; nrepsx = dim[1];
120
121 if (ntimes != dim[2])
122 err("length of 'times' and 3rd dimension of 'x' do not agree.");
123
124 PROTECT(params = as_matrix(params));
125 dim = INTEGER(GET_DIM(params));
126 npars = dim[0]; nrepsp = dim[1];
127
128 nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
129
130 if ((nreps % nrepsp != 0) || (nreps % nrepsx != 0))
131 err("larger number of replicates is not a multiple of smaller.");
132
133 PROTECT(pompfun = GET_SLOT(object,install("emeasure")));
134
135 PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x)));
136 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
137 PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
138 PROTECT(Onames = GET_SLOT(pompfun,install("obsnames")));
139
140 // set up the covariate table
141 covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
142 PROTECT(cvec = NEW_NUMERIC(ncovars));
143 cov = REAL(cvec);
144
145 // extract the user-defined function
146 PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,Onames,Cnames));
147
148 // extract 'userdata' as pairlist
149 PROTECT(args = GET_SLOT(object,install("userdata")));
150
151 int nprotect = 11;
152 int first = 1;
153
154 // first do setup
155 switch (mode) {
156
157 case Rfun: {
158 double *ys, *yt = 0;
159 double *time = REAL(times), *xs = REAL(x), *ps = REAL(params);
160 SEXP ans;
161 int j, k;
162
163 PROTECT(args = add_args(args,Snames,Pnames,Cnames)); nprotect++;
164
165 for (k = 0; k < ntimes; k++, time++) { // loop over times
166
167 R_CheckUserInterrupt(); // check for user interrupt
168
169 table_lookup(&covariate_table,*time,cov); // interpolate the covariates
170
171 for (j = 0; j < nreps; j++) { // loop over replicates
172
173 if (first) {
174
175 PROTECT(
176 ans = eval_call(
177 fn,args,
178 time,
179 xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
180 ps+npars*(j%nrepsp),npars,
181 cov,ncovars
182 )
183 );
184
185 nobs = LENGTH(ans);
186
187 PROTECT(Onames = GET_NAMES(ans));
188 if (invalid_names(Onames))
189 err("'emeasure' must return a named numeric vector.");
190
191 PROTECT(Y = ret_array(nobs,nreps,ntimes,Onames));
192
193 nprotect += 3;
194
195 yt = REAL(Y);
196 ys = REAL(AS_NUMERIC(ans));
197
198 memcpy(yt,ys,nobs*sizeof(double));
199 yt += nobs;
200
201 first = 0;
202
203 } else {
204
205 PROTECT(
206 ans = eval_call(
207 fn,args,
208 time,
209 xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
210 ps+npars*(j%nrepsp),npars,
211 cov,ncovars
212 )
213 );
214
215 if (LENGTH(ans) != nobs)
216 err("'emeasure' returns variable-length results.");
217
218 ys = REAL(AS_NUMERIC(ans));
219
220 memcpy(yt,ys,nobs*sizeof(double));
221 yt += nobs;
222
223 UNPROTECT(1);
224
225 }
226
227 }
228 }
229
230 }
231
232 break;
233
234 case native: case regNative: {
235 double *yt = 0, *xp, *pp;
236 double *time = REAL(times), *xs = REAL(x), *ps = REAL(params);
237 int *oidx, *sidx, *pidx, *cidx;
238 pomp_emeasure *ff = NULL;
239 int j, k;
240
241 nobs = LENGTH(Onames);
242 // extract observable, state, parameter covariate indices
243 sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
244 pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
245 oidx = INTEGER(GET_SLOT(pompfun,install("obsindex")));
246 cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
247
248 // address of native routine
249 *((void **) (&ff)) = R_ExternalPtrAddr(fn);
250
251 PROTECT(Y = ret_array(nobs,nreps,ntimes,Onames)); nprotect++;
252 yt = REAL(Y);
253
254 for (k = 0; k < ntimes; k++, time++) { // loop over times
255
256 R_CheckUserInterrupt(); // check for user interrupt
257
258 // interpolate the covar functions for the covariates
259 table_lookup(&covariate_table,*time,cov);
260
261 for (j = 0; j < nreps; j++, yt += nobs) { // loop over replicates
262
263 xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
264 pp = &ps[npars*(j%nrepsp)];
265
266 (*ff)(yt,xp,pp,oidx,sidx,pidx,cidx,cov,*time);
267
268 }
269 }
270
271 }
272
273 break;
274
275 default: {
276 nobs = LENGTH(Onames);
277 int dim[3] = {nobs, nreps, ntimes};
278 const char *dimnm[3] = {"name",".id","time"};
279 double *yt = 0;
280 int i, n = nobs*nreps*ntimes;
281
282 PROTECT(Y = makearray(3,dim)); nprotect++;
283 setrownames(Y,Onames,3);
284 fixdimnames(Y,dimnm,3);
285
286 for (i = 0, yt = REAL(Y); i < n; i++, yt++) *yt = R_NaReal;
287
288 warn("'emeasure' unspecified: NAs generated.");
289 }
290
291 }
292
293 UNPROTECT(nprotect);
294 return Y;
295}
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
eval_call
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition emeasure.c:60
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition emeasure.c:10
ret_array
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition emeasure.c:84
Y
#define Y
Definition gompertz.c:13
pomp_emeasure
void pomp_emeasure(double *f, const double *x, const double *p, const int *obsindex, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition pomp.h:93
warn
#define warn(...)
Definition pomp.h:22
err
#define err(...)
Definition pomp.h:21
fixdimnames
static R_INLINE void fixdimnames(SEXP x, const char **names, int n)
Definition pomp_defines.h:129
setrownames
static R_INLINE void setrownames(SEXP x, SEXP names, int rank)
Definition pomp_defines.h:110
makearray
static R_INLINE SEXP makearray(int rank, const int *dim)
Definition pomp_defines.h:40
as_state_array
static R_INLINE SEXP as_state_array(SEXP x)
Definition pomp_defines.h:183
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
lookup_table_t
Definition pomp_defines.h:20
npars
int npars
Definition trajectory.c:132
ncovars
int ncovars
Definition trajectory.c:133
params
SEXP params
Definition trajectory.c:128
mode
pompfunmode mode
Definition trajectory.c:126
nreps
int nreps
Definition trajectory.c:134
cov
SEXP cov
Definition trajectory.c:129
fn
SEXP fn
Definition trajectory.c:138
nvars
int nvars
Definition trajectory.c:131
Here is the call graph for this function:

◆ eval_call()

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

Definition at line 60 of file emeasure.c.

66{
67
68 SEXP var = args, ans, ob;
69 int v;
70
71 *(REAL(CAR(var))) = *t; var = CDR(var);
72 for (v = 0; v < nvar; v++, x++, var=CDR(var)) *(REAL(CAR(var))) = *x;
73 for (v = 0; v < npar; v++, p++, var=CDR(var)) *(REAL(CAR(var))) = *p;
74 for (v = 0; v < ncov; v++, c++, var=CDR(var)) *(REAL(CAR(var))) = *c;
75
76 PROTECT(ob = LCONS(fn,args));
77 PROTECT(ans = eval(ob,CLOENV(fn)));
78
79 UNPROTECT(2);
80 return ans;
81
82}
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 84 of file emeasure.c.

84 {
85 int dim[3] = {n, nreps, ntimes};
86 const char *dimnm[3] = {"name", ".id", "time"};
87 SEXP Y;
88
89 PROTECT(Y = makearray(3,dim));
90 setrownames(Y,names,3);
91 fixdimnames(Y,dimnm,3);
92
93 UNPROTECT(1);
94 return Y;
95
96}
Here is the call graph for this function:
Here is the caller graph for this function: