pomp
Inference for partially observed Markov processes
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
dinit.c File Reference
#include <R.h>
#include <Rmath.h>
#include <Rdefines.h>
#include "internal.h"
Include dependency graph for dinit.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 *t0, double *x, int nvar, double *p, int npar, double *c, int ncov)
 
static R_INLINE SEXP ret_array (int nreps)
 
static SEXP init_density (SEXP func, SEXP X, SEXP t0, SEXP params, SEXP covar, SEXP log, SEXP args, SEXP gnsi)
 
SEXP do_dinit (SEXP object, SEXP t0, SEXP x, 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 10 of file dinit.c.

13 {
14
15 SEXP var;
16 int v;
17
18 PROTECT(args = VectorToPairList(args));
19
20 // Covariates
21 for (v = LENGTH(Cnames)-1; v >= 0; v--) {
22 var = NEW_NUMERIC(1);
23 args = LCONS(var,args);
24 UNPROTECT(1);
25 PROTECT(args);
26 SET_TAG(args,installChar(STRING_ELT(Cnames,v)));
27 }
28
29 // Parameters
30 for (v = LENGTH(Pnames)-1; v >= 0; v--) {
31 var = NEW_NUMERIC(1);
32 args = LCONS(var,args);
33 UNPROTECT(1);
34 PROTECT(args);
35 SET_TAG(args,installChar(STRING_ELT(Pnames,v)));
36 }
37
38 // Latent state variables
39 for (v = LENGTH(Snames)-1; v >= 0; v--) {
40 var = NEW_NUMERIC(1);
41 args = LCONS(var,args);
42 UNPROTECT(1);
43 PROTECT(args);
44 SET_TAG(args,installChar(STRING_ELT(Snames,v)));
45 }
46
47 // Time
48 var = NEW_NUMERIC(1);
49 args = LCONS(var,args);
50 UNPROTECT(1);
51 PROTECT(args);
52 SET_TAG(args,install("t0"));
53
54 UNPROTECT(1);
55 return args;
56
57}
SEXP Snames
Definition trajectory.c:140
SEXP args
Definition trajectory.c:139
Here is the caller graph for this function:

◆ do_dinit()

SEXP do_dinit ( SEXP  object,
SEXP  t0,
SEXP  x,
SEXP  params,
SEXP  log,
SEXP  gnsi 
)

Definition at line 220 of file dinit.c.

223 {
224 SEXP F, fn, args, covar;
225 PROTECT(t0=AS_NUMERIC(t0));
226 PROTECT(x = as_matrix(x));
227 PROTECT(params = as_matrix(params));
228 // extract the process function
229 PROTECT(fn = GET_SLOT(object,install("dinit")));
230 // extract other arguments
231 PROTECT(args = GET_SLOT(object,install("userdata")));
232 PROTECT(covar = GET_SLOT(object,install("covar")));
233 // evaluate the density
234 PROTECT(F = init_density(fn,x,t0,params,covar,log,args,gnsi));
235 UNPROTECT(7);
236 return F;
237}
static SEXP init_density(SEXP func, SEXP X, SEXP t0, SEXP params, SEXP covar, SEXP log, SEXP args, SEXP gnsi)
Definition dinit.c:98
static R_INLINE SEXP as_matrix(SEXP x)
SEXP params
Definition trajectory.c:128
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 *  t0,
double *  x,
int  nvar,
double *  p,
int  npar,
double *  c,
int  ncov 
)
static

Definition at line 59 of file dinit.c.

66 {
67
68 SEXP var = args, ans, ob;
69 int v;
70
71 *(REAL(CAR(var))) = *t0; var = CDR(var);
72 for (v = 0; v < nvar; v++, x++) {
73 *(REAL(CAR(var))) = *x; var = CDR(var);
74 }
75 for (v = 0; v < npar; v++, p++, var=CDR(var)) *(REAL(CAR(var))) = *p;
76 for (v = 0; v < ncov; v++, c++, var=CDR(var)) *(REAL(CAR(var))) = *c;
77
78 PROTECT(ob = LCONS(fn,args));
79 PROTECT(ans = eval(ob,CLOENV(fn)));
80
81 UNPROTECT(2);
82 return ans;
83
84}
Here is the caller graph for this function:

◆ init_density()

static SEXP init_density ( SEXP  func,
SEXP  X,
SEXP  t0,
SEXP  params,
SEXP  covar,
SEXP  log,
SEXP  args,
SEXP  gnsi 
)
static

Definition at line 97 of file dinit.c.

101 {
102
104 int give_log;
105 int nvars, npars, nrepsx, nrepsp, nreps, ncovars;
106 SEXP Snames, Pnames, Cnames;
107 SEXP fn, ans;
108 SEXP F, cvec;
109 double *cov;
110 int *dim;
111
112 dim = INTEGER(GET_DIM(X)); nvars = dim[0]; nrepsx = dim[1];
113 dim = INTEGER(GET_DIM(params)); npars = dim[0]; nrepsp = dim[1];
114
115 give_log = *(INTEGER(AS_INTEGER(log)));
116
117 // handle case with different numbers of states and parameters
118 if (nrepsx != nrepsp && nrepsx % nrepsp != 0 && nrepsp % nrepsx != 0) {
119 err("the larger number of replicates is not a multiple of smaller.");
120 } else {
121 nreps = (nrepsx > nrepsp) ? nrepsx : nrepsp;
122 }
123
124 PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(X)));
125 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
126 PROTECT(Cnames = get_covariate_names(covar));
127 PROTECT(F = ret_array(nreps));
128
129 // set up the covariate table
130 lookup_table_t covariate_table = make_covariate_table(covar,&ncovars);
131 PROTECT(cvec = NEW_NUMERIC(ncovars));
132 cov = REAL(cvec);
133
134 // process the pomp_fun
135 PROTECT(fn = pomp_fun_handler(func,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
136
137 int nprotect = 6;
138
139 switch (mode) {
140
141 case Rfun: {
142
143 double *t = REAL(t0);
144 double *ft = REAL(F);
145
146 PROTECT(args = add_args(args,Snames,Pnames,Cnames)); nprotect++;
147
148 // interpolate the covariates at time t1
149 table_lookup(&covariate_table,*t,cov);
150
151 for (int j = 0; j < nreps; j++, ft++) {
152
153 double *xs = REAL(X)+nvars*(j%nrepsx);
154 double *ps = REAL(params)+npars*(j%nrepsp);
155
156 PROTECT(ans = eval_call(fn,args,t,xs,nvars,ps,npars,cov,ncovars));
157
158 *ft = *REAL(AS_NUMERIC(ans));
159
160 UNPROTECT(1);
161
162 if (!give_log) *ft = exp(*ft);
163
164 }
165
166
167 }
168
169 break;
170
171 case native: case regNative: {
172
173 int *sidx, *pidx, *cidx;
174 double *t = REAL(t0);
175 double *ft = REAL(F);
176 pomp_dinit *ff = NULL;
177
178 sidx = INTEGER(GET_SLOT(func,install("stateindex")));
179 pidx = INTEGER(GET_SLOT(func,install("paramindex")));
180 cidx = INTEGER(GET_SLOT(func,install("covarindex")));
181
182 *((void **) (&ff)) = R_ExternalPtrAddr(fn);
183
184 // interpolate the covariates
185 table_lookup(&covariate_table,*t,cov);
186
187 for (int j = 0; j < nreps; j++, ft++) {
188
189 double *xs = REAL(X)+nvars*(j%nrepsx);
190 double *ps = REAL(params)+npars*(j%nrepsp);
191
192 (*ff)(ft,xs,ps,*t,sidx,pidx,cidx,cov);
193
194 if (!give_log) *ft = exp(*ft);
195
196 }
197
198 }
199
200 break;
201
202 default: {
203 double *ft = REAL(F);
204 int j;
205
206 for (j = 0; j < nreps; j++, ft++) { // loop over replicates
207 *ft = R_NaReal;
208 }
209
210 warn("'dinit' unspecified: likelihood undefined.");
211
212 }
213
214 }
215
216 UNPROTECT(nprotect);
217 return F;
218}
lookup_table_t make_covariate_table(SEXP, int *)
SEXP pomp_fun_handler(SEXP, SEXP, pompfunmode *, SEXP, SEXP, SEXP, SEXP)
Definition pomp_fun.c:30
void table_lookup(lookup_table_t *, double, double *)
SEXP get_covariate_names(SEXP)
Definition lookup_table.c:7
static R_INLINE SEXP ret_array(int nreps)
Definition dinit.c:86
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t0, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition dinit.c:60
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition dinit.c:11
#define X
Definition gompertz.c:14
void pomp_dinit(double *lik, const double *x, const double *p, double t0, const int *stateindex, const int *parindex, const int *covindex, const double *covars)
Definition pomp.h:46
#define warn(...)
Definition pomp.h:22
#define err(...)
Definition pomp.h:21
pompfunmode
@ Rfun
@ native
@ undef
@ regNative
int npars
Definition trajectory.c:132
int ncovars
Definition trajectory.c:133
pompfunmode mode
Definition trajectory.c:126
int nreps
Definition trajectory.c:134
SEXP cov
Definition trajectory.c:129
int nvars
Definition trajectory.c:131
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  nreps)
static

Definition at line 86 of file dinit.c.

87{
88 int dim[1] = {nreps};
89 const char *dimnm[1] = {".id"};
90 SEXP F;
91 PROTECT(F = makearray(1,dim));
92 fixdimnames(F,dimnm,1);
93 UNPROTECT(1);
94 return F;
95}
static R_INLINE void fixdimnames(SEXP x, const char **names, int n)
static R_INLINE SEXP makearray(int rank, const int *dim)
Here is the call graph for this function:
Here is the caller graph for this function: