[Pomp-commits] r372 - in pkg: . R inst inst/doc man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 7 15:05:16 CEST 2010
Author: kingaa
Date: 2010-10-07 15:05:16 +0200 (Thu, 07 Oct 2010)
New Revision: 372
Added:
pkg/src/simulate.c
pkg/tests/fhn.R
pkg/tests/fhn.Rout.save
Modified:
pkg/NAMESPACE
pkg/R/pomp-methods.R
pkg/R/rprocess-pomp.R
pkg/R/simulate-pomp.R
pkg/inst/ChangeLog
pkg/inst/NEWS
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/ou2-first-mif.rda
pkg/inst/doc/ou2-trajmatch.rda
pkg/man/pomp-methods.Rd
pkg/man/rprocess-pomp.Rd
pkg/src/pomp_internal.h
pkg/src/rprocess.c
pkg/tests/gillespie.Rout.save
pkg/tests/logistic.Rout.save
pkg/tests/ou2-bsmc.Rout.save
pkg/tests/ou2-forecast.Rout.save
pkg/tests/ou2-kalman.Rout.save
pkg/tests/ou2-mif.Rout.save
pkg/tests/ou2-nlf.Rout.save
pkg/tests/ou2-pmcmc.Rout.save
pkg/tests/ou2-probe.Rout.save
pkg/tests/ou2-simulate.Rout.save
pkg/tests/ou2-trajmatch.Rout.save
pkg/tests/pfilter.Rout.save
pkg/tests/ricker-probe.Rout.save
pkg/tests/ricker.Rout.save
pkg/tests/rw2.Rout.save
pkg/tests/sir.Rout.save
pkg/tests/skeleton.Rout.save
Log:
- add 'offset' argument to 'rprocess'
- implement 'simulate' in C for speed
- add 'obs' method as synonym for 'data.array'
- add FitzHugh-Nagumo model tests
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/NAMESPACE 2010-10-07 13:05:16 UTC (rev 372)
@@ -10,6 +10,7 @@
R_Euler_Multinom,
D_Euler_Multinom,
pfilter_computations,
+ simulation_computations,
apply_probe_data,
apply_probe_sim,
probe_marginal_setup,
@@ -36,9 +37,10 @@
exportMethods(
'plot','show','print','coerce',
'dprocess','rprocess','rmeasure','dmeasure','init.state','skeleton',
- 'data.array','summary','coef','logLik','time','time<-','timezero','timezero<-','window',
+ 'data.array','obs','coef','time','time<-','timezero','timezero<-',
'simulate','pfilter',
'particles','mif','continue','coef<-','states','trajectory',
+ 'summary','logLik','window',
'pred.mean','pred.var','filter.mean','conv.rec',
'bsmc',
'pmcmc','dprior',
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/R/pomp-methods.R 2010-10-07 13:05:16 UTC (rev 372)
@@ -3,6 +3,8 @@
## functions to extract or call the components of a "pomp" object
setGeneric("data.array",function(object,...)standardGeneric("data.array"))
+setGeneric("obs",function(object,...)standardGeneric("obsns"))
+
setGeneric("time<-",function(object,...,value)standardGeneric("time<-"))
setGeneric("coef<-",function(object,...,value)standardGeneric("coef<-"))
@@ -45,6 +47,20 @@
}
)
+## a simple method to extract the data array
+setMethod(
+ "obs",
+ "pomp",
+ function (object, vars, ...) {
+ varnames <- rownames(object at data)
+ if (missing(vars))
+ vars <- varnames
+ else if (!all(vars%in%varnames))
+ stop("some elements of ",sQuote("vars")," correspond to no observed variable")
+ object at data[vars,,drop=FALSE]
+ }
+ )
+
## a simple method to extract the array of states
setMethod(
"states",
Modified: pkg/R/rprocess-pomp.R
===================================================================
--- pkg/R/rprocess-pomp.R 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/R/rprocess-pomp.R 2010-10-07 13:05:16 UTC (rev 372)
@@ -6,13 +6,6 @@
setMethod(
'rprocess',
'pomp',
- function (object, xstart, times, params, ...) { # the package algorithms will only use these arguments
- x <- try(
- .Call(do_rprocess,object,xstart,times,params),
- silent=FALSE
- )
- if (inherits(x,'try-error'))
- stop("rprocess error: error in user ",sQuote("rprocess"),call.=FALSE)
- x
- }
+ function (object, xstart, times, params, offset = 0, ...)
+ .Call(do_rprocess,object,xstart,times,params,offset)
)
Modified: pkg/R/simulate-pomp.R
===================================================================
--- pkg/R/simulate-pomp.R 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/R/simulate-pomp.R 2010-10-07 13:05:16 UTC (rev 372)
@@ -13,97 +13,42 @@
call.=FALSE
)
- if (missing(times)) {
+ if (missing(times))
times <- time(object,t0=FALSE)
- } else {
+ else
times <- as.numeric(times)
- }
- if (length(times)==0)
- stop("if ",sQuote("times")," is empty, there is no work to do",call.=FALSE)
-
- if (any(diff(times)<0))
- stop(sQuote("times")," must be a nondecreasing sequence of times",call.=FALSE)
-
- if (missing(t0)) {
+ if (missing(t0))
t0 <- timezero(object)
- } else {
+ else
t0 <- as.numeric(t0)
- }
- if (t0>times[1])
- stop("the zero-time ",sQuote("t0")," must occur no later than the first observation",call.=FALSE)
- ntimes <- length(times)
+ if (missing(params))
+ params <- coef(object)
- if (missing(params)) {
- if (length(object at params)>0)
- params <- object at params
- else
- stop("no ",sQuote("params")," specified",call.=FALSE)
- }
- if (is.null(dim(params)))
- params <- matrix(params,ncol=1,dimnames=list(names(params),NULL))
- npars <- ncol(params)
+ if (length(params)==0)
+ stop("no ",sQuote("params")," specified",call.=FALSE)
+ params <- as.matrix(params)
+
if (!is.null(seed)) { # set the random seed (be very careful about this)
if (!exists('.Random.seed',envir=.GlobalEnv)) runif(1)
save.seed <- get('.Random.seed',envir=.GlobalEnv)
set.seed(seed)
}
- if (!is.numeric(nsim)||(length(nsim)>1)||nsim<1)
- stop(sQuote("nsim")," must be a positive integer")
- nsim <- as.integer(nsim)
+ retval <- try(
+ .Call(simulation_computations,object,params,times,t0,nsim,obs,states),
+ silent=FALSE
+ )
- xstart <- init.state(object,params=params,t0=t0)
- nreps <- npars*nsim # total number of replicates
- ## we will do the process model simulations with single calls to the user functions
- if (nsim > 1) { # make nsim copies of the IC and parameter matrices
- xstart <- array(rep(xstart,nsim),dim=c(nrow(xstart),nreps),dimnames=list(rownames(xstart),NULL))
- params <- array(rep(params,nsim),dim=c(nrow(params),nreps),dimnames=list(rownames(params),NULL))
- }
+ if (inherits(retval,'try-error'))
+ stop(sQuote("simulate")," error",call.=FALSE)
- x <- rprocess(object,xstart=xstart,times=c(t0,times),params=params) # simulate the process model
- x <- x[,,-1,drop=FALSE]
- ## rprocess returns a rank-3 matrix (nvars x nreps x (ntimes+1))
-
- if (obs || !states) { # we need only simulate the observation process if obs=T or states=F
- y <- rmeasure(object,x=x,times=times,params=params)
- ## rmeasure returns a rank-3 matrix (nobs x nreps x ntimes)
- }
-
if (!is.null(seed)) { # restore the RNG state
assign('.Random.seed',save.seed,envir=.GlobalEnv)
}
- if (!obs && !states) { # both obs=F and states=F, return a list of pomp objects
- nm.y <- rownames(y)
- nobs <- nrow(y)
- retval <- lapply(
- seq_len(nreps),
- function (rep) {
- po <- as(object,'pomp')
- po at t0 <- t0
- po at times <- times
- po at params <- params[,rep]
- po at data <- array(0,dim=c(nobs,ntimes),dimnames=list(nm.y,NULL))
- po at data[,] <- y[,rep,]
- po at states <- array(dim=c(nrow(x),ntimes),dimnames=list(rownames(x),NULL))
- po at states[,] <- x[,rep,]
- po
- }
- )
- if (nreps==1) retval <- retval[[1]]
- } else if (!obs && states) # return states only
- retval <- x
- else if (obs && !states) # return only observations
- retval <- y
- else # return both states and observations
- retval <- list(
- states = x,
- obs = y
- )
-
retval
}
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/inst/ChangeLog 2010-10-07 13:05:16 UTC (rev 372)
@@ -1,5 +1,11 @@
+2010-10-07 kingaa
+
+ * [r371] R/aaa.R: - remove banner
+
2010-10-06 kingaa
+ * [r370] inst/ChangeLog, inst/doc/advanced_topics_in_pomp.pdf,
+ inst/doc/intro_to_pomp.pdf:
* [r369] src/probe_ccf.c, src/probe_marginal.c, src/probe_nlar.c: -
remove unused variables
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/inst/NEWS 2010-10-07 13:05:16 UTC (rev 372)
@@ -5,10 +5,20 @@
The default for 't0' is 'timezero(object)'.
The default for 'times' is now 'time(object,t0=FALSE)' (it was 'time(object,t0=TRUE)').
This change eliminates an infelicity in which 'simulate' and 'trajectory' returned undesired values at the initial time.
+ Users can ensure that their codes will continue to function as intended by specifying the values of the 'times' and 't0' arguments to these functions, thus removing dependence on the defaults.
To aid users in modifying their codes, a warning is now displayed whenever a potentially affected call is made.
+ These warnings will be removed in a future release.
+ See the documentation ("pomp?simulate", "pomp?trajectory") for a description of the new default behaviors.
+ o The 'states' slot now holds values of the state process only at the times given in the 'times' slot.
+ This is a change from the earlier behavior in which the 'states' slot sometimes included the value of the state process at 't0'.
+
+ o 'simulate' has been re-implemented in C for speed.
+
+ o A new method 'obs' for pomp objects returns the contents of the 'data' slot.
+ It is synonymous with 'data.array'.
+
o A banner is now displayed on package attachment.
- This will be used to report important changes.
o 'coef<-' now behaves like ordinary vector assignment.
Modified: pkg/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-10-07 13:05:16 UTC (rev 372)
@@ -64,12 +64,101 @@
This document serves to give some examples of the use of native (C or FORTRAN) codes in \code{pomp}
and to introduce the low-level interface to \code{pomp} objects.
-\section{Acceleration using native codes: writing \code{rprocess} and \code{dprocess} from scratch.}
+\section{Acceleration using native codes: using plug-ins with native code}
Since many of the methods we will use require us to simulate the process and/or measurement models many times, it is a good idea to use native (compiled) codes for the computational heavy lifting.
This can result in many-fold speedup.
-The \code{pomp} package includes some examples that use C codes.
-First we'll have a look at how the discrete-time bivariate AR(1) process with normal measurement error is implemented.
+\code{pomp} provides ``plug-in'' facilities to make it easier to define certain kinds of models.
+These plug-ins can be used with native codes as well, as we'll see in the following examples.
+The \code{pomp} package includes other examples that use C codes: these can be loaded using the \code{data} command.
+
+In the ``intro\_to\_pomp'' vignette, we looked at the SIR model, which we implemented using an Euler-multinomial approximation to the continuous-time Markov process.
+Here is the same model implemented using native C codes:
+<<sir-def,keep.source=T>>=
+pomp(
+ data=data.frame(
+ time=seq(from=1/52,to=4,by=1/52),
+ reports=NA
+ ),
+ times="time",
+ t0=0,
+ ## native routine for the process simulator:
+ rprocess=euler.sim(
+ step.fun="_sir_euler_simulator",
+ delta.t=1/52/20
+ ),
+ ## native routine for the skeleton:
+ skeleton.vectorfield="_sir_ODE",
+ ## binomial measurement model:
+ rmeasure="_sir_binom_rmeasure",
+ dmeasure="_sir_binom_dmeasure",
+ ## name of the shared-object library containing the native routines:
+ PACKAGE="pomp",
+ ## the order of the observable assumed in the native routines:
+ obsnames=c("reports"),
+ ## the order of the state variables assumed in the native routines:
+ statenames=c("S","I","R","cases","W"),
+ ## the order of the parameters assumed in the native routines:
+ paramnames=c(
+ "gamma","mu","iota","beta1","beta.sd",
+ "pop","rho","nbasis","degree","period"
+ ),
+ ## reset cases to zero at each new observation:
+ zeronames=c("cases"),
+ initializer=function(params,t0,...){
+ p <- exp(params)
+ with(
+ as.list(p),
+ {
+ fracs <- c(S.0,I.0,R.0)
+ x0 <- round(c(pop*fracs/sum(fracs),0,0))
+ names(x0) <- c("S","I","R","cases","W")
+ x0
+ }
+ )
+ }
+ ) -> sir
+@
+The source code for the native routines \verb+_sir_euler_simulator+, \verb+_sir_ODE+, \verb+_sir_binom_rmeasure+, and \verb+_sir_binom_dmeasure+ is provided with the package (in the \code{examples} directory).
+To see the source code, do
+<<view-sir-source,eval=F>>=
+edit(file=system.file("examples/sir.c",package="pomp"))
+@
+Also in the \code{examples} directory is an \R\ script that shows how to compile \verb+sir.c+ into a shared-object library and link it with \R.
+Note that the native routines for this model are included in the package, which is why we give the \verb+PACKAGE="pomp"+ argument to \pomp.
+When you write your own model using native routines, you'll compile them into a dynamically-loadable library.
+In this case, you'll want to specify the name of that library using the \code{PACKAGE} argument.
+Again, refer to the SIR example included in the \code{examples} directory to see how this is done.
+
+Let's specify some parameters and simulate:
+<<sir-sim>>=
+params <- c(
+ gamma=26,mu=0.02,iota=0.01,
+ beta1=1200,beta2=1800,beta3=600,
+ beta.sd=1e-3,
+ pop=2.1e6,
+ rho=0.6,
+ S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
+ )
+
+sir <- simulate(sir,params=c(log(params),nbasis=3,degree=3,period=1),seed=3493885L)
+
+tic <- Sys.time()
+sims <- simulate(sir,nsim=10)
+toc <- Sys.time()
+print(toc-tic)
+
+tic <- Sys.time()
+traj <- trajectory(sir,hmax=1/52)
+toc <- Sys.time()
+print(toc-tic)
+@
+
+\section{Acceleration using native codes: writing \code{rprocess} and \code{dprocess} from scratch.}
+
+In the preceding example, we used ``plug-ins'' provided by the package, in conjunction with native C routines, to specify our model.
+One can also write simulators and density functions ``from scratch''.
+Here, we'll have a look at how the discrete-time bivariate AR(1) process with normal measurement error is implemented.
You can load a \code{pomp} object for this model and have a look at its structure with the commands
<<eval=F>>=
require(pomp)
@@ -190,94 +279,6 @@
print(toc-tic)
@
-\section{Acceleration using native codes: using plug-ins with native code}
-
-In the preceding example, we've written our simulators and density functions ``from scratch''.
-\code{pomp} provides ``plug-in'' facilities to make it easier to define certain kinds of models.
-These plug-ins can be used with native codes as well, as we'll see in the next example.
-
-In the ``intro\_to\_pomp'' vignette, we looked at the SIR model, which we implemented using an Euler-multinomial approximation to the continuous-time Markov process.
-Here is the same model implemented using native C codes:
-<<sir-def,keep.source=T>>=
-pomp(
- data=data.frame(
- time=seq(from=1/52,to=4,by=1/52),
- reports=NA
- ),
- times="time",
- t0=0,
- ## native routine for the process simulator:
- rprocess=euler.sim(
- step.fun="_sir_euler_simulator",
- delta.t=1/52/20
- ),
- ## native routine for the skeleton:
- skeleton.vectorfield="_sir_ODE",
- ## binomial measurement model:
- rmeasure="_sir_binom_rmeasure",
- dmeasure="_sir_binom_dmeasure",
- ## name of the shared-object library containing the native routines:
- PACKAGE="pomp",
- ## the order of the observable assumed in the native routines:
- obsnames=c("reports"),
- ## the order of the state variables assumed in the native routines:
- statenames=c("S","I","R","cases","W"),
- ## the order of the parameters assumed in the native routines:
- paramnames=c(
- "gamma","mu","iota","beta1","beta.sd",
- "pop","rho","nbasis","degree","period"
- ),
- ## reset cases to zero at each new observation:
- zeronames=c("cases"),
- initializer=function(params,t0,...){
- p <- exp(params)
- with(
- as.list(p),
- {
- fracs <- c(S.0,I.0,R.0)
- x0 <- round(c(pop*fracs/sum(fracs),0,0))
- names(x0) <- c("S","I","R","cases","W")
- x0
- }
- )
- }
- ) -> sir
-@
-The source code for the native routines \verb+_sir_euler_simulator+, \verb+_sir_ODE+, \verb+_sir_binom_rmeasure+, and \verb+_sir_binom_dmeasure+ is provided with the package (in the \code{examples} directory).
-To see the source code, do
-<<view-sir-source,eval=F>>=
-edit(file=system.file("examples/sir.c",package="pomp"))
-@
-Also in the \code{examples} directory is an \R\ script that shows how to compile \verb+sir.c+ into a shared-object library and link it with \R.
-Note that the native routines for this model are included in the package, which is why we give the \verb+PACKAGE="pomp"+ argument to \pomp.
-When you write your own model using native routines, you'll compile them into a dynamically-loadable library.
-In this case, you'll want to specify the name of that library using the \code{PACKAGE} argument.
-Again, refer to the SIR example included in the \code{examples} directory to see how this is done.
-
-Let's specify some parameters and simulate:
-<<sir-sim>>=
-params <- c(
- gamma=26,mu=0.02,iota=0.01,
- beta1=1200,beta2=1800,beta3=600,
- beta.sd=1e-3,
- pop=2.1e6,
- rho=0.6,
- S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
- )
-
-sir <- simulate(sir,params=c(log(params),nbasis=3,degree=3,period=1),seed=3493885L)
-
-tic <- Sys.time()
-sims <- simulate(sir,nsim=10)
-toc <- Sys.time()
-print(toc-tic)
-
-tic <- Sys.time()
-traj <- trajectory(sir,hmax=1/52)
-toc <- Sys.time()
-print(toc-tic)
-@
-
\section{The low-level interface}
There is a low-level interface to \code{pomp} objects, primarily designed for package developers.
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ou2-first-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ou2-trajmatch.rda
===================================================================
(Binary files differ)
Modified: pkg/man/pomp-methods.Rd
===================================================================
--- pkg/man/pomp-methods.Rd 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/man/pomp-methods.Rd 2010-10-07 13:05:16 UTC (rev 372)
@@ -20,6 +20,9 @@
\alias{timezero<--pomp}
\alias{window,pomp-method}
\alias{window--pomp}
+\alias{obs,pomp-method}
+\alias{obs-pomp}
+\alias{obs}
\alias{data.array,pomp-method}
\alias{data.array-pomp}
\alias{data.array}
@@ -39,6 +42,7 @@
\usage{
\S4method{coef}{pomp}(object, pars, \dots)
\S4method{coef}{pomp}(object, pars, \dots) <- value
+\S4method{obs}{pomp}(object, vars, \dots)
\S4method{data.array}{pomp}(object, vars, \dots)
\S4method{states}{pomp}(object, vars, \dots)
\S4method{time}{pomp}(x, t0 = FALSE, \dots)
@@ -134,9 +138,10 @@
If \code{coef(object)} does not exist, then \code{coef(object,pars) <- value} assigns \code{value} to the parameters of \code{object};
in this case, the names of \code{object} will be \code{pars} and the names of \code{value} will be ignored.
}
- \item{data.array}{
- \code{data.array(object)} returns the array of observations.
- \code{data.array(object,vars)} gives just the observations of variables named \code{vars}.
+ \item{obs, data.array}{
+ These functions are synonymous.
+ \code{obs(object)} returns the array of observations.
+ \code{obs(object,vars)} gives just the observations of variables named \code{vars}.
\code{vars} may specify the variables by position or by name.
}
\item{states}{
Modified: pkg/man/rprocess-pomp.Rd
===================================================================
--- pkg/man/rprocess-pomp.Rd 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/man/rprocess-pomp.Rd 2010-10-07 13:05:16 UTC (rev 372)
@@ -12,7 +12,7 @@
}
\usage{
rprocess(object, xstart, times, params, \dots)
-\S4method{rprocess}{pomp}(object, xstart, times, params, \dots)
+\S4method{rprocess}{pomp}(object, xstart, times, params, offset, \dots)
}
\arguments{
\item{object}{an object of class \code{pomp}.}
@@ -30,14 +30,16 @@
\item{params}{
a rank-2 array of parameters with the parameters corresponding to the columns of \code{xstart}.
}
+ \item{offset}{
+ throw away the first \code{offset} times.
+ }
\item{\dots}{at present, these are ignored.}
}
\value{
\code{rprocess} returns a rank-3 array with rownames.
Suppose \code{x} is the array returned.
- Then \code{dim(x)=c(nvars,nreps,ntimes)}, where \code{nvars} is the number of state variables (=\code{nrow(xstart)}), \code{nreps} is the number of independent realizations simulated (=\code{ncol(xstart)}), and \code{ntimes} is the length of the vector \code{times}.
- \code{x[,j,k]} is the value of the state process in the \code{j}-th realization at time \code{times[k]}.
- In particular, \code{x[,,1]} must be identical to \code{xstart}.
+ Then \code{dim(x)=c(nvars,nreps,ntimes-offset)}, where \code{nvars} is the number of state variables (=\code{nrow(xstart)}), \code{nreps} is the number of independent realizations simulated (=\code{ncol(xstart)}), and \code{ntimes} is the length of the vector \code{times}.
+ \code{x[,j,k]} is the value of the state process in the \code{j}-th realization at time \code{times[k+offset]}.
The rownames of \code{x} must correspond to those of \code{xstart}.
}
\details{
Modified: pkg/src/pomp_internal.h
===================================================================
--- pkg/src/pomp_internal.h 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/src/pomp_internal.h 2010-10-07 13:05:16 UTC (rev 372)
@@ -26,22 +26,32 @@
// setting dydt = 0 in the call to 'table_lookup' will bypass computation of the derivative
void table_lookup (struct lookup_table *tab, double x, double *y, double *dydt);
-/* bspline.c */
+// bspline.c
SEXP bspline_basis(SEXP x, SEXP degree, SEXP knots);
SEXP bspline_basis_function(SEXP x, SEXP i, SEXP degree, SEXP knots);
-/* dsobol.c */
+// dsobol.c
SEXP sobol_sequence(SEXP dim);
-/* pomp_fun.c */
+// pomp_fun.c
SEXP pomp_fun_handler (SEXP pfun, int *use_native);
-/* lookup_table.c */
+// lookup_table.c
SEXP lookup_in_table (SEXP ttable, SEXP xtable, SEXP t, int *index);
-/* resample.c */
+// resample.c
SEXP systematic_resampling(SEXP weights);
+// initstate.c
+SEXP do_init_state (SEXP object, SEXP params, SEXP t0);
+
+// rprocess.c
+SEXP do_rprocess (SEXP object, SEXP xstart, SEXP times, SEXP params, SEXP offset);
+
+// rmeasure.c
+SEXP do_rmeasure (SEXP object, SEXP x, SEXP times, SEXP params);
+
+
static R_INLINE SEXP makearray (int rank, int *dim) {
int nprotect = 0;
int *dimp, k;
Modified: pkg/src/rprocess.c
===================================================================
--- pkg/src/rprocess.c 2010-10-07 10:41:01 UTC (rev 371)
+++ pkg/src/rprocess.c 2010-10-07 13:05:16 UTC (rev 372)
@@ -5,31 +5,32 @@
#include <Rdefines.h>
#include <Rinternals.h>
-SEXP do_rprocess (SEXP object, SEXP xstart, SEXP times, SEXP params)
+#include "pomp_internal.h"
+
+SEXP do_rprocess (SEXP object, SEXP xstart, SEXP times, SEXP params, SEXP offset)
{
int nprotect = 0;
- int *xdim, nvars, nreps, ntimes;
- SEXP X, fn, fcall, rho;
+ int *xdim, nvars, nreps, ntimes, off;
+ SEXP X, Xoff, fn, fcall, rho;
SEXP dimXstart, dimP, dimX;
ntimes = length(times);
if (ntimes < 2) {
- UNPROTECT(nprotect);
error("rprocess error: no transitions, no work to do");
}
+ off = *(INTEGER(AS_INTEGER(offset)));
+ if ((off < 0)||(off>=ntimes))
+ error("illegal 'offset' value %d",off);
PROTECT(dimXstart = GET_DIM(xstart)); nprotect++;
if ((isNull(dimXstart)) || (length(dimXstart)!=2)) {
- UNPROTECT(nprotect);
error("rprocess error: 'xstart' must be a rank-2 array");
}
PROTECT(dimP = GET_DIM(params)); nprotect++;
if ((isNull(dimP)) || (length(dimP)!=2)) {
- UNPROTECT(nprotect);
error("rprocess error: 'params' must be a rank-2 array");
}
xdim = INTEGER(dimXstart);
nvars = xdim[0]; nreps = xdim[1];
if (nreps != INTEGER(dimP)[1]) {
- UNPROTECT(nprotect);
error("rprocess error: number of columns of 'params' and 'xstart' do not agree");
}
// extract the process function
@@ -59,18 +60,28 @@
PROTECT(X = eval(fcall,rho)); nprotect++; // do the call
PROTECT(dimX = GET_DIM(X)); nprotect++;
if ((isNull(dimX)) || (length(dimX) != 3)) {
- UNPROTECT(nprotect);
error("rprocess error: user 'rprocess' must return a rank-3 array");
}
xdim = INTEGER(dimX);
if ((xdim[0] != nvars) || (xdim[1] != nreps) || (xdim[2] != ntimes)) {
- UNPROTECT(nprotect);
error("rprocess error: user 'rprocess' must return a %d x %d x %d array",nvars,nreps,ntimes);
}
if (isNull(GET_ROWNAMES(GET_DIMNAMES(X)))) {
- UNPROTECT(nprotect);
error("rprocess error: user 'rprocess' must return an array with rownames");
}
- UNPROTECT(nprotect);
- return X;
+ if (off > 0) {
+ int i, n;
+ double *s, *t;
+ xdim[2] -= off;
+ PROTECT(Xoff = makearray(3,xdim)); nprotect++;
+ setrownames(Xoff,GET_ROWNAMES(GET_DIMNAMES(X)),3);
+ s = REAL(X)+off*nvars*nreps;
+ t = REAL(Xoff);
+ for (i = 0, n = nvars*nreps*(ntimes-off); i < n; i++, s++, t++) *t = *s;
+ UNPROTECT(nprotect);
+ return Xoff;
+ } else {
+ UNPROTECT(nprotect);
+ return X;
+ }
}
Added: pkg/src/simulate.c
===================================================================
--- pkg/src/simulate.c (rev 0)
+++ pkg/src/simulate.c 2010-10-07 13:05:16 UTC (rev 372)
@@ -0,0 +1,217 @@
+// -*- C++ -*-
+
+#include "pomp_internal.h"
+#include <Rdefines.h>
+
+SEXP simulation_computations (SEXP object, SEXP params, SEXP times, SEXP t0, SEXP nsim, SEXP obs, SEXP states)
+{
+ int nprotect = 0;
+ SEXP xstart, x, xx, x0, y, p, alltimes, coef, yy, offset;
+ SEXP ans, ans_names;
+ SEXP po, popo;
+ SEXP statenames, paramnames, obsnames, statedim, obsdim;
+ int nsims, nparsets, nreps, npars, nvars, ntimes, nobs;
+ int qobs, qstates;
+ int *dim, dims[3];
+ double *s, *t, *xs, *xt, *ys, *yt, *ps, *pt, tt;
+ int i, j, k;
+
+ PROTECT(offset = NEW_INTEGER(1)); nprotect++;
+ INTEGER(offset)[0] = 1;
+
+ nsims = INTEGER(AS_INTEGER(nsim))[0]; // number of simulations per parameter set
+ if (LENGTH(nsim)>1)
+ warning("only the first number in 'nsim' is significant");
+ if (nsims < 1)
+ return R_NilValue; // no work to do
+
+ qobs = LOGICAL(AS_LOGICAL(obs))[0]; // 'obs' flag set?
+ qstates = LOGICAL(AS_LOGICAL(states))[0]; // 'states' flag set?
+
+ PROTECT(paramnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
+ dim = INTEGER(GET_DIM(params));
+ npars = dim[0]; nparsets = dim[1];
+
+ nreps = nsims*nparsets;
+
+ // initialize the simulations
+ PROTECT(xstart = do_init_state(object,params,t0)); nprotect++;
+ PROTECT(statenames = GET_ROWNAMES(GET_DIMNAMES(xstart))); nprotect++;
+ dim = INTEGER(GET_DIM(xstart));
+ nvars = dim[0];
+
+ // augment the 'times' vector with 't0'
+ ntimes = LENGTH(times);
+
+ if (ntimes < 1)
+ error("if 'times' is empty, there is no work to do");
+
+ PROTECT(alltimes = NEW_NUMERIC(ntimes+1)); nprotect++;
+ tt = REAL(t0)[0];
+ s = REAL(times);
+ t = REAL(alltimes);
+
+ if (tt > *s)
+ error("the zero-time 't0' must occur no later than the first observation 'times[1]'");
+
+ *(t++) = tt; // copy t0 into alltimes[1]
+ tt = *(t++) = *(s++); // copy times[1] into alltimes[2]
+
+ for (j = 1; j < ntimes; j++) { // copy times[2:ntimes] into alltimes[3:(ntimes+1)]
+ if (tt >= *s)
+ error("'times' must be a nondecreasing sequence of times");
+ tt = *(t++) = *(s++);
+ }
+
+ // call 'rprocess' to simulate state process
+ if (nsims == 1) { // only one simulation per parameter set
+
+ PROTECT(x = do_rprocess(object,xstart,alltimes,params,offset)); nprotect++;
+
+ } else { // nsim > 1
+
+ dims[0] = npars; dims[1] = nreps;
+ PROTECT(p = makearray(2,dims)); nprotect++;
+ setrownames(p,paramnames,2);
+
+ dims[0] = nvars;
+ PROTECT(x0 = makearray(2,dims)); nprotect++;
+ setrownames(x0,statenames,2);
+
+ // make 'nsims' copies of the parameters and initial states
+ for (k = 0, pt = REAL(p), xt = REAL(x0); k < nsims; k++) {
+ for (j = 0, ps = REAL(params), xs = REAL(xstart); j < nparsets; j++) {
+ for (i = 0; i < npars; i++, pt++, ps++) *pt = *ps;
+ for (i = 0; i < nvars; i++, xt++, xs++) *xt = *xs;
+ }
+ }
+
+ PROTECT(x = do_rprocess(object,x0,alltimes,p,offset)); nprotect++;
+
+ }
+
+ if (!qobs && qstates) { // obs=F,states=T: return states only
+
+ UNPROTECT(nprotect);
+ return x;
+
+ } else {
+
+ if (nsims == 1) {
+ PROTECT(y = do_rmeasure(object,x,times,params)); nprotect++;
+ } else {
+ PROTECT(y = do_rmeasure(object,x,times,p)); nprotect++;
+ }
+
+ if (qobs) {
+
+ if (qstates) { // obs=T,states=T: return a list with states and obs'ns
+
+ PROTECT(ans = NEW_LIST(2)); nprotect++;
+ PROTECT(ans_names = NEW_CHARACTER(2)); nprotect++;
+ SET_NAMES(ans,ans_names);
+ SET_STRING_ELT(ans_names,0,mkChar("states"));
+ SET_STRING_ELT(ans_names,1,mkChar("obs"));
+ SET_ELEMENT(ans,0,x);
+ SET_ELEMENT(ans,1,y);
+ UNPROTECT(nprotect);
+ return ans;
+
+ } else { // obs=T,states=F: return observations only
+
+ UNPROTECT(nprotect);
+ return y;
+
+ }
+
+ } else { // obs=F,states=F: return one or more pomp objects
+
+ PROTECT(obsnames = GET_ROWNAMES(GET_DIMNAMES(y))); nprotect++;
+ nobs = INTEGER(GET_DIM(y))[0];
+
+ PROTECT(obsdim = NEW_INTEGER(2)); nprotect++;
+ INTEGER(obsdim)[0] = nobs;
+ INTEGER(obsdim)[1] = ntimes;
+
+ PROTECT(statedim = NEW_INTEGER(2)); nprotect++;
+ INTEGER(statedim)[0] = nvars;
+ INTEGER(statedim)[1] = ntimes;
+
+ PROTECT(coef = NEW_NUMERIC(npars)); nprotect++;
+ SET_NAMES(coef,paramnames);
+
+ PROTECT(po = duplicate(object)); nprotect++;
+ SET_SLOT(po,install("t0"),t0);
+ SET_SLOT(po,install("times"),times);
+ SET_SLOT(po,install("params"),coef);
+
+ if (nreps == 1) {
+
+ SET_DIM(y,obsdim);
+ setrownames(y,obsnames,2);
+ SET_SLOT(po,install("data"),y);
+
+ SET_DIM(x,statedim);
+ setrownames(x,statenames,2);
+ SET_SLOT(po,install("states"),x);
+
+ ps = REAL(params);
+ pt = REAL(GET_SLOT(po,install("params")));
+ for (i = 0; i < npars; i++, ps++, pt++) *pt = *ps;
+
+ UNPROTECT(nprotect);
+ return po;
+
+ } else {
+
+ // a list to hold the pomp objects
+ PROTECT(ans = NEW_LIST(nreps)); nprotect++;
+
+ // create an array for the 'states' slot
+ PROTECT(xx = makearray(2,INTEGER(statedim))); nprotect++;
+ setrownames(xx,statenames,2);
+ SET_SLOT(po,install("states"),xx);
+
+ // create an array for the 'data' slot
+ PROTECT(yy = makearray(2,INTEGER(obsdim))); nprotect++;
+ setrownames(yy,obsnames,2);
+ SET_SLOT(po,install("data"),yy);
+
+ xs = REAL(x);
+ ys = REAL(y);
+ if (nsims == 1)
+ ps = REAL(params);
+ else
+ ps = REAL(p);
+
+ for (k = 0; k < nreps; k++) { // loop over replicates
+
+ PROTECT(popo = duplicate(po));
+
+ // copy x[,k,] and y[,k,] into popo
+ xt = REAL(GET_SLOT(popo,install("states")));
+ yt = REAL(GET_SLOT(popo,install("data")));
+ pt = REAL(GET_SLOT(popo,install("params")));
+
+ for (i = 0; i < npars; i++, pt++, ps++) *pt = *ps;
+ for (j = 0; j < ntimes; j++) {
+ for (i = 0; i < nvars; i++, xt++) *xt = xs[i+nvars*(k+nreps*j)];
+ for (i = 0; i < nobs; i++, yt++) *yt = ys[i+nobs*(k+nreps*j)];
+ }
+
+ SET_ELEMENT(ans,k,popo);
+ UNPROTECT(1);
+
+ }
+
+ UNPROTECT(nprotect);
+ return ans;
+
+ }
+
+ }
+
+ }
+
+}
+
Property changes on: pkg/src/simulate.c
___________________________________________________________________
Added: svn:eol-style
+ native
Added: pkg/tests/fhn.R
===================================================================
--- pkg/tests/fhn.R (rev 0)
+++ pkg/tests/fhn.R 2010-10-07 13:05:16 UTC (rev 372)
@@ -0,0 +1,55 @@
+library(pomp)
+
+pdf.options(useDingbats=FALSE)
+pdf(file="fhn.pdf")
+
+## autonomous case
+fhn <- pomp(
+ data=data.frame(time=seq(0,60,by=0.1),Vobs=NA),
+ times="time",
+ t0=-0.1,
+ skeleton.vectorfield=function(x,t,params,...) {
+ with(
+ as.list(c(x,params)),
+ c(
+ V=c*(V-V^3/3-R+i),
+ R=(V+a-b*R)/c
+ )
+ )
+ }
+ )
+
+x <- array(c(0,1,1,2,1,1,0,-1),dim=c(2,2,2),dimnames=list(c("V","R"),NULL,NULL))
+params <- rbind(a=c(0.7,0.5),b=c(0.8,0.5),c=c(2,5),V.0=c(1,1),R.0=c(0,0),i=c(0.8,0))
+skeleton(fhn,x,t=c(0,3),params=params)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 372
More information about the pomp-commits
mailing list