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

Go to the source code of this file.

Macros

#define NVAR   (__ssa_nvar)
 
#define NPAR   (__ssa_npar)
 
#define NCOV   (__ssa_ncov)
 
#define ARGS   (__ssa_args)
 
#define RATEFN   (__ssa_ratefn)
 
#define FIRST   (__ssa_first)
 
#define RXR   (__ssa_rxrate)
 

Functions

static R_INLINE SEXP add_args (SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
 
static void gillespie (double *t, double tmax, double *f, double *y, const double *v, int nvar, int nevent, Rboolean hasvname, const int *ivmat)
 
static void SSA (pomp_ssa_rate_fn *ratefun, int irep, int nvar, int nevent, int npar, int nrep, int ntimes, double *xstart, double t, const double *times, const double *params, double *xout, const double *v, int nzero, const int *izero, const int *istate, const int *ipar, int ncovar, const int *icovar, Rboolean hasvname, const int *ivmat, int mcov, lookup_table_t *tab, const double *hmax)
 
static double __pomp_Rfun_ssa_ratefn (int j, double t, const double *x, const double *p, int *stateindex, int *parindex, int *covindex, double *c)
 
SEXP SSA_simulator (SEXP func, SEXP xstart, SEXP tstart, SEXP times, SEXP params, SEXP vmatrix, SEXP covar, SEXP accumvars, SEXP hmax, SEXP args, SEXP gnsi)
 

Variables

static int __ssa_nvar
 
static int __ssa_npar
 
static int __ssa_ncov
 
static SEXP __ssa_args
 
static SEXP __ssa_ratefn
 
static int __ssa_first
 
static pomp_ssa_rate_fn * __ssa_rxrate
 

Macro Definition Documentation

◆ ARGS

#define ARGS   (__ssa_args)

Definition at line 194 of file ssa.c.

◆ FIRST

#define FIRST   (__ssa_first)

Definition at line 196 of file ssa.c.

◆ NCOV

#define NCOV   (__ssa_ncov)

Definition at line 193 of file ssa.c.

◆ NPAR

#define NPAR   (__ssa_npar)

Definition at line 192 of file ssa.c.

◆ NVAR

#define NVAR   (__ssa_nvar)

Definition at line 191 of file ssa.c.

◆ RATEFN

#define RATEFN   (__ssa_ratefn)

Definition at line 195 of file ssa.c.

◆ RXR

#define RXR   (__ssa_rxrate)

Definition at line 197 of file ssa.c.

Function Documentation

◆ __pomp_Rfun_ssa_ratefn()

static double __pomp_Rfun_ssa_ratefn ( int  j,
double  t,
const double *  x,
const double *  p,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  c 
)
static

Definition at line 199 of file ssa.c.

201{
202 SEXP var = ARGS, ans, ob;
203 int v;
204 double rate;
205
206 *(INTEGER(CAR(var))) = j; var = CDR(var);
207 *(REAL(CAR(var))) = t; var = CDR(var);
208 for (v = 0; v < NVAR; v++, x++, var=CDR(var)) *(REAL(CAR(var))) = *x;
209 for (v = 0; v < NPAR; v++, p++, var=CDR(var)) *(REAL(CAR(var))) = *p;
210 for (v = 0; v < NCOV; v++, c++, var=CDR(var)) *(REAL(CAR(var))) = *c;
211
212 PROTECT(ob = LCONS(RATEFN,ARGS));
213 PROTECT(ans = eval(ob,CLOENV(RATEFN)));
214
215 if (FIRST) {
216 if (LENGTH(ans) != 1)
217 err("'rate.fun' must return a single numeric rate.");
218 FIRST = 0;
219 }
220
221 rate = *(REAL(AS_NUMERIC(ans)));
222 UNPROTECT(2);
223 return rate;
224}
err
#define err(...)
Definition pomp.h:21
FIRST
#define FIRST
Definition ssa.c:196
NPAR
#define NPAR
Definition ssa.c:192
NCOV
#define NCOV
Definition ssa.c:193
NVAR
#define NVAR
Definition ssa.c:191
RATEFN
#define RATEFN
Definition ssa.c:195
ARGS
#define ARGS
Definition ssa.c:194
Here is the caller graph for this function:

◆ add_args()

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

Definition at line 9 of file ssa.c.

10{
11
12 SEXP var;
13 int v;
14
15 PROTECT(args = VectorToPairList(args));
16
17 // we construct the call from end to beginning
18 // covariates, parameter, states, then time and 'j'
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("t"));
53
54 // 'j'
55 var = NEW_INTEGER(1);
56 args = LCONS(var,args);
57 UNPROTECT(1);
58 PROTECT(args);
59 SET_TAG(args,install("j"));
60
61 UNPROTECT(1);
62 return args;
63
64}
Snames
SEXP Snames
Definition trajectory.c:140
args
SEXP args
Definition trajectory.c:139
Here is the caller graph for this function:

◆ gillespie()

static void gillespie ( double *  t,
double  tmax,
double *  f,
double *  y,
const double *  v,
int  nvar,
int  nevent,
Rboolean  hasvname,
const int *  ivmat 
)
static

Definition at line 66 of file ssa.c.

69{
70 double tstep, p;
71 double vv;
72 int i, j;
73
74 // Determine time interval and update time
75 double fsum = 0;
76 for (j = 0; j < nevent; j++) fsum += f[j];
77
78 if (fsum > 0.0) {
79 tstep = exp_rand()/fsum;
80 *t = *t+tstep;
81 } else {
82 *t = tmax;
83 return;
84 }
85
86 if (*t >= tmax) {
87
88 *t = tmax;
89
90 } else {
91
92 // Determine event, update pops & events
93 p = fsum*unif_rand();
94 int jevent = nevent-1;
95 for (j = 0; j < jevent; j++) {
96 if (p > f[j]) {
97 p -= f[j];
98 } else {
99 jevent = j;
100 break;
101 }
102 }
103
104 for (i = 0; i < nvar; i++) {
105 vv = v[i+nvar*jevent];
106 if (vv != 0) {
107 if (hasvname) {
108 y[ivmat[i]] += vv;
109 } else {
110 y[i] += vv;
111 }
112 }
113 }
114
115 }
116}
Here is the caller graph for this function:

◆ SSA()

static void SSA ( pomp_ssa_rate_fn *  ratefun,
int  irep,
int  nvar,
int  nevent,
int  npar,
int  nrep,
int  ntimes,
double *  xstart,
double  t,
const double *  times,
const double *  params,
double *  xout,
const double *  v,
int  nzero,
const int *  izero,
const int *  istate,
const int *  ipar,
int  ncovar,
const int *  icovar,
Rboolean  hasvname,
const int *  ivmat,
int  mcov,
lookup_table_t *  tab,
const double *  hmax 
)
static

Definition at line 118 of file ssa.c.

125{
126 double tmax;
127 double *f = NULL;
128 double *covars = NULL;
129 const double *par;
130 double y[nvar];
131 int i, j;
132
133 if (mcov > 0) covars = (double *) R_Calloc(mcov,double);
134 if (nevent > 0) f = (double *) R_Calloc(nevent,double);
135
136 par = params+npar*irep;
137 // Copy state variables
138 memcpy(y,xstart+nvar*irep,nvar*sizeof(double));
139 // Set appropriate states to zero
140 for (i = 0; i < nzero; i++) y[izero[i]] = 0.0;
141 // Initialize the covariate vector
142 if (mcov > 0) table_lookup(tab,t,covars);
143 // Initialise propensity functions & tree
144 for (j = 0; j < nevent; j++) {
145 f[j] = ratefun(j+1,t,y,par,istate,ipar,icovar,covars);
146 if (f[j] < 0.0)
147 err("'rate.fun' returns a negative rate");
148 }
149
150 int icount = 0;
151 while (icount < ntimes) {
152
153 R_CheckUserInterrupt();
154 tmax = t + *hmax;
155 tmax = (tmax > times[icount]) ? times[icount] : tmax;
156 gillespie(&t,tmax,f,y,v,nvar,nevent,hasvname,ivmat);
157
158 if (mcov > 0) table_lookup(tab,t,covars);
159
160 for (j = 0; j < nevent; j++) {
161 f[j] = ratefun(j+1,t,y,par,istate,ipar,icovar,covars);
162 if (f[j] < 0.0)
163 err("'rate.fun' returns a negative rate");
164 }
165
166 // Record output at required time points
167 if (t >= times[icount]) {
168 memcpy(xout+nvar*(irep+nrep*icount),y,nvar*sizeof(double));
169 // Set appropriate states to zero at time of last observation
170 for (i = 0; i < nzero; i++) y[izero[i]] = 0;
171 icount++;
172 }
173
174 }
175
176 if (mcov > 0) R_Free(covars);
177 if (nevent > 0) R_Free(f);
178
179}
table_lookup
void table_lookup(lookup_table_t *, double, double *)
Definition lookup_table.c:53
gillespie
static void gillespie(double *t, double tmax, double *f, double *y, const double *v, int nvar, int nevent, Rboolean hasvname, const int *ivmat)
Definition ssa.c:66
params
SEXP params
Definition trajectory.c:128
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SSA_simulator()

SEXP SSA_simulator ( SEXP  func,
SEXP  xstart,
SEXP  tstart,
SEXP  times,
SEXP  params,
SEXP  vmatrix,
SEXP  covar,
SEXP  accumvars,
SEXP  hmax,
SEXP  args,
SEXP  gnsi 
)

Definition at line 226 of file ssa.c.

228{
229
230 int *dim, xdim[3];
231 int nvar, nvarv, nevent, npar, nrep, ntimes;
232 SEXP statenames, paramnames, covarnames;
233 int ncovars, covdim;
234 int nzeros = LENGTH(accumvars);
235 pompfunmode mode = undef;
236 SEXP X, zindex, vindex;
237 int *sidx, *pidx, *cidx, *zidx, *vidx;
238 SEXP fn, Snames, Pnames, Cnames, Vnames;
239 lookup_table_t covariate_table;
240 Rboolean hasvnames;
241 double t0;
242
243 dim = INTEGER(GET_DIM(xstart)); nvar = dim[0]; nrep = dim[1];
244 dim = INTEGER(GET_DIM(params)); npar = dim[0];
245 dim = INTEGER(GET_DIM(vmatrix)); nvarv = dim[0]; nevent = dim[1];
246
247 if (nvarv != nvar)
248 err("number of state variables must equal the number of rows in 'v'.");
249
250 ntimes = LENGTH(times);
251
252 PROTECT(tstart = AS_NUMERIC(tstart));
253 PROTECT(times = AS_NUMERIC(times));
254 t0 = *(REAL(tstart));
255 if (t0 > *(REAL(times))) err("'t0' must be no later than 'times[1]'.");
256
257 PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(xstart)));
258 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
259 PROTECT(Cnames = get_covariate_names(covar));
260 PROTECT(Vnames = GET_ROWNAMES(GET_DIMNAMES(vmatrix)));
261
262 covariate_table = make_covariate_table(covar,&covdim);
263
264 PROTECT(statenames = GET_SLOT(func,install("statenames")));
265 PROTECT(paramnames = GET_SLOT(func,install("paramnames")));
266 PROTECT(covarnames = GET_SLOT(func,install("covarnames")));
267
268 ncovars = LENGTH(covarnames);
269 hasvnames = !isNull(Vnames);
270
271 PROTECT(hmax = AS_NUMERIC(hmax));
272
273 PROTECT(fn = pomp_fun_handler(func,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
274
275 int nprotect = 11;
276
277 switch (mode) {
278
279 case undef: default: // #nocov
280
281 err("unrecognized 'mode' %d",mode); // #nocov
282
283 case native: case regNative:
284
285 *((void **) (&(RXR))) = R_ExternalPtrAddr(fn);
286
287 break;
288
289 case Rfun:
290
291 RXR = (pomp_ssa_rate_fn *) __pomp_Rfun_ssa_ratefn;
292 NVAR = nvar;
293 NPAR = npar;
294 NCOV = covdim;
295 PROTECT(ARGS = add_args(args,Snames,Pnames,Cnames));
296 PROTECT(RATEFN = fn);
297 nprotect += 2;
298 FIRST = 1;
299
300 break;
301
302 }
303
304 xdim[0] = nvar; xdim[1] = nrep; xdim[2] = ntimes;
305 PROTECT(X = makearray(3,xdim)); nprotect++;
306 setrownames(X,Snames,3);
307
308 sidx = INTEGER(GET_SLOT(func,install("stateindex")));
309 pidx = INTEGER(GET_SLOT(func,install("paramindex")));
310 cidx = INTEGER(GET_SLOT(func,install("covarindex")));
311
312 if (nzeros>0) {
313 PROTECT(zindex = MATCHROWNAMES(xstart,accumvars,"state variables")); nprotect++;
314 zidx = INTEGER(zindex);
315 } else {
316 zidx = 0;
317 }
318
319 if (hasvnames) {
320 PROTECT(vindex = MATCHROWNAMES(xstart,Vnames,"state variables")); nprotect++;
321 vidx = INTEGER(vindex);
322 } else {
323 vidx = 0;
324 }
325
326 GetRNGstate();
327 {
328 int i;
329 for (i = 0; i < nrep; i++) {
330 SSA(RXR,i,nvar,nevent,npar,nrep,ntimes,
331 REAL(xstart),t0,REAL(times),REAL(params),
332 REAL(X),REAL(vmatrix),
333 nzeros,zidx,sidx,pidx,ncovars,cidx,hasvnames,vidx,
334 covdim,&covariate_table,REAL(hmax));
335 }
336 }
337 PutRNGstate();
338
339 UNPROTECT(nprotect);
340 return X;
341}
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
get_covariate_names
SEXP get_covariate_names(SEXP)
Definition lookup_table.c:7
X
#define X
Definition gompertz.c:14
pomp_ssa_rate_fn
double pomp_ssa_rate_fn(int event, double t, const double *x, const double *p, const int *stateindex, const int *parindex, const int *covindex, const double *covars)
Definition pomp.h:52
MATCHROWNAMES
#define MATCHROWNAMES(X, N, W)
Definition pomp_defines.h:13
setrownames
static R_INLINE void setrownames(SEXP x, SEXP names, int rank)
Definition pomp_defines.h:110
makearray
static R_INLINE SEXP makearray(int rank, const int *dim)
Definition pomp_defines.h:40
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
SSA
static void SSA(pomp_ssa_rate_fn *ratefun, int irep, int nvar, int nevent, int npar, int nrep, int ntimes, double *xstart, double t, const double *times, const double *params, double *xout, const double *v, int nzero, const int *izero, const int *istate, const int *ipar, int ncovar, const int *icovar, Rboolean hasvname, const int *ivmat, int mcov, lookup_table_t *tab, const double *hmax)
Definition ssa.c:118
__pomp_Rfun_ssa_ratefn
static double __pomp_Rfun_ssa_ratefn(int j, double t, const double *x, const double *p, int *stateindex, int *parindex, int *covindex, double *c)
Definition ssa.c:199
RXR
#define RXR
Definition ssa.c:197
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition ssa.c:9
lookup_table_t
Definition pomp_defines.h:20
ncovars
int ncovars
Definition trajectory.c:133
mode
pompfunmode mode
Definition trajectory.c:126
fn
SEXP fn
Definition trajectory.c:138
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ __ssa_args

SEXP __ssa_args
static

Definition at line 186 of file ssa.c.

◆ __ssa_first

int __ssa_first
static

Definition at line 188 of file ssa.c.

◆ __ssa_ncov

int __ssa_ncov
static

Definition at line 185 of file ssa.c.

◆ __ssa_npar

int __ssa_npar
static

Definition at line 184 of file ssa.c.

◆ __ssa_nvar

int __ssa_nvar
static

Definition at line 183 of file ssa.c.

◆ __ssa_ratefn

SEXP __ssa_ratefn
static

Definition at line 187 of file ssa.c.

◆ __ssa_rxrate

pomp_ssa_rate_fn* __ssa_rxrate
static

Definition at line 189 of file ssa.c.