R/flow.R

In kingaa/pomp: Statistical Inference for Partially Observed Markov Processes

##' flow workhorse
##'
##' Compute the flow generated by a deterministic vectorfield or map.
##'
##' In the case of a discrete-time system (map), \code{flow} iterates the map to yield trajectories of the system.
##' In the case of a continuous-time system (vectorfield), \code{flow} uses the numerical solvers in \pkg{\link[deSolve]{deSolve}} to integrate the vectorfield starting from given initial conditions.
##'
##' @name flow
##' @rdname flow
##' @include pomp_class.R skeleton_spec.R workhorses.R
##' @aliases flow,missing-method flow,ANY-method
##' @family pomp workhorses
##' @family deterministic methods
##' @concept extending the pomp package
##' @concept low-level interface
##' @importFrom deSolve ode diagnostics
##' @importFrom stats setNames
##' @inheritParams dmeasure
##' @inheritParams pomp
##' @param x0 an array with dimensions \code{nvar} x \code{nrep} giving the initial conditions of the trajectories to be computed.
##' @param t0 the time at which the initial conditions are assumed to hold.
##' By default, this is the zero-time (see \code{\link{timezero}}).
##' @param times a numeric vector (length \code{ntimes}) containing times at which the itineraries are desired.
##' These must be in non-decreasing order with \code{times[1]>t0}.
##' By default, this is the full set of observation times (see \code{\link{time}}).
##' @param ... Additional arguments are passed to the ODE integrator (if the skeleton is a vectorfield) and are ignored if it is a map.
##' See \code{\link[deSolve]{ode}} for a description of the additional arguments accepted by the ODE integrator.
##' By default, this is the parameter vector stored in \code{object} (see \code{\link[pomp]{coef}}).
##' @return
##' \code{flow} returns an array of dimensions \code{nvar} x \code{nrep} x \code{ntimes}.
##' If \code{x} is the returned matrix, \code{x[i,j,k]} is the i-th component of the state vector at time \code{times[k]} given parameters \code{params[,j]}.
##'
NULL

setGeneric(
  "flow",
  function (object, ...)
    standardGeneric("flow")
)

setMethod(
  "flow",
  signature=signature(object="missing"),
  definition=function (...) {
    reqd_arg("flow","object")
  }
)

setMethod(
  "flow",
  signature=signature(object="ANY"),
  definition=function (object, ...) {
    undef_method("flow",object)
  }
)

##' @rdname flow
##' @export
setMethod(
  "flow",
  signature=signature(object="pomp"),
  definition=function (
    object, x0,
    t0 = timezero(object),
    times = time(object),
    params = coef(object),
    ...,
    verbose = getOption("verbose", FALSE)) {

    tryCatch(
      flow_internal(object=object,x0=x0,t0=t0,times=times,params=params,
        ...,verbose=verbose),
      error = function (e) pStop(who="flow",conditionMessage(e))
    )
    
  }
)

flow_internal <- function (object, x0, t0, times, params, ...,
  .gnsi = TRUE, verbose) {

  verbose <- as.logical(verbose)

  t0 <- as.numeric(t0)
  times <- as.numeric(times)

  if (length(times)==0)
    pStop_(sQuote("times")," is empty, there is no work to do.")

  if (any(diff(times)<0))
    pStop_(sQuote("times")," must be a non-decreasing numeric sequence.")

  if (times[1]<t0)
    pStop_(sQuote("times[1]")," must be no later than ",sQuote("t0"),".")

  params <- as.matrix(params)
  nrep <- ncol(params)
  storage.mode(params) <- "double"

  x0 <- as.matrix(x0)
  nvar <- nrow(x0)
  storage.mode(x0) <- "double"

  statenames <- rownames(x0)
  repnames <- colnames(x0)
  dim(x0) <- c(nvar,nrep,1)
  dimnames(x0) <- list(statenames,repnames,NULL)

  ntimes <- length(times)

  type <- object@skeleton@type          # map or vectorfield?

  pompLoad(object)
  on.exit(pompUnload(object))

  if (type == skeletontype$map) {                  ## MAP

    x <- .Call(P_iterate_map,object,times,t0,x0,params,.gnsi)
    .gnsi <- FALSE

  } else if (type == skeletontype$vectorfield) {   ## VECTORFIELD

    znames <- object@accumvars
    if (length(znames)>0) {
      if (!all(znames %in% statenames)) {
        mz <- znames[which(!(znames %in% statenames))[1L]]
        pStop_("accumulator variable ",sQuote(mz)," not found among the state variables.")
      }
      x0[znames,,] <- 0
    }

    .Call(P_pomp_desolve_setup,object,x0,params,.gnsi)
    .gnsi <- FALSE

    X <- tryCatch(
      ode(
        y=x0,
        times=c(t0,times),
        func="pomp_vf_eval",
        dllname="pomp",
        initfunc=NULL,
        parms=NULL,
        ...
      ),
      error = function (e) {
        pStop_("error in ODE integrator: ",conditionMessage(e))
      }
    )

    .Call(P_pomp_desolve_takedown)

    if (attr(X,"istate")[1L] != 2)
      pWarn(who="flow",
        "abnormal exit from ODE integrator, istate = ",attr(X,'istate')[1L])

    if (verbose) diagnostics(X)

    x <- array(data=t(X[-1L,-1L]),dim=c(nvar,nrep,ntimes),
      dimnames=list(statenames,repnames,NULL))

    for (z in znames)
      for (r in seq_len(ncol(x)))
        x[z,r,-1] <- diff(x[z,r,])

  } else {                  ## DEFAULT SKELETON

    x <- array(data=NA_real_,dim=c(nrow(x0),ncol(x0),length(times)),
      dimnames=list(rownames(x0),colnames(x0),NULL))

  }

  dimnames(x) <- setNames(dimnames(x),c("name",".id",object@timename))

  x
}