R/dynemu_pred.R

In dynemu: Emulation of Dynamic Simulators via One-Step-Ahead Approach

#' Predictive posterior computation for dynamic systems using one-steap-ahead approach by Gaussian process (GP) emulators
#'
#' @description The function computes the predictive posterior distribution (mean and variance)
#' at all time steps for dynamic systems modeled by GP emulators.
#' It can use either the exact computation or Monte Carlo (MC) approximation to compute the predictive posterior.
#'
#' @details Given a trained set of GP emulators `dynGP` fitted using \code{dynemu_GP}, this function:
#' 1. Computes the predictive posterior mean and variance for each state variable at all time steps.
#' 2. The function can either:
#'    - Use the exact computation for a closed-form solution by \code{method="Exact"}.
#'    - Use the MC approximation method to generate samples and estimate the posterior by \code{method="MC"}.
#'
#' @param dynGP list of trained GP emulators fitted by \code{dynemu_GP}, each corresponding to a state variable.
#' @param times numeric vector specifying the time sequence at which predictions are to be made.
#' @param xnew numeric vector or single-row matrix specifying the initial state values at time 0. The number of columns must match the input dimension of the GP emulators.
#' @param method character specifying the method for prediction, to be chosen between \code{"Exact"}(exact closed-form solution), or \code{"MC"}(MC approximation). Default is \code{"Exact"}.
#' @param mc.sample a number of mc samples generated for MC approximation. Required only if \code{method="MC"}.
#' @param trace logical indicating whether to print the progress bar. If \code{trace=TRUE}, the progress bar is printed. Default is \code{TRUE}.
#'
#' @return A list containing:
#' \itemize{
#'   \item \code{times}: The time sequence vector.
#'   \item \code{mu}: A matrix of predictive mean values for each state variable. Each row corresponds to a time step, and each column corresponds to a state variable.
#'   \item \code{sig2}: A matrix of predictive variance values for each state variable. Each row corresponds to a time step, and each column corresponds to a state variable.
#' }
#'
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @usage dynemu_pred(dynGP, times, xnew, method="Exact", mc.sample=NULL, trace=TRUE)
#' @export
#'
#' @examples
#' \donttest{
#' library(lhs)
#' ### Lotka-Volterra equations ###
#' LVmod0D <- function(Time, State, Pars) {
#'   with(as.list(c(State, Pars)), {
#'     IngestC <- rI * P * C
#'     GrowthP <- rG * P * (1 - P/K)
#'     MortC <- rM * C
#'
#'     dP <- GrowthP - IngestC
#'     dC <- IngestC * AE - MortC
#'     return(list(c(dP, dC)))
#'   })
#' }
#'
#' ### Define parameters ###
#' pars <- c(rI = 1.5, rG = 1.5, rM = 2, AE = 1, K = 10)
#'
#' ### Define time sequence ###
#' times <- seq(0, 30, by = 0.01)
#'
#' ### Initial conditions ###
#' set.seed(1)
#' d <- 2
#' n <- 12*d
#' Design <- maximinLHS(n, d) * 5 # Generate LHS samples in range [0,5]
#' colnames(Design) <- c("P", "C")
#'
#' ### Fit GP emulators ###
#' fit.dyn <- dynemu_GP(LVmod0D, times, pars, Design)
#'
#' ### Define initial test values ###
#' state <- c(P = 1, C = 2)
#'
#' ### Generate test set ###
#' test <- dyn_solve(LVmod0D, times, pars, state)
#'
#' ### Define number of MC samples ###
#' mc.sample <- 1000
#'
#' ### Prediction ###
#' pred.exact.dyn <- dynemu_pred(fit.dyn, times, state, method="Exact")
#' pred.exact.mu <- pred.exact.dyn$mu
#' pred.exact.sig2 <- pred.exact.dyn$sig2
#'
#' pred.mc.dyn <- dynemu_pred(fit.dyn, times, state, method="MC", mc.sample=mc.sample, trace=TRUE)
#' pred.mc.mu <- pred.mc.dyn$mu
#' pred.mc.sig2 <- pred.mc.dyn$sig2
#'
#' ### Compute RMSE ###
#' sqrt(mean((pred.exact.mu[,1]-test$P)^2)) # 0.0007686341
#' sqrt(mean((pred.exact.mu[,2]-test$C)^2)) # 0.0005255164
#'
#' sqrt(mean((pred.mc.mu[,1]-test$P)^2)) # 0.03050849
#' sqrt(mean((pred.mc.mu[,2]-test$C)^2)) # 0.02330542
#' }
#'

dynemu_pred <- function(dynGP, times, xnew, method="Exact", mc.sample=NULL, trace=TRUE){

  if (is.null(dim(xnew))) xnew <- matrix(xnew, nrow = 1)
  colnames(xnew) <- colnames(dynGP[[1]]$X)

  d <- length(dynGP)

  predy.out <- matrix(0, ncol=d, nrow=length(times))
  predsig2.out <- matrix(0, ncol=d, nrow=length(times))
  colnames(predy.out) <- colnames(predsig2.out) <- colnames(dynGP[[1]]$X)

  predy.out[1,] <- xnew
  predsig2.out[1,] <- matrix(rep(0, d), nrow=1)

  if(method=="MC" && is.null(mc.sample)){
    mc.sample <- 1000
    message("The number of MC samples was not specified. Using the default value of 1000. \n")
  }
  if(trace) pb <- txtProgressBar(min = 0, max = length(times), style = 3)

  warned <- FALSE
  old_warn <- getOption("warn")
  on.exit(options(warn = old_warn), add = TRUE)
  options(warn = 1)  # immediate warning display

  ### Emulating dynamic simulator ###
  if(method=="Exact"){
    for (i in 2:length(times)) {
      pred <- dynemu_exact(predy.out[i-1,], predsig2.out[i-1,], dynGP)

      if (!warned && any(pred$predy > 1e3)) {
        warning(sprintf("Prediction at time %s exceeds 1e3: max = %.2f", times[i], max(pred$predy)))
        warned <- TRUE
      }

      predy.out[i, ] <- pred$predy
      predsig2.out[i, ] <- pred$predsig2
      if(trace) setTxtProgressBar(pb, i) # print(i)
    }
  }else if(method=="MC"){
    for (i in 2:length(times)) {
      pred <- dynemu_mc(predy.out[i-1,], predsig2.out[i-1,], dynGP, mc.sample)

      if (!warned && any(pred$predy > 1e3)) {
        message(sprintf("Prediction at time %s exceeds 1e3: max = %.2f", times[i], max(pred$predy)))
        # warning(sprintf("Prediction at time %s exceeds 1e3: max = %.2f", times[i], max(pred$predy)))
        warned <- TRUE
      }

      predy.out[i, ] <- pred$predy
      predsig2.out[i, ] <- pred$predsig2
      if(trace) setTxtProgressBar(pb, i) # print(i)
    }
  }else{
    stop("method should be Exact or MC. \n")
  }

  return(list(times=times, mu = predy.out, sig2 = predsig2.out))
}