pomp
Inference for partially observed Markov processes
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
decls.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

SEXP bspline_basis (SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP periodic_bspline_basis (SEXP, SEXP, SEXP, SEXP, SEXP)
 
void bspline_basis_eval_deriv (double, double *, int, int, int, double *)
 
void periodic_bspline_basis_eval_deriv (double, double, int, int, int, double *)
 
SEXP do_dinit (SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP R_Euler_Multinom (SEXP, SEXP, SEXP, SEXP)
 
SEXP D_Euler_Multinom (SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP R_GammaWN (SEXP, SEXP, SEXP)
 
SEXP R_BetaBinom (SEXP, SEXP, SEXP, SEXP)
 
SEXP D_BetaBinom (SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP do_dmeasure (SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP do_dprior (SEXP, SEXP, SEXP, SEXP)
 
SEXP do_dprocess (SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP do_emeasure (SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP euler_simulator (SEXP, SEXP, SEXP, SEXP, SEXP, double, rprocmode, SEXP, SEXP, SEXP, SEXP)
 
int num_euler_steps (double, double, double *)
 
int num_map_steps (double, double, double)
 
void _gompertz_normal_dmeasure (double *, double *, double *, double *, int, int *, int *, int *, int *, double *, double)
 
void _gompertz_normal_rmeasure (double *, double *, double *, int *, int *, int *, int *, double *, double)
 
void _gompertz_normal_emeasure (double *, double *, double *, int *, int *, int *, int *, double *, double)
 
void _gompertz_normal_vmeasure (double *, double *, double *, int *, int *, int *, int *, double *, double)
 
void _gompertz_step (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void _gompertz_skeleton (double *, double *, const double *, const int *, const int *, const int *, const double *, double)
 
void _gompertz_to_trans (double *, const double *, const int *)
 
void _gompertz_from_trans (double *, const double *, const int *)
 
void R_init_pomp (DllInfo *)
 
SEXP logmeanexp (const SEXP, const SEXP)
 
SEXP get_covariate_names (SEXP)
 
lookup_table_t make_covariate_table (SEXP, int *)
 
SEXP lookup_in_table (SEXP, SEXP)
 
void table_lookup (lookup_table_t *, double, double *)
 
SEXP randwalk_perturbation (SEXP, SEXP)
 
void _ou2_step (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void _ou2_pdf (double *, double *, double *, double, double, const double *, const int *, const int *, const int *, const double *)
 
void _ou2_skel (double *, double *, double *, int *, int *, int *, double *, double)
 
void _ou2_dmeasure (double *, double *, double *, double *, int, int *, int *, int *, int *, double *, double)
 
void _ou2_rmeasure (double *, double *, double *, int *, int *, int *, int *, double *, double)
 
void _ou2_emeasure (double *, double *, double *, int *, int *, int *, int *, double *, double)
 
void _ou2_vmeasure (double *, double *, double *, int *, int *, int *, int *, double *, double)
 
SEXP do_partrans (SEXP, SEXP, SEXP, SEXP)
 
SEXP pfilter (SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP pomp_fun_handler (SEXP, SEXP, pompfunmode *, SEXP, SEXP, SEXP, SEXP)
 
SEXP load_stack_incr (SEXP)
 
SEXP load_stack_decr (SEXP)
 
SEXP apply_probe_data (SEXP, SEXP)
 
SEXP apply_probe_sim (SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP probe_acf (SEXP, SEXP, SEXP)
 
SEXP probe_ccf (SEXP, SEXP, SEXP, SEXP)
 
SEXP probe_marginal_setup (SEXP, SEXP, SEXP)
 
SEXP probe_marginal_solve (SEXP, SEXP, SEXP)
 
SEXP probe_nlar (SEXP, SEXP, SEXP)
 
SEXP systematic_resampling (SEXP, SEXP)
 
void nosort_resamp (int, double *, int, int *, int)
 
SEXP do_rinit (SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP do_rmeasure (SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP do_rprior (SEXP, SEXP, SEXP)
 
SEXP do_rprocess (SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP do_simulate (SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP add_skel_args (SEXP, SEXP, SEXP, SEXP)
 
void eval_skeleton_R (double *, double *, double *, double *, SEXP, SEXP, SEXP, int, int, int, int, int, int, int, lookup_table_t *, double *)
 
void iterate_skeleton_R (double *, double, double, double *, double *, double *, SEXP, SEXP, SEXP, int, int, int, int, int, int, int, lookup_table_t *, int *, double *)
 
void eval_skeleton_native (double *, double *, double *, double *, int, int, int, int, int, int, int, int *, int *, int *, lookup_table_t *, pomp_skeleton *, SEXP, double *)
 
void iterate_skeleton_native (double *, double, double, double *, double *, double *, int, int, int, int, int, int, int, int *, int *, int *, lookup_table_t *, int *, pomp_skeleton *, SEXP, double *)
 
SEXP do_skeleton (SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP sobol_sequence (SEXP, SEXP)
 
SEXP SSA_simulator (SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP synth_loglik (SEXP, SEXP)
 
SEXP iterate_map (SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP pomp_desolve_setup (SEXP, SEXP, SEXP, SEXP)
 
void pomp_vf_eval (int *, double *, double *, double *, double *, int *)
 
SEXP pomp_desolve_takedown (void)
 
SEXP LogitTransform (SEXP)
 
SEXP ExpitTransform (SEXP)
 
SEXP LogBarycentricTransform (SEXP)
 
SEXP InverseLogBarycentricTransform (SEXP)
 
SEXP set_pomp_userdata (SEXP)
 
const SEXP get_userdata (const char *)
 
const int * get_userdata_int (const char *)
 
const double * get_userdata_double (const char *)
 
SEXP do_vmeasure (SEXP, SEXP, SEXP, SEXP, SEXP)
 
SEXP wpfilter (SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
 

Function Documentation

◆ _gompertz_from_trans()

void _gompertz_from_trans ( double *  __p,
const double *  __pt,
const int *  __parindex 
)

Definition at line 113 of file gompertz.c.

116  {
117  r = exp(T_r);
118  K = exp(T_K);
119  sigma = exp(T_sigma);
120  tau = exp(T_tau);
121  X_0 = exp(T_X_0);
122 }
#define X_0
Definition: gompertz.c:95
#define r
Definition: gompertz.c:91
#define sigma
Definition: gompertz.c:93
#define T_sigma
Definition: gompertz.c:98
#define T_tau
Definition: gompertz.c:99
#define K
Definition: gompertz.c:92
#define tau
Definition: gompertz.c:94
#define T_K
Definition: gompertz.c:97
#define T_r
Definition: gompertz.c:96
#define T_X_0
Definition: gompertz.c:100

◆ _gompertz_normal_dmeasure()

void _gompertz_normal_dmeasure ( double *  lik,
double *  y,
double *  x,
double *  p,
int  give_log,
int *  obsindex,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  covars,
double  t 
)

Definition at line 20 of file gompertz.c.

25  {
26  *lik = dlnorm(Y,log(X),TAU,give_log);
27 }
#define X
Definition: gompertz.c:14
#define TAU
Definition: gompertz.c:11
#define Y
Definition: gompertz.c:13

◆ _gompertz_normal_emeasure()

void _gompertz_normal_emeasure ( double *  f,
double *  x,
double *  p,
int *  obsindex,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  covars,
double  t 
)

Definition at line 40 of file gompertz.c.

45  {
46  EY = X*exp(TAU*TAU/2);
47 }
#define EY
Definition: gompertz.c:16

◆ _gompertz_normal_rmeasure()

void _gompertz_normal_rmeasure ( double *  y,
double *  x,
double *  p,
int *  obsindex,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  covars,
double  t 
)

Definition at line 30 of file gompertz.c.

35  {
36  Y = rlnorm(log(X),TAU);
37 }

◆ _gompertz_normal_vmeasure()

void _gompertz_normal_vmeasure ( double *  f,
double *  x,
double *  p,
int *  vmatindex,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  covars,
double  t 
)

Definition at line 50 of file gompertz.c.

52  {
53  double et = exp(TAU*TAU);
54  VY = X*X*et*(et-1);
55 }
#define VY
Definition: gompertz.c:17

◆ _gompertz_skeleton()

void _gompertz_skeleton ( double *  f,
double *  x,
const double *  p,
const int *  stateindex,
const int *  parindex,
const int *  covindex,
const double *  covar,
double  t 
)

Definition at line 70 of file gompertz.c.

75  {
76  double deltat = 1.0;
77  double S = exp(-R*deltat);
78  XPRIME = R_pow(K,1-S)*R_pow(X,S); // X is not over-written in the skeleton function
79 }
#define XPRIME
Definition: gompertz.c:15
#define R
Definition: gompertz.c:8

◆ _gompertz_step()

void _gompertz_step ( double *  x,
const double *  p,
const int *  stateindex,
const int *  parindex,
const int *  covindex,
const double *  covar,
double  t,
double  deltat 
)

Definition at line 58 of file gompertz.c.

63  {
64  double S = exp(-R*deltat);
65  double eps = (SIGMA > 0.0) ? exp(rnorm(0,SIGMA)) : 1.0;
66  X = R_pow(K,1-S)*R_pow(X,S)*eps; // note X is over-written by this line
67 }
#define SIGMA
Definition: gompertz.c:10

◆ _gompertz_to_trans()

void _gompertz_to_trans ( double *  __pt,
const double *  __p,
const int *  __parindex 
)

Definition at line 102 of file gompertz.c.

105  {
106  T_r = log(r);
107  T_K = log(K);
108  T_sigma = log(sigma);
109  T_tau = log(tau);
110  T_X_0 = log(X_0);
111 }

◆ _ou2_dmeasure()

void _ou2_dmeasure ( double *  lik,
double *  y,
double *  x,
double *  p,
int  give_log,
int *  obsindex,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  covar,
double  t 
)

Definition at line 106 of file ou2.c.

111  {
112  double sd = fabs(TAU);
113  double f = 0.0;
114  f += (ISNA(Y1)) ? 0.0 : dnorm(Y1,x[X1],sd,1);
115  f += (ISNA(Y2)) ? 0.0 : dnorm(Y2,x[X2],sd,1);
116  *lik = (give_log) ? f : exp(f);
117 }
#define X1
Definition: ou2.c:59
#define Y1
Definition: ou2.c:62
#define Y2
Definition: ou2.c:63
#define TAU
Definition: ou2.c:57
#define X2
Definition: ou2.c:60

◆ _ou2_emeasure()

void _ou2_emeasure ( double *  y,
double *  x,
double *  p,
int *  obsindex,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  covar,
double  t 
)

Definition at line 132 of file ou2.c.

137  {
138  Y1 = x[X1];
139  Y2 = x[X2];
140 }

◆ _ou2_pdf()

void _ou2_pdf ( double *  f,
double *  x,
double *  z,
double  t1,
double  t2,
const double *  p,
const int *  stateindex,
const int *  parindex,
const int *  covindex,
const double *  covars 
)

Definition at line 82 of file ou2.c.

88  {
89  if (t2-t1 != 1)
90  err("ou2_pdf error: transitions must be consecutive");
91  f[0] = dens_ou2(x[X1],x[X2],z[X1],z[X2],ALPHA1,ALPHA2,ALPHA3,ALPHA4,
92  SIGMA1,SIGMA2,SIGMA3,1);
93 }
#define SIGMA3
Definition: ou2.c:56
#define ALPHA2
Definition: ou2.c:51
static double dens_ou2(double x1, double x2, double z1, double z2, double alpha1, double alpha2, double alpha3, double alpha4, double sigma1, double sigma2, double sigma3, int give_log)
Definition: ou2.c:30
#define SIGMA1
Definition: ou2.c:54
#define SIGMA2
Definition: ou2.c:55
#define ALPHA3
Definition: ou2.c:52
#define ALPHA4
Definition: ou2.c:53
#define ALPHA1
Definition: ou2.c:50
#define err(...)
Definition: pomp.h:21
Here is the call graph for this function:

◆ _ou2_rmeasure()

void _ou2_rmeasure ( double *  y,
double *  x,
double *  p,
int *  obsindex,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  covar,
double  t 
)

Definition at line 120 of file ou2.c.

125  {
126  double sd = fabs(TAU);
127  Y1 = rnorm(x[X1],sd);
128  Y2 = rnorm(x[X2],sd);
129 }

◆ _ou2_skel()

void _ou2_skel ( double *  f,
double *  x,
double *  p,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  covars,
double  t 
)

Definition at line 95 of file ou2.c.

100  {
101  f[X1] = ALPHA1*x[X1]+ALPHA3*x[X2];
102  f[X2] = ALPHA2*x[X1]+ALPHA4*x[X2];
103 }

◆ _ou2_step()

void _ou2_step ( double *  x,
const double *  p,
const int *  stateindex,
const int *  parindex,
const int *  covindex,
const double *  covars,
double  t,
double  deltat 
)

Definition at line 71 of file ou2.c.

76  {
78 }
static void sim_ou2(double *x1, double *x2, double alpha1, double alpha2, double alpha3, double alpha4, double sigma1, double sigma2, double sigma3)
Definition: ou2.c:10
Here is the call graph for this function:

◆ _ou2_vmeasure()

void _ou2_vmeasure ( double *  f,
double *  x,
double *  p,
int *  vmatindex,
int *  stateindex,
int *  parindex,
int *  covindex,
double *  covar,
double  t 
)

Definition at line 143 of file ou2.c.

148  {
149  double sd = fabs(TAU);
150  V11 = V22 = sd*sd;
151  V12 = V21 = 0;
152 }
#define V12
Definition: ou2.c:67
#define V21
Definition: ou2.c:66
#define V22
Definition: ou2.c:68
#define V11
Definition: ou2.c:65

◆ add_skel_args()

SEXP add_skel_args ( SEXP  args,
SEXP  Snames,
SEXP  Pnames,
SEXP  Cnames 
)

Definition at line 10 of file skeleton.c.

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

◆ apply_probe_data()

SEXP apply_probe_data ( SEXP  object,
SEXP  probes 
)

Definition at line 4 of file probe.c.

4  {
5  SEXP retval, data, vals;
6  int nprobe;
7  int i;
8 
9  nprobe = LENGTH(probes);
10  PROTECT(data = GET_SLOT(object,install("data")));
11  PROTECT(vals = NEW_LIST(nprobe));
12  SET_NAMES(vals,GET_NAMES(probes));
13 
14  for (i = 0; i < nprobe; i++) {
15  SET_ELEMENT(vals,i,eval(PROTECT(lang2(VECTOR_ELT(probes,i),data)),
16  CLOENV(VECTOR_ELT(probes,i))));
17  if (!IS_NUMERIC(VECTOR_ELT(vals,i))) {
18  err("probe %d returns a non-numeric result",i+1);
19  }
20  UNPROTECT(1);
21  }
22  PROTECT(vals = VectorToPairList(vals));
23  PROTECT(retval = eval(PROTECT(LCONS(install("c"),vals)),R_BaseEnv));
24 
25  UNPROTECT(5);
26  return retval;
27 }
Here is the caller graph for this function:

◆ apply_probe_sim()

SEXP apply_probe_sim ( SEXP  object,
SEXP  nsim,
SEXP  params,
SEXP  probes,
SEXP  datval,
SEXP  gnsi 
)

Definition at line 29 of file probe.c.

30  {
31  SEXP x, y, names, sims;
32  SEXP returntype, retval, val, valnames;
33  int nprobe, nsims, nobs, ntimes, nvals;
34  int xdim[2];
35  double *xp, *yp;
36  int p, s, i, j, k, len0 = 0, len = 0;
37 
38  PROTECT(nsim = AS_INTEGER(nsim));
39  if ((LENGTH(nsim)!=1) || (INTEGER(nsim)[0]<=0))
40  err("'nsim' must be a positive integer."); // #nocov
41 
42  PROTECT(gnsi = duplicate(gnsi));
43 
44  // 'names' holds the names of the probe values
45  // we get these from a previous call to 'apply_probe_data'
46  nprobe = LENGTH(probes);
47  nvals = LENGTH(datval);
48  PROTECT(names = GET_NAMES(datval));
49  PROTECT(returntype = NEW_INTEGER(1));
50  *(INTEGER(returntype)) = 0;
51  PROTECT(sims = do_simulate(object,params,nsim,returntype,gnsi));
52  PROTECT(y = VECTOR_ELT(sims,1));
53  *(INTEGER(gnsi)) = 0;
54 
55  nobs = INTEGER(GET_DIM(y))[0];
56  nsims = INTEGER(GET_DIM(y))[1];
57  ntimes = INTEGER(GET_DIM(y))[2];
58  // set up temporary storage
59  xdim[0] = nobs; xdim[1] = ntimes;
60  PROTECT(x = makearray(2,xdim));
61  setrownames(x,GET_ROWNAMES(GET_DIMNAMES(y)),2);
62 
63  // set up matrix to hold results
64  xdim[0] = nsims; xdim[1] = nvals;
65  PROTECT(retval = makearray(2,xdim));
66  PROTECT(valnames = NEW_LIST(2));
67  SET_ELEMENT(valnames,1,names); // set column names
68  SET_DIMNAMES(retval,valnames);
69 
70  for (p = 0, k = 0; p < nprobe; p++, k += len) { // loop over probes
71 
72  R_CheckUserInterrupt();
73 
74  for (s = 0; s < nsims; s++) { // loop over simulations
75 
76  // copy the data from y[,s,] to x[,]
77  xp = REAL(x);
78  yp = REAL(y)+nobs*s;
79  for (j = 0; j < ntimes; j++, yp += nobs*nsims) {
80  for (i = 0; i < nobs; i++, xp++) *xp = yp[i];
81  }
82 
83  // evaluate the probe on the simulated data
84  PROTECT(val = eval(PROTECT(lang2(VECTOR_ELT(probes,p),x)),
85  CLOENV(VECTOR_ELT(probes,p))));
86  if (!IS_NUMERIC(val)) {
87  err("probe %d returns a non-numeric result.",p+1);
88  }
89 
90  len = LENGTH(val);
91  if (s == 0)
92  len0 = len;
93  else if (len != len0) {
94  err("variable-sized results returned by probe %d.",p+1);
95  }
96  if (k+len > nvals)
97  err("probes return different number of values on different datasets.");
98 
99  xp = REAL(retval); yp = REAL(val);
100  for (i = 0; i < len; i++) xp[s+nsims*(i+k)] = yp[i];
101 
102  UNPROTECT(2);
103  }
104 
105  }
106  if (k != nvals)
107  err("probes return different number of values on different datasets."); // #nocov
108 
109  UNPROTECT(9);
110  return retval;
111 }
SEXP do_simulate(SEXP, SEXP, SEXP, SEXP, SEXP)
Definition: simulate.c:6
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
SEXP params
Definition: trajectory.c:128
Here is the call graph for this function:
Here is the caller graph for this function:

◆ bspline_basis()

SEXP bspline_basis ( SEXP  range,
SEXP  x,
SEXP  nbasis,
SEXP  degree,
SEXP  deriv 
)

Definition at line 51 of file bspline.c.

51  {
52  SEXP y, xr;
53  int nx = LENGTH(x);
54  int nb = INTEGER_VALUE(nbasis);
55  int deg = INTEGER_VALUE(degree);
56  int d = INTEGER_VALUE(deriv);
57  int nk = nb+deg+1;
58  double dx, minx, maxx;
59  double *knots = NULL;
60  double *xdata, *ydata;
61  int i;
62  if (deg < 0) err("must have degree > 0");
63  if (nb <= deg) err("must have nbasis > degree");
64  if (d < 0) err("must have deriv >= 0");
65  knots = (double *) R_Calloc(nk,double);
66  PROTECT(xr = AS_NUMERIC(x));
67  PROTECT(y = allocMatrix(REALSXP,nx,nb));
68  xdata = REAL(xr);
69  ydata = REAL(y);
70  minx = REAL(range)[0];
71  maxx = REAL(range)[1];
72  if (minx >= maxx) err("improper range 'rg'");
73  dx = (maxx-minx)/((double) (nb-deg));
74  knots[0] = minx-deg*dx;
75  for (i = 1; i < nk; i++) knots[i] = knots[i-1]+dx;
76  for (i = 0; i < nb; i++) {
77  bspline_eval(ydata,xdata,nx,i,deg,d,&knots[0]);
78  ydata += nx;
79  }
80  R_Free(knots);
81  UNPROTECT(2);
82  return(y);
83 }
static void bspline_eval(double *y, const double *x, int nx, int i, int degree, int deriv, const double *knots)
Definition: bspline.c:10
Here is the call graph for this function:

◆ bspline_basis_eval_deriv()

void bspline_basis_eval_deriv ( double  x,
double *  knots,
int  degree,
int  nbasis,
int  deriv,
double *  y 
)

Definition at line 115 of file bspline.c.

116  {
117  for (int i = 0; i < nbasis; i++) {
118  bspline_eval(&y[i],&x,1,i,degree,deriv,knots);
119  }
120 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ D_BetaBinom()

SEXP D_BetaBinom ( SEXP  x,
SEXP  size,
SEXP  prob,
SEXP  theta,
SEXP  log 
)

Definition at line 121 of file distributions.c.

121  {
122  int k, n, nx, ns, np, nt;
123  double *F, *X, *S, *P, *T;
124  SEXP f;
125  PROTECT(x = AS_NUMERIC(x)); nx = LENGTH(x); X = REAL(x);
126  PROTECT(size = AS_NUMERIC(size)); ns = LENGTH(size); S = REAL(size);
127  PROTECT(prob = AS_NUMERIC(prob)); np = LENGTH(prob); P = REAL(prob);
128  PROTECT(theta = AS_NUMERIC(theta)); nt = LENGTH(theta); T = REAL(theta);
129  PROTECT(log = AS_INTEGER(log));
130  n = (nx > ns) ? nx : ns;
131  n = (n > np) ? n : np;
132  n = (n > nt) ? n : nt;
133  PROTECT(f = NEW_NUMERIC(n)); F = REAL(f);
134  for (k = 0; k < n; k++) {
135  F[k] = dbetabinom(X[k%nx],S[k%ns],P[k%np],T[k%nt],INTEGER(log)[0]);
136  }
137  UNPROTECT(6);
138  return f;
139 }
static R_INLINE double dbetabinom(double x, double size, double prob, double theta, int give_log)
Definition: pomp.h:293
Here is the call graph for this function:

◆ D_Euler_Multinom()

SEXP D_Euler_Multinom ( SEXP  x,
SEXP  size,
SEXP  rate,
SEXP  deltat,
SEXP  log 
)

Definition at line 49 of file distributions.c.

49  {
50  int ntrans = length(rate);
51  int *dim, n;
52  SEXP f;
53  dim = INTEGER(GET_DIM(x));
54  if (dim[0] != ntrans)
55  err("NROW('x') should match length of 'rate'");
56  n = dim[1];
57  if (length(size)>1)
58  warn("in 'deulermultinom': only the first element of 'size' is meaningful");
59  if (length(deltat)>1)
60  warn("in 'deulermultinom': only the first element of 'dt' is meaningful");
61  PROTECT(f = NEW_NUMERIC(n));
62  PROTECT(size = AS_NUMERIC(size));
63  PROTECT(rate = AS_NUMERIC(rate));
64  PROTECT(deltat = AS_NUMERIC(deltat));
65  PROTECT(log = AS_LOGICAL(log));
66  deulermultinom_multi(ntrans,n,REAL(size),REAL(rate),REAL(deltat),REAL(x),INTEGER(log),REAL(f));
67  UNPROTECT(5);
68  return f;
69 }
static void deulermultinom_multi(int m, int n, double *size, double *rate, double *deltat, double *trans, int *give_log, double *f)
Definition: distributions.c:15
#define warn(...)
Definition: pomp.h:22
Here is the call 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 fn
Definition: trajectory.c:138
Here is the call graph for this function:

◆ do_dmeasure()

SEXP do_dmeasure ( SEXP  object,
SEXP  y,
SEXP  x,
SEXP  times,
SEXP  params,
SEXP  log,
SEXP  gnsi 
)

Definition at line 111 of file dmeasure.c.

112 {
113 
115  int ntimes, nvars, npars, ncovars, nreps, nrepsx, nrepsp, nobs;
116  SEXP Snames, Pnames, Cnames, Onames;
117  SEXP cvec, pompfun;
118  SEXP fn, args, ans;
119  SEXP F;
120  int *dim;
121  lookup_table_t covariate_table;
122  double *cov;
123 
124  PROTECT(times = AS_NUMERIC(times));
125  ntimes = length(times);
126  if (ntimes < 1) err("length('times') = 0, no work to do.");
127 
128  PROTECT(y = as_matrix(y));
129  dim = INTEGER(GET_DIM(y));
130  nobs = dim[0];
131 
132  if (ntimes != dim[1])
133  err("length of 'times' and 2nd dimension of 'y' do not agree.");
134 
135  PROTECT(x = as_state_array(x));
136  dim = INTEGER(GET_DIM(x));
137  nvars = dim[0]; nrepsx = dim[1];
138 
139  if (ntimes != dim[2])
140  err("length of 'times' and 3rd dimension of 'x' do not agree.");
141 
142  PROTECT(params = as_matrix(params));
143  dim = INTEGER(GET_DIM(params));
144  npars = dim[0]; nrepsp = dim[1];
145 
146  nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
147 
148  if ((nreps % nrepsp != 0) || (nreps % nrepsx != 0))
149  err("larger number of replicates is not a multiple of smaller.");
150 
151  PROTECT(Onames = GET_ROWNAMES(GET_DIMNAMES(y)));
152  PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x)));
153  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
154  PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
155 
156  // set up the covariate table
157  covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
158  PROTECT(cvec = NEW_NUMERIC(ncovars));
159  cov = REAL(cvec);
160 
161  // extract the user-defined function
162  PROTECT(pompfun = GET_SLOT(object,install("dmeasure")));
163  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,Onames,Cnames));
164 
165  // extract 'userdata' as pairlist
166  PROTECT(args = GET_SLOT(object,install("userdata")));
167 
168  // create array to store results
169  PROTECT(F = ret_array(nreps,ntimes));
170 
171  int nprotect = 13;
172 
173  switch (mode) {
174 
175  case Rfun: {
176 
177  double *ys = REAL(y), *xs = REAL(x), *ps = REAL(params), *time = REAL(times);
178  double *ft = REAL(F);
179  int j, k;
180 
181  // build argument list
182  PROTECT(args = add_args(args,Onames,Snames,Pnames,Cnames,log)); nprotect++;
183 
184  for (k = 0; k < ntimes; k++, time++, ys += nobs) { // loop over times
185 
186  R_CheckUserInterrupt(); // check for user interrupt
187 
188  table_lookup(&covariate_table,*time,cov); // interpolate the covariates
189 
190  for (j = 0; j < nreps; j++, ft++) { // loop over replicates
191 
192  // evaluate the call
193  PROTECT(
194  ans = eval_call(
195  fn,args,
196  time,
197  ys,nobs,
198  xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
199  ps+npars*(j%nrepsp),npars,
200  cov,ncovars
201  )
202  );
203 
204  if (k == 0 && j == 0 && LENGTH(ans) != 1)
205  err("user 'dmeasure' returns a vector of length %d when it should return a scalar.",LENGTH(ans));
206 
207  *ft = *(REAL(AS_NUMERIC(ans)));
208 
209  UNPROTECT(1);
210 
211  }
212  }
213  }
214 
215  break;
216 
217  case native: case regNative: {
218  int *oidx, *sidx, *pidx, *cidx;
219  int give_log;
220  pomp_dmeasure *ff = NULL;
221  double *yp = REAL(y), *xs = REAL(x), *ps = REAL(params), *time = REAL(times);
222  double *ft = REAL(F);
223  double *xp, *pp;
224  int j, k;
225 
226  // extract state, parameter, covariate, observable indices
227  sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
228  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
229  oidx = INTEGER(GET_SLOT(pompfun,install("obsindex")));
230  cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
231 
232  give_log = *(INTEGER(AS_INTEGER(log)));
233 
234  // address of native routine
235  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
236 
237  for (k = 0; k < ntimes; k++, time++, yp += nobs) { // loop over times
238 
239  R_CheckUserInterrupt(); // check for user interrupt
240 
241  // interpolate the covar functions for the covariates
242  table_lookup(&covariate_table,*time,cov);
243 
244  for (j = 0; j < nreps; j++, ft++) { // loop over replicates
245 
246  xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
247  pp = &ps[npars*(j%nrepsp)];
248 
249  (*ff)(ft,yp,xp,pp,give_log,oidx,sidx,pidx,cidx,cov,*time);
250 
251  }
252  }
253 
254  }
255 
256  break;
257 
258  default: {
259  double *ft = REAL(F);
260  int j, k;
261 
262  for (k = 0; k < ntimes; k++) { // loop over times
263  for (j = 0; j < nreps; j++, ft++) { // loop over replicates
264  *ft = R_NaReal;
265  }
266  }
267 
268  warn("'dmeasure' unspecified: likelihood undefined.");
269 
270  }
271 
272  }
273 
274  UNPROTECT(nprotect);
275  return F;
276 }
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 add_args(SEXP args, SEXP Onames, SEXP Snames, SEXP Pnames, SEXP Cnames, SEXP log)
Definition: dmeasure.c:11
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *y, int nobs, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition: dmeasure.c:75
static R_INLINE SEXP ret_array(int nreps, int ntimes)
Definition: dmeasure.c:101
void pomp_dmeasure(double *lik, const double *y, const double *x, const double *p, int give_log, const int *obsindex, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition: pomp.h:86
static R_INLINE SEXP as_state_array(SEXP x)
Definition: pomp_defines.h:183
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:

◆ do_dprior()

SEXP do_dprior ( SEXP  object,
SEXP  params,
SEXP  log,
SEXP  gnsi 
)

Definition at line 51 of file dprior.c.

52 {
53 
55  int npars, nreps;
56  SEXP Pnames, pompfun, fn, args, F;
57  int *dim;
58 
59  PROTECT(params = as_matrix(params));
60  dim = INTEGER(GET_DIM(params));
61  npars = dim[0]; nreps = dim[1];
62 
63  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
64 
65  // extract the user-defined function
66  PROTECT(pompfun = GET_SLOT(object,install("dprior")));
67  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,NA_STRING,Pnames,NA_STRING,NA_STRING));
68 
69  // extract 'userdata' as pairlist
70  PROTECT(args = GET_SLOT(object,install("userdata")));
71 
72  // to store results
73  PROTECT(F = NEW_NUMERIC(nreps));
74 
75  int nprotect = 6;
76 
77  switch (mode) {
78  case Rfun: {
79  SEXP ans;
80  double *ps, *pt;
81  int j;
82 
83  PROTECT(args = add_args(Pnames,log,args)); nprotect++;
84 
85  for (j = 0, ps = REAL(params), pt = REAL(F); j < nreps; j++, ps += npars, pt++) {
86 
87  PROTECT(ans = eval_call(fn,args,ps,npars));
88  *pt = *(REAL(AS_NUMERIC(ans)));
89  UNPROTECT(1);
90 
91  }
92  }
93 
94  break;
95 
96  case native: case regNative: {
97  int give_log, *pidx = 0;
98  pomp_dprior *ff = NULL;
99  double *ps, *pt;
100  int j;
101 
102  // construct state, parameter, covariate, observable indices
103  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
104 
105  // address of native routine
106  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
107 
108  give_log = *(INTEGER(AS_INTEGER(log)));
109 
110  R_CheckUserInterrupt(); // check for user interrupt
111 
112  // loop over replicates
113  for (j = 0, pt = REAL(F), ps = REAL(params); j < nreps; j++, ps += npars, pt++)
114  (*ff)(pt,ps,give_log,pidx);
115 
116  }
117 
118  break;
119 
120  default: {
121  int give_log, j;
122  double *pt;
123 
124  give_log = *(INTEGER(AS_INTEGER(log)));
125 
126  // loop over replicates
127  for (j = 0, pt = REAL(F); j < nreps; j++, pt++)
128  *pt = (give_log) ? 0.0 : 1.0;
129 
130  }
131 
132  }
133 
134  UNPROTECT(nprotect);
135  return F;
136 }
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *p, int n)
Definition: dprior.c:34
static R_INLINE SEXP add_args(SEXP names, SEXP log, SEXP args)
Definition: dprior.c:10
void pomp_dprior(double *lik, const double *p, int give_log, const int *parindex)
Definition: pomp.h:111
Here is the call 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 }
static SEXP onestep_density(SEXP func, SEXP x, SEXP times, SEXP params, SEXP covar, SEXP log, SEXP args, SEXP gnsi)
Definition: dprocess.c:128
Here is the call graph for this function:

◆ do_emeasure()

SEXP do_emeasure ( SEXP  object,
SEXP  x,
SEXP  times,
SEXP  params,
SEXP  gnsi 
)

Definition at line 98 of file emeasure.c.

99 {
101  int ntimes, nvars, npars, ncovars, nreps, nrepsx, nrepsp;
102  int nobs = 0;
103  SEXP Snames, Pnames, Cnames, Onames = R_NilValue;
104  SEXP fn, args;
105  SEXP pompfun;
106  SEXP Y = R_NilValue;
107  int *dim;
108  lookup_table_t covariate_table;
109  SEXP cvec;
110  double *cov;
111 
112  PROTECT(times = AS_NUMERIC(times));
113  ntimes = length(times);
114  if (ntimes < 1)
115  err("length('times') = 0, no work to do.");
116 
117  PROTECT(x = as_state_array(x));
118  dim = INTEGER(GET_DIM(x));
119  nvars = dim[0]; nrepsx = dim[1];
120 
121  if (ntimes != dim[2])
122  err("length of 'times' and 3rd dimension of 'x' do not agree.");
123 
124  PROTECT(params = as_matrix(params));
125  dim = INTEGER(GET_DIM(params));
126  npars = dim[0]; nrepsp = dim[1];
127 
128  nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
129 
130  if ((nreps % nrepsp != 0) || (nreps % nrepsx != 0))
131  err("larger number of replicates is not a multiple of smaller.");
132 
133  PROTECT(pompfun = GET_SLOT(object,install("emeasure")));
134 
135  PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x)));
136  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
137  PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
138  PROTECT(Onames = GET_SLOT(pompfun,install("obsnames")));
139 
140  // set up the covariate table
141  covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
142  PROTECT(cvec = NEW_NUMERIC(ncovars));
143  cov = REAL(cvec);
144 
145  // extract the user-defined function
146  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,Onames,Cnames));
147 
148  // extract 'userdata' as pairlist
149  PROTECT(args = GET_SLOT(object,install("userdata")));
150 
151  int nprotect = 11;
152  int first = 1;
153 
154  // first do setup
155  switch (mode) {
156 
157  case Rfun: {
158  double *ys, *yt = 0;
159  double *time = REAL(times), *xs = REAL(x), *ps = REAL(params);
160  SEXP ans;
161  int j, k;
162 
163  PROTECT(args = add_args(args,Snames,Pnames,Cnames)); nprotect++;
164 
165  for (k = 0; k < ntimes; k++, time++) { // loop over times
166 
167  R_CheckUserInterrupt(); // check for user interrupt
168 
169  table_lookup(&covariate_table,*time,cov); // interpolate the covariates
170 
171  for (j = 0; j < nreps; j++) { // loop over replicates
172 
173  if (first) {
174 
175  PROTECT(
176  ans = eval_call(
177  fn,args,
178  time,
179  xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
180  ps+npars*(j%nrepsp),npars,
181  cov,ncovars
182  )
183  );
184 
185  nobs = LENGTH(ans);
186 
187  PROTECT(Onames = GET_NAMES(ans));
188  if (invalid_names(Onames))
189  err("'emeasure' must return a named numeric vector.");
190 
191  PROTECT(Y = ret_array(nobs,nreps,ntimes,Onames));
192 
193  nprotect += 3;
194 
195  yt = REAL(Y);
196  ys = REAL(AS_NUMERIC(ans));
197 
198  memcpy(yt,ys,nobs*sizeof(double));
199  yt += nobs;
200 
201  first = 0;
202 
203  } else {
204 
205  PROTECT(
206  ans = eval_call(
207  fn,args,
208  time,
209  xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
210  ps+npars*(j%nrepsp),npars,
211  cov,ncovars
212  )
213  );
214 
215  if (LENGTH(ans) != nobs)
216  err("'emeasure' returns variable-length results.");
217 
218  ys = REAL(AS_NUMERIC(ans));
219 
220  memcpy(yt,ys,nobs*sizeof(double));
221  yt += nobs;
222 
223  UNPROTECT(1);
224 
225  }
226 
227  }
228  }
229 
230  }
231 
232  break;
233 
234  case native: case regNative: {
235  double *yt = 0, *xp, *pp;
236  double *time = REAL(times), *xs = REAL(x), *ps = REAL(params);
237  int *oidx, *sidx, *pidx, *cidx;
238  pomp_emeasure *ff = NULL;
239  int j, k;
240 
241  nobs = LENGTH(Onames);
242  // extract observable, state, parameter covariate indices
243  sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
244  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
245  oidx = INTEGER(GET_SLOT(pompfun,install("obsindex")));
246  cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
247 
248  // address of native routine
249  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
250 
251  PROTECT(Y = ret_array(nobs,nreps,ntimes,Onames)); nprotect++;
252  yt = REAL(Y);
253 
254  for (k = 0; k < ntimes; k++, time++) { // loop over times
255 
256  R_CheckUserInterrupt(); // check for user interrupt
257 
258  // interpolate the covar functions for the covariates
259  table_lookup(&covariate_table,*time,cov);
260 
261  for (j = 0; j < nreps; j++, yt += nobs) { // loop over replicates
262 
263  xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
264  pp = &ps[npars*(j%nrepsp)];
265 
266  (*ff)(yt,xp,pp,oidx,sidx,pidx,cidx,cov,*time);
267 
268  }
269  }
270 
271  }
272 
273  break;
274 
275  default: {
276  nobs = LENGTH(Onames);
277  int dim[3] = {nobs, nreps, ntimes};
278  const char *dimnm[3] = {"name",".id","time"};
279  double *yt = 0;
280  int i, n = nobs*nreps*ntimes;
281 
282  PROTECT(Y = makearray(3,dim)); nprotect++;
283  setrownames(Y,Onames,3);
284  fixdimnames(Y,dimnm,3);
285 
286  for (i = 0, yt = REAL(Y); i < n; i++, yt++) *yt = R_NaReal;
287 
288  warn("'emeasure' unspecified: NAs generated.");
289  }
290 
291  }
292 
293  UNPROTECT(nprotect);
294  return Y;
295 }
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition: emeasure.c:60
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition: emeasure.c:10
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition: emeasure.c:84
void pomp_emeasure(double *f, const double *x, const double *p, const int *obsindex, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition: pomp.h:93
static R_INLINE void fixdimnames(SEXP x, const char **names, int n)
Definition: pomp_defines.h:129
static R_INLINE int invalid_names(SEXP names)
Definition: pomp_defines.h:55
Here is the call graph for this function:

◆ do_partrans()

SEXP do_partrans ( SEXP  object,
SEXP  params,
SEXP  dir,
SEXP  gnsi 
)

Definition at line 49 of file partrans.c.

50 {
51 
52  SEXP Pnames, tparams, pompfun, fn, args, ob;
54  direction_t direc;
55  int qvec, npars, nreps;
56  int *dim;
57 
58  qvec = isNull(GET_DIM(params)); // is 'params' a vector?
59 
60  PROTECT(tparams = duplicate(params));
61 
62  // coerce 'params' to matrix
63  PROTECT(tparams = as_matrix(tparams));
64  dim = INTEGER(GET_DIM(tparams));
65  npars = dim[0]; nreps = dim[1];
66 
67  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(tparams)));
68 
69  // determine direction of transformation and extract corresponding pomp_fun
70  direc = (direction_t) *(INTEGER(dir));
71  PROTECT(ob = GET_SLOT(object,install("partrans")));
72  switch (direc) {
73  case from: default: // from estimation scale
74  PROTECT(pompfun = GET_SLOT(ob,install("from")));
75  break;
76  case to: // to estimation scale
77  PROTECT(pompfun = GET_SLOT(ob,install("to")));
78  break;
79  }
80 
81  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,NA_STRING,Pnames,NA_STRING,NA_STRING));
82 
83  // extract 'userdata' as pairlist
84  PROTECT(args = GET_SLOT(object,install("userdata")));
85 
86  int nprotect = 7;
87 
88  switch (mode) {
89 
90  case Rfun: {
91 
92  SEXP ans, nm;
93  double *pa, *ps = REAL(tparams);
94  int *posn;
95  int i, j;
96 
97  PROTECT(args = add_args(args,Pnames));
98  PROTECT(ans = eval_call(fn,args,ps,npars));
99 
100  PROTECT(nm = GET_NAMES(ans));
101  if (invalid_names(nm))
102  err("user transformation functions must return named numeric vectors.");
103  posn = INTEGER(PROTECT(matchnames(Pnames,nm,"parameters")));
104 
105  nprotect += 4;
106 
107  pa = REAL(AS_NUMERIC(ans));
108 
109  for (i = 0; i < LENGTH(ans); i++) ps[posn[i]] = pa[i];
110 
111  for (j = 1, ps += npars; j < nreps; j++, ps += npars) {
112 
113  PROTECT(ans = eval_call(fn,args,ps,npars));
114  pa = REAL(AS_NUMERIC(ans));
115  for (i = 0; i < LENGTH(ans); i++) ps[posn[i]] = pa[i];
116  UNPROTECT(1);
117 
118  }
119 
120  }
121 
122  break;
123 
124  case native: case regNative: {
125 
126  pomp_transform *ff;
127  double *ps, *pt;
128  int *idx;
129  int j;
130 
131  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
132 
133  R_CheckUserInterrupt();
134 
135  idx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
136 
137  for (j = 0, ps = REAL(params), pt = REAL(tparams); j < nreps; j++, ps += npars, pt += npars)
138  (*ff)(pt,ps,idx);
139 
140  }
141 
142  break;
143 
144  default: // #nocov
145 
146  break; // #nocov
147 
148  }
149 
150  if (qvec) {
151  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(tparams))); nprotect++;
152  SET_DIM(tparams,R_NilValue);
153  SET_NAMES(tparams,Pnames);
154  }
155 
156  UNPROTECT(nprotect);
157  return tparams;
158 
159 }
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *p, int n)
Definition: partrans.c:33
static R_INLINE SEXP add_args(SEXP args, SEXP names)
Definition: partrans.c:12
direction_t
Definition: partrans.c:10
@ from
Definition: partrans.c:10
@ to
Definition: partrans.c:10
void pomp_transform(double *pt, const double *p, const int *parindex)
Definition: pomp.h:115
static R_INLINE SEXP matchnames(SEXP provided, SEXP needed, const char *where)
Definition: pomp_defines.h:60
Here is the call graph for this function:

◆ do_rinit()

SEXP do_rinit ( SEXP  object,
SEXP  params,
SEXP  t0,
SEXP  nsim,
SEXP  gnsi 
)

Definition at line 91 of file rinit.c.

92 {
93 
94  SEXP Pnames, Cnames, Snames, pcnames;
95  SEXP x = R_NilValue;
96  SEXP pompfun, fn, args;
98  lookup_table_t covariate_table;
99  SEXP cvec;
100  double *cov;
101  int *dim;
102  int npar, nrep, nvar, ncovars, nsims, ns;
103 
104  nsims = *(INTEGER(AS_INTEGER(nsim)));
105  PROTECT(params = as_matrix(params));
106  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
107  PROTECT(pcnames = GET_COLNAMES(GET_DIMNAMES(params)));
108 
109  dim = INTEGER(GET_DIM(params));
110  npar = dim[0]; nrep = dim[1];
111  ns = nsims*nrep;
112 
113  // set up the covariate table
114  covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
115  PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
116  PROTECT(cvec = NEW_NUMERIC(ncovars));
117  cov = REAL(cvec);
118 
119  table_lookup(&covariate_table,*(REAL(t0)),cov);
120 
121  // extract userdata
122  PROTECT(args = GET_SLOT(object,install("userdata")));
123 
124  PROTECT(pompfun = GET_SLOT(object,install("rinit")));
125  PROTECT(Snames = GET_SLOT(pompfun,install("statenames")));
126 
127  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
128 
129  int nprotect = 9;
130 
131  switch (mode) {
132  case Rfun: {
133 
134  SEXP ans;
135  double *time = REAL(AS_NUMERIC(t0));
136  double *ps = REAL(params);
137  double *xs, *xt = NULL;
138  int *midx;
139  int j;
140 
141  PROTECT(args = add_args(args,Pnames,Cnames));
142  PROTECT(ans = eval_call(fn,args,time,ps,npar,cov,ncovars));
143  PROTECT(ans = AS_NUMERIC(ans));
144  PROTECT(Snames = GET_NAMES(ans));
145 
146  if (invalid_names(Snames))
147  err("user 'rinit' must return a named numeric vector.");
148 
149  nvar = LENGTH(ans);
150  xs = REAL(ans);
151  midx = INTEGER(PROTECT(match(Pnames,Snames,0)));
152 
153  for (j = 0; j < nvar; j++) {
154  if (midx[j] != 0)
155  err("a state variable and a parameter share the name: '%s'.",CHAR(STRING_ELT(Snames,j)));
156  }
157 
158  PROTECT(x = ret_array(nvar,ns,Snames));
159  xt = REAL(x);
160 
161  memcpy(xt,xs,nvar*sizeof(double));
162 
163  nprotect += 6;
164 
165  for (j = 1, xt += nvar; j < ns; j++, xt += nvar) {
166  PROTECT(ans = eval_call(fn,args,time,ps+npar*(j%nrep),npar,cov,ncovars));
167  xs = REAL(ans);
168  if (LENGTH(ans) != nvar)
169  err("user 'rinit' returns vectors of variable length.");
170  memcpy(xt,xs,nvar*sizeof(double));
171  UNPROTECT(1);
172  }
173 
174  }
175 
176  break;
177 
178  case native: case regNative: {
179 
180  int *sidx, *pidx, *cidx;
181  double *xt, *ps, time;
182  pomp_rinit *ff = NULL;
183  int j;
184 
185  nvar = *INTEGER(GET_SLOT(object,install("nstatevars")));
186  PROTECT(x = ret_array(nvar,ns,Snames)); nprotect++;
187 
188  sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
189  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
190  cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
191 
192  // address of native routine
193  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
194 
195  GetRNGstate();
196 
197  time = *(REAL(t0));
198 
199  // loop over replicates
200  for (j = 0, xt = REAL(x), ps = REAL(params); j < ns; j++, xt += nvar)
201  (*ff)(xt,ps+npar*(j%nrep),time,sidx,pidx,cidx,cov);
202 
203  PutRNGstate();
204 
205  }
206 
207  break;
208 
209  default: {
210 
211  PROTECT(x = pomp_default_rinit(params,Pnames,npar,nrep,ns)); nprotect++;
212 
213  }
214 
215  break;
216 
217  }
218 
219  // now add column names
220  if (nrep > 1) {
221  SEXP dn, xn;
222  int k, *p;
223 
224  if (isNull(pcnames)) {
225  PROTECT(pcnames = NEW_INTEGER(nrep)); nprotect++;
226  for (k = 0, p = INTEGER(pcnames); k < nrep; k++, p++) *p = k+1;
227  }
228 
229  if (nsims > 1) {
230  int k, *sp;
231  SEXP us;
232  PROTECT(us = mkString("_"));
233  PROTECT(xn = NEW_INTEGER(ns));
234  for (k = 0, sp = INTEGER(xn); k < ns; k++, sp++) *sp = (k/nrep)+1;
235  PROTECT(xn = paste0(pcnames,us,xn));
236  PROTECT(dn = GET_DIMNAMES(x));
237  nprotect += 4;
238  SET_ELEMENT(dn,1,xn);
239  SET_DIMNAMES(x,dn);
240 
241  } else {
242 
243  PROTECT(dn = GET_DIMNAMES(x)); nprotect++;
244  SET_ELEMENT(dn,1,pcnames);
245  SET_DIMNAMES(x,dn);
246 
247  }
248  }
249 
250  UNPROTECT(nprotect);
251  return x;
252 }
void pomp_rinit(double *x, const double *p, double t0, const int *stateindex, const int *parindex, const int *covindex, const double *covars)
Definition: pomp.h:40
static R_INLINE SEXP ret_array(int m, int n, SEXP names)
Definition: rinit.c:79
static R_INLINE SEXP paste0(SEXP a, SEXP b, SEXP c)
Definition: rinit.c:10
static SEXP pomp_default_rinit(SEXP params, SEXP Pnames, int npar, int nrep, int nsim)
Definition: rinit.c:254
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t0, double *p, int npar, double *c, int ncov)
Definition: rinit.c:59
static R_INLINE SEXP add_args(SEXP args, SEXP Pnames, SEXP Cnames)
Definition: rinit.c:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ do_rmeasure()

SEXP do_rmeasure ( SEXP  object,
SEXP  x,
SEXP  times,
SEXP  params,
SEXP  gnsi 
)

Definition at line 98 of file rmeasure.c.

99 {
101  int ntimes, nvars, npars, ncovars, nreps, nrepsx, nrepsp;
102  int nobs = 0;
103  SEXP Snames, Pnames, Cnames, Onames = R_NilValue;
104  SEXP fn, args;
105  SEXP pompfun;
106  SEXP Y = R_NilValue;
107  int *dim;
108  lookup_table_t covariate_table;
109  SEXP cvec;
110  double *cov;
111 
112  PROTECT(times = AS_NUMERIC(times));
113  ntimes = length(times);
114  if (ntimes < 1)
115  err("length('times') = 0, no work to do.");
116 
117  PROTECT(x = as_state_array(x));
118  dim = INTEGER(GET_DIM(x));
119  nvars = dim[0]; nrepsx = dim[1];
120 
121  if (ntimes != dim[2])
122  err("length of 'times' and 3rd dimension of 'x' do not agree.");
123 
124  PROTECT(params = as_matrix(params));
125  dim = INTEGER(GET_DIM(params));
126  npars = dim[0]; nrepsp = dim[1];
127 
128  nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
129 
130  if ((nreps % nrepsp != 0) || (nreps % nrepsx != 0))
131  err("larger number of replicates is not a multiple of smaller.");
132 
133  PROTECT(pompfun = GET_SLOT(object,install("rmeasure")));
134 
135  PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x)));
136  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
137  PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
138  PROTECT(Onames = GET_SLOT(pompfun,install("obsnames")));
139 
140  // set up the covariate table
141  covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
142  PROTECT(cvec = NEW_NUMERIC(ncovars));
143  cov = REAL(cvec);
144 
145  // extract the user-defined function
146  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,Onames,Cnames));
147 
148  // extract 'userdata' as pairlist
149  PROTECT(args = GET_SLOT(object,install("userdata")));
150 
151  int nprotect = 11;
152  int first = 1;
153 
154  // first do setup
155  switch (mode) {
156 
157  case Rfun: {
158  double *ys, *yt = 0;
159  double *time = REAL(times), *xs = REAL(x), *ps = REAL(params);
160  SEXP ans;
161  int j, k;
162 
163  PROTECT(args = add_args(args,Snames,Pnames,Cnames)); nprotect++;
164 
165  for (k = 0; k < ntimes; k++, time++) { // loop over times
166 
167  R_CheckUserInterrupt(); // check for user interrupt
168 
169  table_lookup(&covariate_table,*time,cov); // interpolate the covariates
170 
171  for (j = 0; j < nreps; j++) { // loop over replicates
172 
173  if (first) {
174 
175  PROTECT(
176  ans = eval_call(
177  fn,args,
178  time,
179  xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
180  ps+npars*(j%nrepsp),npars,
181  cov,ncovars
182  )
183  );
184 
185  nobs = LENGTH(ans);
186 
187  PROTECT(Onames = GET_NAMES(ans));
188  if (invalid_names(Onames))
189  err("'rmeasure' must return a named numeric vector.");
190 
191  PROTECT(Y = ret_array(nobs,nreps,ntimes,Onames));
192 
193  nprotect += 3;
194 
195  yt = REAL(Y);
196  ys = REAL(AS_NUMERIC(ans));
197 
198  memcpy(yt,ys,nobs*sizeof(double));
199  yt += nobs;
200 
201  first = 0;
202 
203  } else {
204 
205  PROTECT(
206  ans = eval_call(
207  fn,args,
208  time,
209  xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
210  ps+npars*(j%nrepsp),npars,
211  cov,ncovars
212  )
213  );
214 
215  if (LENGTH(ans) != nobs)
216  err("'rmeasure' returns variable-length results.");
217 
218  ys = REAL(AS_NUMERIC(ans));
219 
220  memcpy(yt,ys,nobs*sizeof(double));
221  yt += nobs;
222 
223  UNPROTECT(1);
224 
225  }
226 
227  }
228  }
229 
230  }
231 
232  break;
233 
234  case native: case regNative: {
235  double *yt = 0, *xp, *pp;
236  double *time = REAL(times), *xs = REAL(x), *ps = REAL(params);
237  int *oidx, *sidx, *pidx, *cidx;
238  pomp_rmeasure *ff = NULL;
239  int j, k;
240 
241  nobs = LENGTH(Onames);
242  // extract observable, state, parameter covariate indices
243  sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
244  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
245  oidx = INTEGER(GET_SLOT(pompfun,install("obsindex")));
246  cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
247 
248  // address of native routine
249  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
250 
251  PROTECT(Y = ret_array(nobs,nreps,ntimes,Onames)); nprotect++;
252  yt = REAL(Y);
253 
254  GetRNGstate();
255 
256  for (k = 0; k < ntimes; k++, time++) { // loop over times
257 
258  R_CheckUserInterrupt(); // check for user interrupt
259 
260  // interpolate the covar functions for the covariates
261  table_lookup(&covariate_table,*time,cov);
262 
263  for (j = 0; j < nreps; j++, yt += nobs) { // loop over replicates
264 
265  xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
266  pp = &ps[npars*(j%nrepsp)];
267 
268  (*ff)(yt,xp,pp,oidx,sidx,pidx,cidx,cov,*time);
269 
270  }
271  }
272 
273  PutRNGstate();
274 
275  }
276 
277  break;
278 
279  default: {
280  nobs = LENGTH(Onames);
281  int dim[3] = {nobs, nreps, ntimes};
282  const char *dimnm[3] = {"name",".id","time"};
283  double *yt = 0;
284  int i, n = nobs*nreps*ntimes;
285 
286  PROTECT(Y = makearray(3,dim)); nprotect++;
287  setrownames(Y,Onames,3);
288  fixdimnames(Y,dimnm,3);
289 
290  for (i = 0, yt = REAL(Y); i < n; i++, yt++) *yt = R_NaReal;
291 
292  warn("'rmeasure' unspecified: NAs generated.");
293  }
294 
295  }
296 
297  UNPROTECT(nprotect);
298  return Y;
299 }
void pomp_rmeasure(double *y, const double *x, const double *p, const int *obsindex, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition: pomp.h:79
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition: rmeasure.c:60
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition: rmeasure.c:10
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition: rmeasure.c:84
Here is the call graph for this function:
Here is the caller graph for this function:

◆ do_rprior()

SEXP do_rprior ( SEXP  object,
SEXP  params,
SEXP  gnsi 
)

Definition at line 58 of file rprior.c.

59 {
60 
62  int npars, nreps;
63  SEXP Pnames, pompfun, fn, args;
64  int *dim;
65 
66  PROTECT(params = ret_array(params));
67  dim = INTEGER(GET_DIM(params));
68  npars = dim[0]; nreps = dim[1];
69 
70  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
71 
72  // extract the user-defined function
73  PROTECT(pompfun = GET_SLOT(object,install("rprior")));
74  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,NA_STRING,Pnames,NA_STRING,NA_STRING));
75 
76  // extract 'userdata' as pairlist
77  PROTECT(args = GET_SLOT(object,install("userdata")));
78 
79  int nprotect = 5;
80  int first = 1;
81 
82  switch (mode) {
83 
84  case Rfun: {
85 
86  SEXP ans, nm;
87  double *pa, *p = REAL(params);
88  int *posn = NULL;
89  int i, j;
90 
91  // set up the function call
92  PROTECT(args = add_args(args,Pnames)); nprotect++;
93 
94  for (j = 0; j < nreps; j++, p += npars) {
95 
96  if (first) {
97 
98  PROTECT(ans = eval_call(fn,args,p,npars));
99  PROTECT(ans = AS_NUMERIC(ans));
100 
101  PROTECT(nm = GET_NAMES(ans));
102  if (invalid_names(nm))
103  err("'rprior' must return a named numeric vector.");
104  posn = INTEGER(PROTECT(matchnames(Pnames,nm,"parameters")));
105 
106  nprotect += 4;
107 
108  pa = REAL(ans);
109  for (i = 0; i < LENGTH(ans); i++) p[posn[i]] = pa[i];
110 
111  first = 0;
112 
113  } else {
114 
115  PROTECT(ans = eval_call(fn,args,p,npars));
116  PROTECT(ans = AS_NUMERIC(ans));
117 
118  pa = REAL(ans);
119  for (i = 0; i < LENGTH(ans); i++) p[posn[i]] = pa[i];
120 
121  UNPROTECT(2);
122 
123  }
124  }
125  }
126 
127  break;
128 
129  case native: case regNative: {
130 
131  double *p;
132  int *pidx = 0;
133  pomp_rprior *ff = NULL;
134  int j;
135 
136  // extract parameter index
137  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
138 
139  // address of native routine
140  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
141 
142  R_CheckUserInterrupt(); // check for user interrupt
143 
144  GetRNGstate();
145 
146  // loop over replicates
147  for (j = 0, p = REAL(params); j < nreps; j++, p += npars)
148  (*ff)(p,pidx);
149 
150  PutRNGstate();
151 
152  }
153 
154  break;
155 
156  default: // just duplicate
157 
158  warn("'rprior' unspecified: duplicating parameters.");
159 
160  }
161 
162  UNPROTECT(nprotect);
163  return params;
164 }
void pomp_rprior(double *p, const int *parindex)
Definition: pomp.h:107
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *p, int n)
Definition: rprior.c:31
static R_INLINE SEXP ret_array(SEXP params)
Definition: rprior.c:47
static R_INLINE SEXP add_args(SEXP args, SEXP names)
Definition: rprior.c:10
Here is the call graph for this function:

◆ do_rprocess()

SEXP do_rprocess ( SEXP  object,
SEXP  xstart,
SEXP  tstart,
SEXP  times,
SEXP  params,
SEXP  gnsi 
)

Definition at line 25 of file rprocess.c.

26 {
27 
28  int *xdim, type, nvars, npars, nreps, nrepsx, ntimes;
29  SEXP X, copy, rproc, args, accumvars, covar;
30  SEXP dimXstart, dimP;
31  const char *dimnm[3] = {"name",".id","time"};
32 
33  PROTECT(gnsi = duplicate(gnsi));
34 
35  PROTECT(tstart = AS_NUMERIC(tstart));
36 
37  PROTECT(times = AS_NUMERIC(times));
38  ntimes = length(times);
39  if (ntimes < 1) {
40  err("length(times) < 1: no work to do.");
41  }
42 
43  PROTECT(xstart = as_matrix(xstart));
44  PROTECT(dimXstart = GET_DIM(xstart));
45  xdim = INTEGER(dimXstart);
46  nvars = xdim[0]; nrepsx = xdim[1];
47 
48  PROTECT(params = as_matrix(params));
49  PROTECT(dimP = GET_DIM(params));
50  xdim = INTEGER(dimP);
51  npars = xdim[0]; nreps = xdim[1];
52 
53  int nprotect = 7;
54 
55  if (nrepsx > nreps) { // more ICs than parameters
56  if (nrepsx % nreps != 0) {
57  err("the larger number of replicates is not a multiple of smaller.");
58  } else {
59  double *src, *tgt;
60  int dims[2];
61  int j, k;
62  dims[0] = npars; dims[1] = nrepsx;
63  PROTECT(copy = duplicate(params));
64  PROTECT(params = makearray(2,dims));
65  nprotect += 2;
66  setrownames(params,GET_ROWNAMES(GET_DIMNAMES(copy)),2);
67  src = REAL(copy);
68  tgt = REAL(params);
69  for (j = 0; j < nrepsx; j++) {
70  for (k = 0; k < npars; k++, tgt++) {
71  *tgt = src[k+npars*(j%nreps)];
72  }
73  }
74  }
75  nreps = nrepsx;
76  } else if (nrepsx < nreps) { // more parameters than ICs
77  if (nreps % nrepsx != 0) {
78  err("the larger number of replicates is not a multiple of smaller.");
79  } else {
80  double *src, *tgt;
81  int dims[2];
82  int j, k;
83  dims[0] = nvars; dims[1] = nreps;
84  PROTECT(copy = duplicate(xstart));
85  PROTECT(xstart = makearray(2,dims));
86  nprotect += 2;
87  setrownames(xstart,GET_ROWNAMES(GET_DIMNAMES(copy)),2);
88  src = REAL(copy);
89  tgt = REAL(xstart);
90  for (j = 0; j < nreps; j++) {
91  for (k = 0; k < nvars; k++, tgt++) {
92  *tgt = src[k+nvars*(j%nrepsx)];
93  }
94  }
95  }
96  }
97 
98  PROTECT(rproc = GET_SLOT(object,install("rprocess")));
99  PROTECT(args = GET_SLOT(object,install("userdata")));
100  PROTECT(accumvars = GET_SLOT(object,install("accumvars")));
101  PROTECT(covar = GET_SLOT(object,install("covar")));
102 
103  nprotect += 4;
104 
105  // extract the process function
106  type = *(INTEGER(GET_SLOT(rproc,install("type"))));
107  switch (type) {
108  case onestep: // one-step simulator
109  {
110  SEXP fn;
111  double deltat = 1.0;
112  PROTECT(fn = GET_SLOT(rproc,install("step.fn")));
113  PROTECT(X = euler_simulator(fn,xstart,tstart,times,params,deltat,type,
114  accumvars,covar,args,gnsi));
115  nprotect += 2;
116  }
117  break;
118  case discrete: case euler: // discrete-time and Euler
119  {
120  SEXP fn;
121  double deltat;
122  PROTECT(fn = GET_SLOT(rproc,install("step.fn")));
123  deltat = *(REAL(AS_NUMERIC(GET_SLOT(rproc,install("delta.t")))));
124  PROTECT(X = euler_simulator(fn,xstart,tstart,times,params,deltat,type,
125  accumvars,covar,args,gnsi));
126  nprotect += 2;
127  }
128  break;
129  case gill: // Gillespie's method
130  {
131  SEXP fn, vmatrix, hmax;
132  PROTECT(fn = GET_SLOT(rproc,install("rate.fn")));
133  PROTECT(vmatrix = GET_SLOT(rproc,install("v")));
134  PROTECT(hmax = GET_SLOT(rproc,install("hmax")));
135  PROTECT(X = SSA_simulator(fn,xstart,tstart,times,params,vmatrix,covar,
136  accumvars,hmax,args,gnsi));
137  nprotect += 4;
138  }
139  break;
140  case dflt: default:
141  PROTECT(X = pomp_default_rprocess(xstart,nvars,nreps,ntimes));
142  nprotect++;
143  break;
144  }
145 
146  fixdimnames(X,dimnm,3);
147  UNPROTECT(nprotect);
148  return X;
149 
150 }
SEXP euler_simulator(SEXP, SEXP, SEXP, SEXP, SEXP, double, rprocmode, SEXP, SEXP, SEXP, SEXP)
Definition: euler.c:108
SEXP SSA_simulator(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
Definition: ssa.c:228
@ gill
Definition: pomp_defines.h:17
@ discrete
Definition: pomp_defines.h:17
@ dflt
Definition: pomp_defines.h:17
@ onestep
Definition: pomp_defines.h:17
@ euler
Definition: pomp_defines.h:17
static SEXP pomp_default_rprocess(SEXP xstart, int nvars, int nreps, int ntimes)
Definition: rprocess.c:10
Here is the call graph for this function:
Here is the caller graph for this function:

◆ do_simulate()

SEXP do_simulate ( SEXP  object,
SEXP  params,
SEXP  nsim,
SEXP  rettype,
SEXP  gnsi 
)

Definition at line 6 of file simulate.c.

7 {
8 
9  SEXP t0, times, x0, x, y;
10  SEXP ans = R_NilValue, ans_names = R_NilValue;
11  SEXP simnames;
12  int return_type = *(INTEGER(rettype)); // 0 = array, 1 = pomps
13 
14  if (LENGTH(nsim) != 1) err("'nsim' must be a single integer"); // #nocov
15 
16  PROTECT(params = as_matrix(params));
17 
18  PROTECT(t0 = GET_SLOT(object,install("t0")));
19  PROTECT(times = GET_SLOT(object,install("times")));
20 
21  // initialize the simulations
22  PROTECT(x0 = do_rinit(object,params,t0,nsim,gnsi));
23  PROTECT(simnames = GET_COLNAMES(GET_DIMNAMES(x0)));
24 
25  // call 'rprocess' to simulate state process
26  PROTECT(x = do_rprocess(object,x0,t0,times,params,gnsi));
27 
28  // call 'rmeasure' to simulate the measurement process
29  PROTECT(y = do_rmeasure(object,x,times,params,gnsi));
30 
31  setcolnames(x,simnames);
32  setcolnames(y,simnames);
33 
34  int nprotect = 7;
35 
36  switch (return_type) {
37 
38  case 0: // return a list of arrays
39 
40  PROTECT(ans = NEW_LIST(2));
41  PROTECT(ans_names = NEW_CHARACTER(2));
42  nprotect += 2;
43  SET_STRING_ELT(ans_names,0,mkChar("states"));
44  SET_STRING_ELT(ans_names,1,mkChar("obs"));
45  SET_NAMES(ans,ans_names);
46  SET_ELEMENT(ans,0,x);
47  SET_ELEMENT(ans,1,y);
48 
49  break;
50 
51  case 1: default:
52 
53  // a list to hold the pomp objects
54  {
55  SEXP pp, xx, yy, po;
56  const int *xdim;
57  int nvar, npar, nobs, nsim, ntim, nparsets;
58  int dim[2], i, j, k;
59 
60  PROTECT(po = duplicate(object));
61  SET_SLOT(po,install("t0"),t0);
62  SET_SLOT(po,install("times"),times);
63 
64  xdim = INTEGER(GET_DIM(x));
65  nvar = xdim[0]; nsim = xdim[1]; ntim = xdim[2];
66 
67  xdim = INTEGER(GET_DIM(y));
68  nobs = xdim[0]; // second dimensions of 'x' and 'y' must agree
69 
70  xdim = INTEGER(GET_DIM(params));
71  npar = xdim[0]; nparsets = xdim[1];
72 
73  dim[0] = nvar; dim[1] = ntim;
74  PROTECT(xx = makearray(2,dim));
75  setrownames(xx,GET_ROWNAMES(GET_DIMNAMES(x)),2);
76  SET_SLOT(po,install("states"),xx);
77 
78  dim[0] = nobs; dim[1] = ntim;
79  PROTECT(yy = makearray(2,dim));
80  setrownames(yy,GET_ROWNAMES(GET_DIMNAMES(y)),2);
81  SET_SLOT(po,install("data"),yy);
82 
83  PROTECT(pp = NEW_NUMERIC(npar));
84  SET_NAMES(pp,GET_ROWNAMES(GET_DIMNAMES(params)));
85  SET_SLOT(po,install("params"),pp);
86 
87  PROTECT(ans = NEW_LIST(nsim));
88  SET_NAMES(ans,simnames);
89 
90  nprotect += 5;
91 
92  for (k = 0; k < nsim; k++) {
93 
94  SEXP po2;
95  double *xs = REAL(x), *ys = REAL(y), *ps = REAL(params);
96  double *xt, *yt, *pt;
97 
98  PROTECT(po2 = duplicate(po));
99  xt = REAL(GET_SLOT(po2,install("states")));
100  yt = REAL(GET_SLOT(po2,install("data")));
101  pt = REAL(GET_SLOT(po2,install("params")));
102 
103  memcpy(pt,ps+npar*(k%nparsets),npar*sizeof(double));
104 
105  // copy x[,k,] and y[,k,] into po2
106  for (j = 0; j < ntim; j++) {
107  for (i = 0; i < nvar; i++, xt++) *xt = xs[i+nvar*(k+nsim*j)];
108  for (i = 0; i < nobs; i++, yt++) *yt = ys[i+nobs*(k+nsim*j)];
109  }
110 
111  SET_ELEMENT(ans,k,po2);
112  UNPROTECT(1);
113 
114  }
115 
116  break;
117 
118  }
119 
120  }
121 
122  UNPROTECT(nprotect);
123  return ans;
124 
125 }
SEXP do_rmeasure(SEXP, SEXP, SEXP, SEXP, SEXP)
Definition: rmeasure.c:98
SEXP do_rinit(SEXP, SEXP, SEXP, SEXP, SEXP)
Definition: rinit.c:91
SEXP do_rprocess(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
Definition: rprocess.c:25
static R_INLINE void setcolnames(SEXP x, SEXP names)
Definition: pomp_defines.h:121
Here is the call graph for this function:
Here is the caller graph for this function:

◆ do_skeleton()

SEXP do_skeleton ( SEXP  object,
SEXP  x,
SEXP  t,
SEXP  params,
SEXP  gnsi 
)

Definition at line 326 of file skeleton.c.

327 {
328 
329  int nvars, npars, nrepp, nrepx, nreps, ntimes, ncovars;
331  int *dim;
332  SEXP Snames, Cnames, Pnames;
333  SEXP pompfun, cov, ob;
334  SEXP fn, args, F;
335  lookup_table_t covariate_table;
336 
337  PROTECT(t = AS_NUMERIC(t));
338  ntimes = LENGTH(t);
339 
340  PROTECT(x = as_state_array(x));
341  dim = INTEGER(GET_DIM(x));
342  nvars = dim[0]; nrepx = dim[1];
343  if (ntimes != dim[2])
344  err("length of 't' and 3rd dimension of 'x' do not agree.");
345 
346  PROTECT(params = as_matrix(params));
347  dim = INTEGER(GET_DIM(params));
348  npars = dim[0]; nrepp = dim[1];
349 
350  // 2nd dimension of 'x' and 'params' need not agree
351  nreps = (nrepp > nrepx) ? nrepp : nrepx;
352  if ((nreps % nrepp != 0) || (nreps % nrepx != 0))
353  err("2nd dimensions of 'x' and 'params' are incompatible");
354 
355  PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x)));
356  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
357  PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
358 
359  // set up the covariate table
360  covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
361  PROTECT(cov = NEW_NUMERIC(ncovars));
362 
363  // extract the user-defined function
364  PROTECT(ob = GET_SLOT(object,install("skeleton")));
365  PROTECT(pompfun = GET_SLOT(ob,install("skel.fn")));
366  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
367 
368  // extract 'userdata' as pairlist
369  PROTECT(args = GET_SLOT(object,install("userdata")));
370 
371  PROTECT(F = ret_array(nvars,nreps,ntimes,Snames));
372 
373  int nprotect = 12;
374 
375  switch (mode) {
376 
377  case Rfun: {
378 
379  PROTECT(args = add_skel_args(args,Snames,Pnames,Cnames)); nprotect++;
380 
381  eval_skeleton_R(REAL(F),REAL(t),REAL(x),REAL(params),
382  fn,args,Snames,
383  nvars,npars,ncovars,ntimes,nrepx,nrepp,nreps,
384  &covariate_table,REAL(cov));
385 
386  }
387 
388  break;
389 
390  case native: case regNative: {
391  int *sidx, *pidx, *cidx;
392  pomp_skeleton *ff = NULL;
393 
394  sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
395  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
396  cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
397 
398  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
399 
401  REAL(F),REAL(t),REAL(x),REAL(params),
402  nvars,npars,ncovars,ntimes,nrepx,nrepp,nreps,
403  sidx,pidx,cidx,&covariate_table,ff,args,REAL(cov));
404 
405  }
406 
407  break;
408 
409  default: {
410 
411  double *ft = REAL(F);
412  int i, n = nvars*nreps*ntimes;
413  for (i = 0; i < n; i++, ft++) *ft = R_NaReal;
414  warn("'skeleton' unspecified: NAs generated.");
415 
416  }
417 
418  }
419 
420  UNPROTECT(nprotect);
421  return F;
422 }
void pomp_skeleton(double *f, const double *x, const double *p, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition: pomp.h:73
void eval_skeleton_native(double *f, double *time, double *x, double *p, int nvars, int npars, int ncovars, int ntimes, int nrepx, int nrepp, int nreps, int *sidx, int *pidx, int *cidx, lookup_table_t *covar_table, pomp_skeleton *fun, SEXP args, double *cov)
Definition: skeleton.c:245
static R_INLINE SEXP ret_array(int nvars, int nreps, int ntimes, SEXP Snames)
Definition: skeleton.c:81
SEXP add_skel_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition: skeleton.c:10
void eval_skeleton_R(double *f, double *time, double *x, double *p, SEXP fn, SEXP args, SEXP Snames, int nvars, int npars, int ncovars, int ntimes, int nrepx, int nrepp, int nreps, lookup_table_t *covar_table, double *cov)
Definition: skeleton.c:94
Here is the call graph for this function:

◆ do_vmeasure()

SEXP do_vmeasure ( SEXP  object,
SEXP  x,
SEXP  times,
SEXP  params,
SEXP  gnsi 
)

Definition at line 99 of file vmeasure.c.

100 {
102  int ntimes, nvars, npars, ncovars, nreps, nrepsx, nrepsp;
103  int nobs = 0, nans = 0;
104  SEXP Snames, Pnames, Cnames, Onames = R_NilValue;
105  SEXP fn, args;
106  SEXP pompfun;
107  SEXP Y = R_NilValue;
108  int *dim;
109  lookup_table_t covariate_table;
110  SEXP cvec;
111  double *cov;
112 
113  PROTECT(times = AS_NUMERIC(times));
114  ntimes = length(times);
115  if (ntimes < 1)
116  err("length('times') = 0, no work to do.");
117 
118  PROTECT(x = as_state_array(x));
119  dim = INTEGER(GET_DIM(x));
120  nvars = dim[0]; nrepsx = dim[1];
121 
122  if (ntimes != dim[2])
123  err("length of 'times' and 3rd dimension of 'x' do not agree.");
124 
125  PROTECT(params = as_matrix(params));
126  dim = INTEGER(GET_DIM(params));
127  npars = dim[0]; nrepsp = dim[1];
128 
129  nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
130 
131  if ((nreps % nrepsp != 0) || (nreps % nrepsx != 0))
132  err("larger number of replicates is not a multiple of smaller.");
133 
134  PROTECT(pompfun = GET_SLOT(object,install("vmeasure")));
135 
136  PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x)));
137  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
138  PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
139  PROTECT(Onames = GET_SLOT(pompfun,install("obsnames")));
140 
141  // set up the covariate table
142  covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
143  PROTECT(cvec = NEW_NUMERIC(ncovars));
144  cov = REAL(cvec);
145 
146  // extract the user-defined function
147  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,Onames,Cnames));
148 
149  // extract 'userdata' as pairlist
150  PROTECT(args = GET_SLOT(object,install("userdata")));
151 
152  int nprotect = 11;
153  int first = 1;
154 
155  // first do setup
156  switch (mode) {
157 
158  case Rfun: {
159  double *ys, *yt = 0;
160  double *time = REAL(times), *xs = REAL(x), *ps = REAL(params);
161  SEXP ans;
162  int j, k;
163 
164  PROTECT(args = add_args(args,Snames,Pnames,Cnames)); nprotect++;
165 
166  for (k = 0; k < ntimes; k++, time++) { // loop over times
167 
168  R_CheckUserInterrupt(); // check for user interrupt
169 
170  table_lookup(&covariate_table,*time,cov); // interpolate the covariates
171 
172  for (j = 0; j < nreps; j++) { // loop over replicates
173 
174  if (first) {
175 
176  PROTECT(
177  ans = eval_call(
178  fn,args,
179  time,
180  xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
181  ps+npars*(j%nrepsp),npars,
182  cov,ncovars
183  )
184  );
185 
186  nans = LENGTH(ans);
187 
188  nobs = (int) sqrt(nans);
189  if (nans != nobs*nobs)
190  err("'vmeasure' must return a symmetric square matrix.");
191 
192  PROTECT(Onames = GET_ROWNAMES(GET_DIMNAMES(ans)));
193  if (invalid_names(Onames))
194  err("'vmeasure' must return a matrix with row-names.");
195 
196  PROTECT(Y = ret_array(nobs,nreps,ntimes,Onames));
197 
198  nprotect += 3;
199 
200  yt = REAL(Y);
201  ys = REAL(AS_NUMERIC(ans));
202 
203  memcpy(yt,ys,nans*sizeof(double));
204  yt += nans;
205 
206  first = 0;
207 
208  } else {
209 
210  PROTECT(
211  ans = eval_call(
212  fn,args,
213  time,
214  xs+nvars*((j%nrepsx)+nrepsx*k),nvars,
215  ps+npars*(j%nrepsp),npars,
216  cov,ncovars
217  )
218  );
219 
220  if (LENGTH(ans) != nans)
221  err("'vmeasure' returns variable-length results.");
222 
223  ys = REAL(AS_NUMERIC(ans));
224 
225  memcpy(yt,ys,nans*sizeof(double));
226  yt += nans;
227 
228  UNPROTECT(1);
229 
230  }
231 
232  }
233  }
234 
235  }
236 
237  break;
238 
239  case native: case regNative: {
240  double *yt = 0, *xp, *pp;
241  double *time = REAL(times), *xs = REAL(x), *ps = REAL(params);
242  SEXP vmatindex;
243  int *oidx, *sidx, *pidx, *cidx, *vmidx;
244  pomp_vmeasure *ff = NULL;
245  int j, k;
246 
247  nobs = LENGTH(Onames);
248  nans = nobs*nobs;
249  // extract observable, state, parameter covariate indices
250  sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
251  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
252  oidx = INTEGER(GET_SLOT(pompfun,install("obsindex")));
253  cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
254 
255  PROTECT(vmatindex = NEW_INTEGER(nans)); nprotect++;
256  vmidx = INTEGER(vmatindex);
257  for (k = 0; k < nobs; k++) {
258  for (j = 0; j < nobs; j++) {
259  vmidx[j+nobs*k] = oidx[j]+nobs*oidx[k];
260  }
261  }
262 
263  // address of native routine
264  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
265 
266  PROTECT(Y = ret_array(nobs,nreps,ntimes,Onames)); nprotect++;
267  yt = REAL(Y);
268 
269  for (k = 0; k < ntimes; k++, time++) { // loop over times
270 
271  R_CheckUserInterrupt(); // check for user interrupt
272 
273  // interpolate the covar functions for the covariates
274  table_lookup(&covariate_table,*time,cov);
275 
276  for (j = 0; j < nreps; j++, yt += nans) { // loop over replicates
277 
278  xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
279  pp = &ps[npars*(j%nrepsp)];
280 
281  (*ff)(yt,xp,pp,vmidx,sidx,pidx,cidx,cov,*time);
282 
283  }
284  }
285 
286  }
287 
288  break;
289 
290  default: {
291  nobs = LENGTH(Onames);
292  int dim[4] = {nobs, nobs, nreps, ntimes};
293  const char *dimnm[4] = {"var1","var2",".id","time"};
294  double *yt = 0;
295  int i, n = nobs*nobs*nreps*ntimes;
296 
297  PROTECT(Y = makearray(4,dim)); nprotect++;
298  setrownames(Y,Onames,4);
299  setcolnames(Y,Onames);
300  fixdimnames(Y,dimnm,4);
301 
302  for (i = 0, yt = REAL(Y); i < n; i++, yt++) *yt = R_NaReal;
303 
304  warn("'vmeasure' unspecified: NAs generated.");
305  }
306 
307  }
308 
309  UNPROTECT(nprotect);
310  return Y;
311 }
void pomp_vmeasure(double *f, const double *x, const double *p, const int *vmatindex, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition: pomp.h:100
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition: vmeasure.c:60
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition: vmeasure.c:10
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition: vmeasure.c:84
Here is the call graph for this function:

◆ euler_simulator()

SEXP euler_simulator ( SEXP  func,
SEXP  xstart,
SEXP  tstart,
SEXP  times,
SEXP  params,
double  deltat,
rprocmode  method,
SEXP  accumvars,
SEXP  covar,
SEXP  args,
SEXP  gnsi 
)

Definition at line 107 of file euler.c.

112  {
113 
115  int nvars, npars, nreps, ntimes, nzeros, ncovars;
116  double *cov, t0;
117  SEXP cvec, X, fn;
118 
119  if (deltat <= 0) err("'delta.t' should be a positive number."); // #nocov
120 
121  int *dim;
122  dim = INTEGER(GET_DIM(xstart)); nvars = dim[0]; nreps = dim[1];
123  dim = INTEGER(GET_DIM(params)); npars = dim[0];
124  ntimes = LENGTH(times);
125 
126  PROTECT(tstart = AS_NUMERIC(tstart));
127  PROTECT(times = AS_NUMERIC(times));
128  t0 = *(REAL(tstart));
129  if (t0 > *(REAL(times))) err("'t0' must be no later than 'times[1]'.");
130 
131  SEXP Snames, Pnames, Cnames;
132  PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(xstart)));
133  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
134  PROTECT(Cnames = get_covariate_names(covar));
135 
136  // set up the covariate table
137  lookup_table_t covariate_table = make_covariate_table(covar,&ncovars);
138  PROTECT(cvec = NEW_NUMERIC(ncovars));
139  cov = REAL(cvec);
140 
141  // indices of accumulator variables
142  nzeros = LENGTH(accumvars);
143  int *zidx = INTEGER(PROTECT(matchnames(Snames,accumvars,"state variables")));
144 
145  // extract user function
146  PROTECT(fn = pomp_fun_handler(func,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
147 
148  // array to hold results
149  PROTECT(X = ret_array(nvars,nreps,ntimes,Snames));
150 
151  // copy the start values into the result array
152  memcpy(REAL(X),REAL(xstart),nvars*nreps*sizeof(double));
153 
154  // set up
155 
156  int nprotect = 9;
157  int *pidx = 0, *sidx = 0, *cidx = 0;
158  pomp_onestep_sim *ff = NULL;
159 
160  switch (mode) {
161 
162  case Rfun:
163  // construct list of all arguments
164  PROTECT(args = add_args(args,Snames,Pnames,Cnames)); nprotect++;
165  break;
166 
167  case native: case regNative:
168  // construct state, parameter, covariate indices
169  sidx = INTEGER(GET_SLOT(func,install("stateindex")));
170  pidx = INTEGER(GET_SLOT(func,install("paramindex")));
171  cidx = INTEGER(GET_SLOT(func,install("covarindex")));
172 
173  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
174 
175  GetRNGstate();
176  break;
177 
178  default: // #nocov
179  err("unrecognized 'mode' %d",mode); // #nocov
180  break; // #nocov
181  }
182 
183  // main computation loop
184  int step;
185  double *xt, *time, t;
186  int first = 1;
187 
188  for (
189  step = 0, xt = REAL(X), time = REAL(times), t = t0;
190  step < ntimes;
191  step++, xt += nvars*nreps
192  ) {
193 
194  double dt;
195  int nstep = 0;
196  int i, j, k;
197 
198  R_CheckUserInterrupt();
199 
200  if (t > time[step]) err("'times' must be an increasing sequence."); // #nocov
201 
202  // set accumulator variables to zero
203  for (j = 0; j < nreps; j++)
204  for (i = 0; i < nzeros; i++)
205  xt[zidx[i]+nvars*j] = 0.0;
206 
207  // determine size and number of time-steps
208  switch (method) {
209  case onestep: default: // one step
210  dt = time[step]-t;
211  nstep = 1;
212  break;
213  case discrete: // fixed step
214  dt = deltat;
215  nstep = num_map_steps(t,time[step],dt);
216  break;
217  case euler: // Euler method
218  dt = deltat;
219  nstep = num_euler_steps(t,time[step],&dt);
220  break;
221  }
222 
223  // loop over individual steps
224  for (k = 0; k < nstep; k++) {
225 
226  // interpolate the covar functions for the covariates
227  table_lookup(&covariate_table,t,cov);
228 
229  // loop over replicates
230  double *ap, *pm, *xm, *ps = REAL(params);
231 
232  for (j = 0, pm = ps, xm = xt; j < nreps; j++, pm += npars, xm += nvars) {
233 
234  switch (mode) {
235 
236  case Rfun:
237  {
238  SEXP ans, nm;
239 
240  if (first) {
241 
242  PROTECT(ans = eval_call(fn,args,&t,&dt,xm,nvars,pm,npars,cov,ncovars));
243 
244  PROTECT(nm = GET_NAMES(ans));
245  if (invalid_names(nm))
246  err("'rprocess' must return a named numeric vector.");
247  pidx = INTEGER(PROTECT(matchnames(nm,Snames,"state variables")));
248 
249  nprotect += 3;
250 
251  ap = REAL(AS_NUMERIC(ans));
252  for (i = 0; i < nvars; i++) xm[pidx[i]] = ap[i];
253 
254  first = 0;
255 
256  } else {
257 
258  PROTECT(ans = eval_call(fn,args,&t,&dt,xm,nvars,pm,npars,cov,ncovars));
259 
260  ap = REAL(AS_NUMERIC(ans));
261  for (i = 0; i < nvars; i++) xm[pidx[i]] = ap[i];
262 
263  UNPROTECT(1);
264 
265  }
266  }
267  break;
268 
269  case native: case regNative:
270  (*ff)(xm,pm,sidx,pidx,cidx,cov,t,dt);
271  break;
272 
273  default: // #nocov
274  err("unrecognized 'mode' %d",mode); // #nocov
275  break; // #nocov
276 
277  }
278 
279  }
280 
281  t += dt;
282 
283  if ((method == euler) && (k == nstep-2)) { // penultimate step
284  dt = time[step]-t;
285  t = time[step]-dt;
286  }
287 
288  }
289 
290  if (step < ntimes-1)
291  memcpy(xt+nvars*nreps,xt,nreps*nvars*sizeof(double));
292 
293  }
294 
295  // clean up
296  switch (mode) {
297 
298  case native: case regNative: {
299 
300  PutRNGstate();
301 
302  }
303 
304  break;
305 
306  case Rfun: default:
307 
308  break;
309 
310  }
311 
312  UNPROTECT(nprotect);
313  return X;
314 
315 }
int num_map_steps(double t1, double t2, double deltat)
Definition: euler.c:342
int num_euler_steps(double t1, double t2, double *deltat)
Definition: euler.c:317
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *dt, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition: euler.c:67
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition: euler.c:10
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition: euler.c:92
void pomp_onestep_sim(double *x, const double *p, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t, double dt)
Definition: pomp.h:59
Here is the call graph for this function:
Here is the caller graph for this function:

◆ eval_skeleton_native()

void eval_skeleton_native ( double *  f,
double *  time,
double *  x,
double *  p,
int  nvars,
int  npars,
int  ncovars,
int  ntimes,
int  nrepx,
int  nrepp,
int  nreps,
int *  sidx,
int *  pidx,
int *  cidx,
lookup_table_t covar_table,
pomp_skeleton fun,
SEXP  args,
double *  cov 
)

Definition at line 245 of file skeleton.c.

252 {
253  double *xp, *pp;
254  int j, k;
255 
256  for (k = 0; k < ntimes; k++, time++) {
257 
258  R_CheckUserInterrupt();
259 
261 
262  for (j = 0; j < nreps; j++, f += nvars) {
263 
264  xp = &x[nvars*((j%nrepx)+nrepx*k)];
265  pp = &p[npars*(j%nrepp)];
266 
267  (*fun)(f,xp,pp,sidx,pidx,cidx,cov,*time);
268 
269  }
270  }
271 
272 }
lookup_table_t covar_table
Definition: trajectory.c:130
Here is the call graph for this function:
Here is the caller graph for this function:

◆ eval_skeleton_R()

void eval_skeleton_R ( double *  f,
double *  time,
double *  x,
double *  p,
SEXP  fn,
SEXP  args,
SEXP  Snames,
int  nvars,
int  npars,
int  ncovars,
int  ntimes,
int  nrepx,
int  nrepp,
int  nreps,
lookup_table_t covar_table,
double *  cov 
)

Definition at line 93 of file skeleton.c.

100  {
101 
102  SEXP ans, nm;
103  double *fs;
104  int *posn = 0;
105  int i, j, k;
106  int first = 1;
107  int nprotect = 0;
108 
109  for (k = 0; k < ntimes; k++, time++) {
110 
111  R_CheckUserInterrupt();
112 
113  // interpolate the covar functions for the covariates
115 
116  for (j = 0; j < nreps; j++, f += nvars) {
117 
118  if (first) {
119 
120  PROTECT(ans = eval_call(fn,args,time,
121  x+nvars*((j%nrepx)+nrepx*k),nvars,
122  p+npars*(j%nrepp),npars,
123  cov,ncovars));
124 
125  if (LENGTH(ans)!=nvars)
126  err("'skeleton' returns a vector of %d state variables but %d are expected.",LENGTH(ans),nvars);
127 
128  // get name information to fix alignment problems
129  PROTECT(nm = GET_NAMES(ans));
130  if (invalid_names(nm))
131  err("'skeleton' must return a named numeric vector.");
132  posn = INTEGER(PROTECT(matchnames(Snames,nm,"state variables")));
133 
134  fs = REAL(AS_NUMERIC(ans));
135  for (i = 0; i < nvars; i++) f[posn[i]] = fs[i];
136 
137  nprotect += 3;
138  first = 0;
139 
140  } else {
141 
142  PROTECT(ans = eval_call(fn,args,time,
143  x+nvars*((j%nrepx)+nrepx*k),nvars,
144  p+npars*(j%nrepp),npars,
145  cov,ncovars));
146 
147  fs = REAL(AS_NUMERIC(ans));
148  for (i = 0; i < nvars; i++) f[posn[i]] = fs[i];
149 
150  UNPROTECT(1);
151 
152  }
153 
154  }
155  }
156 
157  UNPROTECT(nprotect);
158 
159 }
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition: skeleton.c:57
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ExpitTransform()

SEXP ExpitTransform ( SEXP  X)

Definition at line 26 of file transformations.c.

26  {
27  double *x;
28  int k, n;
29  PROTECT(X = duplicate(AS_NUMERIC(X)));
30  x = REAL(X);
31  n = LENGTH(X);
32  for (k = 0; k < n; k++, x++)
33  *x = 1.0/(1.0+exp(-(*x)));
34  UNPROTECT(1);
35  return X;
36 }

◆ get_covariate_names()

SEXP get_covariate_names ( SEXP  object)

Definition at line 7 of file lookup_table.c.

7  {
8  return GET_ROWNAMES(GET_DIMNAMES(GET_SLOT(object,install("table"))));
9 }
Here is the caller graph for this function:

◆ get_userdata()

const SEXP get_userdata ( const char *  name)

Definition at line 13 of file userdata.c.

13  {
14  SEXP elt = getListElement(USERDATA,name);
15  if (isNull(elt)) err("no user-data element '%s' is found.",name);
16  return elt;
17 }
static R_INLINE SEXP getListElement(SEXP list, const char *str)
Definition: pomp_defines.h:232
#define USERDATA
Definition: userdata.c:6
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_userdata_double()

const double* get_userdata_double ( const char *  name)

Definition at line 26 of file userdata.c.

26  {
27  SEXP elt = getListElement(USERDATA,name);
28  if (isNull(elt)) err("no user-data element '%s' is found.",name);
29  if (!isReal(elt)) err("user-data element '%s' is not a numeric vector.",name);
30  return REAL(elt);
31 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_userdata_int()

const int* get_userdata_int ( const char *  name)

Definition at line 19 of file userdata.c.

19  {
20  SEXP elt = getListElement(USERDATA,name);
21  if (isNull(elt)) err("no user-data element '%s' is found.",name);
22  if (!isInteger(elt)) err("user-data element '%s' is not an integer.",name);
23  return INTEGER(elt);
24 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InverseLogBarycentricTransform()

SEXP InverseLogBarycentricTransform ( SEXP  Y)

Definition at line 47 of file transformations.c.

47  {
48  SEXP X;
49  PROTECT(Y = AS_NUMERIC(Y));
50  PROTECT(X = NEW_NUMERIC(LENGTH(Y)));
51  from_log_barycentric(REAL(X),REAL(Y),LENGTH(Y));
52  UNPROTECT(2);
53  return X;
54 }
static R_INLINE void from_log_barycentric(double *xt, const double *x, int n)
Definition: pomp.h:260
Here is the call graph for this function:

◆ iterate_map()

SEXP iterate_map ( SEXP  object,
SEXP  times,
SEXP  t0,
SEXP  x0,
SEXP  params,
SEXP  gnsi 
)

Definition at line 19 of file trajectory.c.

20 {
21 
23  SEXP fn, args;
24  SEXP X, cov;
25  SEXP Snames, Pnames, Cnames;
26  SEXP skel, pompfun;
27  SEXP accumvars;
28  SEXP repnames;
29  int *zidx = 0;
30  int nvars, npars, nreps, nrepp, ntimes, ncovars, nzeros;
31  int *dim;
32  lookup_table_t covariate_table;
33  double deltat, t;
34 
35  PROTECT(skel = GET_SLOT(object,install("skeleton")));
36  deltat = *(REAL(GET_SLOT(skel,install("delta.t"))));
37  t = *(REAL(AS_NUMERIC(t0)));
38 
39  PROTECT(x0 = duplicate(x0));
40  PROTECT(x0 = as_matrix(x0));
41  dim = INTEGER(GET_DIM(x0));
42  nvars = dim[0]; nreps = dim[1];
43 
44  PROTECT(params = as_matrix(params));
45  dim = INTEGER(GET_DIM(params));
46  npars = dim[0]; nrepp = dim[1];
47  PROTECT(repnames = GET_COLNAMES(GET_DIMNAMES(params)));
48 
49  PROTECT(times = AS_NUMERIC(times));
50  ntimes = LENGTH(times);
51 
52  PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x0)));
53  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
54  PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
55 
56  // set up the covariate table
57  covariate_table = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
58  PROTECT(cov = NEW_NUMERIC(ncovars));
59 
60  // extract user-defined function
61  PROTECT(pompfun = GET_SLOT(skel,install("skel.fn")));
62  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
63 
64  // extract 'userdata' as pairlist
65  PROTECT(args = GET_SLOT(object,install("userdata")));
66 
67  // create array to store results
68  PROTECT(X = ret_array(nvars,nreps,ntimes,Snames,repnames));
69 
70  // get names and indices of accumulator variables
71  PROTECT(accumvars = GET_SLOT(object,install("accumvars")));
72  nzeros = LENGTH(accumvars);
73 
74  int nprotect = 15;
75 
76  if (nzeros > 0) {
77  zidx = INTEGER(PROTECT(matchnames(Snames,accumvars,"state variables"))); nprotect++;
78  }
79 
80  // set up the computations
81  switch (mode) {
82 
83  case Rfun: {
84 
85  PROTECT(args = add_skel_args(args,Snames,Pnames,Cnames)); nprotect++;
86 
87  iterate_skeleton_R(REAL(X),t,deltat,REAL(times),REAL(x0),REAL(params),
88  fn,args,Snames,nvars,npars,ncovars,ntimes,nrepp,nreps,nzeros,
89  &covariate_table,zidx,REAL(cov));
90 
91  }
92 
93  break;
94 
95  case native: case regNative: {
96  int *sidx, *pidx, *cidx;
97  pomp_skeleton *ff;
98 
99  *((void **) (&ff)) = R_ExternalPtrAddr(fn);
100 
101  // construct state, parameter, covariate indices
102  sidx = INTEGER(GET_SLOT(pompfun,install("stateindex")));
103  pidx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
104  cidx = INTEGER(GET_SLOT(pompfun,install("covarindex")));
105 
106  iterate_skeleton_native(REAL(X),t,deltat,REAL(times),REAL(x0),REAL(params),
107  nvars,npars,ncovars,ntimes,nrepp,nreps,nzeros,sidx,pidx,cidx,
108  &covariate_table,zidx,ff,args,REAL(cov));
109 
110  }
111 
112  break;
113 
114  default: // #nocov
115 
116  break; // #nocov
117 
118  }
119 
120  UNPROTECT(nprotect);
121  return X;
122 }
SEXP add_skel_args(SEXP, SEXP, SEXP, SEXP)
Definition: skeleton.c:10
void iterate_skeleton_R(double *, double, double, double *, double *, double *, SEXP, SEXP, SEXP, int, int, int, int, int, int, int, lookup_table_t *, int *, double *)
Definition: skeleton.c:161
void iterate_skeleton_native(double *, double, double, double *, double *, double *, int, int, int, int, int, int, int, int *, int *, int *, lookup_table_t *, int *, pomp_skeleton *, SEXP, double *)
Definition: skeleton.c:274
static R_INLINE SEXP ret_array(int nvars, int nreps, int ntimes, SEXP Snames, SEXP repnames)
Definition: trajectory.c:6
Here is the call graph for this function:

◆ iterate_skeleton_native()

void iterate_skeleton_native ( double *  X,
double  t,
double  deltat,
double *  time,
double *  x,
double *  p,
int  nvars,
int  npars,
int  ncovars,
int  ntimes,
int  nrepp,
int  nreps,
int  nzeros,
int *  sidx,
int *  pidx,
int *  cidx,
lookup_table_t covar_table,
int *  zeroindex,
pomp_skeleton fun,
SEXP  args,
double *  cov 
)

Definition at line 274 of file skeleton.c.

282 {
283  int nsteps;
284  double *xs, *Xs;
285  int h, i, j, k;
286 
287  for (k = 0; k < ntimes; k++, time++, X += nvars*nreps) {
288 
289  R_CheckUserInterrupt();
290 
291  // compute number of steps needed
292  nsteps = num_map_steps(t,*time,deltat);
293 
294  // set accumulator variables to zero
295  for (i = 0; i < nzeros; i++)
296  for (j = 0, xs = &x[zeroindex[i]]; j < nreps; j++, xs += nvars)
297  *xs = 0.0;
298 
299  for (h = 0; h < nsteps; h++) {
300 
301  // interpolate the covariates
303 
304  for (j = 0, Xs = X, xs = x; j < nreps; j++, Xs += nvars, xs += nvars) {
305 
306  (*fun)(Xs,xs,p+npars*(j%nrepp),sidx,pidx,cidx,cov,t);
307 
308  }
309 
310  memcpy(x,X,nvars*nreps*sizeof(double));
311 
312  if (h != nsteps-1) {
313  t += deltat;
314  } else {
315  t = *time;
316  }
317 
318  }
319 
320  if (nsteps == 0) memcpy(X,x,nvars*nreps*sizeof(double));
321 
322  }
323 
324 }
int num_map_steps(double, double, double)
Definition: euler.c:342
Here is the call graph for this function:
Here is the caller graph for this function:

◆ iterate_skeleton_R()

void iterate_skeleton_R ( double *  X,
double  t,
double  deltat,
double *  time,
double *  x,
double *  p,
SEXP  fn,
SEXP  args,
SEXP  Snames,
int  nvars,
int  npars,
int  ncovars,
int  ntimes,
int  nrepp,
int  nreps,
int  nzeros,
lookup_table_t covar_table,
int *  zeroindex,
double *  cov 
)

Definition at line 161 of file skeleton.c.

169 {
170 
171  SEXP ans, nm;
172  double *ap, *xs;
173  int nsteps;
174  int *posn = 0;
175  int h, i, j, k;
176  int first = 1;
177  int nprotect = 0;
178 
179  for (k = 0; k < ntimes; k++, time++, X += nvars*nreps) {
180 
181  R_CheckUserInterrupt();
182 
183  // compute number of steps needed
184  nsteps = num_map_steps(t,*time,deltat);
185 
186  // set accumulator variables to zero
187  for (i = 0; i < nzeros; i++)
188  for (j = 0, xs = &x[zeroindex[i]]; j < nreps; j++, xs += nvars)
189  *xs = 0.0;
190 
191  for (h = 0; h < nsteps; h++) {
192 
193  // interpolate the covariates
195 
196  for (j = 0, xs = x; j < nreps; j++, xs += nvars) {
197 
198  if (first) {
199 
200  PROTECT(ans = eval_call(fn,args,&t,xs,nvars,p+npars*(j%nrepp),npars,cov,ncovars));
201 
202  if (LENGTH(ans) != nvars)
203  err("'skeleton' returns a vector of %d state variables but %d are expected.",LENGTH(ans),nvars);
204 
205  // get name information to fix alignment problems
206  PROTECT(nm = GET_NAMES(ans));
207  if (invalid_names(nm)) err("'skeleton' must return a named numeric vector.");
208  posn = INTEGER(PROTECT(matchnames(Snames,nm,"state variables")));
209 
210  ap = REAL(AS_NUMERIC(ans));
211  for (i = 0; i < nvars; i++) xs[posn[i]] = ap[i];
212 
213  nprotect += 3;
214  first = 0;
215 
216  } else {
217 
218  PROTECT(ans = eval_call(fn,args,&t,xs,nvars,p+npars*(j%nrepp),npars,cov,ncovars));
219 
220  ap = REAL(AS_NUMERIC(ans));
221  for (i = 0; i < nvars; i++) xs[posn[i]] = ap[i];
222 
223  UNPROTECT(1);
224 
225  }
226 
227  }
228 
229  if (h != nsteps-1) {
230  t += deltat;
231  } else {
232  t = *time;
233  }
234 
235  }
236 
237  memcpy(X,x,nvars*nreps*sizeof(double));
238 
239  }
240 
241  UNPROTECT(nprotect);
242 
243 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ load_stack_decr()

SEXP load_stack_decr ( SEXP  pack)

Definition at line 130 of file pomp_fun.c.

130  {
131  SEXP s;
132  const char *pkg;
133  void (*ff)(int *);
134  PROTECT(s = ScalarInteger(NA_INTEGER));
135  pkg = (const char *) CHAR(STRING_ELT(pack,0));
136  ff = (void (*)(int *)) R_GetCCallable(pkg,"__pomp_load_stack_decr");
137  ff(INTEGER(s));
138  if (*(INTEGER(s)) < 0) err("impossible!");
139  UNPROTECT(1);
140  return s;
141 }
Here is the caller graph for this function:

◆ load_stack_incr()

SEXP load_stack_incr ( SEXP  pack)

Definition at line 121 of file pomp_fun.c.

121  {
122  const char *pkg;
123  void (*ff)(void);
124  pkg = (const char *) CHAR(STRING_ELT(pack,0));
125  ff = (void (*)(void)) R_GetCCallable(pkg,"__pomp_load_stack_incr");
126  ff();
127  return R_NilValue;
128 }
Here is the caller graph for this function:

◆ LogBarycentricTransform()

SEXP LogBarycentricTransform ( SEXP  X)

Definition at line 38 of file transformations.c.

38  {
39  SEXP Y;
40  PROTECT(X = AS_NUMERIC(X));
41  PROTECT(Y = NEW_NUMERIC(LENGTH(X)));
42  to_log_barycentric(REAL(Y),REAL(X),LENGTH(X));
43  UNPROTECT(2);
44  return Y;
45 }
static R_INLINE void to_log_barycentric(double *xt, const double *x, int n)
Definition: pomp.h:250
Here is the call graph for this function:

◆ LogitTransform()

SEXP LogitTransform ( SEXP  P)

Definition at line 14 of file transformations.c.

14  {
15  double *p;
16  int k, n;
17  PROTECT(P = duplicate(AS_NUMERIC(P)));
18  p = REAL(P);
19  n = LENGTH(P);
20  for (k = 0; k < n; k++, p++)
21  *p = logit(*p);
22  UNPROTECT(1);
23  return P;
24 }
static R_INLINE double logit(double p)
Definition: pomp.h:118
Here is the call graph for this function:

◆ logmeanexp()

SEXP logmeanexp ( const  SEXP,
const  SEXP 
)

Definition at line 7 of file logmeanexp.c.

7  {
8  int j, n = LENGTH(X);
9  int k = *INTEGER(Drop)-1; // zero-based index
10  double *x = REAL(X);
11  long double m = R_NegInf;
12  long double s = 0;
13  for (j = 0; j < n; j++) {
14  if (j != k)
15  m = (x[j] > m) ? (long double) x[j] : m;
16  }
17  for (j = 0; j < n; j++) {
18  if (j != k)
19  s += expl((long double) x[j] - m);
20  }
21  if (k >= 0 && k < n) n--;
22  return ScalarReal(m + log(s/n));
23 }

◆ lookup_in_table()

SEXP lookup_in_table ( SEXP  covar,
SEXP  t 
)

Definition at line 24 of file lookup_table.c.

24  {
25  int xdim[2], nvar;
26  int j, nt;
27  double *tp, *xp;
28  SEXP Cnames, X;
29 
30  PROTECT(t = AS_NUMERIC(t));
31  nt = LENGTH(t);
32  PROTECT(Cnames = get_covariate_names(covar));
33 
34  lookup_table_t tab = make_covariate_table(covar,&nvar);
35 
36  if (nt > 1) {
37  xdim[0] = nvar; xdim[1] = nt;
38  PROTECT(X = makearray(2,xdim));
39  setrownames(X,Cnames,2);
40  } else {
41  PROTECT(X = NEW_NUMERIC(nvar));
42  SET_NAMES(X,Cnames);
43  }
44 
45  for (j = 0, tp = REAL(t), xp = REAL(X); j < nt; j++, tp++, xp += nvar)
46  table_lookup(&tab,*tp,xp);
47 
48  UNPROTECT(3);
49  return X;
50 }
void table_lookup(lookup_table_t *tab, double x, double *y)
Definition: lookup_table.c:53
lookup_table_t make_covariate_table(SEXP object, int *ncovar)
Definition: lookup_table.c:11
SEXP get_covariate_names(SEXP object)
Definition: lookup_table.c:7
Here is the call graph for this function:
Here is the caller graph for this function:

◆ make_covariate_table()

lookup_table_t make_covariate_table ( SEXP  object,
int *  ncovar 
)

Definition at line 11 of file lookup_table.c.

11  {
12  lookup_table_t tab;
13  int *dim;
14  dim = INTEGER(GET_DIM(GET_SLOT(object,install("table"))));
15  *ncovar = tab.width = dim[0];
16  tab.length = dim[1];
17  tab.index = 0;
18  tab.x = REAL(GET_SLOT(object,install("times")));
19  tab.y = REAL(GET_SLOT(object,install("table")));
20  tab.order = *(INTEGER(GET_SLOT(object,install("order"))));
21  return tab;
22 }
Here is the caller graph for this function:

◆ nosort_resamp()

void nosort_resamp ( int  nw,
double *  w,
int  np,
int *  p,
int  offset 
)

Definition at line 25 of file resample.c.

26 {
27  int i, j;
28  double du, u;
29 
30  for (j = 1; j < nw; j++) w[j] += w[j-1];
31 
32  if (w[nw-1] <= 0.0)
33  err("in 'systematic_resampling': non-positive sum of weights");
34 
35  du = w[nw-1] / ((double) np);
36  u = -du*unif_rand();
37 
38  for (i = 0, j = 0; j < np; j++) {
39  u += du;
40  // In the following line, the second test is needed to correct
41  // the infamous Bug of St. Patrick, 2017-03-17.
42  while ((u > w[i]) && (i < nw-1)) i++;
43  p[j] = i;
44  }
45  if (offset) // add offset if needed
46  for (j = 0; j < np; j++) p[j] += offset;
47 
48 }
Here is the caller graph for this function:

◆ num_euler_steps()

int num_euler_steps ( double  t1,
double  t2,
double *  deltat 
)

Definition at line 317 of file euler.c.

317  {
318  double tol = sqrt(DBL_EPSILON);
319  int nstep;
320  // nstep will be the number of Euler steps to take in going from t1 to t2.
321  // note also that the stepsize changes.
322  // this choice is meant to be conservative
323  // (i.e., so that the actual deltat does not exceed the specified deltat
324  // by more than the relative tolerance 'tol')
325  // and to counteract roundoff error.
326  // It seems to work well, but is not guaranteed:
327  // suggestions would be appreciated.
328 
329  if (t1 >= t2) {
330  *deltat = 0.0;
331  nstep = 0;
332  } else if (t1+*deltat >= t2) {
333  *deltat = t2-t1;
334  nstep = 1;
335  } else {
336  nstep = (int) ceil((t2-t1)/(*deltat)/(1+tol));
337  *deltat = (t2-t1)/((double) nstep);
338  }
339  return nstep;
340 }
Here is the caller graph for this function:

◆ num_map_steps()

int num_map_steps ( double  t1,
double  t2,
double  deltat 
)

Definition at line 342 of file euler.c.

342  {
343  double tol = sqrt(DBL_EPSILON);
344  int nstep;
345  // nstep will be the number of discrete-time steps to take in going from t1 to t2.
346  nstep = (int) floor((t2-t1)/deltat/(1-tol));
347  return (nstep > 0) ? nstep : 0;
348 }
Here is the caller graph for this function:

◆ periodic_bspline_basis()

SEXP periodic_bspline_basis ( SEXP  x,
SEXP  nbasis,
SEXP  degree,
SEXP  period,
SEXP  deriv 
)

Definition at line 85 of file bspline.c.

86  {
87 
88  SEXP y, xr;
89  int nx = LENGTH(x);
90  int nb = INTEGER_VALUE(nbasis);
91  int deg = INTEGER_VALUE(degree);
92  int d = INTEGER_VALUE(deriv);
93  double pd = NUMERIC_VALUE(period);
94  int j, k;
95  double *xrd, *ydata, *val;
96  PROTECT(xr = AS_NUMERIC(x));
97  xrd = REAL(xr);
98  PROTECT(y = allocMatrix(REALSXP,nx,nb));
99  ydata = REAL(y);
100  val = (double *) R_alloc(nb,sizeof(double));
101  for (j = 0; j < nx; j++) {
102  periodic_bspline_basis_eval_deriv(xrd[j],pd,deg,nb,d,val);
103  for (k = 0; k < nb; k++) ydata[j+nx*k] = val[k];
104  }
105  UNPROTECT(2);
106  return y;
107 }
void periodic_bspline_basis_eval_deriv(double x, double period, int degree, int nbasis, int deriv, double *y)
Definition: bspline.c:122
Here is the call graph for this function:

◆ periodic_bspline_basis_eval_deriv()

void periodic_bspline_basis_eval_deriv ( double  x,
double  period,
int  degree,
int  nbasis,
int  deriv,
double *  y 
)

Definition at line 122 of file bspline.c.

124 {
125  int nknots = nbasis+2*degree+1;
126  int shift = (degree-1)/2;
127  double *knots = NULL, *yy = NULL;
128  double dx;
129  int j, k;
130 
131  if (period <= 0.0) err("must have period > 0");
132  if (nbasis <= 0) err("must have nbasis > 0");
133  if (degree < 0) err("must have degree >= 0");
134  if (nbasis < degree) err("must have nbasis >= degree");
135  if (deriv < 0) err("must have deriv >= 0");
136 
137  knots = (double *) R_Calloc(nknots+degree+1,double);
138  yy = (double *) R_Calloc(nknots,double);
139 
140  dx = period/((double) nbasis);
141  for (k = -degree; k <= nbasis+degree; k++) {
142  knots[degree+k] = k*dx;
143  }
144  x = fmod(x,period);
145  if (x < 0.0) x += period;
146  for (k = 0; k < nknots; k++) {
147  bspline_eval(&yy[k],&x,1,k,degree,deriv,knots);
148  }
149  for (k = 0; k < degree; k++) yy[k] += yy[nbasis+k];
150  for (k = 0; k < nbasis; k++) {
151  j = (shift+k)%nbasis;
152  y[k] = yy[j];
153  }
154  R_Free(yy); R_Free(knots);
155 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ pfilter()

SEXP pfilter ( SEXP  x,
SEXP  params,
SEXP  Np,
SEXP  predmean,
SEXP  predvar,
SEXP  filtmean,
SEXP  trackancestry,
SEXP  doparRS,
SEXP  weights,
SEXP  wave 
)

Definition at line 16 of file pfilter.c.

20 {
21 
22  SEXP pm = R_NilValue, pv = R_NilValue, fm = R_NilValue;
23  SEXP anc = R_NilValue, wmean = R_NilValue;
24  SEXP ess, loglik;
25  SEXP newstates = R_NilValue, newparams = R_NilValue;
26  SEXP retval, retvalnames;
27  const char *dimnm[2] = {"name",".id"};
28  double *xw = 0;
29  long double w = 0, ws, maxw, sum;
30  int *xanc = 0;
31  SEXP dimX, dimP, newdim, Xnames, Pnames;
32  int *dim, np;
33  int nvars, npars = 0, nreps;
34  int do_pm, do_pv, do_fm, do_ta, do_pr, do_wave, all_fail = 0;
35  int j, k;
36 
37  PROTECT(dimX = GET_DIM(x));
38  dim = INTEGER(dimX);
39  nvars = dim[0]; nreps = dim[1];
40  PROTECT(Xnames = GET_ROWNAMES(GET_DIMNAMES(x)));
41 
42  PROTECT(params = as_matrix(params));
43  PROTECT(dimP = GET_DIM(params));
44  dim = INTEGER(dimP);
45  npars = dim[0];
46  if (nreps % dim[1] != 0)
47  err("ncol('states') should be a multiple of ncol('params')"); // # nocov
48  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
49 
50  PROTECT(weights = duplicate(weights));
51 
52  np = *(INTEGER(AS_INTEGER(Np))); // number of particles to resample
53 
54  do_pm = *(LOGICAL(AS_LOGICAL(predmean))); // calculate prediction means?
55  do_pv = *(LOGICAL(AS_LOGICAL(predvar))); // calculate prediction variances?
56  do_fm = *(LOGICAL(AS_LOGICAL(filtmean))); // calculate filtering means?
57  do_ta = *(LOGICAL(AS_LOGICAL(trackancestry))); // track ancestry?
58  // Do we need to do parameter resampling?
59  do_pr = *(LOGICAL(AS_LOGICAL(doparRS)));
60  // Do we need to take a weighted average over the parameters?
61  do_wave = *(LOGICAL(AS_LOGICAL(wave)));
62 
63  if (do_pr) {
64  if (dim[1] != nreps)
65  err("ncol('states') should be equal to ncol('params')"); // # nocov
66  }
67 
68  PROTECT(ess = NEW_NUMERIC(1)); // effective sample size
69  PROTECT(loglik = NEW_NUMERIC(1)); // log likelihood
70 
71  int nprotect = 8;
72 
73  xw = REAL(weights);
74 
75  // check and normalize the weights
76  for (k = 0, maxw = R_NegInf; k < nreps; k++) {
77 
78  if (ISNAN(xw[k]) || xw[k] == R_PosInf) { // check the weights
79  PROTECT(retval = NEW_INTEGER(1)); nprotect++;
80  *INTEGER(retval) = k+1; // return the index of the peccant particle
81  UNPROTECT(nprotect);
82  return retval;
83  }
84 
85  if (maxw < xw[k]) maxw = xw[k];
86 
87  }
88 
89  if (maxw == R_NegInf) all_fail = 1;
90 
91  if (all_fail) {
92 
93  *(REAL(loglik)) = R_NegInf;
94  *(REAL(ess)) = 0; // zero effective sample size
95 
96  } else {
97 
98  // compute sum and sum of squares
99  for (k = 0, w = 0, ws = 0; k < nreps; k++) {
100  xw[k] = exp(xw[k]-maxw);
101  if (xw[k] != 0) {
102  w += xw[k];
103  ws += xw[k]*xw[k];
104  }
105  }
106 
107  *(REAL(loglik)) = maxw + log(w/((long double) nreps)); // mean of weights is likelihood
108  *(REAL(ess)) = w*w/ws; // effective sample size
109  }
110 
111  if (do_pm || do_pv) {
112  PROTECT(pm = NEW_NUMERIC(nvars)); nprotect++;
113  }
114 
115  if (do_pv) {
116  PROTECT(pv = NEW_NUMERIC(nvars)); nprotect++;
117  }
118 
119  if (do_fm) {
120  PROTECT(fm = NEW_NUMERIC(nvars)); nprotect++;
121  }
122 
123  if (do_ta) {
124  PROTECT(anc = NEW_INTEGER(np)); nprotect++;
125  xanc = INTEGER(anc);
126  }
127 
128  if (do_wave) {
129  PROTECT(wmean = NEW_NUMERIC(npars)); nprotect++;
130  SET_NAMES(wmean,Pnames);
131  }
132 
133  if (do_pm || do_pv) {
134  double *tmp = (do_pv) ? REAL(pv) : 0;
135  pred_mean_var(nvars,nreps,do_pv,REAL(x),REAL(pm),tmp);
136  }
137 
138  // compute filter mean
139  if (do_fm) {
140  filt_mean(nvars,nreps,all_fail,w,xw,REAL(x),REAL(fm));
141  }
142 
143  // compute weighted average of parameters
144  if (do_wave) {
145  if (all_fail)
146  warn("%s %s","filtering failure at last filter iteration:",
147  "using unweighted mean for point estimate.");
148  double *xwm = REAL(wmean);
149  for (j = 0; j < npars; j++, xwm++) {
150  double *xp = REAL(params)+j;
151  if (all_fail) { // unweighted average
152  for (k = 0, sum = 0; k < nreps; k++, xp += npars) sum += *xp;
153  *xwm = sum/((long double) nreps);
154  } else {
155  for (k = 0, sum = 0; k < nreps; k++, xp += npars) {
156  if (xw[k]!=0) sum += xw[k]*(*xp);
157  }
158  *xwm = sum/w;
159  }
160  }
161  }
162 
163  GetRNGstate();
164 
165  if (!all_fail) { // resample the particles unless we have filtering failure
166  int xdim[2];
167  int sample[np];
168  double *ss = 0, *st = 0, *ps = 0, *pt = 0;
169 
170  // create storage for new states
171  xdim[0] = nvars; xdim[1] = np;
172  PROTECT(newstates = makearray(2,xdim)); nprotect++;
173  setrownames(newstates,Xnames,2);
174  fixdimnames(newstates,dimnm,2);
175  ss = REAL(x);
176  st = REAL(newstates);
177 
178  // create storage for new parameters
179  if (do_pr) {
180  xdim[0] = npars; xdim[1] = np;
181  PROTECT(newparams = makearray(2,xdim)); nprotect++;
182  setrownames(newparams,Pnames,2);
183  fixdimnames(newparams,dimnm,2);
184  ps = REAL(params);
185  pt = REAL(newparams);
186  }
187 
188  // resample
189  nosort_resamp(nreps,REAL(weights),np,sample,0);
190  for (k = 0; k < np; k++) { // copy the particles
191  int sp = sample[k];
192  double *xx, *xp;
193  for (j = 0, xx = ss+nvars*sp; j < nvars; j++, st++, xx++)
194  *st = *xx;
195  if (do_pr) {
196  for (j = 0, xp = ps+npars*sp; j < npars; j++, pt++, xp++)
197  *pt = *xp;
198  }
199  if (do_ta) xanc[k] = sp+1;
200  }
201 
202  } else { // don't resample: just drop 3rd dimension in x prior to return
203 
204  PROTECT(newdim = NEW_INTEGER(2)); nprotect++;
205  dim = INTEGER(newdim);
206  dim[0] = nvars; dim[1] = nreps;
207  SET_DIM(x,newdim);
208  setrownames(x,Xnames,2);
209  fixdimnames(x,dimnm,2);
210 
211  if (do_ta)
212  for (k = 0; k < np; k++) xanc[k] = k+1;
213  }
214 
215  PutRNGstate();
216 
217  PROTECT(retval = NEW_LIST(9));
218  PROTECT(retvalnames = NEW_CHARACTER(9));
219  nprotect += 2;
220  SET_STRING_ELT(retvalnames,0,mkChar("loglik"));
221  SET_STRING_ELT(retvalnames,1,mkChar("ess"));
222  SET_STRING_ELT(retvalnames,2,mkChar("states"));
223  SET_STRING_ELT(retvalnames,3,mkChar("params"));
224  SET_STRING_ELT(retvalnames,4,mkChar("pm"));
225  SET_STRING_ELT(retvalnames,5,mkChar("pv"));
226  SET_STRING_ELT(retvalnames,6,mkChar("fm"));
227  SET_STRING_ELT(retvalnames,7,mkChar("ancestry"));
228  SET_STRING_ELT(retvalnames,8,mkChar("wmean"));
229  SET_NAMES(retval,retvalnames);
230 
231  SET_ELEMENT(retval,0,loglik);
232  SET_ELEMENT(retval,1,ess);
233 
234  if (all_fail) {
235  SET_ELEMENT(retval,2,x);
236  } else {
237  SET_ELEMENT(retval,2,newstates);
238  }
239 
240  if (all_fail || !do_pr) {
241  SET_ELEMENT(retval,3,params);
242  } else {
243  SET_ELEMENT(retval,3,newparams);
244  }
245 
246  if (do_pm) {
247  SET_ELEMENT(retval,4,pm);
248  }
249  if (do_pv) {
250  SET_ELEMENT(retval,5,pv);
251  }
252  if (do_fm) {
253  SET_ELEMENT(retval,6,fm);
254  }
255  if (do_ta) {
256  SET_ELEMENT(retval,7,anc);
257  }
258  if (do_wave) {
259  SET_ELEMENT(retval,8,wmean);
260  }
261 
262  UNPROTECT(nprotect);
263  return(retval);
264 }
void nosort_resamp(int, double *, int, int *, int)
Definition: resample.c:25
static void filt_mean(int, int, int, long double, const double *, const double *, double *)
Definition: pfilter.c:294
static void pred_mean_var(int, int, int, const double *, double *, double *)
Definition: pfilter.c:266
Here is the call graph for this function:

◆ pomp_desolve_setup()

SEXP pomp_desolve_setup ( SEXP  object,
SEXP  x0,
SEXP  params,
SEXP  gnsi 
)

Definition at line 156 of file trajectory.c.

156  {
157 
159  SEXP fn, args;
160  SEXP Snames, Pnames, Cnames;
161  SEXP pompfun, ob;
162  int *dim;
163  int nvars, npars, nreps, ncovars;
164 
165  // extract user-defined skeleton function
166  PROTECT(ob = GET_SLOT(object,install("skeleton")));
167  PROTECT(pompfun = GET_SLOT(ob,install("skel.fn")));
168  // extract 'userdata' as pairlist
169  PROTECT(args = GET_SLOT(object,install("userdata")));
170 
171  COMMON(object) = object;
172  COMMON(params) = params;
173  if (!isNull(COMMON(object))) R_ReleaseObject(COMMON(object));
174  if (!isNull(COMMON(params))) R_ReleaseObject(COMMON(params));
175  R_PreserveObject(COMMON(object));
176  R_PreserveObject(COMMON(params));
177 
178  dim = INTEGER(GET_DIM(x0));
179  nvars = dim[0];
180  COMMON(nvars) = nvars;
181 
182  dim = INTEGER(GET_DIM(params));
183  npars = dim[0]; nreps = dim[1];
184  COMMON(npars) = npars;
185  COMMON(nreps) = nreps;
186 
187  // set up the covariate table
188  COMMON(covar_table) = make_covariate_table(GET_SLOT(object,install("covar")),&ncovars);
190  PROTECT(COMMON(cov) = NEW_NUMERIC(ncovars));
191  R_PreserveObject(COMMON(cov));
192 
193  PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x0)));
194  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
195  PROTECT(Cnames = get_covariate_names(GET_SLOT(object,install("covar"))));
196 
197  PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode,Snames,Pnames,NA_STRING,Cnames));
198 
199  COMMON(mode) = mode;
200 
201  int nprotect = 8;
202 
203  switch (COMMON(mode)) {
204 
205  case Rfun: {
206 
207  PROTECT(RFUN(fn) = fn);
208  PROTECT(RFUN(args) = add_skel_args(args,Snames,Pnames,Cnames));
209  PROTECT(RFUN(Snames) = Snames);
210  nprotect += 3;
211 
212  if (!isNull(RFUN(fn))) R_ReleaseObject(RFUN(fn));
213  if (!isNull(RFUN(args))) R_ReleaseObject(RFUN(args));
214  if (!isNull(RFUN(Snames))) R_ReleaseObject(RFUN(Snames));
215 
216  R_PreserveObject(RFUN(fn));
217  R_PreserveObject(RFUN(args));
218  R_PreserveObject(RFUN(Snames));
219 
220  }
221 
222  break;
223 
224  case native: case regNative: {
225 
226  NAT(args) = args;
227 
228  PROTECT(NAT(sindex) = GET_SLOT(pompfun,install("stateindex")));
229  PROTECT(NAT(pindex) = GET_SLOT(pompfun,install("paramindex")));
230  PROTECT(NAT(cindex) = GET_SLOT(pompfun,install("covarindex")));
231  nprotect += 3;
232 
233  *((void **) (&(NAT(fun)))) = R_ExternalPtrAddr(fn);
234 
235  if (!isNull(NAT(args))) R_ReleaseObject(NAT(args));
236  if (!isNull(NAT(sindex))) R_ReleaseObject(NAT(sindex));
237  if (!isNull(NAT(pindex))) R_ReleaseObject(NAT(pindex));
238  if (!isNull(NAT(cindex))) R_ReleaseObject(NAT(cindex));
239 
240  R_PreserveObject(NAT(args));
241  R_PreserveObject(NAT(sindex));
242  R_PreserveObject(NAT(pindex));
243  R_PreserveObject(NAT(cindex));
244 
245  }
246 
247  break;
248 
249  default: // #nocov
250 
251  err("in 'pomp_desolve_setup': unrecognized 'mode'"); // # nocov
252 
253  }
254 
255  UNPROTECT(nprotect);
256 
257  return R_NilValue;
258 }
#define RFUN(X)
Definition: trajectory.c:153
#define COMMON(X)
Definition: trajectory.c:152
SEXP cindex
Definition: trajectory.c:146
pomp_skeleton * fun
Definition: trajectory.c:147
SEXP sindex
Definition: trajectory.c:144
SEXP pindex
Definition: trajectory.c:145
#define NAT(X)
Definition: trajectory.c:154
SEXP object
Definition: trajectory.c:127
Here is the call graph for this function:

◆ pomp_desolve_takedown()

SEXP pomp_desolve_takedown ( void  )

Definition at line 292 of file trajectory.c.

292  {
293  R_ReleaseObject(COMMON(object));
294  R_ReleaseObject(COMMON(params));
295  R_ReleaseObject(COMMON(cov));
296  COMMON(object) = R_NilValue;
297  COMMON(params) = R_NilValue;
298  COMMON(cov) = R_NilValue;
299  COMMON(nvars) = 0;
300  COMMON(npars) = 0;
301  COMMON(ncovars) = 0;
302  COMMON(nreps) = 0;
303 
304  switch (COMMON(mode)) {
305 
306  case Rfun: {
307 
308  R_ReleaseObject(RFUN(fn));
309  R_ReleaseObject(RFUN(args));
310  R_ReleaseObject(RFUN(Snames));
311  RFUN(fn) = R_NilValue;
312  RFUN(args) = R_NilValue;
313  RFUN(Snames) = R_NilValue;
314 
315  }
316 
317  break;
318 
319  case native: case regNative: {
320 
321  NAT(fun) = 0;
322  R_ReleaseObject(NAT(args));
323  R_ReleaseObject(NAT(sindex));
324  R_ReleaseObject(NAT(pindex));
325  R_ReleaseObject(NAT(cindex));
326  NAT(args) = R_NilValue;
327  NAT(sindex) = R_NilValue;
328  NAT(pindex) = R_NilValue;
329  NAT(cindex) = R_NilValue;
330 
331  }
332 
333  break;
334 
335  default: // #nocov
336 
337  err("in 'pomp_desolve_takedown': unrecognized 'mode'"); // # nocov
338 
339  break;
340 
341  }
342 
343  COMMON(mode) = undef;
344 
345  return R_NilValue;
346 }

◆ pomp_fun_handler()

SEXP pomp_fun_handler ( SEXP  pfun,
SEXP  gnsi,
pompfunmode mode,
SEXP  S,
SEXP  P,
SEXP  O,
SEXP  C 
)

Definition at line 30 of file pomp_fun.c.

32 {
33  int nprotect = 0;
34  SEXP f = R_NilValue;
35  SEXP sidx, pidx, oidx, cidx;
36 
37  *mode = *(INTEGER(GET_SLOT(pfun,install("mode"))));
38 
39  switch (*mode) {
40 
41  case Rfun: // R function
42 
43  PROTECT(f = GET_SLOT(pfun,install("R.fun"))); nprotect++;
44 
45  break;
46 
47  case native: case regNative: // native code
48 
49  if (*(LOGICAL(gnsi))) { // get native symbol information?
50 
51  SEXP nf, pack;
52  PROTECT(nf = GET_SLOT(pfun,install("native.fun")));
53  PROTECT(pack = GET_SLOT(pfun,install("PACKAGE")));
54  nprotect += 2;
55 
56  if (LENGTH(pack) < 1) {
57  PROTECT(pack = mkString("")); nprotect++; // #nocov
58  }
59 
60  if (*mode == native) {
61 
62  SEXP nsi;
63  PROTECT(nsi = eval(PROTECT(lang3(install("getNativeSymbolInfo"),nf,pack)),R_BaseEnv));
64  PROTECT(f = getListElement(nsi,"address"));
65  nprotect += 3;
66 
67  } else if (*mode == regNative) {
68 
69  const char *fname, *pkg;
70  fname = (const char *) CHAR(STRING_ELT(nf,0));
71  pkg = (const char *) CHAR(STRING_ELT(pack,0));
72  DL_FUNC fn;
73  fn = R_GetCCallable(pkg,fname);
74  PROTECT(f = R_MakeExternalPtrFn(fn,R_NilValue,R_NilValue)); nprotect++;
75 
76  }
77 
78  SET_SLOT(pfun,install("address"),f);
79 
80  if (S != NA_STRING) {
81  PROTECT(sidx = name_index(S,pfun,"statenames","state variables")); nprotect++;
82  SET_SLOT(pfun,install("stateindex"),sidx);
83  }
84 
85  if (P != NA_STRING) {
86  PROTECT(pidx = name_index(P,pfun,"paramnames","parameters")); nprotect++;
87  SET_SLOT(pfun,install("paramindex"),pidx);
88  }
89 
90  if (O != NA_STRING) {
91  PROTECT(oidx = name_index(O,pfun,"obsnames","observables")); nprotect++;
92  SET_SLOT(pfun,install("obsindex"),oidx);
93  }
94 
95  if (C != NA_STRING) {
96  PROTECT(cidx = name_index(C,pfun,"covarnames","covariates")); nprotect++;
97  SET_SLOT(pfun,install("covarindex"),cidx);
98  }
99 
100  } else { // native symbol info is stored
101 
102  PROTECT(f = GET_SLOT(pfun,install("address"))); nprotect++;
103 
104  }
105 
106  break;
107 
108  case undef: default:
109 
110  PROTECT(f = R_NilValue); nprotect++;
111  *mode = undef;
112 
113  break;
114 
115  }
116 
117  UNPROTECT(nprotect);
118  return f;
119 }
static R_INLINE SEXP name_index(SEXP provided, SEXP object, const char *slot, const char *humanreadable)
Definition: pomp_fun.c:10
Here is the call graph for this function:
Here is the caller graph for this function:

◆ pomp_vf_eval()

void pomp_vf_eval ( int *  neq,
double *  t,
double *  y,
double *  ydot,
double *  yout,
int *  ip 
)

Definition at line 260 of file trajectory.c.

261 {
262  switch (COMMON(mode)) {
263 
264  case Rfun: // R function
265 
266  eval_skeleton_R(ydot,t,y,REAL(COMMON(params)),
267  RFUN(fn),RFUN(args),RFUN(Snames),
270  &COMMON(covar_table),REAL(COMMON(cov)));
271 
272  break;
273 
274  case native: case regNative: // native code
275  eval_skeleton_native(ydot,t,y,REAL(COMMON(params)),
278  INTEGER(NAT(sindex)),INTEGER(NAT(pindex)),INTEGER(NAT(cindex)),
279  &COMMON(covar_table),NAT(fun),NAT(args),REAL(COMMON(cov)));
280 
281  break;
282 
283  default: // #nocov
284 
285  err("in 'pomp_vf_eval': unrecognized 'mode'"); // # nocov
286 
287  break;
288 
289  }
290 }
void eval_skeleton_R(double *, double *, double *, double *, SEXP, SEXP, SEXP, int, int, int, int, int, int, int, lookup_table_t *, double *)
Definition: skeleton.c:94
void eval_skeleton_native(double *, double *, double *, double *, int, int, int, int, int, int, int, int *, int *, int *, lookup_table_t *, pomp_skeleton *, SEXP, double *)
Definition: skeleton.c:245
Here is the call graph for this function:

◆ probe_acf()

SEXP probe_acf ( SEXP  x,
SEXP  lags,
SEXP  corr 
)

Definition at line 100 of file probe_acf.c.

100  {
101  SEXP ans, ans_names;
102  SEXP X;
103  int nlag, correlation, nvars, n;
104  int j, k, l;
105  double *p, *p1, *cov;
106  int *lag;
107  char tmp[BUFSIZ];
108 
109  nlag = LENGTH(lags); // number of lags
110  PROTECT(lags = AS_INTEGER(lags));
111  lag = INTEGER(lags);
112  correlation = *(INTEGER(AS_INTEGER(corr))); // correlation, or covariance?
113 
114  nvars = INTEGER(GET_DIM(x))[0]; // number of variables
115  n = INTEGER(GET_DIM(x))[1]; // number of observations
116 
117  PROTECT(X = duplicate(AS_NUMERIC(x)));
118 
119  PROTECT(ans = NEW_NUMERIC(nlag*nvars));
120 
121  pomp_acf_compute(REAL(ans),REAL(X),n,nvars,lag,nlag);
122 
123  if (correlation) {
124  l = 0;
125  cov = (double *) R_alloc(nvars,sizeof(double));
126  pomp_acf_compute(cov,REAL(X),n,nvars,&l,1); // compute lag-0 covariance
127  for (j = 0, p = REAL(ans), p1 = cov; j < nvars; j++, p1++)
128  for (k = 0; k < nlag; k++, p++)
129  *p /= *p1;
130  }
131 
132  PROTECT(ans_names = NEW_STRING(nlag*nvars));
133  for (j = 0, l = 0; j < nvars; j++) {
134  for (k = 0; k < nlag; k++, l++) {
135  snprintf(tmp,BUFSIZ,"acf[%d]",lag[k]);
136  SET_STRING_ELT(ans_names,l,mkChar(tmp));
137  }
138  }
139  SET_NAMES(ans,ans_names);
140 
141  UNPROTECT(4);
142  return ans;
143 }
static void pomp_acf_compute(double *acf, double *x, int n, int nvars, int *lags, int nlag)
Definition: probe_acf.c:11
Here is the call graph for this function:

◆ probe_ccf()

SEXP probe_ccf ( SEXP  x,
SEXP  y,
SEXP  lags,
SEXP  corr 
)

Definition at line 145 of file probe_acf.c.

145  {
146  SEXP ans, ans_names;
147  SEXP X, Y;
148  double cov[2], xx, *p;
149  int nlag, n, correlation;
150  int k;
151  char tmp[BUFSIZ];
152 
153  nlag = LENGTH(lags);
154  PROTECT(lags = AS_INTEGER(lags));
155  correlation = *(INTEGER(AS_INTEGER(corr))); // correlation, or covariance?
156 
157  n = LENGTH(x); // n = # of observations
158  if (n != LENGTH(y)) err("'x' and 'y' must have equal lengths"); // #nocov
159 
160  PROTECT(X = duplicate(AS_NUMERIC(x)));
161  PROTECT(Y = duplicate(AS_NUMERIC(y)));
162 
163  PROTECT(ans = NEW_NUMERIC(nlag));
164 
165  pomp_ccf_compute(REAL(ans),REAL(X),REAL(Y),n,INTEGER(lags),nlag);
166 
167  if (correlation) {
168  k = 0;
169  pomp_acf_compute(&cov[0],REAL(X),n,1,&k,1); // compute lag-0 covariance of x
170  pomp_acf_compute(&cov[1],REAL(Y),n,1,&k,1); // compute lag-0 covariance of y
171  xx = sqrt(cov[0]*cov[1]);
172  for (k = 0, p = REAL(ans); k < nlag; k++, p++) *p /= xx; // convert to correlation
173  }
174 
175  PROTECT(ans_names = NEW_STRING(nlag));
176  for (k = 0; k < nlag; k++) {
177  snprintf(tmp,BUFSIZ,"ccf[%d]",INTEGER(lags)[k]);
178  SET_STRING_ELT(ans_names,k,mkChar(tmp));
179  }
180  SET_NAMES(ans,ans_names);
181 
182  UNPROTECT(5);
183  return ans;
184 }
static void pomp_ccf_compute(double *ccf, double *x, double *y, int n, int *lags, int nlag)
Definition: probe_acf.c:47
Here is the call graph for this function:

◆ probe_marginal_setup()

SEXP probe_marginal_setup ( SEXP  ref,
SEXP  order,
SEXP  diff 
)

Definition at line 10 of file probe_marginal.c.

10  {
11  SEXP z, mm, tau, pivot, retval, retvalnames;
12  int n, nx, np, df, dim[2];
13 
14  np = *(INTEGER(AS_INTEGER(order))); // order of polynomial regression
15  df = *(INTEGER(AS_INTEGER(diff))); // order of differencing
16  n = LENGTH(ref); // n = number of observations
17  nx = n - df; // nx = rows in model matrix
18  dim[0] = nx; dim[1] = np; // dimensions of model matrix
19 
20  if (nx < 1) err("must have diff < number of observations");
21 
22  PROTECT(z = duplicate(AS_NUMERIC(ref)));
23  PROTECT(mm = makearray(2,dim));
24  PROTECT(tau = NEW_NUMERIC(np));
25  PROTECT(pivot = NEW_INTEGER(np));
26 
27  PROTECT(retval = NEW_LIST(3));
28  PROTECT(retvalnames = NEW_CHARACTER(3));
29  SET_STRING_ELT(retvalnames,0,mkChar("mm"));
30  SET_STRING_ELT(retvalnames,1,mkChar("tau"));
31  SET_STRING_ELT(retvalnames,2,mkChar("pivot"));
32  SET_ELEMENT(retval,0,mm);
33  SET_ELEMENT(retval,1,tau);
34  SET_ELEMENT(retval,2,pivot);
35  SET_NAMES(retval,retvalnames);
36 
37  order_reg_model_matrix(REAL(z),REAL(mm),REAL(tau),INTEGER(pivot),n,np,df);
38 
39  UNPROTECT(6);
40  return(retval);
41 }
static void order_reg_model_matrix(double *z, double *X, double *tau, int *pivot, int n, int np, int diff)
Here is the call graph for this function:

◆ probe_marginal_solve()

SEXP probe_marginal_solve ( SEXP  x,
SEXP  setup,
SEXP  diff 
)

Definition at line 43 of file probe_marginal.c.

43  {
44  SEXP X, mm, tau, pivot, beta, beta_names;
45  int n, nx, np, df;
46  int i;
47  char tmp[BUFSIZ];
48 
49  df = *(INTEGER(AS_INTEGER(diff))); // order of differencing
50  n = LENGTH(x); // n = number of observations
51 
52  // unpack the setup information
53  PROTECT(mm = VECTOR_ELT(setup,0)); // QR-decomposed model matrix
54  PROTECT(tau = VECTOR_ELT(setup,1)); // diagonals
55  PROTECT(pivot = VECTOR_ELT(setup,2)); // pivots
56 
57  nx = INTEGER(GET_DIM(mm))[0]; // nx = number of rows in model matrix
58  np = INTEGER(GET_DIM(mm))[1]; // np = order of polynomial
59 
60  if (n-df != nx) err("length of 'ref' must equal length of data");
61  PROTECT(X = duplicate(AS_NUMERIC(x)));
62 
63  PROTECT(beta = NEW_NUMERIC(np));
64  PROTECT(beta_names = NEW_STRING(np));
65  for (i = 0; i < np; i++) {
66  snprintf(tmp,BUFSIZ,"marg.%d",i+1);
67  SET_STRING_ELT(beta_names,i,mkChar(tmp));
68  }
69  SET_NAMES(beta,beta_names);
70 
71  order_reg_solve(REAL(beta),REAL(X),REAL(mm),REAL(tau),INTEGER(pivot),n,np,df);
72 
73  UNPROTECT(6);
74  return(beta);
75 }
static void order_reg_solve(double *beta, double *x, double *mm, double *tau, int *pivot, int n, int np, int diff)
Here is the call graph for this function:

◆ probe_nlar()

SEXP probe_nlar ( SEXP  x,
SEXP  lags,
SEXP  powers 
)

Definition at line 9 of file probe_nlar.c.

9  {
10  SEXP y, beta, beta_names;
11  int n, nterms;
12  int k;
13  double *mm;
14  char tmp[BUFSIZ];
15 
16  n = LENGTH(x); // n = # of observations
17  nterms = LENGTH(lags);
18 
19  PROTECT(y = duplicate(AS_NUMERIC(x)));
20  PROTECT(beta = NEW_NUMERIC(nterms));
21 
22  mm = (double *) R_alloc(n*nterms,sizeof(double)); // storage for the model matrix
23  poly_nlar_fit(REAL(beta),REAL(y),n,nterms,INTEGER(lags),INTEGER(powers),mm);
24 
25  PROTECT(beta_names = NEW_STRING(nterms));
26  for (k = 0; k < nterms; k++) {
27  snprintf(tmp,BUFSIZ,"nlar.%d^%d",INTEGER(lags)[k],INTEGER(powers)[k]);
28  SET_STRING_ELT(beta_names,k,mkChar(tmp));
29  }
30  SET_NAMES(beta,beta_names);
31 
32  UNPROTECT(3);
33  return beta;
34 }
static void poly_nlar_fit(double *beta, double *y, int n, int nterms, int *lag, int *power, double *X)
Definition: probe_nlar.c:39
Here is the call graph for this function:

◆ R_BetaBinom()

SEXP R_BetaBinom ( SEXP  n,
SEXP  size,
SEXP  prob,
SEXP  theta 
)

Definition at line 103 of file distributions.c.

103  {
104  int k, nval, ns, np, nt;
105  double *X, *S, *P, *T;
106  SEXP ans;
107  PROTECT(n = AS_INTEGER(n)); nval = INTEGER(n)[0];
108  PROTECT(size = AS_NUMERIC(size)); ns = LENGTH(size); S = REAL(size);
109  PROTECT(prob = AS_NUMERIC(prob)); np = LENGTH(prob); P = REAL(prob);
110  PROTECT(theta = AS_NUMERIC(theta)); nt = LENGTH(theta); T = REAL(theta);
111  PROTECT(ans = NEW_NUMERIC(nval)); X = REAL(ans);
112  GetRNGstate();
113  for (k = 0; k < nval; k++) {
114  X[k] = rbetabinom(S[k%ns],P[k%np],T[k%nt]);
115  }
116  PutRNGstate();
117  UNPROTECT(5);
118  return ans;
119 }
static R_INLINE double rbetabinom(double size, double prob, double theta)
Definition: pomp.h:286
Here is the call graph for this function:

◆ R_Euler_Multinom()

SEXP R_Euler_Multinom ( SEXP  n,
SEXP  size,
SEXP  rate,
SEXP  deltat 
)

Definition at line 23 of file distributions.c.

23  {
24  int ntrans = length(rate);
25  int dim[2];
26  SEXP x, nm;
27  if (length(size)>1)
28  warn("in 'reulermultinom': only the first element of 'size' is meaningful");
29  if (length(deltat)>1)
30  warn("in 'reulermultinom': only the first element of 'dt' is meaningful");
31  PROTECT(n = AS_INTEGER(n));
32  PROTECT(size = AS_NUMERIC(size));
33  PROTECT(rate = AS_NUMERIC(rate));
34  PROTECT(deltat = AS_NUMERIC(deltat));
35  dim[0] = ntrans;
36  dim[1] = *INTEGER(n);
37  if (dim[1] == NA_INTEGER || dim[1] < 0)
38  err("in 'reulermultinom': 'n' must be a non-negative integer.");
39  PROTECT(x = makearray(2,dim));
40  PROTECT(nm = GET_NAMES(rate));
41  setrownames(x,nm,2);
42  GetRNGstate();
43  reulermultinom_multi(dim[0],dim[1],REAL(size),REAL(rate),REAL(deltat),REAL(x));
44  PutRNGstate();
45  UNPROTECT(6);
46  return x;
47 }
static void reulermultinom_multi(int m, int n, double *size, double *rate, double *deltat, double *trans)
Definition: distributions.c:7
Here is the call graph for this function:

◆ R_GammaWN()

SEXP R_GammaWN ( SEXP  n,
SEXP  sigma,
SEXP  deltat 
)

Definition at line 77 of file distributions.c.

77  {
78  int k, nval, nsig, ndt;
79  double *x, *sig, *dt;
80  SEXP ans;
81  PROTECT(n = AS_INTEGER(n));
82  nval = INTEGER(n)[0];
83  PROTECT(sigma = AS_NUMERIC(sigma));
84  nsig = LENGTH(sigma);
85  sig = REAL(sigma);
86  PROTECT(deltat = AS_NUMERIC(deltat));
87  ndt = LENGTH(deltat);
88  dt = REAL(deltat);
89  PROTECT(ans = NEW_NUMERIC(nval));
90  x = REAL(ans);
91  GetRNGstate();
92  for (k = 0; k < nval; k++) {
93  x[k] = rgammawn(sig[k%nsig],dt[k%ndt]);
94  }
95  PutRNGstate();
96  UNPROTECT(4);
97  return ans;
98 }
static R_INLINE double rgammawn(double sigma, double dt)
Definition: pomp.h:129
Here is the call graph for this function:

◆ R_init_pomp()

void R_init_pomp ( DllInfo *  info)

Definition at line 53 of file init.c.

53  {
54  // C functions provided for users
55  R_RegisterCCallable("pomp","bspline_basis_eval_deriv",(DL_FUNC) &bspline_basis_eval_deriv);
56  R_RegisterCCallable("pomp","periodic_bspline_basis_eval_deriv",(DL_FUNC) &periodic_bspline_basis_eval_deriv);
57  R_RegisterCCallable("pomp","get_userdata",(DL_FUNC) &get_userdata);
58  R_RegisterCCallable("pomp","get_userdata_int",(DL_FUNC) &get_userdata_int);
59  R_RegisterCCallable("pomp","get_userdata_double",(DL_FUNC) &get_userdata_double);
60  R_RegisterCCallable("pomp","pomp_fun_handler",(DL_FUNC) &pomp_fun_handler);
61  R_RegisterCCallable("pomp","load_stack_incr",(DL_FUNC) &load_stack_incr);
62  R_RegisterCCallable("pomp","load_stack_decr",(DL_FUNC) &load_stack_decr);
63  R_RegisterCCallable("pomp","make_covariate_table",(DL_FUNC) &make_covariate_table);
64  R_RegisterCCallable("pomp","get_covariate_names",(DL_FUNC) &get_covariate_names);
65  R_RegisterCCallable("pomp","table_lookup",(DL_FUNC) &table_lookup);
66  R_RegisterCCallable("pomp","lookup_in_table",(DL_FUNC) &lookup_in_table);
67  R_RegisterCCallable("pomp","apply_probe_data",(DL_FUNC) &apply_probe_data);
68  R_RegisterCCallable("pomp","apply_probe_sim",(DL_FUNC) &apply_probe_sim);
69  R_RegisterCCallable("pomp","systematic_resampling",(DL_FUNC) &systematic_resampling);
70  R_RegisterCCallable("pomp","randwalk_perturbation", (DL_FUNC) &randwalk_perturbation);
71 
72  // Register routines
73  R_registerRoutines(info,NULL,callMethods,NULL,NULL);
74  R_useDynamicSymbols(info,TRUE);
75  R_forceSymbols(info,FALSE);
76 }
void bspline_basis_eval_deriv(double x, double *knots, int degree, int nbasis, int deriv, double *y)
Definition: bspline.c:115
SEXP randwalk_perturbation(SEXP, SEXP)
Definition: mif2.c:6
const SEXP get_userdata(const char *)
Definition: userdata.c:13
SEXP load_stack_decr(SEXP)
Definition: pomp_fun.c:130
SEXP apply_probe_sim(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
Definition: probe.c:29
SEXP load_stack_incr(SEXP)
Definition: pomp_fun.c:121
SEXP apply_probe_data(SEXP, SEXP)
Definition: probe.c:4
SEXP systematic_resampling(SEXP, SEXP)
Definition: resample.c:9
const double * get_userdata_double(const char *)
Definition: userdata.c:26
const int * get_userdata_int(const char *)
Definition: userdata.c:19
SEXP lookup_in_table(SEXP, SEXP)
Definition: lookup_table.c:24
static const R_CallMethodDef callMethods[]
Definition: init.c:4
Here is the call graph for this function:

◆ randwalk_perturbation()

SEXP randwalk_perturbation ( SEXP  params,
SEXP  rw_sd 
)

Definition at line 6 of file mif2.c.

7 {
8  double *xp = 0, *rw, *xrw, *xs;
9  SEXP Pnames, rwnames, pindex;
10  int *dim, *pidx;
11  int nrw = 0, npars, nreps;
12  int j, k;
13 
14  PROTECT(params = duplicate(params));
15 
16  // unpack parameter matrix
17  xp = REAL(params);
18  dim = INTEGER(GET_DIM(params)); npars = dim[0]; nreps = dim[1];
19  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
20 
21  // names of parameters undergoing random walk
22  PROTECT(rwnames = GET_NAMES(rw_sd));
23  nrw = LENGTH(rwnames); rw = REAL(rw_sd);
24 
25  // indices of parameters undergoing random walk
26  PROTECT(pindex = matchnames(Pnames,rwnames,"parameters"));
27  pidx = INTEGER(pindex);
28 
29  GetRNGstate();
30 
31  for (j = 0, xrw = rw; j < nrw; j++, pidx++, xrw++) {
32  for (k = 0, xs = xp+(*pidx); k < nreps; k++, xs += npars) {
33  *xs += *xrw * norm_rand();
34  }
35  }
36 
37  PutRNGstate();
38 
39  UNPROTECT(4);
40  return(params);
41 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_pomp_userdata()

SEXP set_pomp_userdata ( SEXP  userdata)

Definition at line 8 of file userdata.c.

8  {
9  USERDATA = userdata;
10  return R_NilValue;
11 }

◆ sobol_sequence()

SEXP sobol_sequence ( SEXP  dim,
SEXP  length 
)

Definition at line 234 of file sobolseq.c.

235 {
236  SEXP data;
237  unsigned int d = (unsigned int) INTEGER(dim)[0];
238  unsigned int n = (unsigned int) INTEGER(length)[0];
239  double *dp;
240  unsigned int k;
241  nlopt_sobol s = nlopt_sobol_create((unsigned int) d);
242  if (s==0) err("dimension is too high");
243  PROTECT(data = allocMatrix(REALSXP,d,n)); dp = REAL(data);
244  nlopt_sobol_skip(s,n,dp);
245  for (k = 1; k < n; k++) sobol_gen(s,dp+k*d);
247  UNPROTECT(1);
248  return(data);
249 }
static void nlopt_sobol_destroy(nlopt_sobol s)
Definition: sobolseq.c:196
static nlopt_sobol nlopt_sobol_create(unsigned sdim)
Definition: sobolseq.c:188
static void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x)
Definition: sobolseq.c:225
static int sobol_gen(soboldata *sd, double *x)
Definition: sobolseq.c:96
Here is the call 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 }
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 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
#define FIRST
Definition: ssa.c:198
#define NPAR
Definition: ssa.c:194
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
#define NCOV
Definition: ssa.c:195
#define NVAR
Definition: ssa.c:193
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition: ssa.c:11
#define RATEFN
Definition: ssa.c:197
#define ARGS
Definition: ssa.c:196
Here is the call graph for this function:
Here is the caller graph for this function:

◆ synth_loglik()

SEXP synth_loglik ( SEXP  ysim,
SEXP  ydat 
)

Definition at line 104 of file synth_lik.c.

104  {
105  SEXP loglik, dim, y;
106 
107  PROTECT(y = duplicate(AS_NUMERIC(ydat)));
108  PROTECT(dim = GET_DIM(ysim));
109  PROTECT(ysim = AS_NUMERIC(ysim));
110  PROTECT(loglik = NEW_NUMERIC(1));
111 
112  robust_synth_loglik(REAL(ysim),INTEGER(dim),REAL(y),REAL(loglik));
113 
114  UNPROTECT(4);
115  return loglik;
116 }
static void robust_synth_loglik(double *y, int *dim, double *ydat, double *loglik)
Definition: synth_lik.c:11
Here is the call graph for this function:

◆ systematic_resampling()

SEXP systematic_resampling ( SEXP  weights,
SEXP  np 
)

Definition at line 9 of file resample.c.

10 {
11  int m, n;
12  SEXP perm;
13 
14  m = *(INTEGER(AS_INTEGER(np)));
15  n = LENGTH(weights);
16  PROTECT(perm = NEW_INTEGER(m));
17  PROTECT(weights = AS_NUMERIC(weights));
18  GetRNGstate();
19  nosort_resamp(n,REAL(weights),m,INTEGER(perm),1);
20  PutRNGstate();
21  UNPROTECT(2);
22  return(perm);
23 }
void nosort_resamp(int nw, double *w, int np, int *p, int offset)
Definition: resample.c:25
Here is the call graph for this function:
Here is the caller graph for this function:

◆ table_lookup()

void table_lookup ( lookup_table_t tab,
double  x,
double *  y 
)

Definition at line 53 of file lookup_table.c.

54 {
55  int flag = 0;
56  int j, k, n;
57  double e;
58  if ((tab == 0) || (tab->length < 1) || (tab->width < 1)) return;
59  tab->index = findInterval(tab->x,tab->length,x,1,1,tab->index,&flag);
60  // warn only if we are *outside* the interval
61  if ((x < tab->x[0]) || (x > tab->x[tab->length-1]))
62  warn("in 'table_lookup': extrapolating at %le.", x);
63  switch (tab->order) {
64  case 1: default: // linear interpolation
65  e = (x - tab->x[tab->index-1]) / (tab->x[tab->index] - tab->x[tab->index-1]);
66  for (j = 0, k = tab->index*tab->width, n = k-tab->width; j < tab->width; j++, k++, n++) {
67  y[j] = e*(tab->y[k])+(1-e)*(tab->y[n]);
68  }
69  break;
70  case 0: // piecewise constant
71  if (flag < 0) n = 0;
72  else if (flag > 0) n = tab->index;
73  else n = tab->index-1;
74  for (j = 0, k = n*tab->width; j < tab->width; j++, k++) {
75  y[j] = tab->y[k];
76  }
77  break;
78  }
79 }
Here is the caller graph for this function:

◆ wpfilter()

SEXP wpfilter ( SEXP  X,
SEXP  Params,
SEXP  Weights,
SEXP  W,
SEXP  Trigger,
SEXP  Target,
SEXP  Np 
)

Definition at line 10 of file wpfilter.c.

11 {
12 
13  SEXP dimX, dimP, Xnames, Pnames;
14  int nvars, nreps, np;
15  // REVISIT THIS IF PARAMETER RESAMPLING IS IMPLEMENTED
16  // int npars;
17  // int do_pr = 0; // do parameter resampling? at the moment, we do not
18  int *dim, k;
19  const char *dimnm[] = {"name",".id"};
20 
21  PROTECT(dimX = GET_DIM(X));
22  dim = INTEGER(dimX);
23  nvars = dim[0]; nreps = dim[1];
24  PROTECT(Xnames = GET_ROWNAMES(GET_DIMNAMES(X)));
25 
26  PROTECT(Params = as_matrix(Params));
27  PROTECT(dimP = GET_DIM(Params));
28  dim = INTEGER(dimP);
29  // npars = dim[0];
30  if (dim[1] > 1 && nreps != dim[1]) // #nocov
31  err("ncol('params') does not match ncol('states')"); // #nocov
32  PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(Params)));
33 
34  np = *INTEGER(AS_INTEGER(Np));
35  // REVISIT THIS IF PARAMETER RESAMPLING IS IMPLEMENTED
36  // do_pr = (dim[1] > 1); // resample parameters as well as states
37 
38  PROTECT(Weights = duplicate(AS_NUMERIC(Weights)));
39  PROTECT(W = duplicate(AS_NUMERIC(W)));
40 
41  int nprotect = 7;
42 
43  double trigger, target;
44  trigger = *REAL(AS_NUMERIC(Trigger));
45  target = *REAL(AS_NUMERIC(Target));
46 
47  // add W to Weights. check weights for NAN or +Inf. find max of weights.
48  long double maxlogw = R_NegInf, oldw = 0;
49  double *xW = REAL(Weights);
50  double *xw = REAL(W);
51  for (k = 0; k < nreps; k++, xw++, xW++) {
52 
53  oldw += exp(*xW);
54  *xW += *xw;
55 
56  if (ISNAN(*xW) || *xW == R_PosInf) { // check the weights
57  SEXP rv;
58  PROTECT(rv = ScalarInteger(k+1)); nprotect++;
59  UNPROTECT(nprotect);
60  return rv;
61  }
62 
63  if (maxlogw < *xW) maxlogw = *xW;
64 
65  }
66 
67  // set up return list
68  SEXP retval, retvalnames, newdim;
69 
70  PROTECT(retval = NEW_LIST(5));
71  PROTECT(retvalnames = NEW_CHARACTER(5));
72  PROTECT(newdim = NEW_INTEGER(2));
73  nprotect += 3;
74 
75  dim = INTEGER(newdim);
76  dim[0] = nvars; dim[1] = nreps;
77  SET_DIM(X,newdim);
78  setrownames(X,Xnames,2);
79  fixdimnames(X,dimnm,2);
80 
81  SET_STRING_ELT(retvalnames,0,mkChar("states"));
82  SET_STRING_ELT(retvalnames,1,mkChar("params"));
83  SET_STRING_ELT(retvalnames,2,mkChar("weights"));
84  SET_STRING_ELT(retvalnames,3,mkChar("ess"));
85  SET_STRING_ELT(retvalnames,4,mkChar("loglik"));
86  SET_NAMES(retval,retvalnames);
87 
88  SEXP ess, loglik;
89  PROTECT(ess = NEW_NUMERIC(1)); // effective sample size
90  PROTECT(loglik = NEW_NUMERIC(1)); // conditional log likelihood
91  nprotect += 2;
92 
93  if (maxlogw == R_NegInf) { // all particles have zero likelihood
94  *REAL(ess) = 0;
95  *REAL(loglik) = R_NegInf;
96  } else {
97  // compute sum and sum of squares
98  long double w = 0, ws = 0;
99  xw = REAL(W);
100  xW = REAL(Weights);
101  for (k = 0; k < nreps; k++, xw++, xW++) {
102  *xW -= maxlogw;
103  *xw = exp(*xW);
104  if (*xw != 0) {
105  w += *xw;
106  ws += (*xw)*(*xw);
107  }
108  }
109  *(REAL(ess)) = w*w/ws;
110  *(REAL(loglik)) = maxlogw + log(w) - log(oldw);
111  }
112 
113  SET_ELEMENT(retval,3,ess);
114  SET_ELEMENT(retval,4,loglik);
115 
116  if (maxlogw == R_NegInf || (*REAL(ess) >= nreps*trigger && np == nreps)) { // do not resample
117 
118  SET_ELEMENT(retval,0,X);
119  SET_ELEMENT(retval,1,Params);
120  SET_ELEMENT(retval,2,Weights);
121 
122  } else { // no resampling
123 
124  // create storage for new states
125  SEXP newstates = R_NilValue;
126  double *ss = 0, *st = 0;
127  {
128  int xdim[] = {nvars, np};
129  PROTECT(newstates = makearray(2,xdim)); nprotect++;
130  setrownames(newstates,Xnames,2);
131  fixdimnames(newstates,dimnm,2);
132  ss = REAL(X);
133  st = REAL(newstates);
134  }
135  SET_ELEMENT(retval,0,newstates);
136 
137  // REVISIT THIS IF PARAMETER RESAMPLING IS IMPLEMENTED
138  // create storage for new parameters
139  // SEXP newparams = R_NilValue;
140  // double *ps = 0, *pt = 0;
141  // if (do_pr) {
142  // int xdim[] = {npars, np};
143  // PROTECT(newparams = makearray(2,xdim)); nprotect++;
144  // setrownames(newparams,Pnames,2);
145  // fixdimnames(newparams,dimnm,2);
146  // ps = REAL(Params);
147  // pt = REAL(newparams);
148  // SET_ELEMENT(retval,1,newparams);
149  // } else {
150  SET_ELEMENT(retval,1,Params);
151  // }
152 
153  // create storage for new weights
154  SEXP newweights = R_NilValue;
155  double *ws = 0, *wt = 0;
156  {
157  PROTECT(newweights = NEW_NUMERIC(np)); nprotect++;
158  ws = REAL(Weights);
159  wt = REAL(newweights);
160  }
161  SET_ELEMENT(retval,2,newweights);
162 
163  // compute target resampling weights
164  for (k = 0, xw = REAL(W), xW = REAL(Weights); k < nreps; k++, xw++, xW++) {
165  *xw = *xW;
166  *xW *= target;
167  *xw = exp(*xw - *xW);
168  }
169 
170  // do resampling
171  {
172  int sample[np];
173  GetRNGstate();
174  nosort_resamp(nreps,REAL(W),np,sample,0);
175  PutRNGstate();
176 
177  for (k = 0; k < np; k++) { // copy the particles
178  int sp = sample[k], j;
179  double *xx;
180  // double *xp;
181  for (j = 0, xx = ss+nvars*sp; j < nvars; j++, st++, xx++)
182  *st = *xx;
183  // REVISIT THIS IF PARAMETER RESAMPLING IS IMPLEMENTED
184  // if (do_pr) {
185  // for (j = 0, xp = ps+npars*sp; j < npars; j++, pt++, xp++)
186  // *pt = *xp;
187  // }
188  wt[k] = ws[sp];
189  }
190  }
191  }
192 
193  UNPROTECT(nprotect);
194  return(retval);
195 }
Here is the call graph for this function: