pomp
Inference for partially observed Markov processes
•All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
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 E_Euler_Multinom (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 
)
extern

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}
X_0
#define X_0
Definition gompertz.c:95
r
#define r
Definition gompertz.c:91
sigma
#define sigma
Definition gompertz.c:93
T_sigma
#define T_sigma
Definition gompertz.c:98
T_tau
#define T_tau
Definition gompertz.c:99
K
#define K
Definition gompertz.c:9
tau
#define tau
Definition gompertz.c:94
T_K
#define T_K
Definition gompertz.c:97
T_r
#define T_r
Definition gompertz.c:96
T_X_0
#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 
)
extern

Definition at line 20 of file gompertz.c.

25 {
26 *lik = dlnorm(Y,log(X),TAU,give_log);
27}
X
#define X
Definition gompertz.c:14
TAU
#define TAU
Definition gompertz.c:11
Y
#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 
)
extern

Definition at line 40 of file gompertz.c.

45 {
46 EY = X*exp(TAU*TAU/2);
47}
EY
#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 
)
extern

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 
)
extern

Definition at line 50 of file gompertz.c.

52 {
53 double et = exp(TAU*TAU);
54 VY = X*X*et*(et-1);
55}
VY
#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 
)
extern

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}
XPRIME
#define XPRIME
Definition gompertz.c:15
R
#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 
)
extern

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}
SIGMA
#define SIGMA
Definition gompertz.c:10

◆ _gompertz_to_trans()

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

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 
)
extern

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}
X1
#define X1
Definition ou2.c:59
Y1
#define Y1
Definition ou2.c:62
Y2
#define Y2
Definition ou2.c:63
TAU
#define TAU
Definition ou2.c:57
X2
#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 
)
extern

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 
)
extern

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}
SIGMA3
#define SIGMA3
Definition ou2.c:56
ALPHA2
#define ALPHA2
Definition ou2.c:51
dens_ou2
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
SIGMA1
#define SIGMA1
Definition ou2.c:54
SIGMA2
#define SIGMA2
Definition ou2.c:55
ALPHA3
#define ALPHA3
Definition ou2.c:52
ALPHA4
#define ALPHA4
Definition ou2.c:53
ALPHA1
#define ALPHA1
Definition ou2.c:50
err
#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 
)
extern

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 
)
extern

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 
)
extern

Definition at line 71 of file ou2.c.

76 {
77 sim_ou2(&x[X1],&x[X2],ALPHA1,ALPHA2,ALPHA3,ALPHA4,SIGMA1,SIGMA2,SIGMA3);
78}
sim_ou2
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 
)
extern

Definition at line 143 of file ou2.c.

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

◆ add_skel_args()

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

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}
Snames
SEXP Snames
Definition trajectory.c:140
args
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 
)
extern

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 
)
extern

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}
do_simulate
SEXP do_simulate(SEXP, SEXP, SEXP, SEXP, SEXP)
Definition simulate.c:6
setrownames
static R_INLINE void setrownames(SEXP x, SEXP names, int rank)
Definition pomp_defines.h:110
makearray
static R_INLINE SEXP makearray(int rank, const int *dim)
Definition pomp_defines.h:40
params
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 
)
extern

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}
bspline_eval
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 
)
extern

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 
)
extern

Definition at line 137 of file distributions.c.

137 {
138 int k, n, nx, ns, np, nt;
139 double *F, *X, *S, *P, *T;
140 SEXP f;
141 PROTECT(x = AS_NUMERIC(x)); nx = LENGTH(x); X = REAL(x);
142 PROTECT(size = AS_NUMERIC(size)); ns = LENGTH(size); S = REAL(size);
143 PROTECT(prob = AS_NUMERIC(prob)); np = LENGTH(prob); P = REAL(prob);
144 PROTECT(theta = AS_NUMERIC(theta)); nt = LENGTH(theta); T = REAL(theta);
145 PROTECT(log = AS_INTEGER(log));
146 n = (nx > ns) ? nx : ns;
147 n = (n > np) ? n : np;
148 n = (n > nt) ? n : nt;
149 PROTECT(f = NEW_NUMERIC(n)); F = REAL(f);
150 for (k = 0; k < n; k++) {
151 F[k] = dbetabinom(X[k%nx],S[k%ns],P[k%np],T[k%nt],INTEGER(log)[0]);
152 }
153 UNPROTECT(6);
154 return f;
155}
dbetabinom
static R_INLINE double dbetabinom(double x, double size, double prob, double theta, int give_log)
Definition pomp.h:321
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 
)
extern

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}
deulermultinom_multi
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
warn
#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 
)
extern

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}
init_density
static SEXP init_density(SEXP func, SEXP X, SEXP t0, SEXP params, SEXP covar, SEXP log, SEXP args, SEXP gnsi)
Definition dinit.c:98
as_matrix
static R_INLINE SEXP as_matrix(SEXP x)
Definition pomp_defines.h:145
fn
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 
)
extern

Definition at line 111 of file dmeasure.c.

112{
113
114 pompfunmode mode = undef;
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}
make_covariate_table
lookup_table_t make_covariate_table(SEXP, int *)
Definition lookup_table.c:11
pomp_fun_handler
SEXP pomp_fun_handler(SEXP, SEXP, pompfunmode *, SEXP, SEXP, SEXP, SEXP)
Definition pomp_fun.c:30
table_lookup
void table_lookup(lookup_table_t *, double, double *)
Definition lookup_table.c:53
get_covariate_names
SEXP get_covariate_names(SEXP)
Definition lookup_table.c:7
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Onames, SEXP Snames, SEXP Pnames, SEXP Cnames, SEXP log)
Definition dmeasure.c:11
eval_call
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
ret_array
static R_INLINE SEXP ret_array(int nreps, int ntimes)
Definition dmeasure.c:101
pomp_dmeasure
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
as_state_array
static R_INLINE SEXP as_state_array(SEXP x)
Definition pomp_defines.h:183
pompfunmode
pompfunmode
Definition pomp_defines.h:16
Rfun
@ Rfun
Definition pomp_defines.h:16
native
@ native
Definition pomp_defines.h:16
undef
@ undef
Definition pomp_defines.h:16
regNative
@ regNative
Definition pomp_defines.h:16
lookup_table_t
Definition pomp_defines.h:20
npars
int npars
Definition trajectory.c:132
ncovars
int ncovars
Definition trajectory.c:133
mode
pompfunmode mode
Definition trajectory.c:126
nreps
int nreps
Definition trajectory.c:134
cov
SEXP cov
Definition trajectory.c:129
nvars
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 
)
extern

Definition at line 51 of file dprior.c.

52{
53
54 pompfunmode mode = undef;
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}
eval_call
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *p, int n)
Definition dprior.c:34
add_args
static R_INLINE SEXP add_args(SEXP names, SEXP log, SEXP args)
Definition dprior.c:10
pomp_dprior
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 
)
extern

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}
onestep_density
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 
)
extern

Definition at line 98 of file emeasure.c.

99{
100 pompfunmode mode = undef;
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}
eval_call
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
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition emeasure.c:10
ret_array
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition emeasure.c:84
pomp_emeasure
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
fixdimnames
static R_INLINE void fixdimnames(SEXP x, const char **names, int n)
Definition pomp_defines.h:129
invalid_names
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 
)
extern

Definition at line 49 of file partrans.c.

50{
51
52 SEXP Pnames, tparams, pompfun, fn, args, ob;
53 pompfunmode mode = undef;
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}
eval_call
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *p, int n)
Definition partrans.c:33
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP names)
Definition partrans.c:12
direction_t
direction_t
Definition partrans.c:10
from
@ from
Definition partrans.c:10
to
@ to
Definition partrans.c:10
pomp_transform
void pomp_transform(double *pt, const double *p, const int *parindex)
Definition pomp.h:115
matchnames
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 
)
extern

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;
97 pompfunmode mode = undef;
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}
pomp_rinit
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
ret_array
static R_INLINE SEXP ret_array(int m, int n, SEXP names)
Definition rinit.c:79
paste0
static R_INLINE SEXP paste0(SEXP a, SEXP b, SEXP c)
Definition rinit.c:10
pomp_default_rinit
static SEXP pomp_default_rinit(SEXP params, SEXP Pnames, int npar, int nrep, int nsim)
Definition rinit.c:254
eval_call
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t0, double *p, int npar, double *c, int ncov)
Definition rinit.c:59
add_args
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 
)
extern

Definition at line 98 of file rmeasure.c.

99{
100 pompfunmode mode = undef;
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}
pomp_rmeasure
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
eval_call
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
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition rmeasure.c:10
ret_array
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 
)
extern

Definition at line 58 of file rprior.c.

59{
60
61 pompfunmode mode = undef;
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}
pomp_rprior
void pomp_rprior(double *p, const int *parindex)
Definition pomp.h:107
eval_call
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *p, int n)
Definition rprior.c:31
ret_array
static R_INLINE SEXP ret_array(SEXP params)
Definition rprior.c:47
add_args
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 
)
extern

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}
euler_simulator
SEXP euler_simulator(SEXP, SEXP, SEXP, SEXP, SEXP, double, rprocmode, SEXP, SEXP, SEXP, SEXP)
Definition euler.c:108
SSA_simulator
SEXP SSA_simulator(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
Definition ssa.c:226
gill
@ gill
Definition pomp_defines.h:17
discrete
@ discrete
Definition pomp_defines.h:17
dflt
@ dflt
Definition pomp_defines.h:17
onestep
@ onestep
Definition pomp_defines.h:17
euler
@ euler
Definition pomp_defines.h:17
pomp_default_rprocess
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 
)
extern

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}
do_rmeasure
SEXP do_rmeasure(SEXP, SEXP, SEXP, SEXP, SEXP)
Definition rmeasure.c:98
do_rinit
SEXP do_rinit(SEXP, SEXP, SEXP, SEXP, SEXP)
Definition rinit.c:91
do_rprocess
SEXP do_rprocess(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
Definition rprocess.c:25
setcolnames
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 
)
extern

Definition at line 326 of file skeleton.c.

327{
328
329 int nvars, npars, nrepp, nrepx, nreps, ntimes, ncovars;
330 pompfunmode mode = undef;
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 'times' 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
400 eval_skeleton_native(
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}
pomp_skeleton
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
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 skeleton.c:245
ret_array
static R_INLINE SEXP ret_array(int nvars, int nreps, int ntimes, SEXP Snames)
Definition skeleton.c:81
add_skel_args
SEXP add_skel_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition skeleton.c:10
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 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 
)
extern

Definition at line 99 of file vmeasure.c.

100{
101 pompfunmode mode = undef;
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}
pomp_vmeasure
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
eval_call
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
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition vmeasure.c:10
ret_array
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:

◆ E_Euler_Multinom()

SEXP E_Euler_Multinom ( SEXP  size,
SEXP  rate,
SEXP  deltat 
)
extern

Definition at line 71 of file distributions.c.

71 {
72 int ntrans = length(rate);
73 SEXP x;
74 if (length(size)>1)
75 warn("in 'eeulermultinom': only the first element of 'size' is meaningful");
76 if (length(deltat)>1)
77 warn("in 'eeulermultinom': only the first element of 'dt' is meaningful");
78 PROTECT(size = AS_NUMERIC(size));
79 PROTECT(rate = AS_NUMERIC(rate));
80 PROTECT(deltat = AS_NUMERIC(deltat));
81 PROTECT(x = duplicate(rate));
82 eeulermultinom(ntrans,*REAL(size),REAL(rate),*REAL(deltat),REAL(x));
83 UNPROTECT(4);
84 return x;
85}
eeulermultinom
static R_INLINE void eeulermultinom(int m, double size, const double *rate, double dt, double *trans)
Definition pomp.h:138
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 
)
extern

Definition at line 107 of file euler.c.

112 {
113
114 pompfunmode mode = undef;
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}
num_map_steps
int num_map_steps(double t1, double t2, double deltat)
Definition euler.c:342
num_euler_steps
int num_euler_steps(double t1, double t2, double *deltat)
Definition euler.c:317
eval_call
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
add_args
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition euler.c:10
ret_array
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition euler.c:92
pomp_onestep_sim
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 
)
extern

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
260 table_lookup(covar_table,*time,cov);
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}
covar_table
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 
)
extern

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
114 table_lookup(covar_table,*time,cov);
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}
eval_call
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)
extern

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)
extern

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)
extern

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}
getListElement
static R_INLINE SEXP getListElement(SEXP list, const char *str)
Definition pomp_defines.h:232
USERDATA
#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)
extern

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)
extern

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)
extern

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}
from_log_barycentric
static R_INLINE void from_log_barycentric(double *xt, const double *x, int n)
Definition pomp.h:288
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 
)
extern

Definition at line 19 of file trajectory.c.

20{
21
22 pompfunmode mode = undef;
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}
add_skel_args
SEXP add_skel_args(SEXP, SEXP, SEXP, SEXP)
Definition skeleton.c:10
iterate_skeleton_R
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
iterate_skeleton_native
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
ret_array
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 
)
extern

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
302 table_lookup(covar_table,t,cov);
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}
num_map_steps
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 
)
extern

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
194 table_lookup(covar_table,t,cov);
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)
extern

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)
extern

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)
extern

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}
to_log_barycentric
static R_INLINE void to_log_barycentric(double *xt, const double *x, int n)
Definition pomp.h:278
Here is the call graph for this function:

◆ LogitTransform()

SEXP LogitTransform ( SEXP  P)
extern

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}
logit
static R_INLINE double logit(double p)
Definition pomp.h:118
Here is the call graph for this function:

◆ logmeanexp()

SEXP logmeanexp ( const SEXP  X,
const SEXP  Drop 
)
extern

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 
)
extern

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}
table_lookup
void table_lookup(lookup_table_t *tab, double x, double *y)
Definition lookup_table.c:53
make_covariate_table
lookup_table_t make_covariate_table(SEXP object, int *ncovar)
Definition lookup_table.c:11
get_covariate_names
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 
)
extern

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}
lookup_table_t::index
int index
Definition pomp_defines.h:22
lookup_table_t::width
int width
Definition pomp_defines.h:21
lookup_table_t::order
int order
Definition pomp_defines.h:23
lookup_table_t::length
int length
Definition pomp_defines.h:21
lookup_table_t::x
double * x
Definition pomp_defines.h:24
lookup_table_t::y
double * y
Definition pomp_defines.h:25
Here is the caller graph for this function:

◆ nosort_resamp()

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

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 
)
extern

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 
)
extern

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 
)
extern

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}
periodic_bspline_basis_eval_deriv
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 
)
extern

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 
)
extern

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}
nosort_resamp
void nosort_resamp(int, double *, int, int *, int)
Definition resample.c:25
filt_mean
static void filt_mean(int, int, int, long double, const double *, const double *, double *)
Definition pfilter.c:294
pred_mean_var
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 
)
extern

Definition at line 156 of file trajectory.c.

156 {
157
158 pompfunmode mode = undef;
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);
189 COMMON(ncovars) = 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}
RFUN
#define RFUN(X)
Definition trajectory.c:153
COMMON
#define COMMON(X)
Definition trajectory.c:152
cindex
SEXP cindex
Definition trajectory.c:146
fun
pomp_skeleton * fun
Definition trajectory.c:147
sindex
SEXP sindex
Definition trajectory.c:144
pindex
SEXP pindex
Definition trajectory.c:145
NAT
#define NAT(X)
Definition trajectory.c:154
object
SEXP object
Definition trajectory.c:127
Here is the call graph for this function:

◆ pomp_desolve_takedown()

SEXP pomp_desolve_takedown ( void  )
extern

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 
)
extern

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}
name_index
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 
)
extern

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),
268 COMMON(nvars),COMMON(npars),COMMON(ncovars),1,
269 COMMON(nreps),COMMON(nreps),COMMON(nreps),
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)),
276 COMMON(nvars),COMMON(npars),COMMON(ncovars),1,
277 COMMON(nreps),COMMON(nreps),COMMON(nreps),
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}
eval_skeleton_R
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
eval_skeleton_native
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 
)
extern

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}
pomp_acf_compute
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 
)
extern

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}
pomp_ccf_compute
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 
)
extern

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}
order_reg_model_matrix
static void order_reg_model_matrix(double *z, double *X, double *tau, int *pivot, int n, int np, int diff)
Definition probe_marginal.c:78
Here is the call graph for this function:

◆ probe_marginal_solve()

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

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}
order_reg_solve
static void order_reg_solve(double *beta, double *x, double *mm, double *tau, int *pivot, int n, int np, int diff)
Definition probe_marginal.c:117
Here is the call graph for this function:

◆ probe_nlar()

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

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}
poly_nlar_fit
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 
)
extern

Definition at line 119 of file distributions.c.

119 {
120 int k, nval, ns, np, nt;
121 double *X, *S, *P, *T;
122 SEXP ans;
123 PROTECT(n = AS_INTEGER(n)); nval = INTEGER(n)[0];
124 PROTECT(size = AS_NUMERIC(size)); ns = LENGTH(size); S = REAL(size);
125 PROTECT(prob = AS_NUMERIC(prob)); np = LENGTH(prob); P = REAL(prob);
126 PROTECT(theta = AS_NUMERIC(theta)); nt = LENGTH(theta); T = REAL(theta);
127 PROTECT(ans = NEW_NUMERIC(nval)); X = REAL(ans);
128 GetRNGstate();
129 for (k = 0; k < nval; k++) {
130 X[k] = rbetabinom(S[k%ns],P[k%np],T[k%nt]);
131 }
132 PutRNGstate();
133 UNPROTECT(5);
134 return ans;
135}
rbetabinom
static R_INLINE double rbetabinom(double size, double prob, double theta)
Definition pomp.h:314
Here is the call graph for this function:

◆ R_Euler_Multinom()

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

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}
reulermultinom_multi
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 
)
extern

Definition at line 93 of file distributions.c.

93 {
94 int k, nval, nsig, ndt;
95 double *x, *sig, *dt;
96 SEXP ans;
97 PROTECT(n = AS_INTEGER(n));
98 nval = INTEGER(n)[0];
99 PROTECT(sigma = AS_NUMERIC(sigma));
100 nsig = LENGTH(sigma);
101 sig = REAL(sigma);
102 PROTECT(deltat = AS_NUMERIC(deltat));
103 ndt = LENGTH(deltat);
104 dt = REAL(deltat);
105 PROTECT(ans = NEW_NUMERIC(nval));
106 x = REAL(ans);
107 GetRNGstate();
108 for (k = 0; k < nval; k++) {
109 x[k] = rgammawn(sig[k%nsig],dt[k%ndt]);
110 }
111 PutRNGstate();
112 UNPROTECT(4);
113 return ans;
114}
rgammawn
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)
extern

Definition at line 54 of file init.c.

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

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)
extern

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 
)
extern

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);
246 nlopt_sobol_destroy(s);
247 UNPROTECT(1);
248 return(data);
249}
nlopt_sobol_destroy
static void nlopt_sobol_destroy(nlopt_sobol s)
Definition sobolseq.c:196
nlopt_sobol_create
static nlopt_sobol nlopt_sobol_create(unsigned sdim)
Definition sobolseq.c:188
nlopt_sobol_skip
static void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x)
Definition sobolseq.c:225
sobol_gen
static int sobol_gen(soboldata *sd, double *x)
Definition sobolseq.c:96
nlopt_soboldata_s
Definition sobolseq.c:65
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 
)
extern

Definition at line 226 of file ssa.c.

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

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}
robust_synth_loglik
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 
)
extern

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}
nosort_resamp
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 
)
extern

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 
)
extern

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: