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

Go to the source code of this file.

Functions

static R_INLINE SEXP paste0 (SEXP a, SEXP b)
 
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 *t1, double *t2, double *x1, double *x2, int nvar, double *p, int npar, double *c, int ncov)
 
static R_INLINE SEXP ret_array (int nreps, int ntimes)
 
static SEXP onestep_density (SEXP func, SEXP x, SEXP times, SEXP params, SEXP covar, SEXP log, SEXP args, SEXP gnsi)
 
SEXP do_dprocess (SEXP object, SEXP x, SEXP times, SEXP params, SEXP log, SEXP gnsi)
 

Function Documentation

◆ add_args()

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

Definition at line 18 of file dprocess.c.

19{
20
21 SEXP S1names, S2names;
22 SEXP s1, s2;
23 SEXP var;
24 int v;
25
26 PROTECT(s1 = mkString("_1"));
27 PROTECT(s2 = mkString("_2"));
28 PROTECT(S1names = paste0(Snames,s1));
29 PROTECT(S2names = paste0(Snames,s2));
30
31 PROTECT(args = VectorToPairList(args));
32
33 // Covariates
34 for (v = LENGTH(Cnames)-1; v >= 0; v--) {
35 var = NEW_NUMERIC(1);
36 args = LCONS(var,args);
37 UNPROTECT(1);
38 PROTECT(args);
39 SET_TAG(args,installChar(STRING_ELT(Cnames,v)));
40 }
41
42 // Parameters
43 for (v = LENGTH(Pnames)-1; v >= 0; v--) {
44 var = NEW_NUMERIC(1);
45 args = LCONS(var,args);
46 UNPROTECT(1);
47 PROTECT(args);
48 SET_TAG(args,installChar(STRING_ELT(Pnames,v)));
49 }
50
51 // Latent state variables
52 for (v = LENGTH(Snames)-1; v >= 0; v--) {
53
54 var = NEW_NUMERIC(1);
55 args = LCONS(var,args);
56 UNPROTECT(1);
57 PROTECT(args);
58 SET_TAG(args,installChar(STRING_ELT(S2names,v)));
59
60 var = NEW_NUMERIC(1);
61 args = LCONS(var,args);
62 UNPROTECT(1);
63 PROTECT(args);
64 SET_TAG(args,installChar(STRING_ELT(S1names,v)));
65
66 }
67
68 // Time
69 var = NEW_NUMERIC(1);
70 args = LCONS(var,args);
71 UNPROTECT(1);
72 PROTECT(args);
73 SET_TAG(args,install("t_2"));
74
75 var = NEW_NUMERIC(1);
76 args = LCONS(var,args);
77 UNPROTECT(1);
78 PROTECT(args);
79 SET_TAG(args,install("t_1"));
80
81 UNPROTECT(5);
82 return args;
83
84}
paste0
static R_INLINE SEXP paste0(SEXP a, SEXP b)
Definition dprocess.c:10
Snames
SEXP Snames
Definition trajectory.c:140
args
SEXP args
Definition trajectory.c:139
Here is the call graph for this function:
Here is the caller graph for this function:

◆ do_dprocess()

SEXP do_dprocess ( SEXP  object,
SEXP  x,
SEXP  times,
SEXP  params,
SEXP  log,
SEXP  gnsi 
)

Definition at line 273 of file dprocess.c.

274{
275 SEXP X, fn, args, covar;
276 PROTECT(times=AS_NUMERIC(times));
277 PROTECT(x = as_state_array(x));
278 PROTECT(params = as_matrix(params));
279 // extract the process function
280 PROTECT(fn = GET_SLOT(object,install("dprocess")));
281 // extract other arguments
282 PROTECT(args = GET_SLOT(object,install("userdata")));
283 PROTECT(covar = GET_SLOT(object,install("covar")));
284 // evaluate the density
285 PROTECT(X = onestep_density(fn,x,times,params,covar,log,args,gnsi));
286 UNPROTECT(7);
287 return X;
288}
onestep_density
static SEXP onestep_density(SEXP func, SEXP x, SEXP times, SEXP params, SEXP covar, SEXP log, SEXP args, SEXP gnsi)
Definition dprocess.c:128
X
#define X
Definition gompertz.c:14
as_state_array
static R_INLINE SEXP as_state_array(SEXP x)
Definition pomp_defines.h:183
as_matrix
static R_INLINE SEXP as_matrix(SEXP x)
Definition pomp_defines.h:145
params
SEXP params
Definition trajectory.c:128
fn
SEXP fn
Definition trajectory.c:138
Here is the call graph for this function:

◆ eval_call()

static R_INLINE SEXP eval_call ( SEXP  fn,
SEXP  args,
double *  t1,
double *  t2,
double *  x1,
double *  x2,
int  nvar,
double *  p,
int  npar,
double *  c,
int  ncov 
)
static

Definition at line 86 of file dprocess.c.

93 {
94
95 SEXP var = args, ans, ob;
96 int v;
97
98 *(REAL(CAR(var))) = *t1; var = CDR(var);
99 *(REAL(CAR(var))) = *t2; var = CDR(var);
100 for (v = 0; v < nvar; v++, x1++, x2++) {
101 *(REAL(CAR(var))) = *x1; var = CDR(var);
102 *(REAL(CAR(var))) = *x2; var = CDR(var);
103 }
104 for (v = 0; v < npar; v++, p++, var=CDR(var)) *(REAL(CAR(var))) = *p;
105 for (v = 0; v < ncov; v++, c++, var=CDR(var)) *(REAL(CAR(var))) = *c;
106
107 PROTECT(ob = LCONS(fn,args));
108 PROTECT(ans = eval(ob,CLOENV(fn)));
109
110 UNPROTECT(2);
111 return ans;
112
113}
Here is the caller graph for this function:

◆ onestep_density()

static SEXP onestep_density ( SEXP  func,
SEXP  x,
SEXP  times,
SEXP  params,
SEXP  covar,
SEXP  log,
SEXP  args,
SEXP  gnsi 
)
static

Definition at line 127 of file dprocess.c.

131 {
132
133 pompfunmode mode = undef;
134 int give_log;
135 int nvars, npars, nrepsx, nrepsp, nreps, ntimes, ncovars;
136 SEXP Snames, Pnames, Cnames;
137 SEXP fn, ans;
138 SEXP F, cvec;
139 double *cov;
140 int *dim;
141
142 ntimes = LENGTH(times);
143 dim = INTEGER(GET_DIM(x)); nvars = dim[0]; nrepsx = dim[1];
144 if (ntimes < 2)
145 err("length(times) < 2: with no transitions, there is no work to do.");
146 if (ntimes != dim[2])
147 err("the length of 'times' and 3rd dimension of 'x' do not agree.");
148 dim = INTEGER(GET_DIM(params)); npars = dim[0]; nrepsp = dim[1];
149
150 give_log = *(INTEGER(AS_INTEGER(log)));
151
152 // handle case with different numbers of states and parameters
153 if (nrepsx != nrepsp && nrepsx % nrepsp != 0 && nrepsp % nrepsx != 0) {
154 err("the larger number of replicates is not a multiple of smaller.");
155 } else {
156 nreps = (nrepsx > nrepsp) ? nrepsx : nrepsp;
157 }
158
159 PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x)));
160 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
161 PROTECT(Cnames = get_covariate_names(covar));
162
163 PROTECT(F = ret_array(nreps,ntimes-1));
164
165 // set up the covariate table
166 lookup_table_t covariate_table = make_covariate_table(covar,&ncovars);
167 PROTECT(cvec = NEW_NUMERIC(ncovars));
168 cov = REAL(cvec);
169
170 PROTECT(fn = pomp_fun_handler(func,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
171
172 int nprotect = 6;
173 double *t1 = REAL(times), *t2 = REAL(times)+1;
174 double *x1p = REAL(x);
175 double *x2p = REAL(x)+nrepsx*nvars;
176 double *ft = REAL(F);
177
178 switch (mode) {
179
180 case Rfun: {
181
182 PROTECT(args = add_args(args,Snames,Pnames,Cnames)); nprotect++;
183
184 for (int k = 0; k < ntimes-1; k++, t1++, t2++) { // loop over times
185
186 R_CheckUserInterrupt();
187
188 // interpolate the covariates at time t1
189 table_lookup(&covariate_table,*t1,cov);
190
191 for (int j = 0; j < nreps; j++, ft++) {
192
193 double *p = REAL(params)+npars*(j%nrepsp);
194 double *x1 = x1p+nvars*(j%nrepsx);
195 double *x2 = x2p+nvars*(j%nrepsx);
196
197 PROTECT(ans = eval_call(fn,args,t1,t2,x1,x2,nvars,p,npars,cov,ncovars));
198 *ft = *REAL(AS_NUMERIC(ans));
199 UNPROTECT(1);
200
201 if (!give_log) *ft = exp(*ft);
202
203 }
204
205 x1p = x2p;
206 x2p += nrepsx*nvars;
207
208 }
209
210 }
211
212 break;
213
214 case native: case regNative: {
215
216 int *sidx, *pidx, *cidx;
217 pomp_dprocess *ff = NULL;
218
219 sidx = INTEGER(GET_SLOT(func,install("stateindex")));
220 pidx = INTEGER(GET_SLOT(func,install("paramindex")));
221 cidx = INTEGER(GET_SLOT(func,install("covarindex")));
222
223 *((void **) (&ff)) = R_ExternalPtrAddr(fn);
224
225 for (int k = 0; k < ntimes-1; k++, t1++, t2++) {
226
227 R_CheckUserInterrupt();
228
229 // interpolate the covariates at time t1
230 table_lookup(&covariate_table,*t1,cov);
231
232 for (int j = 0; j < nreps; j++, ft++) {
233
234 double *p = REAL(params)+npars*(j%nrepsp);
235 double *x1 = x1p+nvars*(j%nrepsx);
236 double *x2 = x2p+nvars*(j%nrepsx);
237
238 (*ff)(ft,x1,x2,*t1,*t2,p,sidx,pidx,cidx,cov);
239
240 if (!give_log) *ft = exp(*ft);
241
242 }
243
244 x1p = x2p;
245 x2p += nrepsx*nvars;
246
247 }
248
249 }
250
251 break;
252
253 default: {
254 double *ft = REAL(F);
255 int j, k;
256
257 for (k = 0; k < ntimes-1; k++) { // loop over times
258 for (j = 0; j < nreps; j++, ft++) { // loop over replicates
259 *ft = R_NaReal;
260 }
261 }
262
263 warn("'dprocess' unspecified: likelihood undefined.");
264
265 }
266
267 }
268
269 UNPROTECT(nprotect);
270 return F;
271}
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 *t1, double *t2, double *x1, double *x2, int nvar, double *p, int npar, double *c, int ncov)
Definition dprocess.c:87
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition dprocess.c:18
ret_array
static R_INLINE SEXP ret_array(int nreps, int ntimes)
Definition dprocess.c:115
warn
#define warn(...)
Definition pomp.h:22
pomp_dprocess
void pomp_dprocess(double *loglik, const double *x1, const double *x2, double t1, double t2, const double *p, const int *stateindex, const int *parindex, const int *covindex, const double *covars)
Definition pomp.h:66
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
lookup_table_t
Definition pomp_defines.h:20
npars
int npars
Definition trajectory.c:132
ncovars
int ncovars
Definition trajectory.c:133
mode
pompfunmode mode
Definition trajectory.c:126
nreps
int nreps
Definition trajectory.c:134
cov
SEXP cov
Definition trajectory.c:129
nvars
int nvars
Definition trajectory.c:131
Here is the call graph for this function:
Here is the caller graph for this function:

◆ paste0()

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

Definition at line 10 of file dprocess.c.

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

◆ ret_array()

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

Definition at line 115 of file dprocess.c.

116{
117 int dim[2] = {nreps, ntimes};
118 const char *dimnm[2] = {".id","time"};
119 SEXP F;
120 PROTECT(F = makearray(2,dim));
121 fixdimnames(F,dimnm,2);
122 UNPROTECT(1);
123 return F;
124}
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: