[Pomp-commits] r1200 - in pkg/pomp: . R inst inst/examples src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 11 13:56:30 CEST 2015
Author: kingaa
Date: 2015-06-11 13:56:30 +0200 (Thu, 11 Jun 2015)
New Revision: 1200
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/builder.R
pkg/pomp/R/init-state-pomp.R
pkg/pomp/R/pomp-class.R
pkg/pomp/R/pomp.R
pkg/pomp/inst/NEWS
pkg/pomp/inst/NEWS.Rd
pkg/pomp/inst/examples/bbs.R
pkg/pomp/inst/examples/euler.sir.R
pkg/pomp/src/initstate.c
pkg/pomp/src/pomp_internal.h
pkg/pomp/src/simulate.c
pkg/pomp/src/sir.c
Log:
- modify basic pomp structure and 'init.state' method to support Csnippets for the initializer
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/DESCRIPTION 2015-06-11 11:56:30 UTC (rev 1200)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical Inference for Partially Observed Markov Processes
-Version: 0.66-10
-Date: 2015-06-08
+Version: 0.67-1
+Date: 2015-06-10
Authors at R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="kingaa at umich.edu"),
person(given=c("Edward","L."),family="Ionides",role=c("aut")),
Modified: pkg/pomp/R/builder.R
===================================================================
--- pkg/pomp/R/builder.R 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/R/builder.R 2015-06-11 11:56:30 UTC (rev 1200)
@@ -2,6 +2,7 @@
statenames, paramnames, covarnames, obsnames,
rmeasure, dmeasure, step.fn, skeleton,
fromEstimationScale, toEstimationScale,
+ initializer,
rprior, dprior, globals,
verbose = getOption("verbose",FALSE))
{
@@ -65,6 +66,14 @@
## list of functions to register
registry <- c("load_stack_incr","load_stack_decr")
+ ## initializer function
+ if (!missing(initializer)) {
+ registry <- c(registry,"initializer")
+ cat(file=out,render(header$initializer,name=name))
+ cat(file=out,callable.decl(initializer))
+ cat(file=out,initializer,footer$initializer)
+ }
+
## parameter transformation function
if (!missing(fromEstimationScale)) {
registry <- c(registry,"par_trans")
@@ -317,6 +326,7 @@
dmeasure= "\nvoid __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, int __ncovars, const double *__covars, double t)\n{\n",
step.fn="\nvoid __pomp_stepfn (double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, int __covdim, const double *__covars, double t, double dt)\n{\n",
skeleton="\nvoid __pomp_skelfn (double *__f, const double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, int __ncovars, const double *__covars, double t)\n{\n",
+ initializer="\nvoid __pomp_initializer (double *__x, const double *__p, double t, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars)\n{\n",
fromEstimationScale="\nvoid __pomp_par_trans (double *__pt, const double *__p, const int *__parindex)\n{\n",
toEstimationScale="\nvoid __pomp_par_untrans (double *__pt, const double *__p, const int *__parindex)\n{\n",
rprior="\nvoid __pomp_rprior (double *__p, const int *__parindex)\n{\n",
@@ -328,6 +338,7 @@
dmeasure= "__pomp_dmeasure",
step.fn="__pomp_stepfn",
skeleton="__pomp_skelfn",
+ skeleton="__pomp_initializer",
fromEstimationScale="__pomp_par_trans",
toEstimationScale="__pomp_par_untrans",
rprior="__pomp_rprior",
@@ -348,6 +359,7 @@
dmeasure="\n}\n\n",
step.fn="\n}\n\n",
skeleton="\n}\n\n",
+ initializer="\n}\n\n",
fromEstimationScale="\n}\n\n",
toEstimationScale="\n}\n\n",
rprior="\n}\n\n",
Modified: pkg/pomp/R/init-state-pomp.R
===================================================================
--- pkg/pomp/R/init-state-pomp.R 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/R/init-state-pomp.R 2015-06-11 11:56:30 UTC (rev 1200)
@@ -1,12 +1,13 @@
## initialize the state variables of the process model
-init.state.internal <- function (object, params, t0, ...) {
+init.state.internal <- function (object, params, t0,
+ .getnativesymbolinfo = TRUE, ...) {
if (missing(t0)) t0 <- object at t0
if (missing(params)) params <- coef(object)
pompLoad(object)
- rv <- .Call(do_init_state,object,params,t0)
+ x <- .Call(do_init_state,object,params,t0,.getnativesymbolinfo)
pompUnload(object)
- rv
+ x
}
setMethod('init.state',
Modified: pkg/pomp/R/pomp-class.R
===================================================================
--- pkg/pomp/R/pomp-class.R 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/R/pomp-class.R 2015-06-11 11:56:30 UTC (rev 1200)
@@ -1,14 +1,3 @@
-## this is the initial-condition setting function that is used by default
-## it simply finds all parameters in the vector 'params' that have a name ending in '.0'
-## and returns a vector with their values with names stripped of '.0'
-default.initializer <- function (params, t0, ...) {
- ivpnames <- grep("\\.0$",names(params),value=TRUE)
- if (length(ivpnames)<1)
- stop("default initializer error: no parameter names ending in ",
- sQuote(".0")," found: see ",sQuote("pomp")," documentation")
- setNames(params[ivpnames],sub("\\.0$","",ivpnames))
-}
-
## define the pomp class
setClass(
'pomp',
@@ -25,7 +14,8 @@
skeleton.type = 'character',
skeleton = 'pomp.fun',
skelmap.delta.t = 'numeric',
- initializer = 'function',
+ default.init = 'logical',
+ initializer = 'pomp.fun',
states = 'array',
params = 'numeric',
covar = 'matrix',
@@ -50,7 +40,8 @@
skeleton.type="map",
skeleton=pomp.fun(),
skelmap.delta.t=as.numeric(NA),
- initializer=default.initializer,
+ default.init=TRUE,
+ initializer=pomp.fun(),
states=array(data=numeric(0),dim=c(0,0)),
params=numeric(0),
covar=array(data=numeric(0),dim=c(0,0)),
@@ -97,14 +88,6 @@
sQuote("dprocess(x,times,params,log,...)")
)
)
- if (!all(c('params','t0','...')%in%names(formals(object at initializer))))
- retval <- append(
- retval,
- paste(
- sQuote("initializer"),"must be a function of prototype",
- sQuote("initializer(params,t0,...)")
- )
- )
if (length(object at tcovar)!=nrow(object at covar)) {
retval <- append(
retval,
Modified: pkg/pomp/R/pomp.R
===================================================================
--- pkg/pomp/R/pomp.R 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/R/pomp.R 2015-06-11 11:56:30 UTC (rev 1200)
@@ -142,11 +142,25 @@
storage.mode(tcovar) <- "double"
storage.mode(covar) <- "double"
- ## check initializer
- if (missing(initializer)) initializer <- default.initializer
- if (!is.function(initializer))
- stop("pomp error: ",sQuote("initializer")," must be a function")
+ ## handle initializer
+ default.init <- missing(initializer) || (is(initializer,"pomp.fun") && initializer at mode == pompfunmode$undef)
+ if (default.init) {
+ initializer <- pomp.fun()
+ } else {
+ initializer <- pomp.fun(
+ f=initializer,
+ PACKAGE=PACKAGE,
+ proto=quote(initializer(params,t0,...)),
+ slotname="initializer",
+ libname=libname,
+ statenames=statenames,
+ paramnames=paramnames,
+ obsnames=obsnames,
+ covarnames=covarnames
+ )
+ }
+
## default rprocess & dprocess
if (missing(rprocess))
rprocess <- function (xstart,times,params,...) stop(sQuote("rprocess")," not specified")
@@ -173,6 +187,8 @@
snips <- c(snips,fromEstimationScale=fromEstimationScale at text)
if (is(toEstimationScale,"Csnippet"))
snips <- c(snips,toEstimationScale=toEstimationScale at text)
+ if (is(initializer,"Csnippet"))
+ snips <- c(snips,initializer=initializer at text)
if (length(snips)>0) {
libname <- try(
do.call(
@@ -412,6 +428,7 @@
data = data,
times = times,
t0 = t0,
+ default.init = default.init,
initializer = initializer,
params = params,
covar = covar,
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/inst/NEWS 2015-06-11 11:56:30 UTC (rev 1200)
@@ -1,5 +1,12 @@
_N_e_w_s _f_o_r _p_a_c_k_a_g_e '_p_o_m_p'
+_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _0._6_7-_1:
+
+ • A 'pomp' object's initializer can now be specified as a
+ Csnippet.
+
+ • The default initializer is now considerably faster.
+
_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _0._6_6-_4:
• A full-featured version of IF2, an improved iterated
Modified: pkg/pomp/inst/NEWS.Rd
===================================================================
--- pkg/pomp/inst/NEWS.Rd 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/inst/NEWS.Rd 2015-06-11 11:56:30 UTC (rev 1200)
@@ -1,5 +1,11 @@
\name{NEWS}
\title{News for package `pomp'}
+\section{Changes in \pkg{pomp} version 0.67-1}{
+ \itemize{
+ \item A 'pomp' object's initializer can now be specified as a Csnippet.
+ \item The default initializer is now considerably faster.
+ }
+}
\section{Changes in \pkg{pomp} version 0.66-4}{
\itemize{
\item A full-featured version of IF2, an improved iterated filtering algorithm, is now available as \code{mif2}.
Modified: pkg/pomp/inst/examples/bbs.R
===================================================================
--- pkg/pomp/inst/examples/bbs.R 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/inst/examples/bbs.R 2015-06-11 11:56:30 UTC (rev 1200)
@@ -46,8 +46,6 @@
"S.0","I.0","R.0"
),
zeronames=c("cases"),
- comp.names=c("S","I","R"),
- ic.names=c("S.0","I.0","R.0"),
nbasis=1L,
degree=0L,
period=1.0,
@@ -66,14 +64,7 @@
params[logitvar] <- plogis(params[logitvar])
params
},
- initializer=function(params, t0, comp.names, ic.names, ...) {
- snames <- c("S","I","R","cases","W")
- fracs <- params[ic.names]
- x0 <- numeric(length(snames))
- names(x0) <- snames
- x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
- x0
- }
+ initializer="_sir_init"
) -> bbs
c("bbs")
Modified: pkg/pomp/inst/examples/euler.sir.R
===================================================================
--- pkg/pomp/inst/examples/euler.sir.R 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/inst/examples/euler.sir.R 2015-06-11 11:56:30 UTC (rev 1200)
@@ -240,21 +240,12 @@
"S.0","I.0","R.0"
),
zeronames=c("cases"),
- comp.names=c("S","I","R"),
- ic.names=c("S.0","I.0","R.0"),
fromEstimationScale="_sir_par_trans",
toEstimationScale="_sir_par_untrans",
nbasis=3L,
degree=3L,
period=1.0,
- initializer=function(params, t0, comp.names, ic.names, ...) {
- snames <- c("S","I","R","cases","W")
- fracs <- params[ic.names]
- x0 <- numeric(length(snames))
- names(x0) <- snames
- x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
- x0
- }
+ initializer="_sir_init"
) -> euler.sir
## the following was originally used to generate the data
Modified: pkg/pomp/src/initstate.c
===================================================================
--- pkg/pomp/src/initstate.c 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/src/initstate.c 2015-06-11 11:56:30 UTC (rev 1200)
@@ -7,89 +7,220 @@
#include "pomp_internal.h"
-SEXP do_init_state (SEXP object, SEXP params, SEXP t0)
+// initializer.c
+typedef double pomp_initializer(double *x, const double *p, double t,
+ const int *stateindex, const int *parindex, const int *covindex,
+ const double *covars);
+
+SEXP do_init_state (SEXP object, SEXP params, SEXP t0, SEXP gnsi)
{
int nprotect = 0;
- SEXP x, x1, x2, par;
- SEXP fcall, fn, rho;
- SEXP covar, tcovar, covars;
- SEXP paramnames, statenames, mindex;
- int *dim, *midx;
+ SEXP Pnames, Snames, x;
+ int *dim;
int npar, nrep, nvar;
- int xdim[2], j, k;
- double *p, *pp, *xp, *xpp;
+ int definit;
+ int xdim[2];
const char *dimnms[2] = {"variable","rep"};
PROTECT(params = as_matrix(params)); nprotect++;
+ PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
dim = INTEGER(GET_DIM(params));
npar = dim[0]; nrep = dim[1];
- p = REAL(params);
- PROTECT(par = NEW_NUMERIC(npar)); nprotect++;
- PROTECT(paramnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
- SET_NAMES(par,paramnames);
- pp = REAL(par);
+ definit = *(INTEGER(GET_SLOT(object,install("default.init"))));
- // extract the initializer function and its environment
- PROTECT(fn = GET_SLOT(object,install("initializer"))); nprotect++;
- PROTECT(rho = (CLOENV(fn))); nprotect++;
+ if (definit) { // default initializer
- // begin to construct the call
- PROTECT(fcall = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
+ SEXP fcall, pat, repl, val, ivpnames, statenames;
+ int *pidx, j, k;
+ double *xp, *pp;
+
+ PROTECT(pat = NEW_CHARACTER(1)); nprotect++;
+ SET_STRING_ELT(pat,0,mkChar("\\.0$"));
+ PROTECT(repl = NEW_CHARACTER(1)); nprotect++;
+ SET_STRING_ELT(repl,0,mkChar(""));
+ PROTECT(val = NEW_LOGICAL(1)); nprotect++;
+ *(INTEGER(val)) = 1;
+
+ // extract names of IVPs
+ PROTECT(fcall = LCONS(val,R_NilValue)); nprotect++;
+ SET_TAG(fcall,install("value"));
+ PROTECT(fcall = LCONS(Pnames,fcall)); nprotect++;
+ SET_TAG(fcall,install("x"));
+ PROTECT(fcall = LCONS(pat,fcall)); nprotect++;
+ SET_TAG(fcall,install("pattern"));
+ PROTECT(fcall = LCONS(install("grep"),fcall)); nprotect++;
+ PROTECT(ivpnames = eval(fcall,R_BaseEnv)); nprotect++;
+
+ nvar = LENGTH(ivpnames);
+ if (nvar < 1) {
+ error("default initializer error: no parameter names ending in '.0' found: see 'pomp' documentation");
+ }
+ pidx = INTEGER(PROTECT(match(Pnames,ivpnames,0))); nprotect++;
+ for (k = 0; k < nvar; k++) pidx[k]--;
+
+ // construct names of state variables
+ PROTECT(fcall = LCONS(ivpnames,R_NilValue)); nprotect++;
+ SET_TAG(fcall,install("x"));
+ PROTECT(fcall = LCONS(repl,fcall)); nprotect++;
+ SET_TAG(fcall,install("replacement"));
+ PROTECT(fcall = LCONS(pat,fcall)); nprotect++;
+ SET_TAG(fcall,install("pattern"));
+ PROTECT(fcall = LCONS(install("sub"),fcall)); nprotect++;
+ PROTECT(statenames = eval(fcall,R_BaseEnv)); nprotect++;
- // extract covariates and interpolate
- PROTECT(tcovar = GET_SLOT(object,install("tcovar"))); nprotect++;
- if (LENGTH(tcovar) > 0) { // do table lookup
- PROTECT(covar = GET_SLOT(object,install("covar"))); nprotect++;
- PROTECT(covars = lookup_in_table(tcovar,covar,t0)); nprotect++;
- PROTECT(fcall = LCONS(covars,fcall)); nprotect++;
- SET_TAG(fcall,install("covars"));
- }
+ xdim[0] = nvar; xdim[1] = nrep;
+ PROTECT(x = makearray(2,xdim)); nprotect++;
+ setrownames(x,statenames,2);
+ fixdimnames(x,dimnms,2);
- // finish constructing the call
- PROTECT(fcall = LCONS(t0,fcall)); nprotect++;
- SET_TAG(fcall,install("t0"));
- PROTECT(fcall = LCONS(par,fcall)); nprotect++;
- SET_TAG(fcall,install("params"));
- PROTECT(fcall = LCONS(fn,fcall)); nprotect++;
+ for (j = 0, xp = REAL(x), pp = REAL(params); j < nrep; j++, pp += npar)
+ for (k = 0; k < nvar; k++, xp++)
+ *xp = pp[pidx[k]];
- for (k = 0; k < npar; k++) pp[k] = p[k];
- PROTECT(x1 = eval(fcall,rho)); nprotect++; // do the call
- PROTECT(statenames = GET_NAMES(x1)); nprotect++;
-
- if (!IS_NUMERIC(x1) || isNull(statenames)) {
- UNPROTECT(nprotect);
- error("init.state error: user 'initializer' must return a named numeric vector");
- }
+ } else { // user-supplied initializer
+
+ SEXP pompfun, fcall, fn, tcovar, covar, covars;
+ pompfunmode mode = undef;
+ double *cp = NULL;
- nvar = LENGTH(x1);
- xp = REAL(x1);
- PROTECT(mindex = match(paramnames,statenames,0)); nprotect++;
- midx = INTEGER(mindex);
-
- for (k = 0; k < nvar; k++) {
- if (midx[k]!=0) {
- UNPROTECT(nprotect);
- error("a state variable and a parameter share a single name: '%s'",CHARACTER_DATA(STRING_ELT(statenames,k)));
+ // extract the initializer function and its environment
+ PROTECT(pompfun = GET_SLOT(object,install("initializer"))); nprotect++;
+ PROTECT(fn = pomp_fun_handler(pompfun,gnsi,&mode)); nprotect++;
+
+ // extract covariates and interpolate
+ PROTECT(tcovar = GET_SLOT(object,install("tcovar"))); nprotect++;
+ if (LENGTH(tcovar) > 0) { // do table lookup
+ PROTECT(covar = GET_SLOT(object,install("covar"))); nprotect++;
+ PROTECT(covars = lookup_in_table(tcovar,covar,t0)); nprotect++;
+ cp = REAL(covars);
}
- }
+
+ // extract userdata
+ PROTECT(fcall = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
+
+ switch (mode) {
+ case Rfun: // use R function
- xdim[0] = nvar; xdim[1] = nrep;
- PROTECT(x = makearray(2,xdim)); nprotect++;
- setrownames(x,statenames,2);
- fixdimnames(x,dimnms,2);
- xpp = REAL(x);
+ {
+ SEXP par, rho, x1, x2, mindex;
+ double *p, *pp, *xp, *xt;
+ int j, *midx;
- for (k = 0; k < nvar; k++) xpp[k] = xp[k];
+ // extract covariates and interpolate
+ if (LENGTH(tcovar) > 0) { // add covars to call
+ PROTECT(fcall = LCONS(covars,fcall)); nprotect++;
+ SET_TAG(fcall,install("covars"));
+ }
+
+ // parameter vector
+ PROTECT(par = NEW_NUMERIC(npar)); nprotect++;
+ SET_NAMES(par,Pnames);
+ pp = REAL(par);
+
+ // finish constructing the call
+ PROTECT(fcall = LCONS(t0,fcall)); nprotect++;
+ SET_TAG(fcall,install("t0"));
+ PROTECT(fcall = LCONS(par,fcall)); nprotect++;
+ SET_TAG(fcall,install("params"));
+ PROTECT(fcall = LCONS(fn,fcall)); nprotect++;
+
+ // evaluation environment
+ PROTECT(rho = (CLOENV(fn))); nprotect++;
- for (j = 1; j < nrep; j++) {
- for (k = 0; k < npar; k++) pp[k] = p[j*npar+k];
- PROTECT(x2 = eval(fcall,rho));
- xp = REAL(x2);
- if (LENGTH(x2)!=nvar)
- error("user initializer returns vectors of non-uniform length");
- for (k = 0; k < nvar; k++) xpp[j*nvar+k] = xp[k];
- UNPROTECT(1);
+ p = REAL(params);
+ memcpy(pp,p,npar*sizeof(double)); // copy the parameters
+ PROTECT(x1 = eval(fcall,rho)); nprotect++; // do the call
+ PROTECT(Snames = GET_NAMES(x1)); nprotect++;
+
+ if (!IS_NUMERIC(x1) || isNull(Snames)) {
+ UNPROTECT(nprotect);
+ error("'init.state' error: user 'initializer' must return a named numeric vector");
+ }
+
+ nvar = LENGTH(x1);
+ xp = REAL(x1);
+ midx = INTEGER(PROTECT(match(Pnames,Snames,0))); nprotect++;
+
+ for (j = 0; j < nvar; j++) {
+ if (midx[j]!=0) {
+ UNPROTECT(nprotect);
+ error("a state variable and a parameter share a single name: '%s'",CHARACTER_DATA(STRING_ELT(Snames,j)));
+ }
+ }
+
+ xdim[0] = nvar; xdim[1] = nrep;
+ PROTECT(x = makearray(2,xdim)); nprotect++;
+ setrownames(x,Snames,2);
+ fixdimnames(x,dimnms,2);
+ xt = REAL(x);
+
+ memcpy(xt,xp,nvar*sizeof(double));
+
+ for (j = 1, p += npar, xt += nvar; j < nrep; j++, p += npar, xt += nvar) {
+ memcpy(pp,p,npar*sizeof(double));
+ PROTECT(x2 = eval(fcall,rho));
+ xp = REAL(x2);
+ if (LENGTH(x2)!=nvar)
+ error("user initializer returns vectors of non-uniform length");
+ memcpy(xt,xp,nvar*sizeof(double));
+ UNPROTECT(1);
+ }
+
+ }
+
+ break;
+
+ case native: // use native routine
+
+ {
+
+ SEXP Cnames;
+ int *sidx, *pidx, *cidx;
+ double *xt, *ps, time;
+ pomp_initializer *ff = NULL;
+ int j;
+
+ PROTECT(Snames = GET_SLOT(pompfun,install("statenames"))); nprotect++;
+ PROTECT(Cnames = GET_SLOT(pompfun,install("covarnames"))); nprotect++;
+
+ // construct state, parameter, covariate, observable indices
+ sidx = INTEGER(PROTECT(name_index(Snames,pompfun,"statenames"))); nprotect++;
+ pidx = INTEGER(PROTECT(name_index(Pnames,pompfun,"paramnames"))); nprotect++;
+ cidx = INTEGER(PROTECT(name_index(Cnames,pompfun,"covarnames"))); nprotect++;
+
+ // address of native routine
+ ff = (pomp_initializer *) R_ExternalPtrAddr(fn);
+
+ nvar = LENGTH(Snames);
+ xdim[0] = nvar; xdim[1] = nrep;
+ PROTECT(x = makearray(2,xdim)); nprotect++;
+ setrownames(x,Snames,2);
+ fixdimnames(x,dimnms,2);
+
+ set_pomp_userdata(fcall);
+ GetRNGstate();
+
+ time = *(REAL(t0));
+
+ // loop over replicates
+ for (j = 0, xt = REAL(x), ps = REAL(params); j < nrep; j++, xt += nvar, ps += npar)
+ (*ff)(xt,ps,time,sidx,pidx,cidx,cp);
+
+ PutRNGstate();
+ unset_pomp_userdata();
+
+ }
+
+ break;
+
+ default:
+
+ error("unrecognized 'mode' slot in 'init.state'");
+ break;
+
+ }
+
}
UNPROTECT(nprotect);
Modified: pkg/pomp/src/pomp_internal.h
===================================================================
--- pkg/pomp/src/pomp_internal.h 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/src/pomp_internal.h 2015-06-11 11:56:30 UTC (rev 1200)
@@ -50,7 +50,7 @@
SEXP systematic_resampling(SEXP weights);
// initstate.c
-SEXP do_init_state (SEXP object, SEXP params, SEXP t0);
+SEXP do_init_state (SEXP object, SEXP params, SEXP t0, SEXP gnsi);
// rprocess.c
SEXP do_rprocess (SEXP object, SEXP xstart, SEXP times, SEXP params, SEXP offset, SEXP gnsi);
Modified: pkg/pomp/src/simulate.c
===================================================================
--- pkg/pomp/src/simulate.c 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/src/simulate.c 2015-06-11 11:56:30 UTC (rev 1200)
@@ -41,7 +41,7 @@
nreps = nsims*nparsets;
// initialize the simulations
- PROTECT(xstart = do_init_state(object,params,t0)); nprotect++;
+ PROTECT(xstart = do_init_state(object,params,t0,gnsi)); nprotect++;
PROTECT(statenames = GET_ROWNAMES(GET_DIMNAMES(xstart))); nprotect++;
dim = INTEGER(GET_DIM(xstart));
nvars = dim[0];
Modified: pkg/pomp/src/sir.c
===================================================================
--- pkg/pomp/src/sir.c 2015-06-11 11:56:11 UTC (rev 1199)
+++ pkg/pomp/src/sir.c 2015-06-11 11:56:30 UTC (rev 1200)
@@ -70,6 +70,18 @@
from_log_barycentric(pt+parindex[7],&S0,3);
}
+void _sir_init (double *x, const double *p, double t,
+ const int *stateindex, const int *parindex, const int *covindex,
+ const double *covars) {
+ double m;
+ m = POPSIZE/(S0+I0+R0);
+ SUSC = nearbyint(m*S0);
+ INFD = nearbyint(m*I0);
+ RCVD = nearbyint(m*R0);
+ CASE = 0;
+ W = 0;
+}
+
void _sir_binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t) {
More information about the pomp-commits
mailing list