R/dynemu_GP.R

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

#' Fit Gaussian process (GP) emulators for One-step-ahead approach
#'
#' @description The function fits GP emulators for each state variable
#' by solving an ODE model one step ahead and using the generated outputs as training data.
#'
#' @details Given an initial set of input conditions, this function:
#' 1. Solves the ODE system for a short time interval using \code{dyn_solve}.
#' 2. Fits GP emulators to model the next-step evolution of each state variable.
#'
#' The ODE system is defined by the user through \code{odefun},
#' and the GP emulators are trained using the generated one-step-ahead outputs.
#'
#' @param odefun function defining the ODE system. It should return a list where the first element is a numeric vector of derivatives.
#' @param times numeric vector of time points where the solution is computed. Only the first two time points are used for one-step-ahead prediction.
#' @param param numeric scalar or vector of parameters to be passed to \code{odefun}.
#' @param X A matrix where each row represents an initial condition for the ODE system. Column names must be provided to match the state variables in \code{odefun}.
#' @param ... Additional arguments for compatibility with \code{GP}.
#' @return A list of GP emulators, each corresponding to a state variable.
#'
#' @usage dynemu_GP(odefun, times, param, X, ...)
#' @export
#' @examples
#' 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)
#' print(fit.dyn)
#'

dynemu_GP <- function(odefun, times, param, X, ...){
  if (is.null(dim(X))) X <- matrix(X, nrow = 1)
  if (is.null(colnames(X))){
    stop("X requires column names for ODE function \n")
  }

  ### solving ODEs for output ###
  output <- do.call(rbind, lapply(apply(X, 1, dyn_solve, times = times[1:2], odefun = odefun, param = param), function(res) {
    as.vector(unlist(res))  # Flatten each sublist into a vector
  }))[,-c(1:2)]
  output <- output[, seq(2, ncol(output), by = 2)]

  ### fit GPs ###
  gp_models <- vector("list", ncol(output))
  for (i in 1:ncol(output)) {
    gp_models[[i]] <- GP(X, output[, i], constant = TRUE, ...)
  }

  class(gp_models) <- "dynemu"
  return(dynemuGP=gp_models)
}