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