pomp
Inference for partially observed Markov processes
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 196 of file ssa.c.

◆ FIRST

#define FIRST   (__ssa_first)

Definition at line 198 of file ssa.c.

◆ NCOV

#define NCOV   (__ssa_ncov)

Definition at line 195 of file ssa.c.

◆ NPAR

#define NPAR   (__ssa_npar)

Definition at line 194 of file ssa.c.

◆ NVAR

#define NVAR   (__ssa_nvar)

Definition at line 193 of file ssa.c.

◆ RATEFN

#define RATEFN   (__ssa_ratefn)

Definition at line 197 of file ssa.c.

◆ RXR

#define RXR   (__ssa_rxrate)

Definition at line 199 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 201 of file ssa.c.

203 {
204  SEXP var = ARGS, ans, ob;
205  int v;
206  double rate;
207 
208  *(INTEGER(CAR(var))) = j; var = CDR(var);
209  *(REAL(CAR(var))) = t; var = CDR(var);
210  for (v = 0; v < NVAR; v++, x++, var=CDR(var)) *(REAL(CAR(var))) = *x;
211  for (v = 0; v < NPAR; v++, p++, var=CDR(var)) *(REAL(CAR(var))) = *p;
212  for (v = 0; v < NCOV; v++, c++, var=CDR(var)) *(REAL(CAR(var))) = *c;
213 
214  PROTECT(ob = LCONS(RATEFN,ARGS));
215  PROTECT(ans = eval(ob,CLOENV(RATEFN)));
216 
217  if (FIRST) {
218  if (LENGTH(ans) != 1)
219  err("'rate.fun' must return a single numeric rate.");
220  FIRST = 0;
221  }
222 
223  rate = *(REAL(AS_NUMERIC(ans)));
224  UNPROTECT(2);
225  return rate;
226 }
#define err(...)
Definition: pomp.h:21
#define FIRST
Definition: ssa.c:198
#define NPAR
Definition: ssa.c:194
#define NCOV
Definition: ssa.c:195
#define NVAR
Definition: ssa.c:193
#define RATEFN
Definition: ssa.c:197
#define ARGS
Definition: ssa.c:196
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 11 of file ssa.c.

12 {
13 
14  SEXP var;
15  int v;
16 
17  PROTECT(args = VectorToPairList(args));
18 
19  // we construct the call from end to beginning
20  // covariates, parameter, states, then time and 'j'
21 
22  // Covariates
23  for (v = LENGTH(Cnames)-1; v >= 0; v--) {
24  var = NEW_NUMERIC(1);
25  args = LCONS(var,args);
26  UNPROTECT(1);
27  PROTECT(args);
28  SET_TAG(args,installChar(STRING_ELT(Cnames,v)));
29  }
30 
31  // Parameters
32  for (v = LENGTH(Pnames)-1; v >= 0; v--) {
33  var = NEW_NUMERIC(1);
34  args = LCONS(var,args);
35  UNPROTECT(1);
36  PROTECT(args);
37  SET_TAG(args,installChar(STRING_ELT(Pnames,v)));
38  }
39 
40  // Latent state variables
41  for (v = LENGTH(Snames)-1; v >= 0; v--) {
42  var = NEW_NUMERIC(1);
43  args = LCONS(var,args);
44  UNPROTECT(1);
45  PROTECT(args);
46  SET_TAG(args,installChar(STRING_ELT(Snames,v)));
47  }
48 
49  // Time
50  var = NEW_NUMERIC(1);
51  args = LCONS(var,args);
52  UNPROTECT(1);
53  PROTECT(args);
54  SET_TAG(args,install("t"));
55 
56  // 'j'
57  var = NEW_INTEGER(1);
58  args = LCONS(var,args);
59  UNPROTECT(1);
60  PROTECT(args);
61  SET_TAG(args,install("j"));
62 
63  UNPROTECT(1);
64  return args;
65 
66 }
SEXP Snames
Definition: trajectory.c:140
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 68 of file ssa.c.

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

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

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

◆ __ssa_first

int __ssa_first
static

Definition at line 190 of file ssa.c.

◆ __ssa_ncov

int __ssa_ncov
static

Definition at line 187 of file ssa.c.

◆ __ssa_npar

int __ssa_npar
static

Definition at line 186 of file ssa.c.

◆ __ssa_nvar

int __ssa_nvar
static

Definition at line 185 of file ssa.c.

◆ __ssa_ratefn

SEXP __ssa_ratefn
static

Definition at line 189 of file ssa.c.

◆ __ssa_rxrate

pomp_ssa_rate_fn* __ssa_rxrate
static

Definition at line 191 of file ssa.c.