change 'rprocess' (including interface) to eliminate 'offset' and unn…

…ecessary copy operation

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: pomp2
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
-Version: 2.0.10.0
-Date: 2019-02-28
+Version: 2.0.11.0
+Date: 2019-03-05
 Authors@R: c(person(given=c("Aaron","A."),family="King",
 		role=c("aut","cre"),email="[email protected]"),
 	  person(given=c("Edward","L."),family="Ionides",role=c("aut")),

Makefile

@@ -15,8 +15,9 @@ INSTALL = install
 PKG = $(shell perl -ne 'print $$1 if /Package:\s+((\w+[-\.]?)+)/;' DESCRIPTION)
 VERSION = $(shell perl -ne 'print $$1 if /Version:\s+((\d+[-\.]?)+)/;' DESCRIPTION)
 PKGVERS = $(PKG)_$(VERSION)
-SOURCE=$(shell ls R/*R src/*.c src/*.h data/* tests/*R)
+SOURCE=$(shell ls R/*R src/*.c src/*.h data/*)
 CSOURCE=$(shell ls src/*.c)
+TESTS=$(shell ls tests/*R)
 
 default:
 	@echo $(PKGVERS)
@@ -63,7 +64,7 @@ roxy: $(SOURCE)
 
 dist: NEWS $(PKGVERS).tar.gz
 
-$(PKGVERS).tar.gz: $(SOURCE) includes headers
+$(PKGVERS).tar.gz: $(SOURCE) $(TESTS) includes headers
 	$(RCMD) build --force --no-manual --resave-data --compact-vignettes=both --md5 .
 
 src/pomp_decls.h: $(CSOURCE)
@@ -116,7 +117,7 @@ $(PKG).pdf: $(SOURCE)
 	$(RCMD) Rd2pdf --no-preview --pdf --force -o $(PKG).pdf .
 	$(RSCRIPT) -e "tools::compactPDF(\"$(PKG).pdf\")";
 
-tests: install tests/*.R
+tests: install $(TESTS)
 	export R_LIBS
 	$(MAKE) -C tests
 

R/bsmc2.R

@@ -274,8 +274,8 @@ bsmc2.internal <- function (object, Np, smooth, tol, max.fail,
 
     tparams <- partrans(object,params,dir="fromEst",.gnsi=gnsi)
 
-    xpred <- rprocess(object,xstart=x,times=times[c(nt,nt+1)],params=tparams,
-      offset=1L,.gnsi=gnsi)
+    xpred <- rprocess(object,x0=x,t0=times[nt],times=times[nt+1],
+      params=tparams,.gnsi=gnsi)
 
     ## evaluate likelihood of observation given xpred (from L&W AGM (4))
     weights <- dmeasure(object,y=object@data[,nt,drop=FALSE],x=xpred,

R/flow.R

@@ -63,19 +63,19 @@ setMethod(
 setMethod(
   "flow",
   signature=signature(object="pomp"),
-  definition=function (object, x0, t0, params, times, ...,
+  definition=function (object, x0, t0, times, params, ...,
     verbose = getOption("verbose", FALSE)) {
 
     tryCatch(
-      flow.internal(object=object,x0=x0,t0=t0,params=params,times=times,
+      flow.internal(object=object,x0=x0,t0=t0,times=times,params=params,
         ...,verbose=verbose),
       error = function (e) pStop("flow",conditionMessage(e))
     )
 
   }
 )
 
-flow.internal <- function (object, x0, t0, params, times, ...,
+flow.internal <- function (object, x0, t0, times, params, ...,
   .gnsi = TRUE, verbose) {
 
   verbose <- as.logical(verbose)

R/kalman.R

@@ -281,7 +281,7 @@ enkf.internal <- function (object,
 
   for (k in seq_len(ntimes)) {
     ## advance ensemble according to state process
-    X[,] <- rprocess(object,params=params,xstart=X,times=tt[c(k,k+1)],offset=1)
+    X[,] <- rprocess(object,x0=X,t0=tt[k],times=tt[k+1],params=params)
 
     predMeans[,k] <- pm <- rowMeans(X) # prediction mean
     Y[,] <- apply(X,2,h)               # ensemble of forecasts
@@ -362,7 +362,7 @@ eakf.internal <- function (object,
   for (k in seq_len(ntimes)) {
 
     ## advance ensemble according to state process
-    X[,] <- rprocess(object,xstart=X,params=params,times=tt[c(k,k+1)],offset=1)
+    X[,] <- rprocess(object,x0=X,t0=tt[k],times=tt[k+1],params=params)
 
     predMeans[,k] <- pm <- rowMeans(X) # prediction mean
     X <- X-pm

R/mif2.R

@@ -481,10 +481,10 @@ mif2.pfilter <- function (object, params, Np, mifiter, rw.sd, cooling.fn,
     ## advance the state variables according to the process model
     X <- rprocess(
       object,
-      xstart=x,
-      times=times[c(nt,nt+1)],
+      x0=x,
+      t0=times[nt],
+      times=times[nt+1],
       params=tparams,
-      offset=1,
       .gnsi=gnsi
     )
 

R/pfilter.R

@@ -335,8 +335,8 @@ pfilter.internal <- function (object, Np, tol, max.fail,
   for (nt in seq_len(ntimes)) { ## main loop
 
     ## advance the state variables according to the process model
-    X <- rprocess(object,xstart=x,times=times[c(nt,nt+1)],params=params,
-      offset=1L,.gnsi=gnsi)
+    X <- rprocess(object,x0=x,t0=times[nt],times=times[nt+1],params=params,
+      .gnsi=gnsi)
 
     if (pred.var) { ## check for nonfinite state variables and parameters
       problem.indices <- unique(which(!is.finite(X),arr.ind=TRUE)[,1L])

R/pomp_fun.R

@@ -196,17 +196,13 @@ setMethod(
 
       cat("native function\n  - name: ",sQuote(object@native.fun),"\n",sep="")
       if (length(object@PACKAGE)>0)
-        cat("  - dynamically loaded from: ",
-          sQuote(object@PACKAGE),sep="")
-      cat("\n")
+        cat("  - dynamically loaded from: ",sQuote(object@PACKAGE),sep="")
 
     } else if (mode==pompfunmode$regNative) { # built from C snippets
 
       cat("native function\n  - name: ",sQuote(object@native.fun),"\n",sep="")
       if (length(object@PACKAGE)>0)
-        cat("  - defined by a C snippet in library ",
-          sQuote(object@PACKAGE),sep="")
-      cat("\n")
+        cat("  - defined by a C snippet in library ",sQuote(object@PACKAGE),sep="")
 
     } else {
 

R/workhorses.R

@@ -499,8 +499,8 @@ rprior.internal <- function (object, params, .gnsi = TRUE, ...) {
 ##'
 ##' \code{rprocess} simulates the process-model portion of partially-observed Markov process.
 ##'
-##' When \code{rprocess} is called, the first entry of \code{times} is taken to be the initial time (i.e., that corresponding to \code{xstart}).
-##' Subsequent times are the additional times at which the state of the simulated processes are required.
+##' When \code{rprocess} is called, \code{t0} is taken to be the initial time (i.e., that corresponding to \code{x0}).
+##' The values in \code{times} are the times at which the state of the simulated processes are required.
 ##'
 ##' @name rprocess
 ##' @docType methods
@@ -511,26 +511,26 @@ rprior.internal <- function (object, params, .gnsi = TRUE, ...) {
 ##'
 ##' @inheritParams dmeasure
 ##'
-##' @param xstart an \code{nvar} x \code{nrep} matrix containing the starting state of the system.
-##' Columns of \code{xstart} correspond to states;
+##' @param x0 an \code{nvar} x \code{nrep} matrix containing the starting state of the system.
+##' Columns of \code{x0} correspond to states;
 ##' rows to components of the state vector.
 ##' One independent simulation will be performed for each column.
 ##' Note that in this case, \code{params} must also have \code{nrep} columns.
 ##'
-##' @param params a \code{npar} x \code{nrep} matrix of parameters.
-##' Each column is treated as an independent parameter set, in correspondence with the corresponding column of \code{xstart}.
+##' @param t0 the initial time, i.e., the time corresponding to the state in \code{x0}.
 ##'
-##' @param offset integer; the first \code{offset} times in \code{times} will be discarded.
+##' @param params a \code{npar} x \code{nrep} matrix of parameters.
+##' Each column is treated as an independent parameter set, in correspondence with the corresponding column of \code{x0}.
 ##'
 ##' @return
 ##' \code{rprocess} returns a rank-3 array with rownames.
 ##' Suppose \code{x} is the array returned.
-##' Then \preformatted{dim(x)=c(nvars,nrep,ntimes-offset),}
-##' where \code{nvars} is the number of state variables (=\code{nrow(xstart)}),
-##' \code{nrep} is the number of independent realizations simulated (=\code{ncol(xstart)}), and
+##' Then \preformatted{dim(x)=c(nvars,nrep,ntimes),}
+##' where \code{nvars} is the number of state variables (=\code{nrow(x0)}),
+##' \code{nrep} is the number of independent realizations simulated (=\code{ncol(x0)}), and
 ##' \code{ntimes} is the length of the vector \code{times}.
-##' \code{x[,j,k]} is the value of the state process in the \code{j}-th realization at time \code{times[k+offset]}.
-##' The rownames of \code{x} must correspond to those of \code{xstart}.
+##' \code{x[,j,k]} is the value of the state process in the \code{j}-th realization at time \code{times[k]}.
+##' The rownames of \code{x} will correspond to those of \code{x0}.
 ##'
 NULL
 
@@ -561,23 +561,21 @@ setMethod(
 setMethod(
   "rprocess",
   signature=signature(object="pomp"),
-  definition=function (object, xstart, times, params, offset = 0L, ...) {
+  definition=function (object, x0, t0, times, params, ...) {
     tryCatch(
-      rprocess.internal(object=object,xstart=xstart,times=times,params=params,
-        offset=offset,...),
+      rprocess.internal(object=object,x0=x0,t0=t0,times=times,params=params,...),
       error = function (e) pStop("rprocess",conditionMessage(e))
     )
   }
 )
 
-rprocess.internal <- function (object, xstart, times, params, offset = 0L,
-  .gnsi = TRUE, ...) {
-  storage.mode(xstart) <- "double"
+rprocess.internal <- function (object, x0, t0, times, params, ...,
+  .gnsi = TRUE) {
+  storage.mode(x0) <- "double"
   storage.mode(params) <- "double"
-  if (!is.finite(offset)) pStop_("invalid ",sQuote("offset"),".")
   pompLoad(object)
   on.exit(pompUnload(object))
-  .Call(P_do_rprocess,object,xstart,times,params,offset,.gnsi)
+  .Call(P_do_rprocess,object,x0,t0,times,params,.gnsi)
 }
 
 ##' skeleton

inst/NEWS

@@ -1,11 +1,25 @@
 _N_e_w_s _f_o_r _p_a_c_k_a_g_e '_p_o_m_p'
 
-_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _2._0._1_0:
+_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _2._0._1_1:
 
         • It is now possible to adjust the observation times, zero-time
           in the ‘sir’ and ‘sir2’ examples.  One can also change the
           default RNG seed used in generating the simulated data.
 
+        • The codes underlying the ‘rprocess’ workhorse have been
+          streamlined to obviate an unnecessary copy operation.  As a
+          result, the interface to the ‘rprocess’ workhorse has been
+          changed.
+
+        • The interface to the new ‘flow’ workhorse has been changed to
+          match.
+
+_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _2._0._1_0:
+
+        • The new ‘flow’ workhorse allows integration/iteration of the
+          deterministic skeleton at arbitrary times from an arbitrary
+          initial condition.
+
 _C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _2._0._9:
 
         • When ‘coef’ is applied to a ‘listie’, the result can be more

inst/NEWS.Rd

@@ -1,9 +1,17 @@
 \name{NEWS}
 \title{News for package `pomp'}
-\section{Changes in \pkg{pomp} version 2.0.10}{
+\section{Changes in \pkg{pomp} version 2.0.11}{
   \itemize{
     \item It is now possible to adjust the observation times, zero-time in the \code{sir} and \code{sir2} examples.
     One can also change the default RNG seed used in generating the simulated data.
+    \item The codes underlying the \code{rprocess} workhorse have been streamlined to obviate an unnecessary copy operation.
+    As a result, the interface to the \code{rprocess} workhorse has been changed.
+    \item The interface to the new \code{flow} workhorse has been changed to match.
+  }
+}
+\section{Changes in \pkg{pomp} version 2.0.10}{
+  \itemize{
+    \item The new \code{flow} workhorse allows integration/iteration of the deterministic skeleton at arbitrary times from an arbitrary initial condition.
   }
 }
 \section{Changes in \pkg{pomp} version 2.0.9}{