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

Go to the source code of this file.

Functions

static R_INLINE SEXP add_args (SEXP args, SEXP Onames, SEXP Snames, SEXP Pnames, SEXP Cnames, SEXP log)
 
static R_INLINE SEXP eval_call (SEXP fn, SEXP args, double *t, double *y, int nobs, double *x, int nvar, double *p, int npar, double *c, int ncov)
 
static R_INLINE SEXP ret_array (int nreps, int ntimes)
 
SEXP do_dmeasure (SEXP object, SEXP y, SEXP x, SEXP times, SEXP params, SEXP log, SEXP gnsi)
 

Function Documentation

◆ add_args()

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

Definition at line 10 of file dmeasure.c.

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

◆ do_dmeasure()

SEXP do_dmeasure ( SEXP  object,
SEXP  y,
SEXP  x,
SEXP  times,
SEXP  params,
SEXP  log,
SEXP  gnsi 
)

Definition at line 111 of file dmeasure.c.

112{
113
114 pompfunmode mode = undef;
115 int ntimes, nvars, npars, ncovars, nreps, nrepsx, nrepsp, nobs;
116 SEXP Snames, Pnames, Cnames, Onames;
117 SEXP cvec, pompfun;
118 SEXP fn, args, ans;
119 SEXP F;
120 int *dim;
121 lookup_table_t covariate_table;
122 double *cov;
123
124 PROTECT(times = AS_NUMERIC(times));
125 ntimes = length(times);
126 if (ntimes < 1) err("length('times') = 0, no work to do.");
127
128 PROTECT(y = as_matrix(y));
129 dim = INTEGER(GET_DIM(y));
130 nobs = dim[0];
131
132 if (ntimes != dim[1])
133 err("length of 'times' and 2nd dimension of 'y' do not agree.");
134
135 PROTECT(x = as_state_array(x));
136 dim = INTEGER(GET_DIM(x));
137 nvars = dim[0]; nrepsx = dim[1];
138
139 if (ntimes != dim[2])
140 err("length of 'times' and 3rd dimension of 'x' do not agree.");
141
142 PROTECT(params = as_matrix(params));
143 dim = INTEGER(GET_DIM(params));
144 npars = dim[0]; nrepsp = dim[1];
145
146 nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
147
148 if ((nreps % nrepsp != 0) || (nreps % nrepsx != 0))
149 err("larger number of replicates is not a multiple of smaller.");
150
151 PROTECT(Onames = GET_ROWNAMES(GET_DIMNAMES(y)));
152 PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x)));
153 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
154 PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
155
156 // set up the covariate table
157 covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
158 PROTECT(cvec = NEW_NUMERIC(ncovars));
159 cov = REAL(cvec);
160
161 // extract the user-defined function
162 PROTECT(pompfun = GET_SLOT(object,install("dmeasure")));
163 PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,Onames,Cnames));
164
165 // extract 'userdata' as pairlist
166 PROTECT(args = GET_SLOT(object,install("userdata")));
167
168 // create array to store results
169 PROTECT(F = ret_array(nreps,ntimes));
170
171 int nprotect = 13;
172
173 switch (mode) {
174
175 case Rfun: {
176
177 double *ys = REAL(y), *xs = REAL(x), *ps = REAL(params), *time = REAL(times);
178 double *ft = REAL(F);
179 int j, k;
180
181 // build argument list
182 PROTECT(args = add_args(args,Onames,Snames,Pnames,Cnames,log)); nprotect++;
183
184 for (k = 0; k < ntimes; k++, time++, ys += nobs) { // loop over times
185
186 R_CheckUserInterrupt(); // check for user interrupt
187
188 table_lookup(&covariate_table,*time,cov); // interpolate the covariates
189
190 for (j = 0; j < nreps; j++, ft++) { // loop over replicates
191
192 // evaluate the call
193 PROTECT(
194 ans = eval_call(
195 fn,args,
196 time,
197 ys,nobs,
198 xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
199 ps+npars*(j%nrepsp),npars,
200 cov,ncovars
201 )
202 );
203
204 if (k == 0 && j == 0 && LENGTH(ans) != 1)
205 err("user 'dmeasure' returns a vector of length %d when it should return a scalar.",LENGTH(ans));
206
207 *ft = *(REAL(AS_NUMERIC(ans)));
208
209 UNPROTECT(1);
210
211 }
212 }
213 }
214
215 break;
216
217 case native: case regNative: {
218 int *oidx, *sidx, *pidx, *cidx;
219 int give_log;
220 pomp_dmeasure *ff = NULL;
221 double *yp = REAL(y), *xs = REAL(x), *ps = REAL(params), *time = REAL(times);
222 double *ft = REAL(F);
223 double *xp, *pp;
224 int j, k;
225
226 // extract state, parameter, covariate, observable indices
227 sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
228 pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
229 oidx = INTEGER(GET_SLOT(pompfun,install("obsindex")));
230 cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
231
232 give_log = *(INTEGER(AS_INTEGER(log)));
233
234 // address of native routine
235 *((void **) (&ff)) = R_ExternalPtrAddr(fn);
236
237 for (k = 0; k < ntimes; k++, time++, yp += nobs) { // loop over times
238
239 R_CheckUserInterrupt(); // check for user interrupt
240
241 // interpolate the covar functions for the covariates
242 table_lookup(&covariate_table,*time,cov);
243
244 for (j = 0; j < nreps; j++, ft++) { // loop over replicates
245
246 xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
247 pp = &ps[npars*(j%nrepsp)];
248
249 (*ff)(ft,yp,xp,pp,give_log,oidx,sidx,pidx,cidx,cov,*time);
250
251 }
252 }
253
254 }
255
256 break;
257
258 default: {
259 double *ft = REAL(F);
260 int j, k;
261
262 for (k = 0; k < ntimes; k++) { // loop over times
263 for (j = 0; j < nreps; j++, ft++) { // loop over replicates
264 *ft = R_NaReal;
265 }
266 }
267
268 warn("'dmeasure' unspecified: likelihood undefined.");
269
270 }
271
272 }
273
274 UNPROTECT(nprotect);
275 return F;
276}
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
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Onames, SEXP Snames, SEXP Pnames, SEXP Cnames, SEXP log)
Definition dmeasure.c:11
eval_call
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *y, int nobs, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition dmeasure.c:75
ret_array
static R_INLINE SEXP ret_array(int nreps, int ntimes)
Definition dmeasure.c:101
pomp_dmeasure
void pomp_dmeasure(double *lik, const double *y, const double *x, const double *p, int give_log, const int *obsindex, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition pomp.h:86
warn
#define warn(...)
Definition pomp.h:22
err
#define err(...)
Definition pomp.h:21
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
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 *  y,
int  nobs,
double *  x,
int  nvar,
double *  p,
int  npar,
double *  c,
int  ncov 
)
static

Definition at line 75 of file dmeasure.c.

82{
83
84 SEXP var = args, ans, ob;
85 int v;
86
87 *(REAL(CAR(var))) = *t; var = CDR(var);
88 for (v = 0; v < nobs; v++, y++, var=CDR(var)) *(REAL(CAR(var))) = *y;
89 for (v = 0; v < nvar; v++, x++, var=CDR(var)) *(REAL(CAR(var))) = *x;
90 for (v = 0; v < npar; v++, p++, var=CDR(var)) *(REAL(CAR(var))) = *p;
91 for (v = 0; v < ncov; v++, c++, var=CDR(var)) *(REAL(CAR(var))) = *c;
92
93 PROTECT(ob = LCONS(fn,args));
94 PROTECT(ans = eval(ob,CLOENV(fn)));
95
96 UNPROTECT(2);
97 return ans;
98
99}
Here is the caller graph for this function:

◆ ret_array()

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

Definition at line 101 of file dmeasure.c.

101 {
102 int dim[2] = {nreps, ntimes};
103 const char *dimnm[2] = {".id","time"};
104 SEXP F;
105 PROTECT(F = makearray(2,dim));
106 fixdimnames(F,dimnm,2);
107 UNPROTECT(1);
108 return F;
109}
fixdimnames
static R_INLINE void fixdimnames(SEXP x, const char **names, int n)
Definition pomp_defines.h:129
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: