[Pomp-commits] r484 - in pkg/inst: data-R doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 19 21:22:36 CEST 2011
Author: kingaa
Date: 2011-05-19 21:22:36 +0200 (Thu, 19 May 2011)
New Revision: 484
Modified:
pkg/inst/data-R/ou2.R
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
Log:
- improve the 'advanced topics' vignette
- update the data()-loadable ou2 object
Modified: pkg/inst/data-R/ou2.R
===================================================================
--- pkg/inst/data-R/ou2.R 2011-05-19 19:21:29 UTC (rev 483)
+++ pkg/inst/data-R/ou2.R 2011-05-19 19:22:36 UTC (rev 484)
@@ -39,7 +39,7 @@
ntimes <- length(times)
parindex <- match(paramnames,rownames(params))-1
array(
- .C("_ou2_pdf",
+ .C("ou2_pdf",
d = double(nrep*(ntimes-1)),
X = as.double(x),
par = as.double(params),
@@ -54,8 +54,8 @@
dim=c(nrep,ntimes-1)
)
},
- dmeasure = "_ou2_normal_dmeasure",
- rmeasure = "_ou2_normal_rmeasure",
+ dmeasure = "ou2_normal_dmeasure",
+ rmeasure = "ou2_normal_rmeasure",
skeleton.type="map",
skeleton = function (x, t, params, ...) {
with(
Modified: pkg/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw 2011-05-19 19:21:29 UTC (rev 483)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2011-05-19 19:22:36 UTC (rev 484)
@@ -21,12 +21,12 @@
\newfloat{textbox}{t}{tbx}
\floatname{textbox}{Box}
-\newcommand\code[1]{\texttt{#1}}
+\newcommand{\code}[1]{\texttt{#1}}
+\newcommand{\pkg}[1]{\textsf{#1}}
\newcommand{\R}{\textsf{R}}
-\newcommand{\pomp}{\texttt{pomp}}
\newcommand{\expect}[1]{\mathbb{E}\left[#1\right]}
-\title[Advanced topics in \code{pomp}]{Advanced topics in \code{pomp}}
+\title[Advanced topics in pomp]{Advanced topics in pomp}
\author[A. A. King]{Aaron A. King}
@@ -42,7 +42,7 @@
\urladdr{http:pomp.r-forge.r-project.org}
-\date{\today, \pomp~version~\Sexpr{packageDescription("pomp",fields="Version")}}
+\date{\today, \pkg{pomp}~version~\Sexpr{packageDescription("pomp",fields="Version")}}
\SweaveOpts{echo=T,results=verbatim,print=F,eps=F,pdf=T,keep.source=T}
@@ -61,17 +61,249 @@
set.seed(5384959)
@
-This document gives some examples of the use of native (C or FORTRAN) codes in \code{pomp} and introduces the low-level interface to \code{pomp} objects.
+This document discusses accelerating \pkg{pomp} by vectorizing your codes and/or using native (C or FORTRAN) codes.
+It also introduces \pkg{pomp}'s low-level interface for code developers.
-\section{Acceleration using native codes: using plug-ins with native code}
+\section{Accelerating your codes: vectorizing \code{rprocess} and using native codes}
-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.
-\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 ``Introduction to pomp'' vignette, we used \emph{plug-ins} provided by the package to specify the \code{rprocess} component of partially-observed Markov process models.
+The \code{rprocess} plug-ins require you to write a simulator for a single realization of the process, for a single set of parameters, from one time to another.
+\pkg{pomp} then calls this code many times---using potentially many different parameter values, states, and times---whenever it simulates the process, computes likelihood via Monte Carlo integration, etc.
+The inference methods implemented in \pkg{pomp} are quite computationally intensive, which puts a premium on the speed of your codes.
+Sometimes, you can realize substantial speed-up of your code by vectorizing it.
+This necessitates foregoing the relative simplicity of the plug-in-based implementation and writing \code{rprocess} ``from scratch''.
+Here, we'll develop a vectorized version of \code{rprocess} in \R\ code, then we'll see what the same thing looks like coded in C.
+We'll compare these different versions in terms of their speed at simulation.
-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.
+We'll use a discrete-time bivariate AR(1) process with normal measurement error as our example.
+In this model, the state process $X_{t}\in\mathbb{R}^2$ satisfies
+\begin{equation}\label{eq:ou-proc}
+ X_{t} = \alpha\,X_{t-1}+\sigma\,\varepsilon_{t}.
+\end{equation}
+The measurement process is
+\begin{equation}\label{eq:ou-obs}
+ Y_{t} = \beta\,X_{t}+\tau\,\xi_{t}.
+\end{equation}
+In these equations, $\alpha$ and and $\beta$ are $2\times 2$ constant matrices.
+$\xi_{t}$ and $\varepsilon_{t}$ are mutually-independent families of i.i.d.\ bivariate standard normal random variables.
+$\sigma$ is a lower-triangular matrix such that $\sigma\sigma^T$ is the variance-covariance matrix of $X_{t+1}\vert X_{t}$.
+We'll assume that each component of $X$ is measured independently and with the same error, $\tau$, so that the variance-covariance matrix of $Y_{t}\vert X_{t}$ has $\tau^2$ on the diagonal and zeros elsewhere.
+
+An implementation of this model is included in the package as a \code{pomp} object;
+load it by executing \code{data(ou2)}.
+
+\subsection{An unvectorized implementation using \R\ code only}
+Before we set about vectorizing the codes, let's have a look at what a plug-in based implementation written entirely in \R\ might look like.
+<<plugin-R-code,echo=T>>=
+data(ou2)
+ou2.dat <- as.data.frame(ou2)
+
+pomp(
+ data=ou2.dat[c("time","y1","y2")],
+ times="time",
+ t0=0,
+ rprocess=discrete.time.sim(
+ step.fun=function (x, t, params, ...) {
+ eps <- rnorm(n=2,mean=0,sd=1) # noise terms
+ xnew <- c(
+ x1=params["alpha.1"]*x["x1"]+params["alpha.3"]*x["x2"]+
+ params["sigma.1"]*eps[1],
+ x2=params["alpha.2"]*x["x1"]+params["alpha.4"]*x["x2"]+
+ params["sigma.2"]*eps[1]+params["sigma.3"]*eps[2]
+ )
+ names(xnew) <- c("x1","x2")
+ xnew
+ }
+ )
+ ) -> ou2.Rplug
+@
+<<plugin-R-code-sim,echo=F>>=
+tic <- Sys.time()
+simdat.Rplug <- simulate(ou2.Rplug,params=coef(ou2),nsim=1000,states=T)
+toc <- Sys.time()
+etime.Rplug <- toc-tic
+@
+Notice how we specify the process model simulator using the \code{rprocess} plug-in \code{discrete.time.sim}.
+The latter function's \code{step.fun} argument is itself a function that simulates one realization of the process for one timestep and one set of parameters.
+When we vectorize the code, we'll do many realizations at once.
+
+\subsection{Vectorizing the process simulator using \R\ code only}
+Now, to write a vectorized \code{rprocess} in \R, we must write a function that simulates \code{nrep} realizations of the unobserved process.
+Each of these realizations may start at a different point in state space and each may have a different set of parameters.
+Moreover, this function must be capable of simulating the process over an arbitrary time interval and must be capable of reporting the unobserved states at arbitrary times in that interval.
+We'll accomplish this by writing an \R\ function with arguments \code{xstart}, \code{params}, and \code{times}.
+About these inputs, we must assume:
+\begin{enumerate}
+\item \code{xstart} will be a matrix, each column of which is a vector of initial values of the state process.
+ Each state variable (matrix row) will be named.
+\item \code{params} will be a matrix, the columns of which are parameter vectors.
+ The parameter names will be in the matrix column-names.
+\item \code{times} will be a vector of times at which realizations of the state process are required.
+ We will have \code{times[k] <= times[k+1]} for all indices \code{k}, but we cannot assume that the entries of \code{times} will be unique.
+\item The initial states \code{xstart} are assumed to obtain at time \code{times[1]}.
+\end{enumerate}
+This function must return a rank-3 array, which has the realized values of the state process at the requested times.
+This array must have \code{x} rownames.
+Here is one implementation of such a simulator.
+<<vectorized-R-code>>=
+ou2.Rvect.rprocess <- function (xstart, times, params, ...) {
+ nrep <- ncol(xstart) # number of realizations
+ ntimes <- length(times) # number of timepoints
+ ## unpack the parameters (for legibility only)
+ alpha.1 <- params["alpha.1",]
+ alpha.2 <- params["alpha.2",]
+ alpha.3 <- params["alpha.3",]
+ alpha.4 <- params["alpha.4",]
+ sigma.1 <- params["sigma.1",]
+ sigma.2 <- params["sigma.2",]
+ sigma.3 <- params["sigma.3",]
+ ## x is the array of states to be returned: it must have rownames
+ x <- array(0,dim=c(2,nrep,ntimes))
+ rownames(x) <- rownames(xstart)
+ ## xnow holds the current state values
+ x[,,1] <- xnow <- xstart
+ tnow <- times[1]
+ for (k in seq.int(from=2,to=ntimes,by=1)) {
+ tgoal <- times[k]
+ while (tnow < tgoal) { # take one step at a time
+ eps <- array(rnorm(n=2*nrep,mean=0,sd=1),dim=c(2,nrep))
+ tmp <- alpha.1*xnow['x1',]+alpha.3*xnow['x2',]+
+ sigma.1*eps[1,]
+ xnow['x2',] <- alpha.2*xnow['x1',]+alpha.4*xnow['x2',]+
+ sigma.2*eps[1,]+sigma.3*eps[2,]
+ xnow['x1',] <- tmp
+ tnow <- tnow+1
+ }
+ x[,,k] <- xnow
+ }
+ x
+}
+@
+We can put this into a pomp object that is the same as \code{ou2.Rplug} in every way except in its \code{rprocess} slot by doing
+<<vectorized-R-pomp>>=
+ou2.Rvect <- pomp(ou2.Rplug,rprocess=ou2.Rvect.rprocess)
+@
+
+Let's pick some parameters and simulate some data to see how long it takes this code to run.
+<<theta>>=
+theta <- c(
+ x1.0=-3, x2.0=4,
+ tau=1,
+ alpha.1=0.8, alpha.2=-0.5, alpha.3=0.3, alpha.4=0.9,
+ sigma.1=3, sigma.2=-0.5, sigma.3=2
+ )
+@
+<<vectorized-R-code-sim>>=
+tic <- Sys.time()
+simdat.Rvect <- simulate(ou2.Rvect,params=theta,states=T,nsim=1000)
+toc <- Sys.time()
+etime.Rvect <- toc-tic
+@
+<<echo=F>>=
+units(etime.Rvect) <- units(etime.Rplug)
+@
+Doing \Sexpr{dim(simdat.Rvect)[2]} simulations of \code{ou2.Rvect} took \Sexpr{round(etime.Rvect,2)}~\Sexpr{units(etime.Rvect)}.
+Compared to the \Sexpr{round(etime.Rplug,2)}~\Sexpr{units(etime.Rplug)} it took to run \Sexpr{dim(simdat.Rplug)[2]} simulations of \code{ou2.Rplug}, this is a \Sexpr{round(as.numeric(etime.Rplug)/as.numeric(etime.Rvect))}-fold speed-up.
+
+\subsection{Accelerating the code using C: a plug-in based implementation}
+
+As we've seen, we can usually acheive big accelerations using compiled native code.
+A one-step simulator written in C for use with the \code{discrete.time.sim} plug-in is included with the package and can be viewed by doing
+<<view-ou2-source,eval=F>>=
+file.show(file=system.file("examples/ou2.c",package="pomp"))
+@
+The one-step simulator is in function \code{ou2\_step}.
+We can put this into the \code{pomp} object and simulate as before by doing
+<<plugin-C-code>>=
+ou2.Cplug <- pomp(
+ ou2.Rplug,
+ rprocess=discrete.time.sim("ou2_step"),
+ paramnames=c(
+ "alpha.1","alpha.2","alpha.3","alpha.4",
+ "sigma.1","sigma.2","sigma.3",
+ "tau"
+ ),
+ statenames=c("x1","x2"),
+ obsnames=c("y1","y2")
+ )
+@
+<<plugin-C-sim>>=
+tic <- Sys.time()
+simdat.Cplug <- simulate(ou2.Cplug,params=theta,states=T,nsim=100000)
+toc <- Sys.time()
+etime.Cplug <- toc-tic
+@
+Note that \code{ou2\_step} is written in such a way that we must specify \code{paramnames}, \code{statenames}, and \code{obsnames}.
+These \Sexpr{dim(simdat.Cplug)[2]} simulations of \code{ou2.Cplug} took \Sexpr{round(etime.Cplug,2)}~\Sexpr{units(etime.Cplug)}.
+<<echo=F>>=
+units(etime.Cplug) <- units(etime.Rplug)
+speedup <- as.numeric(etime.Rplug/dim(simdat.Rplug)[2])/as.numeric(etime.Cplug/dim(simdat.Cplug)[2])
+@
+This is a \Sexpr{round(speedup)}-fold speed-up relative to \code{ou2.Rplug}.
+
+
+\subsection{A vectorized C implementation}
+
+The function \code{ou2\_adv} is a fully vectorized version of the simulator written in C.
+View this code by doing
+<<eval=F>>=
+<<view-ou2-source>>
+@
+This function is called in the following \code{rprocess} function.
+Notice that the call to \code{ou2\_adv} uses the \verb+.C+ convention.
+<<vectorized-C-code>>=
+ou2.Cvect.rprocess <- function (xstart, times, params, ...) {
+ nvar <- nrow(xstart)
+ npar <- nrow(params)
+ nrep <- ncol(xstart)
+ ntimes <- length(times)
+ array(
+ .C("ou2_adv",
+ X = double(nvar*nrep*ntimes),
+ xstart = as.double(xstart),
+ par = as.double(params),
+ times = as.double(times),
+ n = as.integer(c(nvar,npar,nrep,ntimes)),
+ DUP = FALSE,
+ NAOK = TRUE
+ )$X,
+ dim=c(nvar,nrep,ntimes),
+ dimnames=list(rownames(xstart),NULL,NULL)
+ )
+}
+@
+
+The call that constructs the \code{pomp} object is:
+<<vectorized-C-code-pomp>>=
+ou2.Cvect <- pomp(
+ ou2.Rplug,
+ rprocess=ou2.Cvect.rprocess
+ )
+@
+<<vectorized-C-code-sim,echo=F>>=
+tic <- Sys.time()
+paramnames <- c(
+ "alpha.1","alpha.2","alpha.3","alpha.4",
+ "sigma.1","sigma.2","sigma.3",
+ "tau",
+ "x1.0","x2.0"
+ )
+simdat.Cvect <- simulate(ou2.Cvect,params=theta[paramnames],nsim=100000,states=T)
+toc <- Sys.time()
+etime.Cvect <- toc-tic
+@
+Note that we've had to rearrange the order of parameters here to ensure that they arrive at the native codes in the right order.
+<<echo=F>>=
+units(etime.Cvect) <- units(etime.Rplug)
+speedup <- as.numeric(etime.Rplug/dim(simdat.Rplug)[2])/as.numeric(etime.Cvect/dim(simdat.Cvect)[2])
+@
+Doing \Sexpr{dim(simdat.Cvect)[2]} simulations of \code{ou2.Cvect} took \Sexpr{round(etime.Cvect,2)}~\Sexpr{units(etime.Cvect)}, a \Sexpr{round(speedup)}-fold speed-up relative to \code{ou2.Rplug}.
+
+
+\subsection{More on native codes and plug-ins}
+
+It's possible to use native codes for \code{dprocess} and for the measurement model portions of the \code{pomp} as well.
+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(
@@ -84,7 +316,8 @@
## native routine for the process simulator:
rprocess=euler.sim(
step.fun="_sir_euler_simulator",
- delta.t=1/52/20
+ delta.t=1/52/20,
+ PACKAGE="pomp"
),
## native routine for the skeleton:
skeleton.type="vectorfield",
@@ -92,7 +325,8 @@
## binomial measurement model:
rmeasure="_sir_binom_rmeasure",
dmeasure="_sir_binom_dmeasure",
- ## name of the shared-object library containing the native routines:
+ ## 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"),
@@ -125,12 +359,25 @@
file.show(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.
+Note that the native routines for this model are included in the package, which is why we give the \verb+PACKAGE="pomp"+ argument to \code{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:
+You can also use the R package \pkg{inline} to put C or FORTRAN codes directly into your R functions.
+
+There is an important issue that arises when using native codes.
+This has to do with the order in which parameters, states, and observables are passed to these codes.
+\pkg{pomp} relies on the names (also row-names and column-names) attributes to identify variables in vectors and arrays.
+When you write a C or FORTRAN version of \code{rprocess} or \code{dmeasure} for example, you write a routine that takes parameters, state variables, and/or observables in the form of a vector.
+However, you have no control over the order in which these are given to you.
+Without some means of knowing which element of each vector corresponds to which variable, you cannot write the codes correctly.
+This is where the \code{paramnames}, \code{statenames}, \code{covarnames}, and \code{obsnames} arguments to \code{pomp} come in: use these arguments to specify the order in which your C code expects to see the parameters, state variables, covariates, and observables (data variables).
+\code{pomp} will match these names against the corresponding names attributes of vectors.
+It will then pass to your native routines index vectors you can use to locate the correct variables.
+See the source code to see how this is done.
+
+Let's specify some parameters, simulate, and compute a deterministic trajectory:
<<sir-sim>>=
params <- c(
gamma=26,mu=0.02,iota=0.01,
@@ -142,143 +389,11 @@
)
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()
+sims <- simulate(sir,nsim=10,obs=T)
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)
-data(ou2)
-show(ou2)
-@
-Here we'll examine how this object is put together.
-The process model simulator and density functions are as follows:
-<<ou2-rprocess>>=
- ou2.rprocess <- function (xstart, times, params, paramnames, ...) {
- nvar <- nrow(xstart)
- npar <- nrow(params)
- nrep <- ncol(xstart)
- ntimes <- length(times)
- ## get indices of the various parameters in the 'params' matrix
- ## C uses zero-based indexing!
- parindex <- match(paramnames,rownames(params))-1
- array(
- .C("_ou2_adv",
- X = double(nvar*nrep*ntimes),
- xstart = as.double(xstart),
- par = as.double(params),
- times = as.double(times),
- n = as.integer(c(nvar,npar,nrep,ntimes)),
- parindex = as.integer(parindex),
- DUP = FALSE,
- NAOK = TRUE,
- PACKAGE = "pomp"
- )$X,
- dim=c(nvar,nrep,ntimes),
- dimnames=list(rownames(xstart),NULL,NULL)
- )
- }
-@
-
-<<ou2-dprocess>>=
- ou2.dprocess <- function (x, times, params, log, paramnames, ...) {
- nvar <- nrow(x)
- npar <- nrow(params)
- nrep <- ncol(x)
- ntimes <- length(times)
- parindex <- match(paramnames,rownames(params))-1
- array(
- .C("_ou2_pdf",
- d = double(nrep*(ntimes-1)),
- X = as.double(x),
- par = as.double(params),
- times = as.double(times),
- n = as.integer(c(nvar,npar,nrep,ntimes)),
- parindex = as.integer(parindex),
- give_log=as.integer(log),
- DUP = FALSE,
- NAOK = TRUE,
- PACKAGE = "pomp"
- )$d,
- dim=c(nrep,ntimes-1)
- )
- }
-@
-
-The call that constructs the \code{pomp} object is:
-<<>>=
-ou2 <- pomp(
- data=data.frame(
- time=seq(1,100),
- y1=NA,
- y2=NA
- ),
- times="time",
- t0=0,
- rprocess = ou2.rprocess,
- dprocess = ou2.dprocess,
- dmeasure = "_ou2_normal_dmeasure",
- rmeasure = "_ou2_normal_rmeasure",
- paramnames=c(
- "alpha.1","alpha.2","alpha.3","alpha.4",
- "sigma.1","sigma.2","sigma.3",
- "tau"
- ),
- statenames = c("x1","x2"),
- obsnames = c("y1","y2"),
- PACKAGE="pomp"
- )
-@
-Notice that the process model is implemented using using \verb+.C+, while the measurement model is specified by giving the names of native C routines.
-Read the source to see the definitions of these functions.
-For convenience, the source codes are provided with the package in the \code{examples} directory.
-Do
-<<view-ou2-source,eval=F>>=
-file.show(file=system.file("examples/ou2.c",package="pomp"))
-@
-to view the source code.
-
-There is an important issue that arises when using native codes.
-This has to do with the order in which parameters, states, and observables are passed into the native codes.
-\pomp\ relies on the names (also row-names and column-names) attributes to identify variables in vectors and arrays.
-When you write a C or FORTRAN version of \code{rprocess} or \code{dmeasure} for example, you write a routine that takes parameters, state variables, and/or observables in the form of a vector.
-However, you have no control over the order in which these are given to you.
-Without some means of knowing which element of each vector corresponds to which variable, you cannot write the codes correctly.
-This is where the \code{paramnames}, \code{statenames}, and \code{obsnames} arguments to \code{pomp} come in.
-When you specifying the names of parameters, state variables, and observables (data variables) here, \code{pomp} matches these names against the corresponding names attributes of vectors and passes to your native routine integer vectors which you can use to identify the correct variables.
-See the source code to see how this is done.
-
-We'll specify some parameters:
-<<>>=
-theta <- c(
- alpha.1=0.8, alpha.2=-0.5, alpha.3=0.3, alpha.4=0.9,
- sigma.1=3, sigma.2=-0.5, sigma.3=2,
- tau=1,
- x1.0=-3, x2.0=4
- )
-@
-
-<<>>=
-tic <- Sys.time()
-x <- simulate(ou2,params=theta,nsim=500,seed=80073088L)
-toc <- Sys.time()
-print(toc-tic)
-@
-
+\clearpage
\section{The low-level interface}
There is a low-level interface to \code{pomp} objects, primarily designed for package developers.
@@ -340,7 +455,7 @@
\section{Other examples}
-There are a number of example \pomp\ objects included with the package.
+There are a number of example \code{pomp} objects included with the package.
These can be found by running
<<all-data-loaable,eval=F>>=
data(package="pomp")
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
More information about the pomp-commits
mailing list