[Pomp-commits] r38 - in pkg: . R data inst/examples man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 13 18:12:12 CEST 2008
Author: kingaa
Date: 2008-08-13 18:12:12 +0200 (Wed, 13 Aug 2008)
New Revision: 38
Added:
pkg/R/trajectory-pomp.R
pkg/man/trajectory-pomp.Rd
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/aaa.R
pkg/R/pomp.R
pkg/data/ou2.rda
pkg/inst/examples/logistic.R
pkg/inst/examples/sir.R
pkg/inst/examples/sir.c
pkg/man/mif-class.Rd
pkg/man/pomp-class.Rd
pkg/man/pomp.Rd
pkg/man/skeleton-pomp.Rd
pkg/src/skeleton.c
pkg/tests/logistic.R
pkg/tests/logistic.Rout.save
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
add support for computing trajectories of the deterministic skeleton
the skeleton must now be specified when the pomp is created: if the pomp is a discrete-time system, 'skeleton.map' is used to specify its deterministic skeleton; if it is a continuous-time system, 'skeleton.vectorfield' is used.
the new method 'trajectory' will, in the first case, iterate the map to produce an itinerary and in the second case, it will call 'lsoda' from the now-required 'odesolve' package to compute the trajectory.
some changes to the basic 'pomp' structure have been made to facilitate all this.
specifically, a new slot, 'skeleton-type', distinguishes between discrete-time and continuous-time systems
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/DESCRIPTION 2008-08-13 16:12:12 UTC (rev 38)
@@ -1,11 +1,11 @@
Package: pomp
Type: Package
Title: Partially-observed Markov processes
-Version: 0.20-8
-Date: 2008-07-31
+Version: 0.21-1
+Date: 2008-08-13
Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
-Depends: R(>= 2.6.1), stats, methods, graphics
+Depends: R(>= 2.6.1), stats, methods, graphics, odesolve
License: GPL (>= 2)
LazyLoad: true
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/NAMESPACE 2008-08-13 16:12:12 UTC (rev 38)
@@ -26,7 +26,7 @@
'simulate','data.array','pfilter',
'coef','logLik',
'pred.mean','pred.var','filter.mean','conv.rec',
- 'particles','mif','continue','coef<-','states'
+ 'particles','mif','continue','coef<-','states','trajectory'
)
export(
Modified: pkg/R/aaa.R
===================================================================
--- pkg/R/aaa.R 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/R/aaa.R 2008-08-13 16:12:12 UTC (rev 38)
@@ -23,6 +23,7 @@
dprocess = 'function',
dmeasure = 'pomp.fun',
rmeasure = 'pomp.fun',
+ skeleton.type = 'character',
skeleton = 'pomp.fun',
initializer = 'function',
states = 'array',
@@ -128,3 +129,7 @@
states <- function (object, ...)
stop("function ",sQuote("states")," is undefined for objects of class ",sQuote(class(object)))
setGeneric('states')
+
+trajectory <- function (object, params, times, ...)
+ stop("function ",sQuote("trajectory")," is undefined for objects of class ",sQuote(class(object)))
+setGeneric('trajectory')
Modified: pkg/R/pomp.R
===================================================================
--- pkg/R/pomp.R 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/R/pomp.R 2008-08-13 16:12:12 UTC (rev 38)
@@ -1,7 +1,7 @@
## constructor of the pomp class
pomp <- function (data, times, t0, ..., rprocess, dprocess,
rmeasure, dmeasure, measurement.model,
- skeleton, initializer, covar, tcovar,
+ skeleton.map, skeleton.vectorfield, initializer, covar, tcovar,
statenames, paramnames, covarnames,
PACKAGE) {
## check the data
@@ -50,9 +50,68 @@
rmeasure <- function(x,t,params,covars,...)stop(sQuote("rmeasure")," not specified")
if (missing(dmeasure))
dmeasure <- function(y,x,t,params,log=FALSE,covars,...)stop(sQuote("dmeasure")," not specified")
- if (missing(skeleton))
- skeleton <- function(x,t,params,covars,...)stop(sQuote("skeleton")," not specified")
-
+
+ if (missing(skeleton.map)) {
+ if (missing(skeleton.vectorfield)) {# skeleton is unspecified
+ skeleton.type <- as.character(NA)
+ skeleton <- new(
+ "pomp.fun",
+ R.fun=function(x,t,params,covars,...)stop(sQuote("skeleton")," not specified"),
+ use=as.integer(1)
+ )
+ } else { # skeleton is a vectorfield (ordinary differential equation)
+ skeleton.type <- "vectorfield"
+ if (is.function(skeleton.vectorfield)) {
+ if (!all(c('x','t','params','...')%in%names(formals(skeleton.vectorfield))))
+ stop(
+ sQuote("skeleton.vectorfield"),
+ " must be a function of prototype ",
+ sQuote("skeleton.vectorfield(x,t,params,...)")
+ )
+ skeleton <- new(
+ "pomp.fun",
+ R.fun=skeleton.vectorfield,
+ use=as.integer(1)
+ )
+ } else if (is.character(skeleton.vectorfield)) {
+ skeleton <- new(
+ "pomp.fun",
+ R.fun=function(x,t,params,...)stop("very bad: skel.fun"),
+ native.fun=skeleton.vectorfield,
+ PACKAGE=PACKAGE,
+ use=as.integer(2)
+ )
+ } else {
+ stop(sQuote("skeleton.vectorfield")," must be either a function or the name of a compiled routine")
+ }
+ }
+ } else {
+ if (missing(skeleton.vectorfield)) { # skeleton is a map (discrete-time system)
+ skeleton.type <- "map"
+ if (is.function(skeleton.map)) {
+ if (!all(c('x','t','params','...')%in%names(formals(skeleton.map))))
+ stop(sQuote("skeleton.map")," must be a function of prototype ",sQuote("skeleton.map(x,t,params,...)"))
+ skeleton <- new(
+ "pomp.fun",
+ R.fun=skeleton.map,
+ use=as.integer(1)
+ )
+ } else if (is.character(skeleton.map)) {
+ skeleton <- new(
+ "pomp.fun",
+ R.fun=function(x,t,params,...)stop("very bad: skel.fun"),
+ native.fun=skeleton.map,
+ PACKAGE=PACKAGE,
+ use=as.integer(2)
+ )
+ } else {
+ stop(sQuote("skeleton.map")," must be either a function or the name of a compiled routine")
+ }
+ } else { # a dynamical system cannot be both a map and a vectorfield
+ stop("it is not permitted to specify both ",sQuote("skeleton.map")," and ",sQuote("skeleton.vectorfield"))
+ }
+ }
+
if (missing(initializer)) {
initializer <- function (params, t0, ...) {
ivpnames <- grep("\\.0$",names(params),val=TRUE)
@@ -111,26 +170,6 @@
stop(sQuote("dmeasure")," must be either a function or the name of a compiled routine")
}
- if (is.function(skeleton)) {
- if (!all(c('x','t','params','...')%in%names(formals(skeleton))))
- stop(sQuote("skeleton")," must be a function of prototype ",sQuote("skeleton(x,t,params,...)"))
- skeleton <- new(
- "pomp.fun",
- R.fun=skeleton,
- use=as.integer(1)
- )
- } else if (is.character(skeleton)) {
- skeleton <- new(
- "pomp.fun",
- R.fun=function(x,t,params,...)stop("very bad: skel.fun"),
- native.fun=skeleton,
- PACKAGE=PACKAGE,
- use=as.integer(2)
- )
- } else {
- stop(sQuote("skeleton")," must be either a function or the name of a compiled routine")
- }
-
if (!is.function(initializer))
stop("pomp error: ",sQuote("initializer")," must be a function")
@@ -207,6 +246,7 @@
dprocess = dprocess,
dmeasure = dmeasure,
rmeasure = rmeasure,
+ skeleton.type = skeleton.type,
skeleton = skeleton,
data = data,
times = times,
Added: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R (rev 0)
+++ pkg/R/trajectory-pomp.R 2008-08-13 16:12:12 UTC (rev 38)
@@ -0,0 +1,63 @@
+setMethod(
+ "trajectory",
+ "pomp",
+ function (object, params, times, ...) {
+ params <- as.matrix(params)
+ if (missing(times))
+ times <- time(object,t0=TRUE)
+ x0 <- init.state(object,params=params,t0=times[1])
+ x <- array(
+ dim=c(nrow(x0),ncol(x0),length(times)),
+ dimnames=list(rownames(x0),NULL,NULL)
+ )
+ switch(
+ object at skeleton.type,
+ map={ # iterate the map
+ x[,,1] <- x0
+ for (k in 2:length(times)) {
+ x[,,k] <- skeleton(
+ object,
+ x=x[,,k-1,drop=FALSE],
+ t=times[k-1],
+ params=params
+ )
+ }
+ },
+ vectorfield={ # integrate the vectorfield
+ for (j in 1:ncol(params)) {
+ X <- try(
+ lsoda(
+ y=x0[,j],
+ times=times,
+ func=function(t,y,parms){
+ list(
+ skeleton(
+ object,
+ x=array(
+ data=y,
+ dim=c(length(y),1,1),
+ dimnames=list(names(y),NULL,NULL)
+ ),
+ t=t,
+ params=as.matrix(parms)
+ ),
+ NULL
+ )
+ },
+ parms=params[,j],
+ ...
+ ),
+ silent=FALSE
+ )
+ if (inherits(X,'try-error'))
+ stop("trajectory error: error in ",sQuote("lsoda"),call.=FALSE)
+ if (attr(X,'istate')!=2)
+ warning("abnormal exit from ",sQuote("lsoda"),", istate = ",attr(X,'istate'),call.=FALSE)
+ x[,j,] <- t(X[,-1])
+ }
+ },
+ unspecified=stop("deterministic skeleton not specified")
+ )
+ x
+ }
+ )
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/examples/logistic.R
===================================================================
--- pkg/inst/examples/logistic.R 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/inst/examples/logistic.R 2008-08-13 16:12:12 UTC (rev 38)
@@ -29,7 +29,7 @@
)
},
measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
- skeleton=function(x,t,params,...){
+ skeleton.vectorfield=function(x,t,params,...){
with(
as.list(c(x,params)),
r*n*(1-n/K)
@@ -41,4 +41,6 @@
params <- c(n.0=10000,K=10000,r=0.9,sigma=0.4,tau=0.1)
set.seed(73658676)
po <- simulate(po,params=params)[[1]]
-plot(po)
+
+params <- cbind(c(n.0=100,K=10000,r=0.2,sigma=0.4,tau=0.1),c(n.0=1000,K=11000,r=0.1,sigma=0.4,tau=0.1))
+x <- trajectory(po,params=params)
Modified: pkg/inst/examples/sir.R
===================================================================
--- pkg/inst/examples/sir.R 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/inst/examples/sir.R 2008-08-13 16:12:12 UTC (rev 38)
@@ -77,7 +77,7 @@
}
)
},
- skeleton=function(x,t,params,covars,...) {
+ skeleton.vectorfield=function(x,t,params,covars,...) {
params <- exp(params)
with(
as.list(c(x,params)),
@@ -92,12 +92,13 @@
mu*I,
mu*R
)
- c(
- terms[1]-terms[2]-terms[3],
- terms[2]-terms[4]-terms[5],
- terms[4]-terms[6],
- terms[4]
- )
+ xdot <- c(
+ terms[1]-terms[2]-terms[3],
+ terms[2]-terms[4]-terms[5],
+ terms[4]-terms[6],
+ terms[4]
+ )
+ ifelse(x>0,xdot,0)
}
)
},
@@ -150,7 +151,7 @@
rprocess=euler.simulate,
dens.fun="sir_euler_density",
dprocess=euler.density,
- skeleton="sir_ODE",
+ skeleton.vectorfield="sir_ODE",
PACKAGE="sir_example", ## name of the shared-object library
measurement.model=measles~binom(size=cases,prob=exp(rho)),
initializer=function(params,t0,...){
Modified: pkg/inst/examples/sir.c
===================================================================
--- pkg/inst/examples/sir.c 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/inst/examples/sir.c 2008-08-13 16:12:12 UTC (rev 38)
@@ -165,6 +165,12 @@
DRDT = term[3]-term[5];
DCDT = term[3]; // cases are cumulative recoveries
+// // trap for negative states
+// DSDT = (SUSC < 0.0) ? DSDT : 0.0;
+// DIDT = (INFD < 0.0) ? DIDT : 0.0;
+// DRDT = (RCVD < 0.0) ? DRDT : 0.0;
+// DCDT = (CASE < 0.0) ? DCDT : 0.0;
+
}
#undef DSDT
Modified: pkg/man/mif-class.Rd
===================================================================
--- pkg/man/mif-class.Rd 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/man/mif-class.Rd 2008-08-13 16:12:12 UTC (rev 38)
@@ -77,7 +77,7 @@
A numeric value containing the value of the log likelihood, as evaluated for the random-parameter model.
Note that this will not be equal to the log likelihood for the fixed-parameter model.
}
- \item{data, times, t0, rprocess, dprocess, dmeasure, rmeasure, skeleton, initializer, states, params, statenames, paramnames, covarnames, tcovar, covar, PACKAGE, userdata}{Inherited from the \code{pomp} class.}
+ \item{data, times, t0, rprocess, dprocess, dmeasure, rmeasure, skeleton.type, skeleton, initializer, states, params, statenames, paramnames, covarnames, tcovar, covar, PACKAGE, userdata}{Inherited from the \code{pomp} class.}
}
}
\section{Extends}{
Modified: pkg/man/pomp-class.Rd
===================================================================
--- pkg/man/pomp-class.Rd 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/man/pomp-class.Rd 2008-08-13 16:12:12 UTC (rev 38)
@@ -31,6 +31,9 @@
\item{rmeasure}{
an object of class "pomp.fun" which encodes the measurement model simulator.
}
+ \item{skeleton.type}{
+ a character variable specifying whether the deterministic skeleton is a map or a vectorfield.
+ }
\item{skeleton}{
an object of class "pomp.fun" which encodes the deterministic skeleton.
}
Modified: pkg/man/pomp.Rd
===================================================================
--- pkg/man/pomp.Rd 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/man/pomp.Rd 2008-08-13 16:12:12 UTC (rev 38)
@@ -6,8 +6,8 @@
}
\usage{
pomp(data, times, t0, \dots, rprocess, dprocess, rmeasure, dmeasure,
- measurement.model, skeleton, initializer, covar, tcovar,
- statenames, paramnames, covarnames, PACKAGE)
+ measurement.model, skeleton.map, skeleton.vectorfield, initializer,
+ covar, tcovar, statenames, paramnames, covarnames, PACKAGE)
}
\arguments{
\item{data}{
@@ -55,12 +55,16 @@
See below for an example.
NB: it will typically be possible to acclerate measurement model computations by writing \code{dmeasure} and/or \code{rmeasure} functions directly.
}
- \item{skeleton}{
- The deterministic skeleton can be provided in one of two ways:
+ \item{skeleton.map}{
+ If we are dealing with a discrete-time Markov process, its deterministic skeleton is a map and can be specified using \code{skeleton.map} in one of two ways:
(1) as a function of prototype \code{skeleton(x,t,params,\dots)} which evaluates the deterministic skeleton (whether vectorfield or map) of at state \code{x} and time \code{t} given the parameters \code{params}, or
(2) as the name of a native (compiled) routine with prototype "pomp\_vectorfield\_map" as defined in the header file "pomp.h".
If the deterministic skeleton depends on covariates, the optional argument \code{covars} will be filled with interpolated values of the covariates at the time \code{t}.
}
+ \item{skeleton.vectorfield}{
+ If we are dealing with a continuous-time Markov process, its deterministic skeleton is a vectorfield and can be specified using \code{skeleton.vectorfield} in either of the two ways described above, under \code{skeleton.map}.
+ Note that it makes no sense to specify the deterministic skeleton both as a map and as a vectorfield: an attempt to do so will generate an error.
+ }
\item{initializer}{
Function of prototype \code{initializer(params,t0,\dots)} which yields initial conditions for the state process when given a vector, \code{params}, of parameters.
By default, i.e., if it is unspecified when \code{pomp} is called, the initializer assumes any names in \code{params} that end in ".0" correspond to initial value parameters.
Modified: pkg/man/skeleton-pomp.Rd
===================================================================
--- pkg/man/skeleton-pomp.Rd 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/man/skeleton-pomp.Rd 2008-08-13 16:12:12 UTC (rev 38)
@@ -38,7 +38,7 @@
}
\details{
This function makes repeated calls to the user-supplied \code{skeleton} of the \code{pomp} object.
- For specifications on supplying these, see \code{\link{pomp}}.
+ For specifications on supplying this, see \code{\link{pomp}}.
}
\references{}
\author{Aaron A. King (kingaa at umich dot edu)}
Added: pkg/man/trajectory-pomp.Rd
===================================================================
--- pkg/man/trajectory-pomp.Rd (rev 0)
+++ pkg/man/trajectory-pomp.Rd 2008-08-13 16:12:12 UTC (rev 38)
@@ -0,0 +1,43 @@
+\name{trajectory-pomp}
+\docType{methods}
+\alias{trajectory}
+\alias{trajectory,pomp-method}
+\alias{trajectory-pomp}
+
+\title{Compute trajectories of the determinstic skeleton.}
+\description{
+ The method \code{trajectory} computes a trajectory of the deterministic skeleton of a Markov process.
+ In the case of a discrete-time system, the deterministic skeleton is a map and a trajectory is obtained by iterating the map.
+ In the case of a continuous-time system, the deterministic skeleton is a vector-field; \code{trajectory} integrates the vectorfield to obtain a trajectory.
+}
+\usage{
+trajectory(object, params, times, \dots)
+\S4method{trajectory}{pomp}(object, params, times = time(object), \dots)
+}
+\arguments{
+ \item{object}{an object of class \code{pomp}.}
+ \item{times}{
+ a numeric vector containing the times at which a trajectory is desired.
+ The first of these will be the initial time.
+ }
+ \item{params}{
+ a rank-2 array of parameters.
+ Each column of \code{params} is a distinct parameter vector.
+ }
+ \item{\dots}{at present, these are ignored.}
+}
+\value{
+ Returns an array of dimensions \code{nvar} x \code{nreps} x \code{ntimes}.
+ If \code{x} is the returned matrix, \code{x[i,j,k]} is the i-th component of the state vector at time \code{times[k]} given parameters \code{params[,j]}.
+}
+\details{
+ This function makes repeated calls to the user-supplied \code{skeleton} of the \code{pomp} object.
+ For specifications on supplying this, see \code{\link{pomp}}.
+
+ What is the behavior if \code{type="map"} and \code{times} are non-integral?
+}
+\references{}
+\author{Aaron A. King (kingaa at umich dot edu)}
+\seealso{\code{\link{pomp-class}}, \code{\link{pomp}}}
+\keyword{models}
+\keyword{ts}
Modified: pkg/src/skeleton.c
===================================================================
--- pkg/src/skeleton.c 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/src/skeleton.c 2008-08-13 16:12:12 UTC (rev 38)
@@ -6,12 +6,12 @@
#include <Rdefines.h>
#include <Rinternals.h>
-static void eval_vf (pomp_vectorfield_map *vf,
- double *f,
- double *x, double *times, double *params,
- int *ndim,
- int *stateindex, int *parindex, int *covindex,
- double *time_table, double *covar_table)
+static void eval_skel (pomp_vectorfield_map *vf,
+ double *f,
+ double *x, double *times, double *params,
+ int *ndim,
+ int *stateindex, int *parindex, int *covindex,
+ double *time_table, double *covar_table)
{
double t, *xp, *pp, *fp;
int nvar = ndim[0];
@@ -48,7 +48,7 @@
}
}
-// these global objects will pass the needed information to the user-defined function (see 'default_vf_fn')
+// these global objects will pass the needed information to the user-defined function (see 'default_skel_fn')
// each of these is allocated once, globally, and refilled many times
static SEXP _pomp_skel_Xvec; // state variable vector
static SEXP _pomp_skel_Pvec; // parameter vector
@@ -71,9 +71,9 @@
// this is the vectorfield that is evaluated when the user supplies an R function
// (and not a native routine)
// Note that stateindex, parindex, covindex are ignored.
-static void default_vf_fn (double *f, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
- int covdim, double *covar, double t)
+static void default_skel_fn (double *f, double *x, double *p,
+ int *stateindex, int *parindex, int *covindex,
+ int covdim, double *covar, double t)
{
int nprotect = 0;
int k;
@@ -88,6 +88,10 @@
xp = REAL(TIME);
xp[0] = t;
PROTECT(ans = eval(FCALL,RHO)); nprotect++; // evaluate the call
+ if (LENGTH(ans)!=NVAR) {
+ UNPROTECT(nprotect);
+ error("skeleton error: user 'skeleton' must return a vector of length %d",NVAR);
+ }
xp = REAL(AS_NUMERIC(ans));
for (k = 0; k < NVAR; k++) f[k] = xp[k];
UNPROTECT(nprotect);
@@ -100,12 +104,14 @@
int fdim[3], ndim[6];
int use_native;
int nstates, nparams, ncovars;
+ double *xp;
SEXP dimP, dimX, fn, F;
SEXP tcovar, covar;
SEXP statenames, paramnames, covarnames;
SEXP sindex, pindex, cindex;
SEXP Xnames, Pnames, Cnames;
pomp_vectorfield_map *ff = NULL;
+ int k;
ntimes = LENGTH(t);
@@ -158,12 +164,14 @@
if (use_native) {
ff = (pomp_vectorfield_map *) R_ExternalPtrAddr(fn);
} else { // else construct a call to the R function
- ff = (pomp_vectorfield_map *) default_vf_fn;
+ ff = (pomp_vectorfield_map *) default_skel_fn;
PROTECT(RHO = (CLOENV(fn))); nprotect++;
NVAR = nvars; // for internal use
NPAR = npars; // for internal use
PROTECT(TIME = NEW_NUMERIC(1)); nprotect++; // for internal use
PROTECT(XVEC = NEW_NUMERIC(nvars)); nprotect++; // for internal use
+ xp = REAL(XVEC);
+ for (k = 0; k < nvars; k++) xp[k] = 0.0;
PROTECT(PVEC = NEW_NUMERIC(npars)); nprotect++; // for internal use
PROTECT(CVEC = NEW_NUMERIC(covdim)); nprotect++; // for internal use
SET_NAMES(XVEC,Xnames); // make sure the names attribute is copied
@@ -185,6 +193,9 @@
fdim[0] = nvars; fdim[1] = nreps; fdim[2] = ntimes;
PROTECT(F = makearray(3,fdim)); nprotect++;
setrownames(F,Xnames,3);
+ xp = REAL(F);
+ for (k = 0; k < nvars*nreps*ntimes; k++) xp[k] = 0.0;
+
if (nstates > 0) {
PROTECT(sindex = MATCHROWNAMES(x,statenames)); nprotect++;
} else {
@@ -200,12 +211,14 @@
} else {
PROTECT(cindex = NEW_INTEGER(0)); nprotect++;
}
+
ndim[0] = nvars; ndim[1] = npars; ndim[2] = nreps; ndim[3] = ntimes;
ndim[4] = covlen; ndim[5] = covdim;
- eval_vf(ff,REAL(F),REAL(x),REAL(t),REAL(params),
- ndim,INTEGER(sindex),INTEGER(pindex),INTEGER(cindex),
- REAL(tcovar),REAL(covar));
+ eval_skel(ff,REAL(F),REAL(x),REAL(t),REAL(params),
+ ndim,INTEGER(sindex),INTEGER(pindex),INTEGER(cindex),
+ REAL(tcovar),REAL(covar));
+
UNPROTECT(nprotect);
return F;
}
Modified: pkg/tests/logistic.R
===================================================================
--- pkg/tests/logistic.R 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/tests/logistic.R 2008-08-13 16:12:12 UTC (rev 38)
@@ -28,7 +28,7 @@
)
},
measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
- skeleton=function(x,t,params,...){
+ skeleton.vectorfield=function(x,t,params,...){
with(
as.list(c(x,params)),
r*n*(1-n/K)
@@ -83,3 +83,8 @@
digits=4
)
+params <- cbind(c(n.0=100,K=10000,r=0.2,sigma=0.4,tau=0.1),c(n.0=1000,K=11000,r=0.1,sigma=0.4,tau=0.1))
+x <- trajectory(po,params=params)
+pdf(file='logistic.pdf')
+matplot(time(po,t0=TRUE),t(x['n',,]),type='l')
+dev.off()
Modified: pkg/tests/logistic.Rout.save
===================================================================
--- pkg/tests/logistic.Rout.save 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/tests/logistic.Rout.save 2008-08-13 16:12:12 UTC (rev 38)
@@ -1,6 +1,6 @@
-R version 2.6.1 (2007-11-26)
-Copyright (C) 2007 The R Foundation for Statistical Computing
+R version 2.7.1 (2008-06-23)
+Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -16,6 +16,7 @@
Type 'q()' to quit R.
> library(pomp)
+Loading required package: odesolve
>
> po <- pomp(
+ data=rbind(obs=rep(0,1000)),
@@ -45,7 +46,7 @@
+ )
+ },
+ measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
-+ skeleton=function(x,t,params,...){
++ skeleton.vectorfield=function(x,t,params,...){
+ with(
+ as.list(c(x,params)),
+ r*n*(1-n/K)
@@ -106,4 +107,11 @@
[1] 0 810 1440 1890 2160 2250 2160 1890 1440 810 0 -990
[13] -2160
>
+> params <- cbind(c(n.0=100,K=10000,r=0.2,sigma=0.4,tau=0.1),c(n.0=1000,K=11000,r=0.1,sigma=0.4,tau=0.1))
+> x <- trajectory(po,params=params)
+> pdf(file='logistic.pdf')
+> matplot(time(po,t0=TRUE),t(x['n',,]),type='l')
+> dev.off()
+null device
+ 1
>
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/tests/sir.R 2008-08-13 16:12:12 UTC (rev 38)
@@ -77,7 +77,8 @@
}
)
},
- skeleton=function(x,t,params,covars,...) {
+ skeleton.vectorfield=function(x,t,params,covars,...) {
+ xdot <- rep(0,length(x))
params <- exp(params)
with(
as.list(c(x,params)),
@@ -92,12 +93,13 @@
mu*I,
mu*R
)
- c(
- terms[1]-terms[2]-terms[3],
- terms[2]-terms[4]-terms[5],
- terms[4]-terms[6],
- terms[4]
- )
+ xdot[1:4] <- c(
+ terms[1]-terms[2]-terms[3],
+ terms[2]-terms[4]-terms[5],
+ terms[4]-terms[6],
+ terms[4]
+ )
+ xdot
}
)
},
@@ -127,6 +129,8 @@
x <- simulate(po,params=log(params),nsim=3)
toc <- Sys.time()
print(toc-tic)
+pdf(file='sir.pdf')
+plot(x[[1]],variables=c("S","I","R","cases","W"))
t <- seq(0,4/52,by=1/52/25)
X <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t)
@@ -171,6 +175,13 @@
)
print(h[c("S","I","R"),,],digits=4)
+t <- seq(0,20,by=1/52)
+tic <- Sys.time()
+X <- trajectory(po,params=log(params),times=t,hmax=1/52)
+toc <- Sys.time()
+print(toc-tic)
+plot(t,X['I',1,],type='l')
+
po <- pomp(
times=seq(1/52,4,by=1/52),
data=rbind(measles=numeric(52*4)),
@@ -186,7 +197,7 @@
rprocess=euler.simulate,
dens.fun="sir_euler_density",
dprocess=euler.density,
- skeleton="sir_ODE",
+ skeleton.vectorfield="sir_ODE",
PACKAGE="pomp",
measurement.model=measles~binom(size=cases,prob=exp(rho)),
initializer=function(params,t0,...){
@@ -212,6 +223,7 @@
x <- simulate(po,params=log(params),nsim=3)
toc <- Sys.time()
print(toc-tic)
+plot(x[[1]],variables=c("S","I","R","cases","W"))
t <- seq(0,4/52,by=1/52/25)
X <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t)
@@ -255,3 +267,11 @@
params=as.matrix(log(params))
)
print(h[c("S","I","R"),,],digits=4)
+
+t <- seq(0,20,by=1/52)
+tic <- Sys.time()
+X <- trajectory(po,params=log(params),times=t,hmax=1/52)
+toc <- Sys.time()
+print(toc-tic)
+plot(t,X['I',1,],type='l')
+dev.off()
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2008-08-01 15:38:26 UTC (rev 37)
+++ pkg/tests/sir.Rout.save 2008-08-13 16:12:12 UTC (rev 38)
@@ -1,6 +1,6 @@
-R version 2.6.1 (2007-11-26)
-Copyright (C) 2007 The R Foundation for Statistical Computing
+R version 2.7.1 (2008-06-23)
+Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -16,6 +16,7 @@
Type 'q()' to quit R.
> library(pomp)
+Loading required package: odesolve
>
> tbasis <- seq(0,25,by=1/52)
> basis <- periodic.bspline.basis(tbasis,nbasis=3)
@@ -94,7 +95,8 @@
+ }
+ )
+ },
-+ skeleton=function(x,t,params,covars,...) {
++ skeleton.vectorfield=function(x,t,params,covars,...) {
++ xdot <- rep(0,length(x))
+ params <- exp(params)
+ with(
+ as.list(c(x,params)),
@@ -109,12 +111,13 @@
+ mu*I,
+ mu*R
+ )
-+ c(
-+ terms[1]-terms[2]-terms[3],
-+ terms[2]-terms[4]-terms[5],
-+ terms[4]-terms[6],
-+ terms[4]
-+ )
++ xdot[1:4] <- c(
++ terms[1]-terms[2]-terms[3],
++ terms[2]-terms[4]-terms[5],
++ terms[4]-terms[6],
++ terms[4]
++ )
++ xdot
+ }
+ )
+ },
@@ -144,7 +147,9 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 3.671906 secs
+Time difference of 4.001635 secs
+> pdf(file='sir.pdf')
+> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t <- seq(0,4/52,by=1/52/25)
> X <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t)
@@ -199,6 +204,14 @@
I 2847 4893 7554 9945 14452 20004 27020.1
R -36462 -33640 -30013 -26751 -20293 -11629 593.4
>
+> t <- seq(0,20,by=1/52)
+> tic <- Sys.time()
+> X <- trajectory(po,params=log(params),times=t,hmax=1/52)
+> toc <- Sys.time()
+> print(toc-tic)
+Time difference of 9.436516 secs
+> plot(t,X['I',1,],type='l')
+>
> po <- pomp(
+ times=seq(1/52,4,by=1/52),
+ data=rbind(measles=numeric(52*4)),
@@ -214,7 +227,7 @@
+ rprocess=euler.simulate,
+ dens.fun="sir_euler_density",
+ dprocess=euler.density,
-+ skeleton="sir_ODE",
++ skeleton.vectorfield="sir_ODE",
+ PACKAGE="pomp",
+ measurement.model=measles~binom(size=cases,prob=exp(rho)),
+ initializer=function(params,t0,...){
@@ -240,7 +253,8 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.3779790 secs
+Time difference of 0.3902411 secs
+> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t <- seq(0,4/52,by=1/52/25)
> X <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t)
@@ -295,3 +309,14 @@
I 17033 23456 31943 39251 46266 49391
R -12929 -2160 13830 31164 54730 82862
>
+> t <- seq(0,20,by=1/52)
+> tic <- Sys.time()
+> X <- trajectory(po,params=log(params),times=t,hmax=1/52)
+> toc <- Sys.time()
+> print(toc-tic)
+Time difference of 7.409947 secs
+> plot(t,X['I',1,],type='l')
+> dev.off()
+null device
+ 1
+>
More information about the pomp-commits
mailing list