R/simulators.R

In readsdr: Translate Models from System Dynamics Software into 'R'

#' Simulate a System Dynamics model
#'
#' @param ds_inputs A list of deSolve inputs generated by read_xmile
#' @param start_time A number indicating the time at which the simulation begins.
#' @param stop_time A number indicating the time at which the simulation ends.
#' @param timestep A number indicating the time interval for the simulation.
#'  Also known as \code{dt}.
#' @param integ_method A string indicating the integration method. It can be
#'  either "euler" or "rk4"
#'
#' @return a data frame
#' @export
#'
#' @examples
#' path      <- system.file("models", "SIR.stmx", package = "readsdr")
#' ds_inputs <- xmile_to_deSolve(path)
#' sd_simulate(ds_inputs, 0, 1, 0.25, "rk4")
sd_simulate <- function(ds_inputs, start_time = NULL, stop_time = NULL,
                        timestep = NULL, integ_method = "euler") {

  .env <- new.env()

  if(!(integ_method %in% c("euler", "rk4"))) stop("Invalid integration method")

  if(is.null(start_time)) start_time  <- ds_inputs$sim_params$start
  if(is.null(stop_time))  stop_time   <- ds_inputs$sim_params$stop
  if(is.null(timestep))   timestep    <- ds_inputs$sim_params$dt

  # Create time vector
  simtime  <- seq(start_time, stop_time, by = timestep)

  ds_func              <- ds_inputs$func
  environment(ds_func) <- .env

  if(!is.null(ds_inputs$delayed_vars)) {
    n_points  <- length(simtime)
    n_vars    <- length(ds_inputs$delayed_vars)
    NAs       <- rep(NA, n_points * n_vars)
    memory_df <- matrix(c(simtime, NAs), ncol = n_vars + 1) %>% as.data.frame()
    colnames(memory_df) <- c("time", ds_inputs$delayed_vars)
    rownames(memory_df) <- simtime
    .env$.memory        <- memory_df
  }

  ode_args <- list(y      = ds_inputs$stocks,
                   times  = simtime,
                   func   = ds_func,
                   parms  = ds_inputs$consts,
                   method = integ_method)

  if("graph_funs" %in% names(ds_inputs)) {
    ode_args$graph_funs <- ds_inputs$graph_funs
  }

  results <- do.call(deSolve::ode, ode_args)

  data.frame(results)
}

#' Perform a sensitivity run on a System Dynamics model
#'
#' \code{sd_sensitivity_run} returns a data frame with the simulation of a
#' model for several iterations of different inputs.
#'
#' @param consts_df A data frame that contains the values of constants to
#'   simulate. Each column corresponds to a constant and each row to an
#'   iteration. If \code{stocks_df} is also supplied, both data frames must have
#'   the same number of rows.
#' @param stocks_df A data frame that containts the initial value of stocks to
#'   be explored. Each column corresponds to a stock and each row to an
#'   iteration. If \code{consts_df} is also supplied, both data frames must have
#'   the same number of rows.
#' @param multicore A boolean value that indicates whether the process
#' is parallelised.
#' @param n_cores An integer indicating the number of cores for the parallel run.
#' @param reporting_interval A real number indicating the interval at which the
#'  simulation results are returned. The default is set to \code{1}. For
#'  instance, if the simulation runs from 0 to 10. This function returns the
#'  results at times 0, 1, 2, ..., 10.
#' @inheritParams sd_simulate
#'
#' @return A data frame
#' @export
#'
#' @examples
#' path      <- system.file("models", "SIR.stmx", package = "readsdr")
#' ds_inputs <- xmile_to_deSolve(path)
#' consts_df <- data.frame(i = c(0.25, 0.30))
#' sd_sensitivity_run(ds_inputs, consts_df)
sd_sensitivity_run <- function(ds_inputs, consts_df = NULL, stocks_df = NULL,
                               start_time = NULL, stop_time = NULL,
                               timestep = NULL, integ_method = "euler",
                               multicore = FALSE,
                               n_cores = NULL,
                               reporting_interval = 1) {

  if(!(integ_method %in% c("euler", "rk4"))) stop("Invalid integration method")

  if(is.null(start_time)) start_time  <- ds_inputs$sim_params$start
  if(is.null(stop_time))  stop_time   <- ds_inputs$sim_params$stop
  if(is.null(timestep))   timestep    <- ds_inputs$sim_params$dt

  # Create time vectors
  simtime     <- seq(start_time, stop_time, by = timestep)
  report_time <- seq(start_time, stop_time, by = reporting_interval)

  ode_args <- list(y      = NULL,
                   times  = simtime,
                   func   = ds_inputs$func,
                   parms  = NULL,
                   method = integ_method)

  if("graph_funs" %in% names(ds_inputs)) {
    ode_args$graph_funs <- ds_inputs$graph_funs
  }

  consts_names <- names(ds_inputs$consts)
  stocks_names <- names(ds_inputs$stocks)


  #-----------------------------------------------------------------------------
  if(!is.null(consts_df)) {

    sens_consts     <- colnames(consts_df)

    validate_consts(sens_consts, consts_names)

    missing_consts  <- consts_names[!consts_names %in% sens_consts]

    if(length(missing_consts) > 0) {
      consts_df <- fill_df(consts_df, missing_consts, ds_inputs$consts)
    }

    const_sensitivity_list <- do.call(function(...) Map(list,...), consts_df)
  }
  #-----------------------------------------------------------------------------

  #-----------------------------------------------------------------------------
  if(!is.null(stocks_df)) {

    sens_stocks     <- colnames(stocks_df)
    missing_stocks  <- stocks_names[!stocks_names %in% sens_stocks]

    if(length(missing_stocks) > 0) {
      stocks_df <- fill_df(stocks_df, missing_stocks, ds_inputs$stocks)
    }

    stocks_df              <- stocks_df[, c(stocks_names), drop = FALSE]
    stock_sensitivity_list <- do.call(function(...) Map(list,...), stocks_df)
  }

  #-----------------------------------------------------------------------------

  if(!is.null(consts_df) & is.null(stocks_df)) {

    iters                  <- nrow(consts_df)
    stock_sensitivity_list <- rep(list(as.list(ds_inputs$stocks)), iters)
  }

  if(is.null(consts_df) & !is.null(stocks_df)) {

    iters   <- nrow(stocks_df)
    const_sensitivity_list <- rep(list(as.list(ds_inputs$consts)), iters)
  }

  if(!is.null(consts_df) & !is.null(stocks_df)) {
    row_consts <- nrow(consts_df)
    row_stocks <- nrow(stocks_df)

    if(row_consts != row_stocks) {
      stop("the number of rows in both data frames (consts & stocks) must be of equal size", call. = FALSE)
    }

    iters   <- nrow(consts_df)
  }

  df_list <- const_stock_sensitivity(const_sensitivity_list,
                                     stock_sensitivity_list,
                                     ode_args, multicore, n_cores,
                                     report_time = report_time)

  sensitivity_df <- do.call("rbind", df_list)

  sensitivity_df$iter  <- rep(1:iters, each = length(report_time))
  sensitivity_df
}

const_stock_sensitivity <- function(const_sensitivity_list,
                                    stock_sensitivity_list, ode_args,
                                    multicore, n_cores, report_time) {

  sens_list <- Map(list, const_sensitivity_list, stock_sensitivity_list)

  if(multicore) {

    if(is.null(n_cores)) n_cores <- future::availableCores() - 1

    if(n_cores > 1) {

      future::plan(future::multisession, workers = n_cores)

      progressr::with_progress({

        p <- progressr::progressor(along = sens_list)

        df_list <- future.apply::future_lapply(sens_list, function(sens_obj,
                                                                   do_sens) {
          p()
          do_sens(sens_obj, ode_args, report_time)
        }, do_sens = do_sens)
      })


        return(df_list)
    }
  }

  df_list <- lapply(sens_list, do_sens,
                    ode_args = ode_args,
                    report_time = report_time)
}

do_sens <- function(sens_obj, ode_args, report_time) {

  const_list     <- sens_obj[1]
  ode_args$parms <- unlist(const_list)

  stock_list <- sens_obj[2]
  ode_args$y <- unlist(stock_list)

  result_matrix  <- do.call(deSolve::ode, ode_args)
  df             <- data.frame(result_matrix)

  df <- df[df$time %in% report_time, ]

  df
}

fill_df <- function(df, missing, elems) {

  for(ms in missing) df[ms] <- elems[ms]

  df
}

update_sim_params <- function(ds_inputs, start_time, stop_time, timestep) {

  if(!is.null(start_time)) ds_inputs$sim_params$start <- start_time
  if(!is.null(stop_time))  ds_inputs$sim_params$stop  <- stop_time
  if(!is.null(timestep))   ds_inputs$sim_params$dt    <- timestep

  ds_inputs
}

validate_consts <- function(sens_consts, consts_names) {

  invalid_names <- sens_consts[!sens_consts %in% consts_names]


  if(length(invalid_names) > 0L) {

    inv_consts <- paste0("`", invalid_names, "`") # invalid consts
    inv_consts <- paste(inv_consts, sep = ", ")
    msg <- stringr::str_glue("{inv_consts} doesn't/don't exist in the model")
    stop(msg, call. = FALSE)
  }

  invisible(NULL)
}