pomp
Inference for partially observed Markov processes
•All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
vmeasure.c File Reference
#include <R.h>
#include <Rmath.h>
#include <Rdefines.h>
#include "internal.h"
Include dependency graph for vmeasure.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_vmeasure (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 vmeasure.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_vmeasure()

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

Definition at line 99 of file vmeasure.c.

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

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