contactsimulator: Contact model simulator

##' Simulations of a contact process for spatio-temporal model
##'
# ##' \code{simulation_contact} generates simulations of the contact process for spatio-temporal model
# ##' processes.
# ##'
# ##' @name simulation_contact
# ##' @docType methods
# ##' @rdname simulation_contact
# ##' @include contact_class.R contact.R
# ##' @importFrom stats simulation_contact
# ##'
# ##' @author Hola K. Adrakey
# ##' @family elementary contact methods
# ##'
# ##' @inheritParams contact
# ##'
# ##' @param object optional; if present, it should be the output of one of \pkg{contact}'s methods
# ##'
# ##' @param params a named numeric vector or a matrix with rownames
# ##' containing the parameters at which the simulations are to be performed.
# ##'
# ##' @param nsim The number of simulations to perform.
# ##' Note that the number of replicates will be \code{nsim} times \code{ncol(params)}.
# ##'
# ##' @param seed optional;
# ##' if set, the pseudorandom number generator (RNG) will be initialized with \code{seed}.  the random seed to use.
# ##' The RNG will be restored to its original state afterward.
# ##'
# ##' @param format the format in which to return the results.
# ##'
# ##' \code{format = "contacts"} causes the results to be returned as a single \dQuote{contact} object,
# ##' identical to \code{object} except for the latent states and observations,
# ##' which have been replaced by the simulated values.
# ##'
# ##' \code{format = "arrays"} causes the results to be returned as a list of two arrays.
# ##' The \dQuote{states} element will contain the simulated state trajectories in a rank-3  array with dimensions
# ##' \code{nvar} x \code{(ncol(params)*nsim)} x \code{ntimes}.
# ##' Here, \code{nvar} is the number of state variables and \code{ntimes} the length of the argument \code{times}.
# ##' The \dQuote{obs} element will contain the simulated data, returned as a rank-3 array with dimensions
# ##' \code{nobs} x \code{(ncol(params)*nsim)} x \code{ntimes}.
# ##' Here, \code{nobs} is the number of observables.
# ##'
# ##' \code{format = "data.frame"} causes the results to be returned as a single data frame containing
# ##' the time, states, and observations.
# ##' An ordered factor variable, \sQuote{.id}, distinguishes one simulation from another.
# ##'
# ##' @param include.data if \code{TRUE}, the original data are included (with \code{.id = "rep"}).
# ##' This option is ignored unless \code{format = "data.frame"}.
# ##'
# ##' @return
# ##' A single \dQuote{contact} object,
# ##' a \dQuote{contactList} object,
# ##' a named list of two arrays,
# ##' or a data frame, according to the \code{format} option.
# ##'
# ##' If \code{params} is a matrix, each column is treated as a distinct parameter set.
# ##' In this case, if \code{nsim=1},
# ##' then \code{simulate} will return one simulation for each parameter set.
# ##' If \code{nsim>1},
# ##' then \code{simulate} will yield \code{nsim} simulations for each parameter set.
# ##' These will be ordered such that
# ##' the first \code{ncol(params)} simulations represent one simulation
# ##' from each of the distinct parameter sets,
# ##' the second \code{ncol(params)} simulations represent a second simulation from each,
# ##' and so on.
# ##'
# ##' Adding column names to \code{params} can be helpful.
# ##'
# NULL
#
# setGeneric(
#   "simulation_contact",
#   function (object, nsim=1, seed=NULL, ...)
#     standardGeneric("simulation_contact")
# )

# ##' @name simulation-missing
# ##' @aliases simulation,missing-method
# ##' @rdname simulation
# ##' @export
# setMethod(
#   "simulate",
#   signature=signature(object="missing"),
#   definition=function (nsim = 1, seed = NULL,
#     times, t0,
#     params, rinit, rprocess, rmeasure,
#     format = c("contacts", "arrays", "data.frame"),
#     include.data = FALSE,
#     ..., verbose = getOption("verbose", FALSE)) {
#
#     tryCatch(
#       simulate.internal(
#         object=NULL,
#         nsim=nsim,
#         seed=seed,
#         times=times,
#         t0=t0,
#         params=params,
#         rinit=rinit,
#         rprocess=rprocess,
#         rmeasure=rmeasure,
#         format=match.arg(format),
#         include.data=include.data,
#         ...,
#         verbose=verbose
#       ),
#       error = function (e) pStop("simulate",conditionMessage(e))
#     )
#
#   }
# )
#
# ##' @name simulate-data.frame
# ##' @aliases simulate simulate,data.frame-method
# ##' @rdname simulate
# ##' @export
# setMethod(
#   "simulate",
#   signature=signature(object="data.frame"),
#   definition=function (object,
#     nsim = 1, seed = NULL,
#     times, t0,
#     params, rinit, rprocess, rmeasure,
#     format = c("contacts", "arrays", "data.frame"),
#     include.data = FALSE,
#     ..., verbose = getOption("verbose", FALSE)) {
#
#     tryCatch(
#       simulate.internal(
#         object,
#         nsim=nsim,
#         seed=seed,
#         times=times,
#         t0=t0,
#         params=params,
#         rinit=rinit,
#         rprocess=rprocess,
#         rmeasure=rmeasure,
#         format=match.arg(format),
#         include.data=include.data,
#         ...,
#         verbose=verbose
#       ),
#       error = function (e) pStop("simulate",conditionMessage(e))
#     )
#
#   }
# )
#
# ##' @name simulation-contact
# ##' @aliases simulate simulate,pomp-method
# ##' @rdname simulation
# ##' @export
# setMethod(
#   "simulate",
#   signature=signature(object="contact"),
#   definition=function (object,
#     nsim = 1, seed = NULL,
#     format = c("contacts", "arrays", "data.frame"),
#     include.data = FALSE,
#     ..., verbose = getOption("verbose", FALSE)) {
#
#     tryCatch(
#       simulate.internal(
#         object,
#         nsim=nsim,
#         seed=seed,
#         format=match.arg(format),
#         include.data=include.data,
#         ...,
#         verbose=verbose
#       ),
#       error = function (e) pStop("simulate",conditionMessage(e))
#     )
#
#   }
# )
#
# simulate.internal <- function (object, nsim = 1L, seed = NULL, params,
#   format, include.data = FALSE, ..., .gnsi = TRUE, verbose) {
#
#   object <- contact(object,...,verbose=verbose)
#
#   if (undefined(object@rprocess)) pStop_(sQuote("rprocess")," is undefined.")
#
#   include.data <- as.logical(include.data)
#
#   if (length(nsim)!=1 || !is.numeric(nsim) || !is.finite(nsim) || nsim < 1)
#     pStop_(sQuote("nsim")," must be a positive integer.")
#   nsim <- as.integer(nsim)
#
#   if (missing(params)) params <- coef(object)
#   if (is.list(params)) params <- unlist(params)
#   if (is.null(params)) params <- numeric(0)
#   if (!is.numeric(params))
#     pStop_(sQuote("params")," must be named and numeric.")
#   params <- as.matrix(params)
#   storage.mode(params) <- "double"
#
#   if (ncol(params) == 1) coef(object) <- params[,1L]
#
#   pompLoad(object,verbose=verbose)
#   on.exit(pompUnload(object,verbose=verbose))
#
#   return.type <- switch(format,arrays=0L,data.frame=0L,pomps=1L)
#
#   sims <- freeze(
#     .Call(P_do_simulate,object,params,nsim,return.type,.gnsi),
#     seed=seed
#   )
#
#   if (format == "data.frame") {
#
#     nsims <- ncol(sims$states)
#     ntimes <- length(time(object))
#     simnames <- colnames(sims$states)
#     if (is.null(simnames)) simnames <- seq_len(nsims)
#
#     dm <- dim(sims$states)
#     nm <- rownames(sims$states)
#     dim(sims$states) <- c(dm[1L],prod(dm[-1L]))
#     rownames(sims$states) <- nm
#
#     dm <- dim(sims$obs)
#     nm <- rownames(sims$obs)
#     dim(sims$obs) <- c(dm[1L],prod(dm[-1L]))
#     rownames(sims$obs) <- nm
#
#     sims <- cbind(
#       time=rep(time(object),each=length(simnames)),
#       as.data.frame(t(sims$states)),
#       as.data.frame(t(sims$obs)),
#       .id=rep(simnames,times=ntimes)
#     )
#
#     names(sims)[[1]] <- object@timename
#
#     if (include.data) {
#       dat <- as.data.frame(object)
#       dat$.id <- "data"
#       allnm <- union(names(dat),names(sims))
#       for (n in setdiff(allnm,names(dat))) dat[[n]] <- NA
#       for (n in setdiff(allnm,names(sims))) sims[[n]] <- NA
#       sims <- rbind(dat,sims)
#       sims$.id <- ordered(sims$.id,levels=c("data",simnames))
#     } else {
#       sims$.id <- ordered(sims$.id,levels=simnames)
#     }
#     othernm <- setdiff(names(sims),c(object@timename,".id"))
#
#     sims <- sims[,c(object@timename,".id",othernm)]
#
#   } else if (format == "contacts") {
#
#     if (length(sims) == 1) sims <- sims[[1]]
#     else sims <- do.call(concat,sims)
#
#   }
#
#   sims
#
# }