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)
Definition: pomp_defines.h:145
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 *)
Definition: lookup_table.c:11
SEXP pomp_fun_handler(SEXP, SEXP, pompfunmode *, SEXP, SEXP, SEXP, SEXP)
Definition: pomp_fun.c:30
void table_lookup(lookup_table_t *, double, double *)
Definition: lookup_table.c:53
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
Definition: pomp_defines.h:16
@ Rfun
Definition: pomp_defines.h:16
@ native
Definition: pomp_defines.h:16
@ undef
Definition: pomp_defines.h:16
@ regNative
Definition: pomp_defines.h:16
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)
Definition: pomp_defines.h:129
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: