[Pomp-commits] r683 - in pkg/pomp: . R data demo inst inst/data-R inst/doc inst/examples inst/include man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 26 16:41:38 CEST 2012
Author: kingaa
Date: 2012-04-26 16:41:38 +0200 (Thu, 26 Apr 2012)
New Revision: 683
Added:
pkg/pomp/man/builder.Rd
pkg/pomp/src/userdata.c
Removed:
pkg/pomp/inst/examples/gompertz.c
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/NAMESPACE
pkg/pomp/R/builder.R
pkg/pomp/data/bbs.rda
pkg/pomp/data/blowflies.rda
pkg/pomp/data/dacca.rda
pkg/pomp/data/euler.sir.rda
pkg/pomp/data/gillespie.sir.rda
pkg/pomp/data/gompertz.rda
pkg/pomp/data/ou2.rda
pkg/pomp/data/ricker.rda
pkg/pomp/data/rw2.rda
pkg/pomp/data/verhulst.rda
pkg/pomp/demo/gompertz.R
pkg/pomp/demo/logistic.R
pkg/pomp/inst/ChangeLog
pkg/pomp/inst/NEWS
pkg/pomp/inst/data-R/sir.R
pkg/pomp/inst/doc/advanced_topics_in_pomp.Rnw
pkg/pomp/inst/doc/advanced_topics_in_pomp.pdf
pkg/pomp/inst/doc/intro_to_pomp.Rnw
pkg/pomp/inst/doc/intro_to_pomp.pdf
pkg/pomp/inst/include/pomp.h
pkg/pomp/src/R_init_pomp.c
pkg/pomp/src/SSA_wrapper.c
pkg/pomp/src/dmeasure.c
pkg/pomp/src/euler.c
pkg/pomp/src/partrans.c
pkg/pomp/src/pomp.h
pkg/pomp/src/pomp_internal.h
pkg/pomp/src/rmeasure.c
pkg/pomp/src/sir.c
pkg/pomp/src/skeleton.c
pkg/pomp/tests/bbs.Rout.save
pkg/pomp/tests/pfilter.Rout.save
pkg/pomp/tests/skeleton.R
pkg/pomp/tests/skeleton.Rout.save
Log:
- add new facility for accessing 'userdata' from within user-defined native routines: get_pomp_userdata{,_int,_double}
- change some data()-loadable examples to make use of this
- make 'pompBuilder' visible
- modify 'gompertz' demo to make use of 'pompBuilder'
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/DESCRIPTION 2012-04-26 14:41:38 UTC (rev 683)
@@ -1,9 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.41-8
-Date: 2012-04-24
-Revision: $Rev$
+Version: 0.42-1
+Date: 2012-04-26
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>
URL: http://pomp.r-forge.r-project.org
Modified: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/NAMESPACE 2012-04-26 14:41:38 UTC (rev 683)
@@ -90,5 +90,6 @@
probe.marginal,
sannbox,
spect.match,
- traj.match.objfun
+ traj.match.objfun,
+ pompBuilder
)
Modified: pkg/pomp/R/builder.R
===================================================================
--- pkg/pomp/R/builder.R 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/R/builder.R 2012-04-26 14:41:38 UTC (rev 683)
@@ -86,7 +86,8 @@
decl <- list(
periodic_bspline_basis_eval="void (*periodic_bspline_basis_eval)(double,double,int,int,double*);\nperiodic_bspline_basis_eval = (void (*)(double,double,int,int,double*)) R_GetCCallable(\"pomp\",\"periodic_bspline_basis_eval\");\n",
- reulermultinom="void (*reulermultinom)(int,double,double*,double,double*);\nreulermultinom = (void (*)(int,double,double*,double,double*)) R_GetCCallable(\"pomp\",\"reulermultinom\");\n"
+ reulermultinom="void (*reulermultinom)(int,double,double*,double,double*);\nreulermultinom = (void (*)(int,double,double*,double,double*)) R_GetCCallable(\"pomp\",\"reulermultinom\");\n",
+ get_pomp_userdata="const SEXP (*get_pomp_userdata)(const char *);\npomp_get_userdata = (const SEXP (*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata\");\n"
)
footer <- list(
Modified: pkg/pomp/data/bbs.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/blowflies.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/gompertz.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/demo/gompertz.R
===================================================================
--- pkg/pomp/demo/gompertz.R 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/demo/gompertz.R 2012-04-26 14:41:38 UTC (rev 683)
@@ -44,98 +44,82 @@
return(f)
},
parameter.inv.transform=function(params,...){
- params <- log(params[c("X.0","r","K","tau","sigma")])
- names(params) <- c("log.X.0","log.r","log.K","log.tau","log.sigma")
- params
+ log(params)
},
parameter.transform=function(params,...){
- params <- exp(params[c("log.X.0","log.r","log.K","log.tau","log.sigma")])
- names(params) <- c("X.0","r","K","tau","sigma")
- params
+ exp(params)
}
) -> gompertz
## Now code up the Gompertz example using native routines results in much faster computations.
-## The C codes are included in the "examples" directory (file "gompertz.c")
-if (Sys.info()['sysname']=='Linux') { # only run this under linux
+dmeas <- "
+ lik = dlnorm(Y,log(X),tau,give_log);
+"
+rmeas <- "
+ Y = rlnorm(log(X),tau);
+
+"
+step.fn <- "
+ double S = exp(-r*dt);
+ double logeps = (sigma > 0.0) ? rnorm(0,sigma) : 0.0;
+ /* note that X is over-written by the next line */
+ X = pow(K,(1-S))*pow(X,S)*exp(logeps);
+"
+skel <- "
+ double dt = 1.0;
+ double S = exp(-r*dt);
+ /* note that X is not over-written in the skeleton function */
+ DX = pow(K,1-S)*pow(X,S);
+"
- model <- "gompertz"
- ## name of the file holding the native codes for this model:
- modelfile <- system.file("examples",paste(model,".c",sep=""),package="pomp")
- ## name of the shared-object library
- solib <- paste(model,.Platform$dynlib.ext,sep="")
- ## environment variables needed to locate the pomp header file:
- cflags <- paste("PKG_CFLAGS=\"-I",system.file("include/",package="pomp"),"\"")
+pompBuilder(
+ name="Gompertz",
+ data=data.frame(t=1:100,Y=NA),
+ times="t",
+ t0=0,
+ paramnames=c("r","K","sigma","X.0","tau"),
+ statenames=c("X"),
+ zeronames=character(0),
+ dmeasure=dmeas,
+ rmeasure=rmeas,
+ step.fn=step.fn,
+ step.fn.delta.t=1,
+ skeleton=skel,
+ skeleton.type="map",
+ skelmap.delta.t=1
+ ) -> Gompertz
- ## compile the shared-object library containing the model codes:
- rv <- system2(
- command=R.home("bin/R"),
- args=c("CMD","SHLIB","-o",solib,modelfile),
- env=cflags
- )
- ## compile the shared-object library containing the model codes:
- if (rv!=0)
- stop("cannot compile shared-object library ",sQuote(solib))
+## simulate some data
+Gompertz <- simulate(Gompertz,params=c(K=1,r=0.1,sigma=0.1,tau=0.1,X.0=1))
- dyn.load(solib) ## load the shared-object library
+p <- parmat(coef(Gompertz),nrep=4)
+p["X.0",] <- c(0.5,0.9,1.1,1.5)
- po <- pomp(
- data=data.frame(time=seq(0,100,by=1),Y=NA),
- times="time",
- t0=0,
- rprocess=discrete.time.sim(
- step.fun="gompertz_simulator",
- PACKAGE="gompertz"
- ),
- rmeasure="gompertz_normal_rmeasure",
- dmeasure="gompertz_normal_dmeasure",
- skeleton.type="map",
- skeleton="gompertz_skeleton",
- PACKAGE="gompertz",
- paramnames=c("r","K","sigma","tau"),
- statenames=c("X"),
- obsnames=c("Y"),
- parameter.transform=function(params,...){
- params <- log(params[c("X.0","r","K","tau","sigma")])
- names(params) <- c("log.X.0","log.r","log.K","log.tau","log.sigma")
- params
- },
- parameter.inv.transform=function(params,...){
- params <- exp(params[c("log.X.0","log.r","log.K","log.tau","log.sigma")])
- names(params) <- c("X.0","r","K","tau","sigma")
- params
- }
- )
+## compute a trajectory of the deterministic skeleton
+tic <- Sys.time()
+X <- trajectory(Gompertz,params=p,as.data.frame=TRUE)
+toc <- Sys.time()
+print(toc-tic)
+X <- reshape(X,dir="wide",v.names="X",timevar="traj",idvar="time")
+matplot(X$time,X[-1],type='l',lty=1,bty='l',xlab="time",ylab="X",
+ main="Gompertz model\ndeterministic trajectories")
- ## set the parameters
- coef(po) <- c(K=1,r=0.1,sigma=0.1,tau=0.1,X.0=1)
+## simulate from the model
+tic <- Sys.time()
+x <- simulate(Gompertz,params=p,as.data.frame=TRUE)
+toc <- Sys.time()
+print(toc-tic)
- p <- parmat(coef(po),nrep=4)
- p["X.0",] <- c(0.5,0.9,1.1,1.5)
-
- ## compute a trajectory of the deterministic skeleton
- tic <- Sys.time()
- X <- trajectory(po,params=p,as.data.frame=TRUE)
- toc <- Sys.time()
- print(toc-tic)
- X <- reshape(X,dir="wide",v.names="X",timevar="traj",idvar="time")
- matplot(X$time,X[-1],type='l',lty=1,bty='l',xlab="time",ylab="X",
- main="Gompertz model\ndeterministic trajectories")
-
- ## simulate from the model
- tic <- Sys.time()
- x <- simulate(po,params=p,as.data.frame=TRUE)
- toc <- Sys.time()
- print(toc-tic)
+x <- reshape(x,dir="wide",v.names=c("Y","X"),timevar="sim",idvar="time")
+op <- par(mfrow=c(2,1),mgp=c(2,1,0),mar=c(3,3,0,0),bty='l')
+matplot(x$time,x[c("X.1","X.2","X.3")],lty=1,type='l',xlab="time",ylab="X",
+ main="Gompertz model\nstochastic simulations")
+matplot(x$time,x[c("Y.1","Y.2","Y.3")],lty=1,type='l',xlab="time",ylab="Y")
+par(op)
- x <- reshape(x,dir="wide",v.names=c("Y","X"),timevar="sim",idvar="time")
- op <- par(mfrow=c(2,1),mgp=c(2,1,0),mar=c(3,3,0,0),bty='l')
- matplot(x$time,x[c("X.1","X.2","X.3")],lty=1,type='l',xlab="time",ylab="X",
- main="Gompertz model\nstochastic simulations")
- matplot(x$time,x[c("Y.1","Y.2","Y.3")],lty=1,type='l',xlab="time",ylab="Y")
- par(op)
-
- dyn.unload(solib)
-
-}
+## run a particle filter
+tic <- Sys.time()
+pf <- pfilter(Gompertz,Np=500)
+toc <- Sys.time()
+print(toc-tic)
Modified: pkg/pomp/demo/logistic.R
===================================================================
--- pkg/pomp/demo/logistic.R 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/demo/logistic.R 2012-04-26 14:41:38 UTC (rev 683)
@@ -49,12 +49,16 @@
params <- c(n.0=10000,K=10000,r=0.9,sigma=0.4,tau=0.1)
set.seed(73658676)
po <- simulate(po,params=params)
-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,as.data.frame=TRUE)
-x <- reshape(x,dir="wide",idvar="time",timevar="traj")
-matplot(x$time,x[-1],type='l',bty='l',lty=1,xlab="time",ylab="n")
+traj <- trajectory(po,params=params,as.data.frame=TRUE)
+traj <- reshape(traj,dir="wide",idvar="time",timevar="traj")
+sim <- simulate(po,params=params,as.data.frame=TRUE,seed=34597368L)
+sim <- reshape(sim,dir="wide",idvar="time",timevar="sim")
+matplot(range(time(po)),range(c(traj[-1],sim[-1])),type='n',bty='l',lty=1,xlab="time",ylab="n",
+ main="simulations vs trajectories")
+matlines(traj$time,traj[-1],type='l',bty='l',lty=1,xlab="time",ylab="n")
+matlines(sim$time,sim[c("n.1","n.2")],type='l',bty='l',lty=1,xlab="time",ylab="n")
Modified: pkg/pomp/inst/ChangeLog
===================================================================
--- pkg/pomp/inst/ChangeLog 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/inst/ChangeLog 2012-04-26 14:41:38 UTC (rev 683)
@@ -1,5 +1,18 @@
+2012-04-25 kingaa
+
+ * [r682] inst/doc/advanced_topics_in_pomp.Rnw: - minor tweak to
+ pompBuilder example
+ * [r681] NAMESPACE, R/trajectory-pomp.R, src/trajectory.c: - put in
+ a routine to reset external variables after deSolve call
+ * [r680] src/bspline.c: - fix error message
+
2012-04-24 kingaa
+ * [r679] DESCRIPTION, R/builder.R, inst/ChangeLog, inst/NEWS,
+ inst/TODO, inst/doc/Makefile,
+ inst/doc/advanced_topics_in_pomp.Rnw,
+ inst/doc/advanced_topics_in_pomp.pdf: - put 'pompBuilder'
+ demonstration into 'advanced topics' vignette
* [r678] R/aaa.R, R/nlf-funcs.R, R/plot-pomp.R, R/pomp.R,
R/probe-match.R, R/probe.R, R/spect.R: - use 'paste0' where
appropriate. define 'paste0' for use with R <= 2.15.0
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/inst/NEWS 2012-04-26 14:41:38 UTC (rev 683)
@@ -1,4 +1,15 @@
NEWS
+0.42-1
+ o An EXPERIMENTAL facility for constructing pomp objects out of native C routines is now included (pompBuilder).
+ This facility is being actively developed and future changes may not be backward compatible.
+
+ o A new facility allowing access to the 'userdata' slot of a pomp object from compiled 'rprocess', 'rmeasure', 'dprocess', 'dmeasure', 'skeleton', and 'partrans' codes is now available.
+ Calls to the new C routines 'get_pomp_userdata', 'get_pomp_userdata_int', and 'get_pomp_userdata_double' allow retrieval of these elements.
+
+ o Some of the data()-loadable examples have been reworked to make use of the above facility.
+ The parameterization of these examples has changed.
+ These changes are not backward compatible: codes that depend on the data()-loadable examples may be broken.
+
0.41-8
o A demonstration of 'pompBuilder' has been put into the 'advanced topics' vignette.
Modified: pkg/pomp/inst/data-R/sir.R
===================================================================
--- pkg/pomp/inst/data-R/sir.R 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/inst/data-R/sir.R 2012-04-26 14:41:38 UTC (rev 683)
@@ -22,7 +22,6 @@
paramnames=c(
"gamma","mu","iota",
"beta1","beta.sd","pop","rho",
- "nbasis","degree","period",
"S.0","I.0","R.0"
),
zeronames=c("cases"),
@@ -30,6 +29,9 @@
ic.names=c("S.0","I.0","R.0"),
parameter.transform="_sir_par_trans",
parameter.inv.transform="_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]
@@ -42,7 +44,6 @@
coef(po) <- c(
gamma=26,mu=0.02,iota=0.01,
- nbasis=3,degree=3,period=1,
beta1=400,beta2=480,beta3=320,
beta.sd=1e-3,
pop=2.1e6,
@@ -98,7 +99,6 @@
rho=0.1,
S.0=0.05,I.0=1e-4,R.0=0.95,
pop=1000000,
- nbasis=3,degree=3,period=1,
beta.sd=0
)
@@ -149,12 +149,14 @@
paramnames=c(
"gamma","mu","iota",
"beta","beta.sd","pop","rho",
- "nbasis","degree","period",
"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,
logvar=c(
"beta","gamma","mu","iota","sigma","beta.sd",
"S.0","I.0","R.0"
@@ -182,7 +184,6 @@
coef(po) <- c(
gamma=1/3,mu=0.0,iota=0.0,
- nbasis=1,degree=0,period=1,
beta=1.4,
beta.sd=0,
pop=1400,
Modified: pkg/pomp/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/pomp/inst/doc/advanced_topics_in_pomp.Rnw 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/inst/doc/advanced_topics_in_pomp.Rnw 2012-04-26 14:41:38 UTC (rev 683)
@@ -360,54 +360,59 @@
In the ``Introduction 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,
- PACKAGE="pomp"
- ),
- ## native routine for the skeleton:
- skeleton.type="vectorfield",
- skeleton="_sir_ODE",
- ## binomial measurement model:
- rmeasure="_sir_binom_rmeasure",
- dmeasure="_sir_binom_dmeasure",
- ## name of the shared-object library containing the
- ## native measurement-model 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"
- ),
- ## designate 'cases' as an accumulator variable
- ## i.e., set it to zero after each observation
- zeronames=c("cases"),
- comp.names=c("S","I","R"),
- parameter.transform="_sir_par_trans",
- parameter.inv.transform="_sir_par_untrans",
- initializer=function(params, t0, comp.names, ...) {
- ic.names <- paste(comp.names,".0",sep="")
- 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["cases"] <- 0
- x0
- }
- ) -> sir
+ 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,
+ PACKAGE="pomp"
+ ),
+ ## native routine for the skeleton:
+ skeleton.type="vectorfield",
+ skeleton="_sir_ODE",
+ ## native measurement-model routines:
+ rmeasure="_sir_binom_rmeasure",
+ dmeasure="_sir_binom_dmeasure",
+ ## name of the shared-object library containing the
+ 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",
+ "S.0","I.0","R.0"
+ ),
+ nbasis=3L, # three seasonal basis functions
+ degree=3L, # use cubic B-splines
+ period=1.0, # seasonality has period 1yr
+ ## designate 'cases' as an accumulator variable
+ ## i.e., set it to zero after each observation
+ zeronames=c("cases"),
+ ## parameter transformations in native routines:
+ parameter.transform="_sir_par_trans",
+ parameter.inv.transform="_sir_par_untrans",
+ ## some variables to be used in the initializer
+ comp.names=c("S","I","R"),
+ ic.names=c("S.0","I.0","R.0"),
+ ## parameterization of the initial conditions:
+ 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
+ }
+ ) -> 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
@@ -675,7 +680,6 @@
<<pomp-builder-eval,echo=F,eval=T,results=hide>>=
if (Sys.getenv("POMP_BUILD_VIGNETTES")=="yes") {
require(pomp)
- pompBuilder <- pomp:::pompBuilder
<<pomp-builder>>
simulate(sir) -> x
pfilter(sir,Np=500) -> pf
Modified: pkg/pomp/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/pomp/inst/doc/intro_to_pomp.Rnw 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/inst/doc/intro_to_pomp.Rnw 2012-04-26 14:41:38 UTC (rev 683)
@@ -570,15 +570,11 @@
A \code{pomp} object corresponding to the one just created (but with the \code{rprocess}, \code{rmeasure}, and \code{dmeasure} bits coded in C for speed) can be loaded by executing \verb+data(gompertz)+.
For your convenience, codes creating this \code{pomp} object are included with the package.
-To view them, do:
+To view them, run the demo
<<eval=F>>=
-file.show(system.file("demo/gompertz.R",package="pomp"))
-file.show(system.file("examples/gompertz.c",package="pomp"))
-@
-or run the demo
-<<eval=F>>=
demo(gompertz)
@
+
\clearpage
\section{Estimating parameters using iterated filtering: \code{mif}}
Modified: pkg/pomp/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Deleted: pkg/pomp/inst/examples/gompertz.c
===================================================================
--- pkg/pomp/inst/examples/gompertz.c 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/inst/examples/gompertz.c 2012-04-26 14:41:38 UTC (rev 683)
@@ -1,54 +0,0 @@
-// dear emacs, please treat this as -*- C++ -*-
-
-// Gompertz example as described in the "Introduction to pomp" vignette.
-// for a demonstration of how to compile, link, and run this example,
-// do 'demo("gompertz",package="pomp")' at the R prompt.
-
-#include <Rmath.h>
-
-#include "pomp.h" // in R, do 'system.file("include","pomp.h",package="pomp")' to find this header file
-
-// define some macros to make the code easier to read
-#define R (p[parindex[0]]) // growth rate
-#define K (p[parindex[1]]) // carrying capacity
-#define SIGMA (p[parindex[2]]) // process noise level
-#define TAU (p[parindex[3]]) // measurement noise level
-#define Y (y[obsindex[0]]) // observed population size
-#define X (x[stateindex[0]]) // actual population size
-#define XPRIME (f[stateindex[0]]) // new population size (for skeleton function only)
-
-// log-normal measurement model density
-void gompertz_normal_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) {
- *lik = dlnorm(Y,log(X),TAU,give_log);
-}
-
-// log-normal measurement model simulator
-void gompertz_normal_rmeasure (double *y, double *x, double *p,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
- Y = rlnorm(log(X),TAU);
-}
-
-// stochastic Gompertz model with log-normal process noise
-void gompertz_simulator (double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar,
- double t, double dt)
-{
- double S = exp(-R*dt);
- double sigma = SIGMA;
- double logeps = (sigma > 0.0) ? rnorm(0,sigma) : 0.0;
- X = pow(K,(1-S))*pow(X,S)*exp(logeps); // note X is over-written by this line
-}
-
-// the deterministic skeleton
-void gompertz_skeleton (double *f, double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar, double t)
-{
- double dt = 1.0;
- double S = exp(-R*dt);
- XPRIME = pow(K,1-S)*pow(X,S); // X is not over-written in the skeleton function
-}
Modified: pkg/pomp/inst/include/pomp.h
===================================================================
--- pkg/pomp/inst/include/pomp.h 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/inst/include/pomp.h 2012-04-26 14:41:38 UTC (rev 683)
@@ -7,6 +7,11 @@
#include <Rmath.h>
#include <Rdefines.h>
+// facility for extracting R objects from the 'userdata' slot
+const SEXP get_pomp_userdata (const char *name);
+const int *get_pomp_userdata_int (const char *name);
+const double *get_pomp_userdata_double (const char *name);
+
// facility for computing evaluating a basis of periodic bsplines
void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
Added: pkg/pomp/man/builder.Rd
===================================================================
--- pkg/pomp/man/builder.Rd (rev 0)
+++ pkg/pomp/man/builder.Rd 2012-04-26 14:41:38 UTC (rev 683)
@@ -0,0 +1,56 @@
+\name{pompBuilder}
+\alias{pompBuilder}
+\title{Write, compile, link, and build a pomp object using native codes}
+\description{
+ \code{pompBuilder} is an EXPERIMENTAL facility for producing compiled \code{pomp} objects.
+}
+\usage{
+pompBuilder(data, times, t0, name, statenames, paramnames, zeronames,
+ rmeasure, dmeasure, step.fn, step.fn.delta.t,
+ skeleton, skeleton.type, skelmap.delta.t = 1, link = TRUE)
+}
+\arguments{
+ \item{data, times, t0}{
+ The data, times, and zero-time.
+ See \code{\link{pomp}} for more information.
+ \code{data} must be a data-frame.
+ }
+ \item{name}{
+ character; the stem of the name for the files that will be produced.
+ }
+ \item{statenames, paramnames, zeronames}{
+ names of state-variables, parameters, and accumulator variables, respectively
+ }
+ \item{rmeasure, dmeasure}{
+ C codes implementing the measurement model
+ }
+ \item{step.fn, step.fn.delta.t}{
+ \code{step.fn} is a C code that implements an Euler step function.
+ The Euler time-step is \code{step.fn.delta.t}, which should be a positive number.
+ }
+ \item{skeleton, skeleton.type, skelmap.delta.t}{
+ \code{skeleton} is a C code that implements the deterministic skeleton.
+ As in \code{pomp}, \code{skeleton.type} indicates whether the skeleton is a map (discrete-time) or vectorfield (continuous-time).
+ If the former, \code{skelmap.delta.t} is the time-step of the map.
+ }
+ \item{link}{
+ logical; if TRUE, the resulting code will be linked after compilation.
+ }
+}
+\value{
+ The constructed \code{pomp} object.
+ A side-effect is the writing and compilation of a C code into a shared-object library.
+ These files will reside in the current working directory (see \code{\link{getwd}}).
+ If \code{pompBuilder} has been called with \code{link=FALSE}, this shared-object library must be linked (see \code{\link{dyn.load}}) before the \code{pomp} object can be used.
+}
+\details{
+ \code{pompBuilder} assumes that files can be written to the current working directory and that shared-object libraries can be compiled and linked, i.e., that \code{R CMD SHLIB} will work.
+ This will not typically be the case in out-of-the-box Windows installations.
+}
+\seealso{
+ \code{\link{pomp}} and the demos.
+}
+\author{Aaron A. King \email{kingaa at umich dot edu}}
+\examples{
+ ## see the demos
+}
Modified: pkg/pomp/src/R_init_pomp.c
===================================================================
--- pkg/pomp/src/R_init_pomp.c 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/src/R_init_pomp.c 2012-04-26 14:41:38 UTC (rev 683)
@@ -7,4 +7,5 @@
R_RegisterCCallable("pomp","dot_product",(DL_FUNC) &dot_product);
R_RegisterCCallable("pomp","reulermultinom",(DL_FUNC) &reulermultinom);
R_RegisterCCallable("pomp","deulermultinom",(DL_FUNC) &deulermultinom);
+ R_RegisterCCallable("pomp","get_pomp_userdata",(DL_FUNC) &get_pomp_userdata);
}
Modified: pkg/pomp/src/SSA_wrapper.c
===================================================================
--- pkg/pomp/src/SSA_wrapper.c 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/src/SSA_wrapper.c 2012-04-26 14:41:38 UTC (rev 683)
@@ -172,13 +172,23 @@
zidx = 0;
}
+ if (use_native) {
+ set_pomp_userdata(args);
+ }
+
GetRNGstate();
F77_CALL(driverssa)(F77_SUB(reactionrate),&nvar,&nevent,&npar,&nrep,&ntimes,
INTEGER(mflag),REAL(xstart),REAL(times),REAL(params),
REAL(X),REAL(e),REAL(vmatrix),REAL(dmatrix),
&nzeros,zidx,sidx,pidx,&ncovars,cidx,
&covlen,&covdim,REAL(tcovar),REAL(covar));
+
PutRNGstate();
+
+ if (use_native) {
+ unset_pomp_userdata();
+ }
+
UNPROTECT(nprotect);
return X;
}
Modified: pkg/pomp/src/dmeasure.c
===================================================================
--- pkg/pomp/src/dmeasure.c 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/src/dmeasure.c 2012-04-26 14:41:38 UTC (rev 683)
@@ -245,9 +245,18 @@
ndim[0] = nvars; ndim[1] = npars; ndim[2] = nrepsp; ndim[3] = nrepsx; ndim[4] = ntimes;
ndim[5] = covlen; ndim[6] = covdim; ndim[7] = nobs;
+ if (use_native) {
+ PROTECT(FCALL = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
+ set_pomp_userdata(FCALL);
+ }
+
dens_meas(ff,REAL(F),REAL(y),REAL(x),REAL(times),REAL(params),INTEGER(log),ndim,
oidx,sidx,pidx,cidx,REAL(tcovar),REAL(covar));
+ if (use_native) {
+ unset_pomp_userdata();
+ }
+
UNPROTECT(nprotect);
return F;
}
Modified: pkg/pomp/src/euler.c
===================================================================
--- pkg/pomp/src/euler.c 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/src/euler.c 2012-04-26 14:41:38 UTC (rev 683)
@@ -415,7 +415,10 @@
ndim[0] = nvar; ndim[1] = npar; ndim[2] = nrep; ndim[3] = ntimes;
ndim[4] = covlen; ndim[5] = covdim; ndim[6] = nzeros;
- if (use_native) GetRNGstate();
+ if (use_native) {
+ set_pomp_userdata(args);
+ GetRNGstate();
+ }
switch (meth) {
case 0:
@@ -435,7 +438,10 @@
break;
}
- if (use_native) PutRNGstate();
+ if (use_native) {
+ PutRNGstate();
+ unset_pomp_userdata();
+ }
VINDEX = 0;
@@ -647,10 +653,14 @@
ndim[0] = nvar; ndim[1] = npar; ndim[2] = nrep; ndim[3] = ntimes;
ndim[4] = covlen; ndim[5] = covdim;
+ if (use_native) set_pomp_userdata(args);
+
euler_densities(ff,REAL(F),REAL(x),REAL(times),REAL(params),
ndim,sidx,pidx,cidx,
REAL(tcovar),REAL(covar),INTEGER(log));
+ if (use_native) unset_pomp_userdata();
+
UNPROTECT(nprotect);
return F;
}
Modified: pkg/pomp/src/partrans.c
===================================================================
--- pkg/pomp/src/partrans.c 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/src/partrans.c 2012-04-26 14:41:38 UTC (rev 683)
@@ -89,7 +89,7 @@
SEXP do_partrans (SEXP object, SEXP params, SEXP fun)
{
int nprotect = 0;
- SEXP ptrans, fn;
+ SEXP ptrans, fn, userdata;
int use, drop;
PROTECT(fn = VECTOR_ELT(fun,0)); nprotect++;
@@ -130,6 +130,9 @@
idx = 0;
}
+ PROTECT(userdata = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
+ set_pomp_userdata(userdata);
+
PROTECT(ptrans = duplicate(params)); nprotect++;
for (k = 0, ps = REAL(params), pt = REAL(ptrans); k < nrep; k++, ps += npar, pt += npar) {
@@ -137,6 +140,8 @@
(*ff)(pt,ps,idx);
}
+ unset_pomp_userdata();
+
}
break;
default:
Modified: pkg/pomp/src/pomp.h
===================================================================
--- pkg/pomp/src/pomp.h 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/src/pomp.h 2012-04-26 14:41:38 UTC (rev 683)
@@ -7,6 +7,11 @@
#include <Rmath.h>
#include <Rdefines.h>
+// facility for extracting R objects from the 'userdata' slot
+const SEXP get_pomp_userdata (const char *name);
+const int *get_pomp_userdata_int (const char *name);
+const double *get_pomp_userdata_double (const char *name);
+
// facility for computing evaluating a basis of periodic bsplines
void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
Modified: pkg/pomp/src/pomp_internal.h
===================================================================
--- pkg/pomp/src/pomp_internal.h 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/src/pomp_internal.h 2012-04-26 14:41:38 UTC (rev 683)
@@ -59,6 +59,9 @@
// skeleton.c
SEXP do_skeleton (SEXP object, SEXP x, SEXP t, SEXP params, SEXP fun);
+//userdata.c
+void set_pomp_userdata (SEXP object);
+void unset_pomp_userdata (void);
static R_INLINE SEXP makearray (int rank, int *dim) {
int nprotect = 0;
@@ -219,6 +222,27 @@
return log(x/(1-x));
}
+static R_INLINE SEXP getListElement (SEXP list, const char *str)
+{
+ SEXP elmt = R_NilValue;
+ SEXP names = getAttrib(list,R_NamesSymbol);
+ for (R_len_t i = 0; i < length(list); i++)
+ if(strcmp(CHAR(STRING_ELT(names,i)),str) == 0) {
+ elmt = VECTOR_ELT(list,i);
+ break;
+ }
+ return elmt;
+}
+
+static R_INLINE SEXP getPairListElement (SEXP list, const char *name)
+{
+ while (list != R_NilValue) {
+ if (strcmp(STRING_PTR(PRINTNAME(TAG(list))),name)==0) break;
+ list = CDR(list);
+ }
+ return CAR(list);
+}
+
#ifdef __cplusplus
template <class Scalar>
Modified: pkg/pomp/src/rmeasure.c
===================================================================
--- pkg/pomp/src/rmeasure.c 2012-04-25 21:53:47 UTC (rev 682)
+++ pkg/pomp/src/rmeasure.c 2012-04-26 14:41:38 UTC (rev 683)
@@ -261,13 +261,20 @@
ndim[0] = nvars; ndim[1] = npars; ndim[2] = nrepsp; ndim[3] = nrepsx; ndim[4] = ntimes;
ndim[5] = covlen; ndim[6] = covdim; ndim[7] = nobs;
- if (use_native) GetRNGstate();
+ if (use_native) {
+ PROTECT(FCALL = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
+ set_pomp_userdata(FCALL);
+ GetRNGstate();
+ }
simul_meas(ff,REAL(Y),REAL(x),REAL(times),REAL(params),ndim,
oidx,sidx,pidx,cidx,
REAL(tcovar),REAL(covar));
- if (use_native) PutRNGstate();
+ if (use_native) {
+ PutRNGstate();
+ unset_pomp_userdata();
+ }
OIDX = 0;
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 683
More information about the pomp-commits
mailing list