[Pomp-commits] r423 - in pkg: . R data inst inst/doc man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Nov 23 22:02:58 CET 2010
Author: kingaa
Date: 2010-11-23 22:02:58 +0100 (Tue, 23 Nov 2010)
New Revision: 423
Added:
pkg/src/trajectory.c
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/aaa.R
pkg/R/dmeasure-pomp.R
pkg/R/pomp-fun.R
pkg/R/pomp.R
pkg/R/rmeasure-pomp.R
pkg/R/skeleton-pomp.R
pkg/R/traj-match.R
pkg/R/trajectory-pomp.R
pkg/data/dacca.rda
pkg/data/gillespie.sir.rda
pkg/data/ou2.rda
pkg/data/ricker.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/ChangeLog
pkg/inst/NEWS
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/man/traj-match.Rd
pkg/man/trajectory-pomp.Rd
pkg/src/dmeasure.c
pkg/src/pomp_fun.c
pkg/src/pomp_internal.h
pkg/src/rmeasure.c
pkg/src/simulate.c
pkg/src/skeleton.c
Log:
- when PACKAGE is undefined in 'pomp' or 'pomp.fun', the default is now "" (rather than 'character(0)')
- rewrite 'simulate' and 'trajectory' to minimize extraction from 'pomp.fun' objects (and consequent lookup of dynamically-loaded symbols).
- rewrite 'skeleton', 'rmeasure', and 'dmeasure' in keeping with above
- rewrite 'traj.match' internals
- add 'sannbox' box-constrained simulated annealing algorithm to 'traj.match' suite of optimizers
- 'traj.match' generates an error when 'est' is not supplied and 'eval.only=FALSE'
- for continuous-time POMPs, 'trajectory' computation is now done by a single call to 'deSolve::ode' (vs. one call for each parameter set as before)
- new .Call-able function 'get_pomp_fun' to extract from 'pomp.fun' objects
- iteration of discrete-time deterministic skeleton maps is now done in C for speed (new file 'trajectory.c')
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/DESCRIPTION 2010-11-23 21:02:58 UTC (rev 423)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.36-1
-Date: 2010-11-14
+Date: 2010-11-23
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing,
Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
@@ -13,7 +13,7 @@
LazyLoad: true
LazyData: false
Collate: aaa.R eulermultinom.R plugins.R
- slice-design.R profile-design.R sobol.R bsplines.R
+ slice-design.R profile-design.R sobol.R bsplines.R sannbox.R
pomp-fun.R pomp.R pomp-methods.R rmeasure-pomp.R rprocess-pomp.R init-state-pomp.R
dmeasure-pomp.R dprocess-pomp.R skeleton-pomp.R simulate-pomp.R trajectory-pomp.R plot-pomp.R
pfilter.R traj-match.R bsmc.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/NAMESPACE 2010-11-23 21:02:58 UTC (rev 423)
@@ -1,5 +1,6 @@
useDynLib(
pomp,
+ get_pomp_fun,
bspline_basis,
periodic_bspline_basis,
bspline_basis_function,
@@ -7,16 +8,13 @@
euler_model_simulator,
euler_model_density,
SSA_simulator,
- R_Euler_Multinom,
- D_Euler_Multinom,
+ R_Euler_Multinom,D_Euler_Multinom,
pfilter_computations,
simulation_computations,
- apply_probe_data,
- apply_probe_sim,
- probe_marginal_setup,
- probe_marginal_solve,
- probe_acf,
- probe_ccf,
+ iterate_map,traj_transp_and_copy,
+ apply_probe_data,apply_probe_sim,
+ probe_marginal_setup,probe_marginal_solve,
+ probe_acf,probe_ccf,
probe_nlar,
synth_loglik,
do_rprocess,
@@ -31,7 +29,7 @@
importFrom(stats,simulate,time,coef,logLik,window)
importFrom(mvtnorm,dmvnorm,rmvnorm)
importFrom(subplex,subplex)
-importFrom(deSolve,ode,lsoda)
+importFrom(deSolve,ode)
exportClasses(
"pomp","pfilterd.pomp",
@@ -71,6 +69,8 @@
periodic.bspline.basis,
compare.mif,
nlf,
+ skeleton,
+ trajectory,
traj.match,
probe.mean,
probe.median,
Modified: pkg/R/aaa.R
===================================================================
--- pkg/R/aaa.R 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/R/aaa.R 2010-11-23 21:02:58 UTC (rev 423)
@@ -6,13 +6,12 @@
## cat("See the NEWS file for important information\n")
##}
-setGeneric("print",function(x,...)standardGeneric("print"))
-setGeneric("summary",function(object,...)standardGeneric("summary"))
-setGeneric("logLik",function(object,...)standardGeneric("logLik"))
+setGeneric("print")
+setGeneric("summary")
+setGeneric("plot",function(x,y,...)standardGeneric("plot"))
setGeneric("simulate",function(object,nsim=1,seed=NULL,...)standardGeneric("simulate"))
setGeneric("time",function(x,...)standardGeneric("time"))
setGeneric("coef",function(object,...)standardGeneric("coef"))
+setGeneric("logLik",function(object,...)standardGeneric("logLik"))
setGeneric("window",function(x,...)standardGeneric("window"))
-setGeneric("plot",function(x,y,...)standardGeneric("plot"))
-
setGeneric("continue",function(object,...)standardGeneric("continue"))
Modified: pkg/R/dmeasure-pomp.R
===================================================================
--- pkg/R/dmeasure-pomp.R 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/R/dmeasure-pomp.R 2010-11-23 21:02:58 UTC (rev 423)
@@ -4,6 +4,8 @@
setMethod(
'dmeasure',
'pomp',
- function (object, y, x, times, params, log = FALSE, ...)
- .Call(do_dmeasure,object,y,x,times,params,log),
+ function (object, y, x, times, params, log = FALSE, ...) {
+ fun <- get.pomp.fun(object at dmeasure)
+ .Call(do_dmeasure,object,y,x,times,params,log,fun)
+ }
)
Modified: pkg/R/pomp-fun.R
===================================================================
--- pkg/R/pomp-fun.R 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/R/pomp-fun.R 2010-11-23 21:02:58 UTC (rev 423)
@@ -11,8 +11,7 @@
## constructor
pomp.fun <- function (f = NULL, PACKAGE, proto = NULL) {
- if (missing(PACKAGE))
- PACKAGE <- character(0)
+ if (missing(PACKAGE)) PACKAGE <- ""
if (is.function(f)) {
if (!is.null(proto)) {
if (!is.call(proto))
@@ -70,3 +69,15 @@
'pomp.fun',
function (x, ...) show(x)
)
+
+get.pomp.fun <- function (object) {
+ use <- as.integer(object at use-1)
+ if (use==0) {
+ f <- object at R.fun
+ } else if (use==1) {
+ f <- getNativeSymbolInfo(name=object at native.fun,PACKAGE=object at PACKAGE)$address
+ } else {
+ stop("invalid ",sQuote("use")," value")
+ }
+ list(f,use)
+}
Modified: pkg/R/pomp.R
===================================================================
--- pkg/R/pomp.R 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/R/pomp.R 2010-11-23 21:02:58 UTC (rev 423)
@@ -72,6 +72,8 @@
stop("pomp error: the zero-time ",sQuote("t0")," must occur no later than the first observation",call.=TRUE)
storage.mode(t0) <- 'double'
+ if (missing(PACKAGE)) PACKAGE <- ""
+
if (missing(rprocess))
rprocess <- function(xstart,times,params,...)stop(sQuote("rprocess")," not specified")
if (missing(dprocess))
@@ -115,8 +117,6 @@
initializer <- default.initializer
}
- if (missing(PACKAGE)) PACKAGE <- character(0)
-
if (!is.function(rprocess))
stop(
"pomp error: ",sQuote("rprocess")," must be a function",
Modified: pkg/R/rmeasure-pomp.R
===================================================================
--- pkg/R/rmeasure-pomp.R 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/R/rmeasure-pomp.R 2010-11-23 21:02:58 UTC (rev 423)
@@ -4,6 +4,8 @@
setMethod(
'rmeasure',
'pomp',
- function (object, x, times, params, ...)
- .Call(do_rmeasure,object,x,times,params),
+ function (object, x, times, params, ...) {
+ fun <- get.pomp.fun(object at rmeasure)
+ .Call(do_rmeasure,object,x,times,params,fun)
+ }
)
Modified: pkg/R/skeleton-pomp.R
===================================================================
--- pkg/R/skeleton-pomp.R 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/R/skeleton-pomp.R 2010-11-23 21:02:58 UTC (rev 423)
@@ -1,9 +1,14 @@
-setGeneric("skeleton",function(object,...)standardGeneric("skeleton"))
+skeleton <- function (object, x, t, params, ...)
+ stop("function ",sQuote("skeleton")," is undefined for objects of class ",sQuote(class(object)))
+setGeneric('skeleton')
+
## evaluate the measurement model density function
setMethod(
'skeleton',
'pomp',
- function (object, x, t, params, ...)
- .Call(do_skeleton,object,x,t,params),
+ function (object, x, t, params, ...) {
+ skel <- .Call(get_pomp_fun,object at skeleton)
+ .Call(do_skeleton,object,x,t,params,skel)
+ }
)
Modified: pkg/R/traj-match.R
===================================================================
--- pkg/R/traj-match.R 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/R/traj-match.R 2010-11-23 21:02:58 UTC (rev 423)
@@ -10,13 +10,13 @@
)
)
-setMethod("$",signature(x="traj.matched.pomp"),function(x, name) slot(x,name))
+setMethod("$",signature(x="traj.matched.pomp"),function(x, name)slot(x,name))
setMethod("logLik",signature(object="traj.matched.pomp"),function(object,...)object at value)
setMethod(
"summary",
- "traj.matched.pomp",
+ signature=signature(object="traj.matched.pomp"),
function (object, ...) {
c(
list(
@@ -48,45 +48,53 @@
coef(obj,names(start)) <- unname(start)
pmat <- as.matrix(start)
+ obj.fn <- function (x, object, params, t0, ind) {
+ params[ind,] <- x
+ X <- trajectory(object,params=params,t0=t0)
+ d <- dmeasure(
+ object,
+ y=obs(object),
+ x=X,
+ times=time(object),
+ params=params,
+ log=TRUE
+ )
+ -sum(d)
+ }
+
if (eval.only) {
- val <- -sum(
- dmeasure(
- obj,
- y=obs(obj),
- x=trajectory(obj,params=pmat,t0=t0),
- times=time(obj),
- params=pmat,
- log=TRUE
- )
- )
+ val <- obj.fn(numeric(0),object=obj,params=pmat,t0=t0,ind=par.est)
conv <- NA
evals <- c(1,0)
msg <- "no optimization performed"
} else {
- obj.fn <- function (x) {
- pmat[par.est,] <- x
- d <- dmeasure(
- obj,
- y=obs(obj),
- x=trajectory(obj,params=pmat,t0=t0),
- times=time(obj),
- params=pmat,
- log=TRUE
- )
- -sum(d)
- }
-
if (method=="subplex") {
opt <- subplex::subplex(
par=guess,
fn=obj.fn,
- control=list(...)
+ control=list(...),
+ object=obj,
+ params=pmat,
+ t0=t0,
+ ind=par.est
)
+ } else if (method=="sannbox") {
+
+ opt <- sannbox(
+ par=guess,
+ fn=obj.fn,
+ control=list(...),
+ object=obj,
+ params=pmat,
+ t0=t0,
+ ind=par.est
+ )
+
} else {
opt <- optim(
@@ -94,7 +102,11 @@
fn=obj.fn,
gr=gr,
method=method,
- control=list(...)
+ control=list(...),
+ object=obj,
+ params=pmat,
+ t0=t0,
+ ind=par.est
)
}
@@ -123,19 +135,21 @@
)
}
-setGeneric("traj.match",function(object,...)standardGeneric("traj.match"))
+traj.match <- function (object, ...)
+ stop("function ",sQuote("traj.match")," is undefined for objects of class ",sQuote(class(object)))
+setGeneric("traj.match")
+
setMethod(
"traj.match",
signature=signature(object="pomp"),
function (object, start, est,
- method = c("Nelder-Mead","SANN","subplex"),
+ method = c("Nelder-Mead","sannbox","subplex"),
gr = NULL, eval.only = FALSE, ...) {
if (missing(start)) start <- coef(object)
- if (missing(est)) {
- est <- character(0)
- eval.only <- TRUE
- }
+ if (!eval.only && missing(est))
+ stop(sQuote("est")," must be supplied if optimization is to be done")
+ if (eval.only) est <- character(0)
method <- match.arg(method)
traj.match.internal(
object=object,
@@ -153,7 +167,7 @@
"traj.match",
signature=signature(object="traj.matched.pomp"),
function (object, start, est,
- method = c("Nelder-Mead","SANN","subplex"),
+ method = c("Nelder-Mead","sannbox","subplex"),
gr = NULL, eval.only = FALSE, ...) {
if (missing(start)) start <- coef(object)
if (missing(est)) est <- object at est
Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/R/trajectory-pomp.R 2010-11-23 21:02:58 UTC (rev 423)
@@ -1,3 +1,8 @@
+trajectory <- function (object, params, times, t0, ...)
+ stop("function ",sQuote("trajectory")," is undefined for objects of class ",sQuote(class(object)))
+
+setGeneric('trajectory')
+
trajectory.internal <- function (object, params, times, t0, ...) {
warn.condition <- missing(t0)
@@ -10,11 +15,10 @@
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)
@@ -22,11 +26,10 @@
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)
@@ -55,67 +58,65 @@
stop("pfilter error: ",sQuote("params")," must have rownames",call.=FALSE)
params <- as.matrix(params)
- tm <- t0
- x0 <- init.state(object,params=params,t0=tm)
- nm <- rownames(x0)
- dim(x0) <- c(nrow(x0),nrep,1)
- dimnames(x0) <- list(nm,NULL,NULL)
+ x0 <- init.state(object,params=params,t0=t0)
+ nvar <- nrow(x0)
+ statenames <- rownames(x0)
+ dim(x0) <- c(nvar,nrep,1)
+ dimnames(x0) <- list(statenames,NULL,NULL)
- x <- array(
- dim=c(nrow(x0),nrep,length(times)),
- dimnames=list(rownames(x0),NULL,NULL)
+ type <- object at skeleton.type # map or vectorfield?
+
+ if (type=="map") {
+
+ x <- .Call(iterate_map,object,times,t0,x0,params)
+
+ } else if (type=="vectorfield") {
+
+ skel <- get.pomp.fun(object at skeleton)
+
+ ## vectorfield function (RHS of ODE) in 'deSolve' format
+ vf.eval <- function (t, y, ...) {
+ list(
+ .Call(
+ do_skeleton,
+ object,
+ x=array(
+ data=y,
+ dim=c(nvar,nrep,1),
+ dimnames=list(statenames,NULL,NULL)
+ ),
+ t=t,
+ params=params,
+ skel=skel
+ ),
+ NULL
+ )
+ }
+
+ X <- try(
+ ode(
+ y=x0,
+ times=c(t0,times),
+ func=vf.eval,
+ method="lsoda",
+ ...
+ ),
+ silent=FALSE
)
- switch(
- object at skeleton.type,
- map={ # iterate the map
- for (k in seq_along(times)) {
- if (tm < times[k]) {
- while (tm < times[k]) {
- x0[,,] <- skeleton(object,x=x0,t=tm,params=params)
- tm <- tm+1
- }
- }
- x[,,k] <- x0
- }
- },
- vectorfield={ # integrate the vectorfield
- for (j in seq_len(nrep)) {
- X <- try(
- deSolve::lsoda(
- y=x0[,j,1],
- times=c(t0,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')[[1]]!=2)
- warning("abnormal exit from ",sQuote("lsoda"),", istate = ",attr(X,'istate'),call.=FALSE)
- x[,j,] <- t(X[-1,-1])
- }
- },
- unspecified=stop("deterministic skeleton not specified")
- )
+ if (inherits(X,'try-error'))
+ stop("trajectory error: error in ODE integrator",call.=FALSE)
+ if (attr(X,'istate')[1]!=2)
+ warning("abnormal exit from ODE integrator, istate = ",attr(X,'istate'),call.=FALSE)
+
+ x <- array(data=t(X[-1,-1]),dim=c(nvar,nrep,ntimes),dimnames=list(statenames,NULL,NULL))
+
+ } else {
+
+ stop("deterministic skeleton not specified")
+
+ }
+
x
}
-setGeneric('trajectory',function(object,params,times,t0,...)standardGeneric("trajectory"))
-
setMethod("trajectory",signature=signature(object="pomp"),definition=trajectory.internal)
Modified: pkg/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/inst/ChangeLog 2010-11-23 21:02:58 UTC (rev 423)
@@ -1,3 +1,24 @@
+2010-11-17 kingaa
+
+ * [r420] src/dmeasure.c, src/rmeasure.c, src/skeleton.c: - remove
+ some redundant/unnecessary lines
+
+2010-11-16 kingaa
+
+ * [r419] R/skeleton-pomp.R: - remove one layer of try/catch
+ * [r418] src/skeleton.c: - minor tweaks
+
+2010-11-14 kingaa
+
+ * [r417] DESCRIPTION, NAMESPACE, R/aaa.R, R/bsmc.R,
+ R/dmeasure-pomp.R, R/dprocess-pomp.R, R/init-state-pomp.R,
+ R/mif.R, R/rmeasure-pomp.R, R/rprocess-pomp.R, inst/ChangeLog,
+ man/dmeasure-pomp.Rd, man/dprocess-pomp.Rd,
+ man/init.state-pomp.Rd, man/pomp-methods.Rd,
+ man/rmeasure-pomp.Rd, man/rprocess-pomp.Rd, man/skeleton-pomp.Rd,
+ man/traj-match.Rd, man/trajectory-pomp.Rd: - make declarations of
+ generics uniform
+
2010-11-09 kingaa
* [r416] DESCRIPTION, R/traj-match.R,
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/inst/NEWS 2010-11-23 21:02:58 UTC (rev 423)
@@ -4,7 +4,7 @@
o When 'traj.match' is called, the returned 'traj.matched.pomp' object has its 'states' slot filled with the model trajectory at the fitted parameters.
- o More optimization methods are now provided in 'traj.match'.
+ o More optimization methods are now provided in 'traj.match'. These include the new algorithm "sannbox", an optionally box-constrained simulated annealing algorithm.
o There is a change to the interface to 'pfilter' in that 'save.states' results in the filling of the 'last.states' slot. Before version 0.36-1, the 'pfilter' returned a list. The element 'states' of that list corresponds to the slot 'last.states' in version 0.36-1. It was necessary to make this name-change in order to avoid a conflict with the 'states' slot inherited from the 'pomp' class.
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/man/traj-match.Rd
===================================================================
--- pkg/man/traj-match.Rd 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/man/traj-match.Rd 2010-11-23 21:02:58 UTC (rev 423)
@@ -17,10 +17,10 @@
}
\usage{
\S4method{traj.match}{pomp}(object, start, est,
- method = c("Nelder-Mead", "SANN", "subplex"),
+ method = c("Nelder-Mead", "sannbox", "subplex"),
gr = NULL, eval.only = FALSE, \dots)
\S4method{traj.match}{traj.matched.pomp}(object, start, est,
- method = c("Nelder-Mead", "SANN", "subplex"),
+ method = c("Nelder-Mead", "sannbox", "subplex"),
gr = NULL, eval.only = FALSE, \dots)
}
\arguments{
@@ -29,7 +29,7 @@
\item{est}{character vector containing the names of parameters to be estimated.}
\item{method}{
Optimization method.
- Choices are \code{\link[subplex]{subplex}} and any of the methods used by \code{\link{optim}}.
+ Choices are \code{\link[subplex]{subplex}}, \dQuote{sannbox}, and any of the methods used by \code{\link{optim}}.
}
\item{gr}{
Passed to \code{\link{optim}}.
Modified: pkg/man/trajectory-pomp.Rd
===================================================================
--- pkg/man/trajectory-pomp.Rd 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/man/trajectory-pomp.Rd 2010-11-23 21:02:58 UTC (rev 423)
@@ -24,7 +24,10 @@
\code{t0} specifies the start time (the time at which the initial conditions hold).
The default for \code{times} is \code{times=time(object,t0=FALSE)} and \code{t0=timezero(object)}, respectively.
}
- \item{\dots}{at present, these are passed to the ODE integrator if the skeleton is a vectorfield and ignored if it is a map.}
+ \item{\dots}{
+ additional arguments are passed to the ODE integrator if the skeleton is a vectorfield and ignored if it is a map.
+ See \code{\link[deSolve]{ode}} for a description of the additional arguments accepted.
+ }
}
\value{
Returns an array of dimensions \code{nvar} x \code{nreps} x \code{ntimes}.
@@ -34,7 +37,7 @@
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}}.
- When the skeleton is a vectorfield, \code{trajectory} integrates it using \code{\link[deSolve]{lsoda}}.
+ When the skeleton is a vectorfield, \code{trajectory} integrates it using \code{\link[deSolve]{ode}}.
Unresolved issue: What is the behavior if \code{type="map"} and \code{times} are non-integral?
}
@@ -50,6 +53,6 @@
plot(time(euler.sir),x["I",1,],type='l',xlab='time',ylab='I')
lines(time(euler.sir)[-1],diff(x["cases",1,]),col='red')
}
-\seealso{\code{\link{pomp-class}}, \code{\link{pomp}}, \code{\link[deSolve]{lsoda}}}
+\seealso{\code{\link{pomp}}, \code{\link{traj.match}}, \code{\link[deSolve]{ode}}}
\keyword{models}
\keyword{ts}
Modified: pkg/src/dmeasure.c
===================================================================
--- pkg/src/dmeasure.c 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/src/dmeasure.c 2010-11-23 21:02:58 UTC (rev 423)
@@ -104,7 +104,7 @@
UNPROTECT(nprotect);
}
-SEXP do_dmeasure (SEXP object, SEXP y, SEXP x, SEXP times, SEXP params, SEXP log)
+SEXP do_dmeasure (SEXP object, SEXP y, SEXP x, SEXP times, SEXP params, SEXP log, SEXP fun)
{
int nprotect = 0;
int *dim, nvars, npars, nreps, ntimes, covlen, covdim, nobs;
@@ -163,12 +163,14 @@
PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
- PROTECT(
- fn = pomp_fun_handler(
- GET_SLOT(object,install("dmeasure")),
- &use_native
- )
- ); nprotect++;
+ PROTECT(fn = VECTOR_ELT(fun,0)); nprotect++;
+ use_native = INTEGER(VECTOR_ELT(fun,1))[0];
+// PROTECT(
+// fn = pomp_fun_handler(
+// GET_SLOT(object,install("dmeasure")),
+// &use_native
+// )
+// ); nprotect++;
switch (use_native) {
case 0: // use R function
@@ -245,8 +247,7 @@
ndim[4] = covlen; ndim[5] = covdim; ndim[6] = nobs;
dens_meas(ff,REAL(F),REAL(y),REAL(x),REAL(times),REAL(params),INTEGER(log),ndim,
- oidx,sidx,pidx,cidx,
- REAL(tcovar),REAL(covar));
+ oidx,sidx,pidx,cidx,REAL(tcovar),REAL(covar));
UNPROTECT(nprotect);
return F;
Modified: pkg/src/pomp_fun.c
===================================================================
--- pkg/src/pomp_fun.c 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/src/pomp_fun.c 2010-11-23 21:02:58 UTC (rev 423)
@@ -3,42 +3,55 @@
#include <R.h>
#include <Rdefines.h>
#include <Rinternals.h>
+#include <R_ext/Rdynload.h>
#include "pomp_internal.h"
// returns either the R function or the address of the native routine
// on return, use_native tells whether to use the native or the R function
-
SEXP pomp_fun_handler (SEXP pfun, int *use_native)
{
int nprotect = 0;
- SEXP nsi, f = R_NilValue;
- int use = INTEGER(GET_SLOT(pfun,install("use")))[0];
+ SEXP nf, pack, nsi, f = R_NilValue;
+ int use;
+ use = INTEGER(GET_SLOT(pfun,install("use")))[0];
+
switch (use) {
case 1: // use R function
PROTECT(f = GET_SLOT(pfun,install("R.fun"))); nprotect++;
*use_native = 0;
break;
case 2: // use native routine
- PROTECT(nsi = eval(
- lang3(
- install("getNativeSymbolInfo"),
- GET_SLOT(pfun,install("native.fun")),
- GET_SLOT(pfun,install("PACKAGE"))
- ),
- R_GlobalEnv
- )
- ); nprotect++;
+ PROTECT(nf = GET_SLOT(pfun,install("native.fun"))); nprotect++;
+ PROTECT(pack = GET_SLOT(pfun,install("PACKAGE"))); nprotect++;
+ if (LENGTH(pack) < 1) {
+ PROTECT(pack = NEW_CHARACTER(1)); nprotect++;
+ SET_STRING_ELT(pack,0,mkChar(""));
+ }
+ PROTECT(nsi = eval(lang3(install("getNativeSymbolInfo"),nf,pack),R_BaseEnv)); nprotect++;
*use_native = 1;
- PROTECT(f = VECTOR_ELT(nsi,1)); nprotect++;
+ PROTECT(f = VECTOR_ELT(nsi,1)); nprotect++; // return a pointer to the loaded routine
break;
default:
- UNPROTECT(nprotect);
- error("'pomp_fun_handler': invalid value of 'use'");
+ error("'pomp_fun_handler': invalid 'use' value");
break;
}
UNPROTECT(nprotect);
return f;
}
+
+// returns either the R function or the address of the native routine
+// returns a list of length 2
+SEXP get_pomp_fun (SEXP pfun) {
+ int nprotect = 0;
+ SEXP ans, use, f = R_NilValue;
+ PROTECT(use = NEW_INTEGER(1)); nprotect++;
+ PROTECT(f = pomp_fun_handler(pfun,INTEGER(use))); nprotect++;
+ PROTECT(ans = NEW_LIST(2)); nprotect++;
+ SET_ELEMENT(ans,0,f); // R function or pointer to native routine
+ SET_ELEMENT(ans,1,use); // =0 if R function; =1 if native
+ UNPROTECT(nprotect);
+ return ans;
+}
Modified: pkg/src/pomp_internal.h
===================================================================
--- pkg/src/pomp_internal.h 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/src/pomp_internal.h 2010-11-23 21:02:58 UTC (rev 423)
@@ -35,6 +35,7 @@
// pomp_fun.c
SEXP pomp_fun_handler (SEXP pfun, int *use_native);
+SEXP get_pomp_fun (SEXP pfun);
// lookup_table.c
SEXP lookup_in_table (SEXP ttable, SEXP xtable, SEXP t, int *index);
@@ -49,9 +50,12 @@
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);
+SEXP do_rmeasure (SEXP object, SEXP x, SEXP times, SEXP params, SEXP fun);
+// skeleton.c
+SEXP do_skeleton (SEXP object, SEXP x, SEXP t, SEXP params, SEXP fun);
+
static R_INLINE SEXP makearray (int rank, int *dim) {
int nprotect = 0;
int *dimp, k;
Modified: pkg/src/rmeasure.c
===================================================================
--- pkg/src/rmeasure.c 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/src/rmeasure.c 2010-11-23 21:02:58 UTC (rev 423)
@@ -119,7 +119,7 @@
UNPROTECT(nprotect);
}
-SEXP do_rmeasure (SEXP object, SEXP x, SEXP times, SEXP params)
+SEXP do_rmeasure (SEXP object, SEXP x, SEXP times, SEXP params, SEXP fun)
{
int nprotect = 0;
int *dim, nvars, npars, nreps, ntimes, covlen, covdim, nobs;
@@ -174,12 +174,8 @@
PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
PROTECT(OBNM = GET_ROWNAMES(GET_DIMNAMES(GET_SLOT(object,install("data"))))); nprotect++;
- PROTECT(
- fn = pomp_fun_handler(
- GET_SLOT(object,install("rmeasure")),
- &use_native
- )
- ); nprotect++;
+ PROTECT(fn = VECTOR_ELT(fun,0)); nprotect++;
+ use_native = INTEGER(VECTOR_ELT(fun,1))[0];
switch (use_native) {
case 0: // use R function
Modified: pkg/src/simulate.c
===================================================================
--- pkg/src/simulate.c 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/src/simulate.c 2010-11-23 21:02:58 UTC (rev 423)
@@ -13,6 +13,7 @@
SEXP ans, ans_names;
SEXP po, popo;
SEXP statenames, paramnames, obsnames, statedim, obsdim;
+ SEXP rmeas;
int nsims, nparsets, nreps, npars, nvars, ntimes, nobs;
int qobs, qstates;
int *dim, dims[3];
@@ -100,10 +101,11 @@
} else { // we must do 'rmeasure'
+ PROTECT(rmeas = get_pomp_fun((GET_SLOT(object,install("rmeasure"))))); nprotect++;
if (nsims > 1) {
- PROTECT(y = do_rmeasure(object,x,times,p)); nprotect++;
+ PROTECT(y = do_rmeasure(object,x,times,p,rmeas)); nprotect++;
} else {
- PROTECT(y = do_rmeasure(object,x,times,params)); nprotect++;
+ PROTECT(y = do_rmeasure(object,x,times,params,rmeas)); nprotect++;
}
if (qobs) {
Modified: pkg/src/skeleton.c
===================================================================
--- pkg/src/skeleton.c 2010-11-23 20:40:41 UTC (rev 422)
+++ pkg/src/skeleton.c 2010-11-23 21:02:58 UTC (rev 423)
@@ -115,7 +115,7 @@
UNPROTECT(nprotect);
}
-SEXP do_skeleton (SEXP object, SEXP x, SEXP t, SEXP params)
+SEXP do_skeleton (SEXP object, SEXP x, SEXP t, SEXP params, SEXP fun)
{
int nprotect = 0;
int *dim, nvars, npars, nreps, ntimes, covlen, covdim;
@@ -127,7 +127,7 @@
SEXP tcovar, covar;
SEXP statenames, paramnames, covarnames;
SEXP sindex, pindex, cindex;
- SEXP Pnames, Cnames;
+ SEXP Snames, Pnames, Cnames;
pomp_vectorfield_map *ff = NULL;
int k, len;
@@ -148,7 +148,7 @@
npars = dim[0];
if (nreps != dim[1])
error("skeleton error: 2nd dimensions of 'params' and 'x' do not agree");
-
+
PROTECT(tcovar = GET_SLOT(object,install("tcovar"))); nprotect++;
PROTECT(covar = GET_SLOT(object,install("covar"))); nprotect++;
PROTECT(statenames = GET_SLOT(object,install("statenames"))); nprotect++;
@@ -159,17 +159,13 @@
ncovars = LENGTH(covarnames);
dim = INTEGER(GET_DIM(covar)); covlen = dim[0]; covdim = dim[1];
- PROTECT(VNAMES = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
+ PROTECT(Snames = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
- PROTECT(
- fn = pomp_fun_handler(
- GET_SLOT(object,install("skeleton")),
- &use_native
- )
- ); nprotect++;
-
+ PROTECT(fn = VECTOR_ELT(fun,0)); nprotect++;
+ use_native = INTEGER(VECTOR_ELT(fun,1))[0];
+
if (use_native) {
ff = (pomp_vectorfield_map *) R_ExternalPtrAddr(fn);
} else { // else construct a call to the R function
@@ -183,7 +179,7 @@
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,VNAMES); // make sure the names attribute is copied
+ SET_NAMES(XVEC,Snames); // make sure the names attribute is copied
SET_NAMES(PVEC,Pnames); // make sure the names attribute is copied
SET_NAMES(CVEC,Cnames); // make sure the names attribute is copied
// set up the function call
@@ -197,11 +193,12 @@
PROTECT(FCALL = LCONS(XVEC,FCALL)); nprotect++;
SET_TAG(FCALL,install("x"));
PROTECT(FCALL = LCONS(fn,FCALL)); nprotect++;
+ PROTECT(VNAMES = duplicate(Snames)); nprotect++; // for internal use
}
fdim[0] = nvars; fdim[1] = nreps; fdim[2] = ntimes;
PROTECT(F = makearray(3,fdim)); nprotect++;
- setrownames(F,VNAMES,3);
+ setrownames(F,Snames,3);
xp = REAL(F);
for (k = 0, len = nvars*nreps*ntimes; k < len; k++) xp[k] = 0.0;
@@ -232,11 +229,12 @@
return F;
}
-#undef FCALL
#undef XVEC
#undef PVEC
#undef CVEC
#undef TIME
#undef NVAR
-#undef NPAR
-#undef RHO
+#undef NPAR
+#undef RHO
+#undef FCALL
+#undef VNAMES
Added: pkg/src/trajectory.c
===================================================================
--- pkg/src/trajectory.c (rev 0)
+++ pkg/src/trajectory.c 2010-11-23 21:02:58 UTC (rev 423)
@@ -0,0 +1,89 @@
+// -*- C++ -*-
+
+#include <Rdefines.h>
+
+#include "pomp_internal.h"
+
+// copy t(x[-1,-1]) -> y[,rep,]
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 423
More information about the pomp-commits
mailing list