Loading...
Searching...
No Matches
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()
|
extern |
Definition at line 113 of file gompertz.c.
116 {
122}
◆ _gompertz_normal_dmeasure()
|
extern |
Definition at line 20 of file gompertz.c.
25 {
27}
◆ _gompertz_normal_emeasure()
|
extern |
Definition at line 40 of file gompertz.c.
◆ _gompertz_normal_rmeasure()
|
extern |
Definition at line 30 of file gompertz.c.
◆ _gompertz_normal_vmeasure()
|
extern |
Definition at line 50 of file gompertz.c.
◆ _gompertz_skeleton()
|
extern |
Definition at line 70 of file gompertz.c.
75 {
76 double deltat = 1.0;
79}
◆ _gompertz_step()
|
extern |
Definition at line 58 of file gompertz.c.
◆ _gompertz_to_trans()
|
extern |
◆ _ou2_dmeasure()
|
extern |
◆ _ou2_emeasure()
|
extern |
◆ _ou2_pdf()
|
extern |
Definition at line 82 of file ou2.c.
88 {
89 if (t2-t1 != 1)
93}
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
Here is the call graph for this function:
◆ _ou2_rmeasure()
|
extern |
◆ _ou2_skel()
|
extern |
◆ _ou2_step()
|
extern |
◆ _ou2_vmeasure()
|
extern |
◆ add_skel_args()
|
extern |
Definition at line 10 of file skeleton.c.
11{
12
13 SEXP var;
14 int v;
15
17
18 // Covariates
19 for (v = LENGTH(Cnames)-1; v >= 0; v--) {
20 var = NEW_NUMERIC(1);
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);
31 UNPROTECT(1);
32 PROTECT(args);
33 SET_TAG(args,installChar(STRING_ELT(Pnames,v)));
34 }
35
36 // Latent state variables
38 var = NEW_NUMERIC(1);
40 UNPROTECT(1);
41 PROTECT(args);
43 }
44
45 // Time
46 var = NEW_NUMERIC(1);
48 UNPROTECT(1);
49 PROTECT(args);
51
52 UNPROTECT(1);
54
55}
Here is the caller graph for this function:
◆ apply_probe_data()
|
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 R_ClosureEnv(VECTOR_ELT(probes,i))));
17 if (!IS_NUMERIC(VECTOR_ELT(vals,i))) {
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()
|
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))
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;
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 R_ClosureEnv(VECTOR_ELT(probes,p))));
86 if (!IS_NUMERIC(val)) {
88 }
89
90 len = LENGTH(val);
91 if (s == 0)
92 len0 = len;
93 else if (len != len0) {
95 }
96 if (k+len > nvals)
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)
108
109 UNPROTECT(9);
110 return retval;
111}
static R_INLINE void setrownames(SEXP x, SEXP names, int rank)
Definition pomp_defines.h:110
Here is the call graph for this function:
Here is the caller graph for this function:
◆ bspline_basis()
|
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;
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];
73 dx = (maxx-minx)/((double) (nb-deg));
74 knots[0] = minx-deg*dx;
75 for (i = 1; i < nk; i++) knots[i] = knots[i-1]+dx;
76 for (i = 0; i < nb; i++) {
77 bspline_eval(ydata,xdata,nx,i,deg,d,&knots[0]);
78 ydata += nx;
79 }
80 R_Free(knots);
81 UNPROTECT(2);
82 return(y);
83}
static void bspline_eval(double *y, const double *x, int nx, int i, int degree, int deriv, const double *knots)
Definition bspline.c:10
Here is the call graph for this function:
◆ bspline_basis_eval_deriv()
|
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()
|
extern |
Definition at line 137 of file distributions.c.
137 {
138 int k, n, nx, ns, np, nt;
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++) {
152 }
153 UNPROTECT(6);
154 return f;
155}
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()
|
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)
56 n = dim[1];
57 if (length(size)>1)
59 if (length(deltat)>1)
61 PROTECT(f = NEW_NUMERIC(n));
62 PROTECT(size = AS_NUMERIC(size));
63 PROTECT(rate = AS_NUMERIC(rate));
64 PROTECT(deltat = AS_NUMERIC(deltat));
65 PROTECT(log = AS_LOGICAL(log));
66 deulermultinom_multi(ntrans,n,REAL(size),REAL(rate),REAL(deltat),REAL(x),INTEGER(log),REAL(f));
67 UNPROTECT(5);
68 return f;
69}
static void deulermultinom_multi(int m, int n, double *size, double *rate, double *deltat, double *trans, int *give_log, double *f)
Definition distributions.c:15
Here is the call graph for this function:
◆ do_dinit()
|
extern |
Definition at line 220 of file dinit.c.
223 {
225 PROTECT(t0=AS_NUMERIC(t0));
226 PROTECT(x = as_matrix(x));
228 // extract the process function
230 // extract other arguments
232 PROTECT(covar = GET_SLOT(object,install("covar")));
233 // evaluate the density
235 UNPROTECT(7);
236 return F;
237}
static SEXP init_density(SEXP func, SEXP X, SEXP t0, SEXP params, SEXP covar, SEXP log, SEXP args, SEXP gnsi)
Definition dinit.c:98
Here is the call graph for this function:
◆ do_dmeasure()
|
extern |
Definition at line 111 of file dmeasure.c.
112{
113
116 SEXP Snames, Pnames, Cnames, Onames;
117 SEXP cvec, pompfun;
119 SEXP F;
120 int *dim;
121 lookup_table_t covariate_table;
123
124 PROTECT(times = AS_NUMERIC(times));
125 ntimes = length(times);
127
128 PROTECT(y = as_matrix(y));
129 dim = INTEGER(GET_DIM(y));
130 nobs = dim[0];
131
132 if (ntimes != dim[1])
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])
141
143 dim = INTEGER(GET_DIM(params));
144 npars = dim[0]; nrepsp = dim[1];
145
146 nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
147
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)));
155
156 // set up the covariate table
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")));
164
165 // extract 'userdata' as pairlist
167
168 // create array to store results
170
171 int nprotect = 13;
172
174
176
178 double *ft = REAL(F);
179 int j, k;
180
181 // build argument list
183
184 for (k = 0; k < ntimes; k++, time++, ys += nobs) { // loop over times
185
186 R_CheckUserInterrupt(); // check for user interrupt
187
189
191
192 // evaluate the call
193 PROTECT(
194 ans = eval_call(
196 time,
197 ys,nobs,
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
218 int *oidx, *sidx, *pidx, *cidx;
219 int give_log;
220 pomp_dmeasure *ff = NULL;
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
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
243
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
264 *ft = R_NaReal;
265 }
266 }
267
269
270 }
271
272 }
273
274 UNPROTECT(nprotect);
275 return F;
276}
SEXP pomp_fun_handler(SEXP, SEXP, pompfunmode *, SEXP, SEXP, SEXP, SEXP)
Definition pomp_fun.c:30
static R_INLINE SEXP add_args(SEXP args, SEXP Onames, SEXP Snames, SEXP Pnames, SEXP Cnames, SEXP log)
Definition dmeasure.c:11
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *y, int nobs, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition dmeasure.c:75
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
Definition pomp_defines.h:20
Here is the call graph for this function:
◆ do_dprior()
|
extern |
Definition at line 51 of file dprior.c.
52{
53
57 int *dim;
58
60 dim = INTEGER(GET_DIM(params));
62
63 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
64
65 // extract the user-defined function
66 PROTECT(pompfun = GET_SLOT(object,install("dprior")));
68
69 // extract 'userdata' as pairlist
71
72 // to store results
73 PROTECT(F = NEW_NUMERIC(nreps));
74
75 int nprotect = 6;
76
79 SEXP ans;
80 double *ps, *pt;
81 int j;
82
84
86
88 *pt = *(REAL(AS_NUMERIC(ans)));
89 UNPROTECT(1);
90
91 }
92 }
93
94 break;
95
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
107
108 give_log = *(INTEGER(AS_INTEGER(log)));
109
110 R_CheckUserInterrupt(); // check for user interrupt
111
112 // loop over replicates
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
128 *pt = (give_log) ? 0.0 : 1.0;
129
130 }
131
132 }
133
134 UNPROTECT(nprotect);
135 return F;
136}
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()
|
extern |
Definition at line 273 of file dprocess.c.
274{
276 PROTECT(times=AS_NUMERIC(times));
277 PROTECT(x = as_state_array(x));
279 // extract the process function
281 // extract other arguments
283 PROTECT(covar = GET_SLOT(object,install("covar")));
284 // evaluate the density
286 UNPROTECT(7);
288}
static SEXP onestep_density(SEXP func, SEXP x, SEXP times, SEXP params, SEXP covar, SEXP log, SEXP args, SEXP gnsi)
Definition dprocess.c:128
Here is the call graph for this function:
◆ do_emeasure()
|
extern |
Definition at line 98 of file emeasure.c.
99{
102 int nobs = 0;
103 SEXP Snames, Pnames, Cnames, Onames = R_NilValue;
105 SEXP pompfun;
106 SEXP Y = R_NilValue;
107 int *dim;
108 lookup_table_t covariate_table;
109 SEXP cvec;
111
112 PROTECT(times = AS_NUMERIC(times));
113 ntimes = length(times);
114 if (ntimes < 1)
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])
123
125 dim = INTEGER(GET_DIM(params));
126 npars = dim[0]; nrepsp = dim[1];
127
128 nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
129
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)));
138 PROTECT(Onames = GET_SLOT(pompfun,install("obsnames")));
139
140 // set up the covariate table
142 PROTECT(cvec = NEW_NUMERIC(ncovars));
143 cov = REAL(cvec);
144
145 // extract the user-defined function
147
148 // extract 'userdata' as pairlist
150
151 int nprotect = 11;
152 int first = 1;
153
154 // first do setup
156
158 double *ys, *yt = 0;
160 SEXP ans;
161 int j, k;
162
164
165 for (k = 0; k < ntimes; k++, time++) { // loop over times
166
167 R_CheckUserInterrupt(); // check for user interrupt
168
170
172
173 if (first) {
174
175 PROTECT(
176 ans = eval_call(
178 time,
182 )
183 );
184
185 nobs = LENGTH(ans);
186
187 PROTECT(Onames = GET_NAMES(ans));
190
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(
208 time,
212 )
213 );
214
215 if (LENGTH(ans) != nobs)
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
235 double *yt = 0, *xp, *pp;
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
250
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
260
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);
278 const char *dimnm[3] = {"name",".id","time"};
279 double *yt = 0;
281
285
287
289 }
290
291 }
292
293 UNPROTECT(nprotect);
295}
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition emeasure.c:60
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition emeasure.c:10
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition emeasure.c:84
void pomp_emeasure(double *f, const double *x, const double *p, const int *obsindex, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition pomp.h:93
static R_INLINE void fixdimnames(SEXP x, const char **names, int n)
Definition pomp_defines.h:129
Here is the call graph for this function:
◆ do_partrans()
|
extern |
Definition at line 49 of file partrans.c.
50{
51
54 direction_t direc;
56 int *dim;
57
59
60 PROTECT(tparams = duplicate(params));
61
62 // coerce 'params' to matrix
63 PROTECT(tparams = as_matrix(tparams));
64 dim = INTEGER(GET_DIM(tparams));
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) {
74 PROTECT(pompfun = GET_SLOT(ob,install("from")));
75 break;
77 PROTECT(pompfun = GET_SLOT(ob,install("to")));
78 break;
79 }
80
82
83 // extract 'userdata' as pairlist
85
86 int nprotect = 7;
87
89
91
92 SEXP ans, nm;
93 double *pa, *ps = REAL(tparams);
94 int *posn;
95 int i, j;
96
99
100 PROTECT(nm = GET_NAMES(ans));
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
112
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
125
126 pomp_transform *ff;
127 double *ps, *pt;
128 int *idx;
129 int j;
130
132
133 R_CheckUserInterrupt();
134
135 idx = INTEGER(GET_SLOT(pompfun,install("paramindex")));
136
138 (*ff)(pt,ps,idx);
139
140 }
141
142 break;
143
144 default: // #nocov
145
146 break; // #nocov
147
148 }
149
150 if (qvec) {
151 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(tparams))); nprotect++;
152 SET_DIM(tparams,R_NilValue);
153 SET_NAMES(tparams,Pnames);
154 }
155
156 UNPROTECT(nprotect);
157 return tparams;
158
159}
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *p, int n)
Definition partrans.c:33
void pomp_transform(double *pt, const double *p, const int *parindex)
Definition pomp.h:115
static R_INLINE SEXP matchnames(SEXP provided, SEXP needed, const char *where)
Definition pomp_defines.h:60
Here is the call graph for this function:
◆ do_rinit()
|
extern |
Definition at line 91 of file rinit.c.
92{
93
94 SEXP Pnames, Cnames, Snames, pcnames;
95 SEXP x = R_NilValue;
98 lookup_table_t covariate_table;
99 SEXP cvec;
101 int *dim;
103
104 nsims = *(INTEGER(AS_INTEGER(nsim)));
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
116 PROTECT(cvec = NEW_NUMERIC(ncovars));
117 cov = REAL(cvec);
118
120
121 // extract userdata
123
124 PROTECT(pompfun = GET_SLOT(object,install("rinit")));
126
128
129 int nprotect = 9;
130
133
134 SEXP ans;
135 double *time = REAL(AS_NUMERIC(t0));
137 double *xs, *xt = NULL;
138 int *midx;
139 int j;
140
143 PROTECT(ans = AS_NUMERIC(ans));
144 PROTECT(Snames = GET_NAMES(ans));
145
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)
156 }
157
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) {
167 xs = REAL(ans);
168 if (LENGTH(ans) != nvar)
170 memcpy(xt,xs,nvar*sizeof(double));
171 UNPROTECT(1);
172 }
173
174 }
175
176 break;
177
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")));
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
194
195 GetRNGstate();
196
197 time = *(REAL(t0));
198
199 // loop over replicates
201 (*ff)(xt,ps+npar*(j%nrep),time,sidx,pidx,cidx,cov);
202
203 PutRNGstate();
204
205 }
206
207 break;
208
209 default: {
210
212
213 }
214
215 break;
216
217 }
218
219 // now add column names
220 if (nrep > 1) {
221 SEXP dn, xn;
222 int k, *p;
223
224 if (isNull(pcnames)) {
225 PROTECT(pcnames = NEW_INTEGER(nrep)); nprotect++;
226 for (k = 0, p = INTEGER(pcnames); k < nrep; k++, p++) *p = k+1;
227 }
228
229 if (nsims > 1) {
230 int k, *sp;
231 SEXP us;
232 PROTECT(us = mkString("_"));
233 PROTECT(xn = NEW_INTEGER(ns));
234 for (k = 0, sp = INTEGER(xn); k < ns; k++, sp++) *sp = (k/nrep)+1;
235 PROTECT(xn = paste0(pcnames,us,xn));
236 PROTECT(dn = GET_DIMNAMES(x));
237 nprotect += 4;
238 SET_ELEMENT(dn,1,xn);
239 SET_DIMNAMES(x,dn);
240
241 } else {
242
243 PROTECT(dn = GET_DIMNAMES(x)); nprotect++;
244 SET_ELEMENT(dn,1,pcnames);
245 SET_DIMNAMES(x,dn);
246
247 }
248 }
249
250 UNPROTECT(nprotect);
251 return x;
252}
void pomp_rinit(double *x, const double *p, double t0, const int *stateindex, const int *parindex, const int *covindex, const double *covars)
Definition pomp.h:40
static SEXP pomp_default_rinit(SEXP params, SEXP Pnames, int npar, int nrep, int nsim)
Definition rinit.c:254
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t0, double *p, int npar, double *c, int ncov)
Definition rinit.c:59
Here is the call graph for this function:
Here is the caller graph for this function:
◆ do_rmeasure()
|
extern |
Definition at line 98 of file rmeasure.c.
99{
102 int nobs = 0;
103 SEXP Snames, Pnames, Cnames, Onames = R_NilValue;
105 SEXP pompfun;
106 SEXP Y = R_NilValue;
107 int *dim;
108 lookup_table_t covariate_table;
109 SEXP cvec;
111
112 PROTECT(times = AS_NUMERIC(times));
113 ntimes = length(times);
114 if (ntimes < 1)
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])
123
125 dim = INTEGER(GET_DIM(params));
126 npars = dim[0]; nrepsp = dim[1];
127
128 nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
129
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)));
138 PROTECT(Onames = GET_SLOT(pompfun,install("obsnames")));
139
140 // set up the covariate table
142 PROTECT(cvec = NEW_NUMERIC(ncovars));
143 cov = REAL(cvec);
144
145 // extract the user-defined function
147
148 // extract 'userdata' as pairlist
150
151 int nprotect = 11;
152 int first = 1;
153
154 // first do setup
156
158 double *ys, *yt = 0;
160 SEXP ans;
161 int j, k;
162
164
165 for (k = 0; k < ntimes; k++, time++) { // loop over times
166
167 R_CheckUserInterrupt(); // check for user interrupt
168
170
172
173 if (first) {
174
175 PROTECT(
176 ans = eval_call(
178 time,
182 )
183 );
184
185 nobs = LENGTH(ans);
186
187 PROTECT(Onames = GET_NAMES(ans));
190
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(
208 time,
212 )
213 );
214
215 if (LENGTH(ans) != nobs)
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
235 double *yt = 0, *xp, *pp;
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
250
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
262
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);
282 const char *dimnm[3] = {"name",".id","time"};
283 double *yt = 0;
285
289
291
293 }
294
295 }
296
297 UNPROTECT(nprotect);
299}
void pomp_rmeasure(double *y, const double *x, const double *p, const int *obsindex, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition pomp.h:79
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition rmeasure.c:60
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition rmeasure.c:10
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition rmeasure.c:84
Here is the call graph for this function:
Here is the caller graph for this function:
◆ do_rprior()
|
extern |
Definition at line 58 of file rprior.c.
59{
60
64 int *dim;
65
67 dim = INTEGER(GET_DIM(params));
69
70 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
71
72 // extract the user-defined function
73 PROTECT(pompfun = GET_SLOT(object,install("rprior")));
75
76 // extract 'userdata' as pairlist
78
79 int nprotect = 5;
80 int first = 1;
81
83
85
86 SEXP ans, nm;
88 int *posn = NULL;
89 int i, j;
90
91 // set up the function call
93
95
96 if (first) {
97
99 PROTECT(ans = AS_NUMERIC(ans));
100
101 PROTECT(nm = GET_NAMES(ans));
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
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
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
141
142 R_CheckUserInterrupt(); // check for user interrupt
143
144 GetRNGstate();
145
146 // loop over replicates
148 (*ff)(p,pidx);
149
150 PutRNGstate();
151
152 }
153
154 break;
155
156 default: // just duplicate
157
159
160 }
161
162 UNPROTECT(nprotect);
164}
Here is the call graph for this function:
◆ do_rprocess()
|
extern |
Definition at line 25 of file rprocess.c.
26{
27
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) {
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
49 PROTECT(dimP = GET_DIM(params));
50 xdim = INTEGER(dimP);
52
53 int nprotect = 7;
54
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));
65 nprotect += 2;
67 src = REAL(copy);
68 tgt = REAL(params);
69 for (j = 0; j < nrepsx; j++) {
72 }
73 }
74 }
75 nreps = nrepsx;
79 } else {
80 double *src, *tgt;
81 int dims[2];
82 int j, k;
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);
92 *tgt = src[k+nvars*(j%nrepsx)];
93 }
94 }
95 }
96 }
97
98 PROTECT(rproc = GET_SLOT(object,install("rprocess")));
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) {
109 {
110 SEXP fn;
111 double deltat = 1.0;
114 accumvars,covar,args,gnsi));
115 nprotect += 2;
116 }
117 break;
119 {
120 SEXP fn;
121 double deltat;
123 deltat = *(REAL(AS_NUMERIC(GET_SLOT(rproc,install("delta.t")))));
125 accumvars,covar,args,gnsi));
126 nprotect += 2;
127 }
128 break;
130 {
131 SEXP fn, vmatrix, hmax;
133 PROTECT(vmatrix = GET_SLOT(rproc,install("v")));
134 PROTECT(hmax = GET_SLOT(rproc,install("hmax")));
136 accumvars,hmax,args,gnsi));
137 nprotect += 4;
138 }
139 break;
142 nprotect++;
143 break;
144 }
145
147 UNPROTECT(nprotect);
149
150}
SEXP euler_simulator(SEXP, SEXP, SEXP, SEXP, SEXP, double, rprocmode, SEXP, SEXP, SEXP, SEXP)
Definition euler.c:108
SEXP SSA_simulator(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP)
Definition ssa.c:226
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()
|
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
15
17
18 PROTECT(t0 = GET_SLOT(object,install("t0")));
19 PROTECT(times = GET_SLOT(object,install("times")));
20
21 // initialize the simulations
23 PROTECT(simnames = GET_COLNAMES(GET_DIMNAMES(x0)));
24
25 // call 'rprocess' to simulate state process
27
28 // call 'rmeasure' to simulate the measurement process
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;
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}
Here is the call graph for this function:
Here is the caller graph for this function:
◆ do_skeleton()
|
extern |
Definition at line 326 of file skeleton.c.
327{
328
331 int *dim;
332 SEXP Snames, Cnames, Pnames;
333 SEXP pompfun, cov, ob;
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])
345
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;
354
355 PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x)));
356 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
358
359 // set up the covariate table
362
363 // extract the user-defined function
364 PROTECT(ob = GET_SLOT(object,install("skeleton")));
365 PROTECT(pompfun = GET_SLOT(ob,install("skel.fn")));
367
368 // extract 'userdata' as pairlist
370
372
373 int nprotect = 12;
374
376
378
380
384 &covariate_table,REAL(cov));
385
386 }
387
388 break;
389
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
399
400 eval_skeleton_native(
401 REAL(F),REAL(t),REAL(x),REAL(params),
404
405 }
406
407 break;
408
409 default: {
410
411 double *ft = REAL(F);
413 for (i = 0; i < n; i++, ft++) *ft = R_NaReal;
415
416 }
417
418 }
419
420 UNPROTECT(nprotect);
421 return F;
422}
void pomp_skeleton(double *f, const double *x, const double *p, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition pomp.h:73
void eval_skeleton_native(double *f, double *time, double *x, double *p, int nvars, int npars, int ncovars, int ntimes, int nrepx, int nrepp, int nreps, int *sidx, int *pidx, int *cidx, lookup_table_t *covar_table, pomp_skeleton *fun, SEXP args, double *cov)
Definition skeleton.c:245
static R_INLINE SEXP ret_array(int nvars, int nreps, int ntimes, SEXP Snames)
Definition skeleton.c:81
SEXP add_skel_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition skeleton.c:10
void eval_skeleton_R(double *f, double *time, double *x, double *p, SEXP fn, SEXP args, SEXP Snames, int nvars, int npars, int ncovars, int ntimes, int nrepx, int nrepp, int nreps, lookup_table_t *covar_table, double *cov)
Definition skeleton.c:94
Here is the call graph for this function:
◆ do_vmeasure()
|
extern |
Definition at line 99 of file vmeasure.c.
100{
103 int nobs = 0, nans = 0;
104 SEXP Snames, Pnames, Cnames, Onames = R_NilValue;
106 SEXP pompfun;
107 SEXP Y = R_NilValue;
108 int *dim;
109 lookup_table_t covariate_table;
110 SEXP cvec;
112
113 PROTECT(times = AS_NUMERIC(times));
114 ntimes = length(times);
115 if (ntimes < 1)
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])
124
126 dim = INTEGER(GET_DIM(params));
127 npars = dim[0]; nrepsp = dim[1];
128
129 nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
130
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)));
139 PROTECT(Onames = GET_SLOT(pompfun,install("obsnames")));
140
141 // set up the covariate table
143 PROTECT(cvec = NEW_NUMERIC(ncovars));
144 cov = REAL(cvec);
145
146 // extract the user-defined function
148
149 // extract 'userdata' as pairlist
151
152 int nprotect = 11;
153 int first = 1;
154
155 // first do setup
157
159 double *ys, *yt = 0;
161 SEXP ans;
162 int j, k;
163
165
166 for (k = 0; k < ntimes; k++, time++) { // loop over times
167
168 R_CheckUserInterrupt(); // check for user interrupt
169
171
173
174 if (first) {
175
176 PROTECT(
177 ans = eval_call(
179 time,
183 )
184 );
185
186 nans = LENGTH(ans);
187
188 nobs = (int) sqrt(nans);
189 if (nans != nobs*nobs)
191
192 PROTECT(Onames = GET_ROWNAMES(GET_DIMNAMES(ans)));
195
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(
213 time,
217 )
218 );
219
220 if (LENGTH(ans) != nans)
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
240 double *yt = 0, *xp, *pp;
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
265
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
275
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);
293 const char *dimnm[4] = {"var1","var2",".id","time"};
294 double *yt = 0;
296
301
303
305 }
306
307 }
308
309 UNPROTECT(nprotect);
311}
void pomp_vmeasure(double *f, const double *x, const double *p, const int *vmatindex, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t)
Definition pomp.h:100
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition vmeasure.c:60
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition vmeasure.c:10
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition vmeasure.c:84
Here is the call graph for this function:
◆ E_Euler_Multinom()
|
extern |
Definition at line 71 of file distributions.c.
71 {
72 int ntrans = length(rate);
73 SEXP x;
74 if (length(size)>1)
76 if (length(deltat)>1)
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}
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()
|
extern |
Definition at line 107 of file euler.c.
112 {
113
118
120
121 int *dim;
124 ntimes = LENGTH(times);
125
126 PROTECT(tstart = AS_NUMERIC(tstart));
127 PROTECT(times = AS_NUMERIC(times));
128 t0 = *(REAL(tstart));
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
138 PROTECT(cvec = NEW_NUMERIC(ncovars));
139 cov = REAL(cvec);
140
141 // indices of accumulator variables
142 nzeros = LENGTH(accumvars);
144
145 // extract user function
147
148 // array to hold results
150
151 // copy the start values into the result array
153
154 // set up
155
156 int nprotect = 9;
157 int *pidx = 0, *sidx = 0, *cidx = 0;
158 pomp_onestep_sim *ff = NULL;
159
161
163 // construct list of all arguments
165 break;
166
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
174
175 GetRNGstate();
176 break;
177
178 default: // #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;
192 ) {
193
194 double dt;
195 int nstep = 0;
196 int i, j, k;
197
198 R_CheckUserInterrupt();
199
201
202 // set accumulator variables to zero
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) {
210 dt = time[step]-t;
211 nstep = 1;
212 break;
214 dt = deltat;
215 nstep = num_map_steps(t,time[step],dt);
216 break;
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
228
229 // loop over replicates
231
233
235
237 {
238 SEXP ans, nm;
239
240 if (first) {
241
243
244 PROTECT(nm = GET_NAMES(ans));
248
249 nprotect += 3;
250
251 ap = REAL(AS_NUMERIC(ans));
253
254 first = 0;
255
256 } else {
257
259
260 ap = REAL(AS_NUMERIC(ans));
262
263 UNPROTECT(1);
264
265 }
266 }
267 break;
268
270 (*ff)(xm,pm,sidx,pidx,cidx,cov,t,dt);
271 break;
272
273 default: // #nocov
275 break; // #nocov
276
277 }
278
279 }
280
281 t += dt;
282
284 dt = time[step]-t;
285 t = time[step]-dt;
286 }
287
288 }
289
290 if (step < ntimes-1)
292
293 }
294
295 // clean up
297
299
300 PutRNGstate();
301
302 }
303
304 break;
305
307
308 break;
309
310 }
311
312 UNPROTECT(nprotect);
314
315}
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *dt, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition euler.c:67
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition euler.c:10
static R_INLINE SEXP ret_array(int n, int nreps, int ntimes, SEXP names)
Definition euler.c:92
void pomp_onestep_sim(double *x, const double *p, const int *stateindex, const int *parindex, const int *covindex, const double *covars, double t, double dt)
Definition pomp.h:59
Here is the call graph for this function:
Here is the caller graph for this function:
◆ eval_skeleton_native()
|
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
261
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}
Here is the call graph for this function:
Here is the caller graph for this function:
◆ eval_skeleton_R()
|
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
115
117
118 if (first) {
119
124
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));
133
134 fs = REAL(AS_NUMERIC(ans));
136
137 nprotect += 3;
138 first = 0;
139
140 } else {
141
146
147 fs = REAL(AS_NUMERIC(ans));
149
150 UNPROTECT(1);
151
152 }
153
154 }
155 }
156
157 UNPROTECT(nprotect);
158
159}
static R_INLINE SEXP eval_call(SEXP fn, SEXP args, double *t, double *x, int nvar, double *p, int npar, double *c, int ncov)
Definition skeleton.c:57
Here is the call graph for this function:
Here is the caller graph for this function:
◆ ExpitTransform()
|
extern |
◆ get_covariate_names()
|
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()
|
extern |
Definition at line 13 of file userdata.c.
13 {
16 return elt;
17}
static R_INLINE SEXP getListElement(SEXP list, const char *str)
Definition pomp_defines.h:232
Here is the call graph for this function:
Here is the caller graph for this function:
◆ get_userdata_double()
|
extern |
Definition at line 26 of file userdata.c.
26 {
30 return REAL(elt);
31}
Here is the call graph for this function:
Here is the caller graph for this function:
◆ get_userdata_int()
|
extern |
Definition at line 19 of file userdata.c.
19 {
23 return INTEGER(elt);
24}
Here is the call graph for this function:
Here is the caller graph for this function:
◆ InverseLogBarycentricTransform()
|
extern |
Definition at line 47 of file transformations.c.
47 {
48 SEXP X;
52 UNPROTECT(2);
54}
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()
|
extern |
Definition at line 19 of file trajectory.c.
20{
21
25 SEXP Snames, Pnames, Cnames;
26 SEXP skel, pompfun;
27 SEXP accumvars;
28 SEXP repnames;
29 int *zidx = 0;
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));
43
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)));
55
56 // set up the covariate table
59
60 // extract user-defined function
61 PROTECT(pompfun = GET_SLOT(skel,install("skel.fn")));
63
64 // extract 'userdata' as pairlist
66
67 // create array to store results
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) {
78 }
79
80 // set up the computations
82
84
86
89 &covariate_table,zidx,REAL(cov));
90
91 }
92
93 break;
94
96 int *sidx, *pidx, *cidx;
97 pomp_skeleton *ff;
98
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
109
110 }
111
112 break;
113
114 default: // #nocov
115
116 break; // #nocov
117
118 }
119
120 UNPROTECT(nprotect);
122}
void iterate_skeleton_R(double *, double, double, double *, double *, double *, SEXP, SEXP, SEXP, int, int, int, int, int, int, int, lookup_table_t *, int *, double *)
Definition skeleton.c:161
void iterate_skeleton_native(double *, double, double, double *, double *, double *, int, int, int, int, int, int, int, int *, int *, int *, lookup_table_t *, int *, pomp_skeleton *, SEXP, double *)
Definition skeleton.c:274
static R_INLINE SEXP ret_array(int nvars, int nreps, int ntimes, SEXP Snames, SEXP repnames)
Definition trajectory.c:6
Here is the call graph for this function:
◆ iterate_skeleton_native()
|
extern |
Definition at line 274 of file skeleton.c.
282{
283 int nsteps;
284 double *xs, *Xs;
285 int h, i, j, k;
286
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++)
297 *xs = 0.0;
298
299 for (h = 0; h < nsteps; h++) {
300
301 // interpolate the covariates
303
305
307
308 }
309
311
312 if (h != nsteps-1) {
313 t += deltat;
314 } else {
315 t = *time;
316 }
317
318 }
319
321
322 }
323
324}
Here is the call graph for this function:
Here is the caller graph for this function:
◆ iterate_skeleton_R()
|
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
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++)
189 *xs = 0.0;
190
191 for (h = 0; h < nsteps; h++) {
192
193 // interpolate the covariates
195
197
198 if (first) {
199
201
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));
209
210 ap = REAL(AS_NUMERIC(ans));
212
213 nprotect += 3;
214 first = 0;
215
216 } else {
217
219
220 ap = REAL(AS_NUMERIC(ans));
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
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()
|
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));
139 UNPROTECT(1);
140 return s;
141}
Here is the caller graph for this function:
◆ load_stack_incr()
|
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()
|
extern |
Definition at line 38 of file transformations.c.
38 {
39 SEXP Y;
43 UNPROTECT(2);
45}
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()
|
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}
Here is the call graph for this function:
◆ logmeanexp()
|
extern |
Definition at line 7 of file logmeanexp.c.
7 {
9 int k = *INTEGER(Drop)-1; // zero-based index
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()
|
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
35
36 if (nt > 1) {
37 xdim[0] = nvar; xdim[1] = nt;
40 } else {
41 PROTECT(X = NEW_NUMERIC(nvar));
42 SET_NAMES(X,Cnames);
43 }
44
46 table_lookup(&tab,*tp,xp);
47
48 UNPROTECT(3);
50}
lookup_table_t make_covariate_table(SEXP object, int *ncovar)
Definition lookup_table.c:11
Here is the call graph for this function:
Here is the caller graph for this function:
◆ make_covariate_table()
|
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;
21 return tab;
22}
Here is the caller graph for this function:
◆ nosort_resamp()
|
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)
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()
|
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()
|
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()
|
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}
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()
|
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
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()
|
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;
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);
40 PROTECT(Xnames = GET_ROWNAMES(GET_DIMNAMES(x)));
41
43 PROTECT(dimP = GET_DIM(params));
44 dim = INTEGER(dimP);
45 npars = dim[0];
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) {
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
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
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
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;
136 }
137
138 // compute filter mean
139 if (do_fm) {
141 }
142
143 // compute weighted average of parameters
144 if (do_wave) {
145 if (all_fail)
147 "using unweighted mean for point estimate.");
148 double *xwm = REAL(wmean);
151 if (all_fail) { // unweighted average
154 } else {
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
190 for (k = 0; k < np; k++) { // copy the particles
191 int sp = sample[k];
192 double *xx, *xp;
194 *st = *xx;
195 if (do_pr) {
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);
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}
static void filt_mean(int, int, int, long double, const double *, const double *, double *)
Definition pfilter.c:294
static void pred_mean_var(int, int, int, const double *, double *, double *)
Definition pfilter.c:266
Here is the call graph for this function:
◆ pomp_desolve_setup()
|
extern |
Definition at line 156 of file trajectory.c.
156 {
157
160 SEXP Snames, Pnames, Cnames;
161 SEXP pompfun, ob;
162 int *dim;
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
170
177
178 dim = INTEGER(GET_DIM(x0));
179 nvars = dim[0];
181
182 dim = INTEGER(GET_DIM(params));
186
187 // set up the covariate table
192
193 PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x0)));
194 PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params)));
196
198
200
201 int nprotect = 8;
202
204
206
210 nprotect += 3;
211
215
219
220 }
221
222 break;
223
225
227
231 nprotect += 3;
232
234
239
244
245 }
246
247 break;
248
249 default: // #nocov
250
252
253 }
254
255 UNPROTECT(nprotect);
256
257 return R_NilValue;
258}
Here is the call graph for this function:
◆ pomp_desolve_takedown()
|
extern |
Definition at line 292 of file trajectory.c.
292 {
303
305
307
314
315 }
316
317 break;
318
320
330
331 }
332
333 break;
334
335 default: // #nocov
336
338
339 break;
340
341 }
342
344
345 return R_NilValue;
346}
◆ pomp_fun_handler()
|
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
38
40
42
43 PROTECT(f = GET_SLOT(pfun,install("R.fun"))); nprotect++;
44
45 break;
46
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
61
62 SEXP nsi;
63 PROTECT(nsi = eval(PROTECT(lang3(install("getNativeSymbolInfo"),nf,pack)),R_BaseEnv));
65 nprotect += 3;
66
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) {
82 SET_SLOT(pfun,install("stateindex"),sidx);
83 }
84
85 if (P != NA_STRING) {
87 SET_SLOT(pfun,install("paramindex"),pidx);
88 }
89
90 if (O != NA_STRING) {
92 SET_SLOT(pfun,install("obsindex"),oidx);
93 }
94
95 if (C != NA_STRING) {
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
109
110 PROTECT(f = R_NilValue); nprotect++;
112
113 break;
114
115 }
116
117 UNPROTECT(nprotect);
118 return f;
119}
static R_INLINE SEXP name_index(SEXP provided, SEXP object, const char *slot, const char *humanreadable)
Definition pomp_fun.c:10
Here is the call graph for this function:
Here is the caller graph for this function:
◆ pomp_vf_eval()
|
extern |
Definition at line 260 of file trajectory.c.
261{
263
265
271
272 break;
273
280
281 break;
282
283 default: // #nocov
284
286
287 break;
288
289 }
290}
void eval_skeleton_R(double *, double *, double *, double *, SEXP, SEXP, SEXP, int, int, int, int, int, int, int, lookup_table_t *, double *)
Definition skeleton.c:94
void eval_skeleton_native(double *, double *, double *, double *, int, int, int, int, int, int, int, int *, int *, int *, lookup_table_t *, pomp_skeleton *, SEXP, double *)
Definition skeleton.c:245
Here is the call graph for this function:
◆ probe_acf()
|
extern |
Definition at line 100 of file probe_acf.c.
100 {
101 SEXP ans, ans_names;
102 SEXP X;
104 int j, k, l;
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
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
122
123 if (correlation) {
124 l = 0;
128 for (k = 0; k < nlag; k++, p++)
129 *p /= *p1;
130 }
131
132 PROTECT(ans_names = NEW_STRING(nlag*nvars));
134 for (k = 0; k < nlag; k++, l++) {
135 snprintf(tmp,BUFSIZ,"acf[%d]",lag[k]);
136 SET_STRING_ELT(ans_names,l,mkChar(tmp));
137 }
138 }
139 SET_NAMES(ans,ans_names);
140
141 UNPROTECT(4);
142 return ans;
143}
static void pomp_acf_compute(double *acf, double *x, int n, int nvars, int *lags, int nlag)
Definition probe_acf.c:11
Here is the call graph for this function:
◆ probe_ccf()
|
extern |
Definition at line 145 of file probe_acf.c.
145 {
146 SEXP ans, ans_names;
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
159
160 PROTECT(X = duplicate(AS_NUMERIC(x)));
161 PROTECT(Y = duplicate(AS_NUMERIC(y)));
162
163 PROTECT(ans = NEW_NUMERIC(nlag));
164
166
167 if (correlation) {
168 k = 0;
172 for (k = 0, p = REAL(ans); k < nlag; k++, p++) *p /= xx; // convert to correlation
173 }
174
175 PROTECT(ans_names = NEW_STRING(nlag));
176 for (k = 0; k < nlag; k++) {
177 snprintf(tmp,BUFSIZ,"ccf[%d]",INTEGER(lags)[k]);
178 SET_STRING_ELT(ans_names,k,mkChar(tmp));
179 }
180 SET_NAMES(ans,ans_names);
181
182 UNPROTECT(5);
183 return ans;
184}
static void pomp_ccf_compute(double *ccf, double *x, double *y, int n, int *lags, int nlag)
Definition probe_acf.c:47
Here is the call graph for this function:
◆ probe_marginal_setup()
|
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
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
38
39 UNPROTECT(6);
40 return(retval);
41}
static void order_reg_model_matrix(double *z, double *X, double *tau, int *pivot, int n, int np, int diff)
Definition probe_marginal.c:78
Here is the call graph for this function:
◆ probe_marginal_solve()
|
extern |
Definition at line 43 of file probe_marginal.c.
43 {
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
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
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
72
73 UNPROTECT(6);
74 return(beta);
75}
static void order_reg_solve(double *beta, double *x, double *mm, double *tau, int *pivot, int n, int np, int diff)
Definition probe_marginal.c:117
Here is the call graph for this function:
◆ probe_nlar()
|
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}
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()
|
extern |
Definition at line 119 of file distributions.c.
119 {
120 int k, nval, ns, np, nt;
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++) {
131 }
132 PutRNGstate();
133 UNPROTECT(5);
134 return ans;
135}
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()
|
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)
29 if (length(deltat)>1)
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)
39 PROTECT(x = makearray(2,dim));
40 PROTECT(nm = GET_NAMES(rate));
41 setrownames(x,nm,2);
42 GetRNGstate();
43 reulermultinom_multi(dim[0],dim[1],REAL(size),REAL(rate),REAL(deltat),REAL(x));
44 PutRNGstate();
45 UNPROTECT(6);
46 return x;
47}
static void reulermultinom_multi(int m, int n, double *size, double *rate, double *deltat, double *trans)
Definition distributions.c:7
Here is the call graph for this function:
◆ R_GammaWN()
|
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];
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}
Here is the call graph for this function:
◆ R_init_pomp()
|
extern |
Definition at line 54 of file init.c.
54 {
55 // C functions provided for users
57 R_RegisterCCallable("pomp","periodic_bspline_basis_eval_deriv",(DL_FUNC) &periodic_bspline_basis_eval_deriv);
72
73 // Register routines
74 R_registerRoutines(info,NULL,callMethods,NULL,NULL);
75 R_useDynamicSymbols(info,TRUE);
76 R_forceSymbols(info,FALSE);
77}
void bspline_basis_eval_deriv(double x, double *knots, int degree, int nbasis, int deriv, double *y)
Definition bspline.c:115
Here is the call graph for this function:
◆ randwalk_perturbation()
|
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;
12 int j, k;
13
15
16 // unpack parameter matrix
17 xp = REAL(params);
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
27 pidx = INTEGER(pindex);
28
29 GetRNGstate();
30
31 for (j = 0, xrw = rw; j < nrw; j++, pidx++, xrw++) {
33 *xs += *xrw * norm_rand();
34 }
35 }
36
37 PutRNGstate();
38
39 UNPROTECT(4);
41}
Here is the call graph for this function:
Here is the caller graph for this function:
◆ set_pomp_userdata()
|
extern |
Definition at line 8 of file userdata.c.
◆ sobol_sequence()
|
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;
243 PROTECT(data = allocMatrix(REALSXP,d,n)); dp = REAL(data);
244 nlopt_sobol_skip(s,n,dp);
246 nlopt_sobol_destroy(s);
247 UNPROTECT(1);
248 return(data);
249}
static void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x)
Definition sobolseq.c:225
Here is the call graph for this function:
◆ SSA_simulator()
|
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;
234 int nzeros = LENGTH(accumvars);
236 SEXP X, zindex, vindex;
237 int *sidx, *pidx, *cidx, *zidx, *vidx;
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)
249
250 ntimes = LENGTH(times);
251
252 PROTECT(tstart = AS_NUMERIC(tstart));
253 PROTECT(times = AS_NUMERIC(times));
254 t0 = *(REAL(tstart));
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
274
275 int nprotect = 11;
276
278
280
282
284
286
287 break;
288
290
292 NVAR = nvar;
293 NPAR = npar;
294 NCOV = covdim;
297 nprotect += 2;
298 FIRST = 1;
299
300 break;
301
302 }
303
304 xdim[0] = nvar; xdim[1] = nrep; xdim[2] = ntimes;
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) {
314 zidx = INTEGER(zindex);
315 } else {
316 zidx = 0;
317 }
318
319 if (hasvnames) {
321 vidx = INTEGER(vindex);
322 } else {
323 vidx = 0;
324 }
325
326 GetRNGstate();
327 {
328 int i;
329 for (i = 0; i < nrep; i++) {
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);
341}
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
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
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
static R_INLINE SEXP add_args(SEXP args, SEXP Snames, SEXP Pnames, SEXP Cnames)
Definition ssa.c:9
Here is the call graph for this function:
Here is the caller graph for this function:
◆ synth_loglik()
|
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}
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()
|
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}
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()
|
extern |
Definition at line 53 of file lookup_table.c.
54{
55 int flag = 0;
56 int j, k, n;
57 double e;
60 // warn only if we are *outside* the interval
64 case 1: default: // linear interpolation
68 }
69 break;
70 case 0: // piecewise constant
71 if (flag < 0) n = 0;
75 y[j] = tab->y[k];
76 }
77 break;
78 }
79}
Here is the caller graph for this function:
◆ wpfilter()
|
extern |
Definition at line 10 of file wpfilter.c.
11{
12
13 SEXP dimX, dimP, Xnames, Pnames;
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);
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];
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);
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);
77 SET_DIM(X,newdim);
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);
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
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 {
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
165 *xw = *xW;
166 *xW *= target;
167 *xw = exp(*xw - *xW);
168 }
169
170 // do resampling
171 {
172 int sample[np];
173 GetRNGstate();
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;
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:
