From f874efbe9d8baa96dfb2b1151aef55d3783ed0a1 Mon Sep 17 00:00:00 2001 From: arindam fadikar Date: Wed, 25 Mar 2026 20:21:10 -0400 Subject: [PATCH 1/5] added stochastic mode handling --- R/classes.R | 19 +++ R/metaRVM.R | 130 ++++++++++++------ R/parse_config.R | 56 ++++++-- README.Rmd | 11 +- README.md | 10 ++ inst/extdata/example_config.yaml | 1 + inst/extdata/example_config_checkpoint.yaml | 1 + inst/extdata/example_config_dist.yaml | 1 + inst/extdata/example_config_stochastic.yaml | 29 ++++ .../extdata/example_config_subgroup_dist.yaml | 1 + man/parse_config.Rd | 5 + tests/testthat/test-metaRVM.R | 112 +++++++++++++++ tests/testthat/test-parse_config.R | 92 ++++++++++++- vignettes/checkpointing.Rmd | 14 +- vignettes/getting-started.Rmd | 11 +- vignettes/running-a-simulation.Rmd | 67 ++++++++- vignettes/yaml-configuration.Rmd | 85 ++++++++++-- 17 files changed, 571 insertions(+), 74 deletions(-) create mode 100644 inst/extdata/example_config_stochastic.yaml diff --git a/R/classes.R b/R/classes.R index 9667639..a2843ba 100644 --- a/R/classes.R +++ b/R/classes.R @@ -342,6 +342,18 @@ MetaRVMResults <- R6::R6Class( cat("Instances:", self$run_info$n_instances, "\n") cat("Populations:", self$run_info$n_populations, "\n") cat("Date range:", paste(self$run_info$date_range, collapse = " to "), "\n") + if (!is.null(self$run_info$nsim)) { + cat("Parameter sets (nsim):", self$run_info$nsim, "\n") + } + if (!is.null(self$run_info$nrep)) { + cat("Replicates per set (nrep):", self$run_info$nrep, "\n") + } + if (!is.null(self$run_info$simulation_mode)) { + cat("Simulation mode:", self$run_info$simulation_mode, "\n") + } + if (!is.null(self$run_info$random_seed)) { + cat("Random seed:", self$run_info$random_seed, "\n") + } cat("Total observations:", nrow(self$results), "\n") cat("Disease states:", paste(unique(self$results$disease_state), collapse = ", "), "\n") invisible(self) @@ -435,6 +447,13 @@ MetaRVMResults <- R6::R6Class( new_run_info <- list( created_at = Sys.time(), original_created_at = self$run_info$created_at, + N_pop = self$run_info$N_pop, + nsim = self$run_info$nsim, + nrep = self$run_info$nrep, + simulation_mode = self$run_info$simulation_mode, + random_seed = self$run_info$random_seed, + delta_t = self$run_info$delta_t, + checkpointing_enabled = self$run_info$checkpointing_enabled, n_instances = length(unique(subset_results$instance)), n_populations = private$calculate_n_populations(subset_results), date_range = if(nrow(subset_results) > 0) range(subset_results$date, na.rm = TRUE) else c(NA, NA), diff --git a/R/metaRVM.R b/R/metaRVM.R index 19c68f7..d08f3d0 100644 --- a/R/metaRVM.R +++ b/R/metaRVM.R @@ -114,6 +114,34 @@ metaRVM <- function(config_input) { # pass inputs to meta_sim nsim <- config_obj$config_data$nsim + nrep <- if ("nrep" %in% names(config_obj$config_data)) { + as.integer(config_obj$config_data$nrep) + } else { + 1L + } + if (is.na(nrep) || nrep < 1) { + stop("nrep must be a positive integer") + } + simulation_mode <- if ("simulation_mode" %in% names(config_obj$config_data)) { + config_obj$config_data$simulation_mode + } else { + "deterministic" + } + is_stoch <- identical(tolower(as.character(simulation_mode)), "stochastic") + run_seed_base <- NA_integer_ + if (is_stoch) { + run_seed <- config_obj$config_data$random_seed + if (is.null(run_seed)) { + run_seed <- sample.int(.Machine$integer.max, 1) + } else { + run_seed <- suppressWarnings(as.integer(run_seed)[1]) + if (is.na(run_seed)) { + stop("random_seed must be coercible to a single integer value") + } + } + config_obj$config_data$random_seed <- run_seed + run_seed_base <- run_seed + } nsteps <- floor(config_obj$config_data$sim_length / config_obj$config_data$delta_t) day_name <- weekdays(config_obj$config_data$start_date) start_day <- dplyr::case_when( @@ -128,51 +156,75 @@ metaRVM <- function(config_input) { ) out <- data.table::data.table() + run_idx <- 0L for (ii in 1:nsim){ + for (rr in 1:nrep){ + run_idx <- run_idx + 1L + if (is_stoch) { + set.seed(run_seed_base + run_idx - 1L) + } - o <- meta_sim(is.stoch = 0, - nsteps = nsteps, - N_pop = config_obj$config_data$N_pop, - S0 = config_obj$config_data$S_ini, - I0 = config_obj$config_data$I_symp_ini, - P0 = config_obj$config_data$P_ini, - V0 = config_obj$config_data$V_ini, - R0 = config_obj$config_data$R_ini, - H0 = config_obj$config_data$H_ini, - D0 = config_obj$config_data$D_ini, - E0 = config_obj$config_data$E_ini, - Ia0 = config_obj$config_data$I_asymp_ini, - Ip0 = config_obj$config_data$I_presymp_ini, - m_weekday_day = config_obj$config_data$m_wd_d, - m_weekday_night = config_obj$config_data$m_wd_n, - m_weekend_day = config_obj$config_data$m_we_d, - m_weekend_night = config_obj$config_data$m_we_n, - start_day = start_day, - delta_t = config_obj$config_data$delta_t, - vac_mat = config_obj$config_data$vac_mat, - ts = config_obj$config_data$ts[ii, ], - dv = config_obj$config_data$dv[ii, ], - de = config_obj$config_data$de[ii, ], - pea = config_obj$config_data$pea[ii, ], - dp = config_obj$config_data$dp[ii, ], - da = config_obj$config_data$da[ii, ], - ds = config_obj$config_data$ds[ii, ], - psr = config_obj$config_data$psr[ii, ], - dh = config_obj$config_data$dh[ii, ], - phr = config_obj$config_data$phr[ii, ], - dr = config_obj$config_data$dr[ii, ], - ve = config_obj$config_data$ve[ii, ], - do_chk = config_obj$config_data$do_chk, - chk_time_steps = config_obj$config_data$chk_time_steps, - chk_file_names = config_obj$config_data$chk_file_names[ii, ]) + o <- meta_sim(is.stoch = is_stoch, + nsteps = nsteps, + N_pop = config_obj$config_data$N_pop, + S0 = config_obj$config_data$S_ini, + I0 = config_obj$config_data$I_symp_ini, + P0 = config_obj$config_data$P_ini, + V0 = config_obj$config_data$V_ini, + R0 = config_obj$config_data$R_ini, + H0 = config_obj$config_data$H_ini, + D0 = config_obj$config_data$D_ini, + E0 = config_obj$config_data$E_ini, + Ia0 = config_obj$config_data$I_asymp_ini, + Ip0 = config_obj$config_data$I_presymp_ini, + m_weekday_day = config_obj$config_data$m_wd_d, + m_weekday_night = config_obj$config_data$m_wd_n, + m_weekend_day = config_obj$config_data$m_we_d, + m_weekend_night = config_obj$config_data$m_we_n, + start_day = start_day, + delta_t = config_obj$config_data$delta_t, + vac_mat = config_obj$config_data$vac_mat, + ts = config_obj$config_data$ts[ii, ], + dv = config_obj$config_data$dv[ii, ], + de = config_obj$config_data$de[ii, ], + pea = config_obj$config_data$pea[ii, ], + dp = config_obj$config_data$dp[ii, ], + da = config_obj$config_data$da[ii, ], + ds = config_obj$config_data$ds[ii, ], + psr = config_obj$config_data$psr[ii, ], + dh = config_obj$config_data$dh[ii, ], + phr = config_obj$config_data$phr[ii, ], + dr = config_obj$config_data$dr[ii, ], + ve = config_obj$config_data$ve[ii, ], + do_chk = config_obj$config_data$do_chk, + chk_time_steps = config_obj$config_data$chk_time_steps, + chk_file_names = config_obj$config_data$chk_file_names[run_idx, ]) - o$instance <- ii - out <- rbind(out, o) + o$instance <- run_idx + out <- rbind(out, o) + } } - + run_info <- list( + created_at = Sys.time(), + N_pop = config_obj$config_data$N_pop, + nsim = nsim, + nrep = nrep, + n_instances = nsim * nrep, + simulation_mode = simulation_mode, + random_seed = config_obj$config_data$random_seed, + delta_t = config_obj$config_data$delta_t, + date_range = if (nrow(out) > 0) { + c(config_obj$config_data$start_date + 1, + config_obj$config_data$start_date + floor(config_obj$config_data$sim_length)) + } else { + c(NA, NA) + }, + checkpointing_enabled = isTRUE(config_obj$config_data$do_chk) + ) + # Create and return MetaRVMResults object - results_obj <- MetaRVMResults$new(out, config_obj) + results_obj <- MetaRVMResults$new(out, config_obj, run_info = run_info) return(results_obj) # return(out) } diff --git a/R/parse_config.R b/R/parse_config.R index bb9f688..55c8f23 100644 --- a/R/parse_config.R +++ b/R/parse_config.R @@ -19,6 +19,9 @@ #' \itemize{ #' \item \code{random_seed}: Optional random seed for reproducibility in case of stochastic simulations or stochastic parameters #' \item \code{nsim}: Number of simulation instances (default: 1) +#' \item \code{nrep}: Number of stochastic replicates per parameter set (default: 1) +#' \item \code{simulation_mode}: Optional simulation mode. Must be one of +#' \code{"deterministic"} or \code{"stochastic"} (default: \code{"deterministic"}). #' \item \code{start_date}: Simulation start date in MM/DD/YYYY format #' \item \code{length}: Simulation length in days #' \item \code{checkpoint_dir}: Optional checkpoint directory for saving intermediate results @@ -72,6 +75,8 @@ #' \item{start_date}{Simulation start date as Date object} #' \item{sim_length}{Simulation length in days} #' \item{nsim}{Number of simulation instances} +#' \item{nrep}{Number of stochastic replicates per parameter set} +#' \item{simulation_mode}{Simulation mode: \code{"deterministic"} or \code{"stochastic"}} #' \item{random_seed}{Random seed used (if any)} #' \item{delta_t}{Time step size (fixed at 0.5)} #' \item{chk_file_names, chk_time_steps, do_chk}{Checkpointing configuration} @@ -142,14 +147,6 @@ parse_config <- function(config_file, return_object = FALSE){ setwd(yaml_file_path) on.exit(setwd(old_wd)) - # check random seed - if(!is.null(yaml_data$simulation_config$random_seed)){ - random_seed <- yaml_data$simulation_config$random_seed - set.seed(random_seed) - } else { - random_seed <- NULL - } - # ===================================================== # read mandatory parameters is_restore <- !is.null(yaml_data$simulation_config$restore_from) @@ -169,6 +166,43 @@ parse_config <- function(config_file, return_object = FALSE){ sim_length <- yaml_data$simulation_config$length nsim <- ifelse(!is.null(yaml_data$simulation_config$nsim), yaml_data$simulation_config$nsim, 1) + nrep <- ifelse(!is.null(yaml_data$simulation_config$nrep), + yaml_data$simulation_config$nrep, 1) + nrep <- suppressWarnings(as.integer(nrep)[1]) + if (is.na(nrep) || nrep < 1) { + setwd(old_wd) + stop("nrep must be a positive integer") + } + simulation_mode <- "deterministic" + if (!is.null(yaml_data$simulation_config$simulation_mode)) { + simulation_mode <- tolower(trimws(as.character(yaml_data$simulation_config$simulation_mode))) + } else if (!is.null(yaml_data$simulation_config$is_stoch)) { + is_stoch <- yaml_data$simulation_config$is_stoch + simulation_mode <- if (isTRUE(is_stoch) || (is.numeric(is_stoch) && is_stoch != 0)) { + "stochastic" + } else { + "deterministic" + } + } + if (!simulation_mode %in% c("deterministic", "stochastic")) { + setwd(old_wd) + stop("simulation_mode must be either 'deterministic' or 'stochastic'") + } + + # check random seed + if(!is.null(yaml_data$simulation_config$random_seed)){ + random_seed <- suppressWarnings(as.integer(yaml_data$simulation_config$random_seed)[1]) + if (is.na(random_seed)) { + setwd(old_wd) + stop("random_seed must be coercible to a single integer value") + } + set.seed(random_seed) + } else if (simulation_mode == "stochastic") { + random_seed <- sample.int(.Machine$integer.max, 1) + set.seed(random_seed) + } else { + random_seed <- NULL + } chk_time_steps <- NULL chk_file_names <- NULL @@ -192,7 +226,7 @@ parse_config <- function(config_file, return_object = FALSE){ # Generate a matrix of file names: rows for instances, columns for dates chk_file_names <- sapply(chk_dates, function(date) { date_str <- format(date, "%Y-%m-%d") - paste0(checkpoint_dir, "/checkpoint_", date_str, "_instance_", 1:nsim, ".Rda") + paste0(checkpoint_dir, "/checkpoint_", date_str, "_instance_", 1:(nsim * nrep), ".Rda") }) if (is.vector(chk_file_names)) { chk_file_names <- matrix(chk_file_names, ncol = 1) @@ -207,7 +241,7 @@ parse_config <- function(config_file, return_object = FALSE){ chk_time_steps <- as.integer(sim_length / delta_t) # Create a 1-column matrix for the single checkpoint date - chk_file_names <- matrix(paste0(checkpoint_dir, "/chk_", date_str, "_", 1:nsim, ".Rda"), ncol = 1) + chk_file_names <- matrix(paste0(checkpoint_dir, "/chk_", date_str, "_", 1:(nsim * nrep), ".Rda"), ncol = 1) } } @@ -576,6 +610,8 @@ parse_config <- function(config_file, return_object = FALSE){ start_date = start_date, sim_length = sim_length, nsim = nsim, + nrep = nrep, + simulation_mode = simulation_mode, random_seed = random_seed, delta_t = delta_t, chk_file_names = chk_file_names, diff --git a/README.Rmd b/README.Rmd index 37f70b0..9a6df47 100644 --- a/README.Rmd +++ b/README.Rmd @@ -88,6 +88,9 @@ simulation_config: start_date: 10/01/2023 # m/d/Y length: 150 nsim: 1 + nrep: 1 + simulation_mode: deterministic + random_seed: 42 ``` @@ -123,6 +126,13 @@ At the start of a simulation, nearly all individuals are in `S`, with optional s Transmission is stratified by user-defined demographic groups (e.g., age, zone, race). Time-varying mixing matrices define how these strata interact (daytime vs. nighttime, weekday vs. weekend), and `MetaRVM` computes stratum-specific forces of infection for susceptible and vaccinated individuals. Hospitalized and deceased individuals are excluded from the “effective” mixing population. The same model can be run in deterministic or stochastic mode, and parameters are supplied through a YAML configuration. +For simulation control: +- `nsim` is the number of parameter sets (rows in sampled parameter matrices when distributions are used) +- `nrep` is the number of simulation replicates per parameter set +- total runs = `nsim * nrep` +- `simulation_mode` chooses `"deterministic"` or `"stochastic"` +- `random_seed` ensures reproducibility for both parameter sampling and stochastic trajectories + --- ### Core disease progression parameters @@ -209,4 +219,3 @@ List of vignettes: For the complete function reference, visit: **https://resume-epi.github.io/MetaRVM/reference/** - diff --git a/README.md b/README.md index 8f6f4e8..666b59f 100644 --- a/README.md +++ b/README.md @@ -95,6 +95,9 @@ simulation_config: start_date: 10/01/2023 # m/d/Y length: 150 nsim: 1 + nrep: 1 + simulation_mode: deterministic + random_seed: 42 ``` ``` r @@ -150,6 +153,13 @@ excluded from the “effective” mixing population. The same model can be run in deterministic or stochastic mode, and parameters are supplied through a YAML configuration. +For simulation control: - `nsim` is the number of parameter sets (rows +in sampled parameter matrices when distributions are used) - `nrep` is +the number of simulation replicates per parameter set - total runs = +`nsim * nrep` - `simulation_mode` chooses `"deterministic"` or +`"stochastic"` - `random_seed` ensures reproducibility for both +parameter sampling and stochastic trajectories + ------------------------------------------------------------------------ ### Core disease progression parameters diff --git a/inst/extdata/example_config.yaml b/inst/extdata/example_config.yaml index 7117481..1ad7947 100644 --- a/inst/extdata/example_config.yaml +++ b/inst/extdata/example_config.yaml @@ -24,3 +24,4 @@ simulation_config: start_date: 10/01/2023 # m/d/Y length: 150 nsim: 1 + simulation_mode: deterministic diff --git a/inst/extdata/example_config_checkpoint.yaml b/inst/extdata/example_config_checkpoint.yaml index c0d1f87..667e99d 100644 --- a/inst/extdata/example_config_checkpoint.yaml +++ b/inst/extdata/example_config_checkpoint.yaml @@ -24,4 +24,5 @@ simulation_config: start_date: 10/01/2024 # m/d/Y length: 90 nsim: 1 + simulation_mode: deterministic # The checkpoint_dir will be set dynamically in the vignette diff --git a/inst/extdata/example_config_dist.yaml b/inst/extdata/example_config_dist.yaml index d3de7d3..0e76626 100644 --- a/inst/extdata/example_config_dist.yaml +++ b/inst/extdata/example_config_dist.yaml @@ -36,3 +36,4 @@ simulation_config: start_date: 01/01/2023 # m/d/Y length: 150 nsim: 20 # Increased nsim for meaningful summary statistics + simulation_mode: deterministic diff --git a/inst/extdata/example_config_stochastic.yaml b/inst/extdata/example_config_stochastic.yaml new file mode 100644 index 0000000..816b50c --- /dev/null +++ b/inst/extdata/example_config_stochastic.yaml @@ -0,0 +1,29 @@ +run_id: ExampleRun_Stochastic_Static +population_data: + initialization: population_init_n24.csv + vaccination: vaccination_n24.csv +mixing_matrix: + weekday_day: m_weekday_day.csv + weekday_night: m_weekday_night.csv + weekend_day: m_weekend_day.csv + weekend_night: m_weekend_night.csv +disease_params: + ts: 0.5 + ve: 0.4 + dv: 180 + dp: 1 + de: 3 + da: 5 + ds: 6 + dh: 8 + dr: 180 + pea: 0.3 + psr: 0.95 + phr: 0.97 +simulation_config: + start_date: 01/01/2023 # m/d/Y + length: 150 + nsim: 1 + nrep: 5 + simulation_mode: stochastic + random_seed: 42 \ No newline at end of file diff --git a/inst/extdata/example_config_subgroup_dist.yaml b/inst/extdata/example_config_subgroup_dist.yaml index 9f60cf3..1be0332 100644 --- a/inst/extdata/example_config_subgroup_dist.yaml +++ b/inst/extdata/example_config_subgroup_dist.yaml @@ -40,3 +40,4 @@ simulation_config: start_date: 01/01/2023 # m/d/Y length: 150 nsim: 20 + simulation_mode: deterministic diff --git a/man/parse_config.Rd b/man/parse_config.Rd index 58df26b..ac9e36a 100644 --- a/man/parse_config.Rd +++ b/man/parse_config.Rd @@ -26,6 +26,8 @@ If \code{return_object = FALSE} (default), returns a named list containing: \item{start_date}{Simulation start date as Date object} \item{sim_length}{Simulation length in days} \item{nsim}{Number of simulation instances} +\item{nrep}{Number of stochastic replicates per parameter set} +\item{simulation_mode}{Simulation mode: \code{"deterministic"} or \code{"stochastic"}} \item{random_seed}{Random seed used (if any)} \item{delta_t}{Time step size (fixed at 0.5)} \item{chk_file_names, chk_time_steps, do_chk}{Checkpointing configuration} @@ -47,6 +49,9 @@ main sections: \itemize{ \item \code{random_seed}: Optional random seed for reproducibility in case of stochastic simulations or stochastic parameters \item \code{nsim}: Number of simulation instances (default: 1) +\item \code{nrep}: Number of stochastic replicates per parameter set (default: 1) +\item \code{simulation_mode}: Optional simulation mode. Must be one of +\code{"deterministic"} or \code{"stochastic"} (default: \code{"deterministic"}). \item \code{start_date}: Simulation start date in MM/DD/YYYY format \item \code{length}: Simulation length in days \item \code{checkpoint_dir}: Optional checkpoint directory for saving intermediate results diff --git a/tests/testthat/test-metaRVM.R b/tests/testthat/test-metaRVM.R index 912859a..dbb1a2d 100644 --- a/tests/testthat/test-metaRVM.R +++ b/tests/testthat/test-metaRVM.R @@ -107,3 +107,115 @@ test_that("subset_data reports valid values when category value is invalid", { regexp = "Invalid values for category 'age': NOT_A_VALID_AGE\\. Valid values are:" ) }) + +test_that("metaRVM uses provided random_seed for stochastic runs", { + ext <- function(x) system.file("extdata", x, package = "MetaRVM") + cfg_path <- tempfile(fileext = ".yaml") + + cfg <- list( + run_id = "StochasticSeedProvided", + population_data = list( + initialization = ext("population_init_n24.csv"), + vaccination = ext("vaccination_n24.csv") + ), + mixing_matrix = list( + weekday_day = ext("m_weekday_day.csv"), + weekday_night = ext("m_weekday_night.csv"), + weekend_day = ext("m_weekend_day.csv"), + weekend_night = ext("m_weekend_night.csv") + ), + disease_params = list( + ts = 0.5, ve = 0.4, dv = 180, dp = 1, de = 3, + da = 5, ds = 6, dh = 8, dr = 180, pea = 0.3, psr = 0.95, phr = 0.97 + ), + simulation_config = list( + start_date = "10/01/2023", + length = 20, + nsim = 1, + simulation_mode = "stochastic", + random_seed = 12345 + ) + ) + + yaml::write_yaml(cfg, cfg_path) + + res1 <- metaRVM(cfg_path) + res2 <- metaRVM(cfg_path) + + expect_equal(res1$config$get("random_seed"), 12345L) + expect_equal(res2$config$get("random_seed"), 12345L) + expect_equal(res1$results$value, res2$results$value) +}) + +test_that("metaRVM generates and stores random_seed for stochastic runs when missing", { + ext <- function(x) system.file("extdata", x, package = "MetaRVM") + cfg_path <- tempfile(fileext = ".yaml") + + cfg <- list( + run_id = "StochasticSeedGenerated", + population_data = list( + initialization = ext("population_init_n24.csv"), + vaccination = ext("vaccination_n24.csv") + ), + mixing_matrix = list( + weekday_day = ext("m_weekday_day.csv"), + weekday_night = ext("m_weekday_night.csv"), + weekend_day = ext("m_weekend_day.csv"), + weekend_night = ext("m_weekend_night.csv") + ), + disease_params = list( + ts = 0.5, ve = 0.4, dv = 180, dp = 1, de = 3, + da = 5, ds = 6, dh = 8, dr = 180, pea = 0.3, psr = 0.95, phr = 0.97 + ), + simulation_config = list( + start_date = "10/01/2023", + length = 20, + nsim = 1, + simulation_mode = "stochastic" + ) + ) + + yaml::write_yaml(cfg, cfg_path) + + res <- metaRVM(cfg_path) + generated_seed <- res$config$get("random_seed") + + expect_true(!is.null(generated_seed)) + expect_true(is.numeric(generated_seed)) + expect_length(generated_seed, 1) +}) + +test_that("metaRVM runs nsim * nrep instances", { + ext <- function(x) system.file("extdata", x, package = "MetaRVM") + cfg_path <- tempfile(fileext = ".yaml") + + cfg <- list( + run_id = "ReplicateCount", + population_data = list( + initialization = ext("population_init_n24.csv"), + vaccination = ext("vaccination_n24.csv") + ), + mixing_matrix = list( + weekday_day = ext("m_weekday_day.csv"), + weekday_night = ext("m_weekday_night.csv"), + weekend_day = ext("m_weekend_day.csv"), + weekend_night = ext("m_weekend_night.csv") + ), + disease_params = list( + ts = 0.5, ve = 0.4, dv = 180, dp = 1, de = 3, + da = 5, ds = 6, dh = 8, dr = 180, pea = 0.3, psr = 0.95, phr = 0.97 + ), + simulation_config = list( + start_date = "10/01/2023", + length = 20, + nsim = 2, + nrep = 3, + simulation_mode = "deterministic" + ) + ) + + yaml::write_yaml(cfg, cfg_path) + + res <- metaRVM(cfg_path) + expect_equal(length(unique(res$results$instance)), 6L) +}) diff --git a/tests/testthat/test-parse_config.R b/tests/testthat/test-parse_config.R index 25c568f..9aaaea0 100644 --- a/tests/testthat/test-parse_config.R +++ b/tests/testthat/test-parse_config.R @@ -15,7 +15,7 @@ test_that("parse_config returns expected list structure for example config", { "m_wd_d", "m_wd_n", "m_we_d", "m_we_n", "ts", "ve", "dv", "de", "dp", "da", "ds", "dh", "dr", "pea", "psr", "phr", - "start_date", "sim_length", "nsim", "random_seed", + "start_date", "sim_length", "nsim", "nrep", "simulation_mode", "random_seed", "delta_t", "chk_file_names", "chk_time_steps", "do_chk" ) @@ -56,6 +56,8 @@ test_that("parse_config returns expected list structure for example config", { expect_true(cfg$sim_length > 0) expect_true(is.numeric(cfg$delta_t)) expect_equal(cfg$delta_t, 0.5) + expect_equal(cfg$nrep, 1L) + expect_true(cfg$simulation_mode %in% c("deterministic", "stochastic")) }) test_that("parse_config expands disease parameters with correct dimensions and ranges", { @@ -158,3 +160,91 @@ test_that("parse_config reports valid category values when sub_disease_params va regexp = "Invalid values for category 'age' in sub_disease_params: INVALID_AGE_GROUP\\. Valid values are:" ) }) + +test_that("parse_config accepts simulation_mode and defaults to deterministic", { + ext <- function(x) system.file("extdata", x, package = "MetaRVM") + + cfg_default_path <- tempfile(fileext = ".yaml") + cfg_stoch_path <- tempfile(fileext = ".yaml") + cfg_invalid_path <- tempfile(fileext = ".yaml") + + base_cfg <- list( + run_id = "SimulationMode", + population_data = list( + initialization = ext("population_init_n24.csv"), + vaccination = ext("vaccination_n24.csv") + ), + mixing_matrix = list( + weekday_day = ext("m_weekday_day.csv"), + weekday_night = ext("m_weekday_night.csv"), + weekend_day = ext("m_weekend_day.csv"), + weekend_night = ext("m_weekend_night.csv") + ), + disease_params = list( + ts = 0.5, ve = 0.4, dv = 180, dp = 1, de = 3, + da = 5, ds = 6, dh = 8, dr = 180, pea = 0.3, psr = 0.95, phr = 0.97 + ), + simulation_config = list( + start_date = "10/01/2023", + length = 30, + nsim = 1 + ) + ) + + yaml::write_yaml(base_cfg, cfg_default_path) + expect_equal(parse_config(cfg_default_path)$simulation_mode, "deterministic") + + stoch_cfg <- base_cfg + stoch_cfg$simulation_config$simulation_mode <- "stochastic" + yaml::write_yaml(stoch_cfg, cfg_stoch_path) + expect_equal(parse_config(cfg_stoch_path)$simulation_mode, "stochastic") + + invalid_cfg <- base_cfg + invalid_cfg$simulation_config$simulation_mode <- "invalid_mode" + yaml::write_yaml(invalid_cfg, cfg_invalid_path) + expect_error( + parse_config(cfg_invalid_path), + regexp = "simulation_mode must be either 'deterministic' or 'stochastic'" + ) +}) + +test_that("parse_config accepts nrep and validates it", { + ext <- function(x) system.file("extdata", x, package = "MetaRVM") + cfg_ok_path <- tempfile(fileext = ".yaml") + cfg_bad_path <- tempfile(fileext = ".yaml") + + base_cfg <- list( + run_id = "Replicates", + population_data = list( + initialization = ext("population_init_n24.csv"), + vaccination = ext("vaccination_n24.csv") + ), + mixing_matrix = list( + weekday_day = ext("m_weekday_day.csv"), + weekday_night = ext("m_weekday_night.csv"), + weekend_day = ext("m_weekend_day.csv"), + weekend_night = ext("m_weekend_night.csv") + ), + disease_params = list( + ts = 0.5, ve = 0.4, dv = 180, dp = 1, de = 3, + da = 5, ds = 6, dh = 8, dr = 180, pea = 0.3, psr = 0.95, phr = 0.97 + ), + simulation_config = list( + start_date = "10/01/2023", + length = 30, + nsim = 2, + nrep = 3 + ) + ) + + yaml::write_yaml(base_cfg, cfg_ok_path) + expect_equal(parse_config(cfg_ok_path)$nrep, 3L) + + bad_cfg <- base_cfg + bad_cfg$simulation_config$nrep <- 0 + yaml::write_yaml(bad_cfg, cfg_bad_path) + expect_error( + parse_config(cfg_bad_path), + regexp = "nrep must be a positive integer" + ) +}) diff --git a/vignettes/checkpointing.Rmd b/vignettes/checkpointing.Rmd index 4d592b6..beea3b9 100644 --- a/vignettes/checkpointing.Rmd +++ b/vignettes/checkpointing.Rmd @@ -32,7 +32,7 @@ The checkpointing system uses three main configuration options in the `simulatio ## Example -For the purpose of this example, we will set the checkpointed directory to be a temporary directory, and will create a temporary yaml configuration file based on a template. +For this example, the checkpoint directory is set to a temporary directory, and a temporary YAML configuration file is created from a template. ### Set Up @@ -63,6 +63,9 @@ yml_tmp <- data.table::copy(yml) # Assign checkpoint directory yml_tmp$simulation_config$checkpoint_dir <- checkpoint_dir +yml_tmp$simulation_config$nrep <- 1 +yml_tmp$simulation_config$simulation_mode <- "deterministic" +yml_tmp$simulation_config$random_seed <- 42 # Create a temporary config temp_config <- tempfile(tmpdir = dirname(example_config), fileext = ".yaml") @@ -122,7 +125,7 @@ print(head(results$results, 10)) ### Resume from a Checkpoint -We'll create a new configuration that restores from the checkpoint. We need to provide the checkpoint file name, and a start date which should be the next day of the end date in the previous run. +A new configuration is created to restore from the checkpoint. The checkpoint file name must be provided, along with a start date that should be the day after the end date in the previous run. ``` r @@ -139,8 +142,11 @@ if (length(checkpoint_files_full) > 0) { new_start_date <- "12/30/2024" yml_tmp <- data.table::copy(yml) yml_tmp$simulation_config$start_date <- new_start_date + yml_tmp$simulation_config$nrep <- 1 + yml_tmp$simulation_config$simulation_mode <- "deterministic" + yml_tmp$simulation_config$random_seed <- 42 yml_tmp$simulation_config$restore_from <- checkpoint_to_restore - yml_tmp$population_data$initialization <- "population_init_n24_only_mapping.csv" # ensure that we don't want to reinitialize the population with the original initial conditions + yml_tmp$population_data$initialization <- "population_init_n24_only_mapping.csv" # ensure that the original initial conditions are not reapplied during restore temp_config_resume <- tempfile(tmpdir = dirname(example_config), fileext = ".yaml") yaml::write_yaml(yml_tmp, temp_config_resume) @@ -197,5 +203,3 @@ Where: For more information on configuring MetaRVM simulations, see the [YAML Configuration](yaml-configuration.html) vignette. - - diff --git a/vignettes/getting-started.Rmd b/vignettes/getting-started.Rmd index 4b75b7c..3cfbfed 100644 --- a/vignettes/getting-started.Rmd +++ b/vignettes/getting-started.Rmd @@ -18,11 +18,11 @@ knitr::opts_chunk$set( ## Introduction -MetaRVM is a comprehensive R package for simulating respiratory virus epidemics using meta-population compartmental models. This vignette will guide you through the basic usage of the package. +MetaRVM is a comprehensive R package for simulating respiratory virus epidemics using meta-population compartmental models. This vignette describes the basic usage of the package. ## Installation -You can install the development version of MetaRVM from GitHub: +The development version of MetaRVM can be installed from GitHub: ```{r eval=FALSE} # install.packages("devtools") @@ -41,7 +41,7 @@ options(odin.verbose = FALSE) This example shows how to run a basic meta-population simulation. -The `metaRVM` package includes a set of example files in its `extdata` directory. To run the example, we first need to locate these files. The `system.file()` function in R is the recommended way to do this, as it will find the files wherever the package is installed. +The `metaRVM` package includes a set of example files in its `extdata` directory. To run the example, these files must first be located. The `system.file()` function in R is the recommended way to do this, as it finds the files wherever the package is installed. ```{r} # Locate the example YAML configuration file @@ -78,13 +78,16 @@ simulation_config: start_date: 01/01/2023 # m/d/Y length: 150 nsim: 1 + nrep: 1 + simulation_mode: deterministic + random_seed: 42 ``` For a detailed explanation of all the configuration options, please see the `yaml-configuration.html` vignette. ## Running the Simulation -Once we have the path to the configuration file, the simulation can be run using the `metaRVM()` function. +Once the path to the configuration file is available, the simulation can be run using the `metaRVM()` function. ```{r, results='hide'} # Load the metaRVM library diff --git a/vignettes/running-a-simulation.Rmd b/vignettes/running-a-simulation.Rmd index 8bd209e..747388c 100644 --- a/vignettes/running-a-simulation.Rmd +++ b/vignettes/running-a-simulation.Rmd @@ -20,7 +20,7 @@ This vignette demonstrates how to run a `metaRVM` simulation using the example c ## Locating the Example Files -The `metaRVM` package includes a set of example files in its `extdata` directory. To run the example, we first need to locate these files. The `system.file()` function in R is the recommended way to do this, as it will find the files wherever the package is installed. +The `metaRVM` package includes a set of example files in its `extdata` directory. To run the example, these files must first be located. The `system.file()` function in R is the recommended way to do this, as it finds the files wherever the package is installed. ```{r} # Locate the example YAML configuration file @@ -57,11 +57,14 @@ simulation_config: start_date: 01/01/2023 # m/d/Y length: 150 nsim: 1 + nrep: 1 + simulation_mode: deterministic + random_seed: 42 ``` ## Running the Simulation -Once we have the path to the configuration file, the simulation can be run using the `metaRVM()` function. +Once the path to the configuration file is available, the simulation can be run using the `metaRVM()` function. ```{r, results='hide'} # Load the metaRVM library @@ -221,18 +224,22 @@ simulation_config: start_date: 01/01/2023 # m/d/Y length: 150 nsim: 20 # Increased nsim for meaningful summary statistics + nrep: 1 + simulation_mode: deterministic + random_seed: 42 ``` -To run a simulation with this configuration, we pass the file path to `metaRVM`. +To run a simulation with this configuration, the file path is passed to `metaRVM`. ```{r, results='hide', message=FALSE, warning=FALSE} # Run the simulation with the new configuration sim_out_dist <- metaRVM(yaml_file_dist) ``` + ## Generating Summary Statistics across Demographics -The `MetaRVMResults` class provides basic summarization functionality across multiple instances of the simulation, when one or more disease parameters are specified via distribution, and there are more than one simulations per configurations. The `summarize` method generates output of class `MetaRVMSummary` which has a `plot` method available to use. Now that we have run a simulation with parameter distributions, we can use the `summarize` method to see the variability in the results. +The `MetaRVMResults` class provides basic summarization functionality across multiple instances of the simulation, when one or more disease parameters are specified via distribution, and there are more than one simulations per configurations. The `summarize` method generates output of class `MetaRVMSummary` which has a `plot` method available. After a simulation is run with parameter distributions, the `summarize` method can be used to inspect variability in the results. ```{r, fig.height = 4, fig.width = 8, fig.align = "center"} library(ggplot2) @@ -265,6 +272,53 @@ hospital_summary hospital_summary$plot() + ggtitle("Daily Hospitalizations by Age and Zone") + theme_bw() ``` +# Running a Stochastic Simulation with Static Parameters + +A stochastic simulation can also be run by setting `simulation_mode: stochastic`. + +An example YAML file with parameter distributions is included in the package, `example_config_stochastic.yaml`. Here is its content: + +```{r} +# Locate the example YAML configuration file with distributions +yaml_file_dist <- system.file("extdata", "example_config_dist.yaml", package = "MetaRVM") +``` + +```yaml +run_id: ExampleRun_Stochastic_Static +population_data: + initialization: population_init_n24.csv + vaccination: vaccination_n24.csv +mixing_matrix: + weekday_day: m_weekday_day.csv + weekday_night: m_weekday_night.csv + weekend_day: m_weekend_day.csv + weekend_night: m_weekend_night.csv +disease_params: + ts: 0.5 + ve: 0.4 + dv: 180 + dp: 1 + de: 3 + da: 5 + ds: 6 + dh: 8 + dr: 180 + pea: 0.3 + psr: 0.95 + phr: 0.97 +simulation_config: + start_date: 01/01/2023 # m/d/Y + length: 150 + nsim: 1 + nrep: 5 + simulation_mode: stochastic + random_seed: 42 +``` + +```{r, results='hide', message=FALSE, warning=FALSE} +sim_out_stoch <- metaRVM(yaml_file_stoch) +``` + # Specifying Disease Parameters by Demographics The disease parameters can also be specified for different demographic subgroups. These subgroup-specific parameters will override the global parameters. For more details, refer to the `yaml-configuration` vignette. An example YAML file is provided, `example_config_subgroup_dist.yaml`, that demonstrates this feature. It also includes parameters defined by distributions. @@ -317,6 +371,9 @@ simulation_config: start_date: 01/01/2023 # m/d/Y length: 150 nsim: 20 + nrep: 1 + simulation_mode: deterministic + random_seed: 42 ``` Now, let's run the simulation with this configuration. @@ -326,7 +383,7 @@ Now, let's run the simulation with this configuration. sim_out_subgroup <- metaRVM(yaml_file_subgroup) ``` -We can now plot the results to see the impact of the subgroup-specific parameters. For example, we can compare the number of hospitalizations in the "65+" age group, which has a `dh` of 10, to other age groups that use the global `dh` drawn from a lognormal distribution. +The results can now be plotted to evaluate the impact of subgroup-specific parameters. For example, the number of hospitalizations in the "65+" age group, which has a `dh` of 10, can be compared to other age groups that use the global `dh` drawn from a lognormal distribution. ```{r, fig.height = 6, fig.width = 8, fig.align = "center"} # Summarize hospitalizations by age group diff --git a/vignettes/yaml-configuration.Rmd b/vignettes/yaml-configuration.Rmd index 4a7ae93..8595171 100644 --- a/vignettes/yaml-configuration.Rmd +++ b/vignettes/yaml-configuration.Rmd @@ -49,12 +49,14 @@ simulation_config: start_date: 01/01/2025 # m/d/Y length: 90 nsim: 1 + nrep: 1 + simulation_mode: deterministic random_seed: 42 ``` ### Configuration Sections -- **`run_id`**: A unique name for your simulation. +- **`run_id`**: A unique name for the simulation. - **`population_data`**: Paths to CSV files for population demographics, initial state, and vaccination schedules. - **`mixing_matrix`**: Paths to CSV files defining contact patterns for different times of the week. - **`disease_params`**: Disease characteristics. In this example, all parameters are single, fixed values. @@ -139,7 +141,7 @@ Below is a list of the disease parameters used in `metaRVM`: ## Defining Parameters with Distributions -Instead of fixed values, you can define disease parameters using statistical distributions. This is useful for capturing uncertainty in the parameters. `metaRVM` supports `uniform` and `lognormal` distributions. +Instead of fixed values, disease parameters can be defined using statistical distributions. This is useful for capturing uncertainty in the parameters. `metaRVM` supports `uniform` and `lognormal` distributions. Here is an example of defining `ve`, `da`, `ds`, and `dh` with distributions: @@ -171,14 +173,14 @@ disease_params: phr: 0.97 ``` -- For a `uniform` distribution, you must specify `min` and `max` values. -- For a `lognormal` distribution, you must specify `mu` and `sd` (mean and standard deviation on the log scale). +- For a `uniform` distribution, `min` and `max` values must be specified. +- For a `lognormal` distribution, `mu` and `sd` (mean and standard deviation on the log scale) must be specified. ## Specifying Subgroup Parameters -`metaRVM` allows you to specify different disease parameters for various demographic subgroups using the `sub_disease_params` section. These subgroup-specific parameters will override the global parameters defined in `disease_params`. +`metaRVM` allows different disease parameters to be specified for demographic subgroups using the `sub_disease_params` section. These subgroup-specific parameters override the global parameters defined in `disease_params`. -The demographic categories used in this section must match the user-defined category column names in the initialization CSV file specified under `population_data`. For example, if your initialization file has columns named `age`, `income_level`, and `occupation`, you can use any of these categories in `sub_disease_params`. The specific values (e.g., `"0-4"`, `"low"`, `"healthcare"`) must exactly match the values in those columns. +The demographic categories used in this section must match the user-defined category column names in the initialization CSV file specified under `population_data`. For example, if the initialization file has columns named `age`, `income_level`, and `occupation`, any of these categories can be used in `sub_disease_params`. The specific values (e.g., `"0-4"`, `"low"`, `"healthcare"`) must exactly match the values in those columns. The following example defines different parameters for different age groups: @@ -204,13 +206,73 @@ sub_disease_params: In this configuration, individuals in the "0-4" age group will have a `dh` (duration of hospitalization) of 4, overriding any global `dh` value. Similarly, the transmission rate `ts` for the "18-49" group is set to 0.01. +## Stochastic Simulation with Distributional Parameters + +When both parameter uncertainty and stochastic disease transitions are represented, +set `simulation_mode: stochastic` and define one or more disease +parameters as distributions. + +- `nsim`: number of sampled parameter sets +- `nrep`: number of stochastic replicates per parameter set +- total runs = `nsim * nrep` + +Example: + +```yaml +run_id: StochasticDistRun +population_data: + initialization: data/population_init.csv + vaccination: data/vaccination.csv +mixing_matrix: + weekday_day: data/m_weekday_day.csv + weekday_night: data/m_weekday_night.csv + weekend_day: data/m_weekend_day.csv + weekend_night: data/m_weekend_night.csv +disease_params: + ts: 0.5 + ve: + dist: uniform + min: 0.29 + max: 0.53 + dv: 158 + dp: 1 + de: 3 + da: + dist: uniform + min: 3 + max: 7 + ds: + dist: uniform + min: 5 + max: 7 + dh: + dist: lognormal + mu: 2.0 + sd: 0.5 + dr: 187 + pea: 0.333 + psr: 0.95 + phr: 0.97 +simulation_config: + start_date: 01/01/2025 + length: 90 + nsim: 20 + nrep: 5 + simulation_mode: stochastic + random_seed: 42 +``` + +For reproducibility, provide `random_seed`. This seed is used to reproduce both +the parameter draws (for distributional parameters) and the stochastic model +replicates. + ## Checkpointing and Restoring Simulations -For long-running simulations, it is useful to save the state of the model at intermediate points. This is known as checkpointing. `metaRVM` allows you to save checkpoints and restore a simulation from a saved state. +For long-running simulations, it is useful to save the state of the model at intermediate points. This is known as checkpointing. `metaRVM` allows checkpoints to be saved and simulations to be restored from a saved state. ### Enabling Checkpointing -To enable checkpointing, you need to add the `checkpoint_dir` and optionally `checkpoint_dates` to the `simulation_config` section of your YAML file. +To enable checkpointing, `checkpoint_dir` and optionally `checkpoint_dates` need to be added to the `simulation_config` section of the YAML file. - `checkpoint_dir`: The directory where checkpoint files will be saved. - `checkpoint_dates`: A list of dates (in `MM/DD/YYYY` format) on which to save a checkpoint. If this is not provided, a single checkpoint will be saved at the end of the simulation. @@ -222,6 +284,8 @@ simulation_config: start_date: 01/01/2025 length: 90 nsim: 10 + nrep: 1 + simulation_mode: deterministic random_seed: 42 checkpoint_dir: "path/to/checkpoints" checkpoint_dates: ["01/15/2025", "01/30/2025"] @@ -229,13 +293,16 @@ simulation_config: ### Restoring from a Checkpoint -To restore a simulation from a checkpoint file, use the `restore_from` parameter in the `simulation_config` section. This will initialize the model with the state saved in the specified checkpoint file. +To restore a simulation from a checkpoint file, the `restore_from` parameter is used in the `simulation_config` section. The model is initialized with the state saved in the specified checkpoint file. ```yaml simulation_config: start_date: 01/30/2025 # Should be the next date of the checkpoint date length: 60 # Remaining simulation length nsim: 10 + nrep: 1 + simulation_mode: deterministic + random_seed: 42 restore_from: "path/to/checkpoints/checkpoint_2025-01-30_instance_1.Rda" ``` From 679ca26d6949f9bca6944cc89ec8677cb92c2078 Mon Sep 17 00:00:00 2001 From: arindam fadikar Date: Thu, 26 Mar 2026 11:14:14 -0400 Subject: [PATCH 2/5] changed matrix indices in distributing one subpopulation to all --- R/classes.R | 184 +++++++++++++++++++++++------------- R/metaODIN.R | 16 +++- man/MetaRVMResults.Rd | 50 ++++++++++ vignettes/checkpointing.Rmd | 20 ++-- 4 files changed, 187 insertions(+), 83 deletions(-) diff --git a/R/classes.R b/R/classes.R index a2843ba..b7b6deb 100644 --- a/R/classes.R +++ b/R/classes.R @@ -359,6 +359,40 @@ MetaRVMResults <- R6::R6Class( invisible(self) }, + #' @description Plot simulation trajectories directly from results + #' @param group_by Vector of demographic category names to group/facet by + #' @param disease_states Optional disease states to include + #' @param date_range Optional date range for filtering + #' @param instances Optional instance IDs to include + #' @param stats Statistics for summary plotting. If NULL, plots raw trajectories. + #' @param quantiles Quantiles for uncertainty bands when summary plotting + #' @param exclude_p_columns Logical, whether to exclude p_ columns (default: TRUE) + #' @param ci_level Confidence level label used in summary plot title + #' @param theme ggplot2 theme function (default: theme_minimal()) + #' @param title Optional custom title + #' @return ggplot object + plot = function(group_by = character(0), disease_states = NULL, date_range = NULL, + instances = NULL, stats = NULL, quantiles = c(0.25, 0.75), + exclude_p_columns = TRUE, ci_level = 0.95, + theme = theme_minimal(), title = NULL) { + + plot_data <- self$subset_data( + disease_states = disease_states, + date_range = date_range, + instances = instances, + exclude_p_columns = exclude_p_columns + ) + + summary_obj <- plot_data$summarize( + group_by = group_by, + stats = stats, + quantiles = quantiles, + exclude_p_columns = FALSE + ) + + summary_obj$plot(ci_level = ci_level, theme = theme, title = title) + }, + #' @description Subset the data based on any combination of parameters #' @param ... Named arguments for category filters (e.g., age = c("0-17"), income = c("low", "high")) #' @param disease_states Vector of disease states to include (default: all, excludes p_ columns) @@ -529,60 +563,45 @@ MetaRVMResults <- R6::R6Class( # Step 3: Calculate statistics across instances for each date/demographic/disease combination final_group_vars <- c("date", group_by, "disease_state") - - # Start with empty result - summary_result <- summed_data[, .(temp = mean(value)), by = final_group_vars][, temp := NULL] - - # Add each statistic individually - if ("mean" %in% stats) { - temp_mean <- summed_data[, .(mean_value = mean(value, na.rm = TRUE)), by = final_group_vars] - summary_result <- merge(summary_result, temp_mean, by = final_group_vars, all = TRUE) - } - - if ("median" %in% stats) { - temp_median <- summed_data[, .(median_value = median(value, na.rm = TRUE)), by = final_group_vars] - summary_result <- merge(summary_result, temp_median, by = final_group_vars, all = TRUE) - } - - if ("sd" %in% stats) { - temp_sd <- summed_data[, .(sd_value = sd(value, na.rm = TRUE)), by = final_group_vars] - summary_result <- merge(summary_result, temp_sd, by = final_group_vars, all = TRUE) - } - - if ("min" %in% stats) { - temp_min <- summed_data[, .(min_value = min(value, na.rm = TRUE)), by = final_group_vars] - summary_result <- merge(summary_result, temp_min, by = final_group_vars, all = TRUE) - } - - if ("max" %in% stats) { - temp_max <- summed_data[, .(max_value = max(value, na.rm = TRUE)), by = final_group_vars] - summary_result <- merge(summary_result, temp_max, by = final_group_vars, all = TRUE) - } - - if ("sum" %in% stats) { - temp_sum <- summed_data[, .(sum_value = sum(value, na.rm = TRUE)), by = final_group_vars] - summary_result <- merge(summary_result, temp_sum, by = final_group_vars, all = TRUE) - } - - # Handle quantiles - if ("quantile" %in% stats) { - for (q in quantiles) { - q_name <- paste0("q", sprintf("%02d", round(q * 100))) - temp_q <- summed_data[, .(temp_quantile = quantile(value, q, na.rm = TRUE)), by = final_group_vars] - setnames(temp_q, "temp_quantile", q_name) - summary_result <- merge(summary_result, temp_q, by = final_group_vars, all = TRUE) + summary_result <- summed_data[, { + out <- list() + + if ("mean" %in% stats) { + out$mean_value <- mean(value, na.rm = TRUE) } - } - - # Sort results for better readability - setorder(summary_result, date) + if ("median" %in% stats) { + out$median_value <- median(value, na.rm = TRUE) + } + if ("sd" %in% stats) { + out$sd_value <- stats::sd(value, na.rm = TRUE) + } + if ("min" %in% stats) { + out$min_value <- min(value, na.rm = TRUE) + } + if ("max" %in% stats) { + out$max_value <- max(value, na.rm = TRUE) + } + if ("sum" %in% stats) { + out$sum_value <- sum(value, na.rm = TRUE) + } + if ("quantile" %in% stats) { + for (q in quantiles) { + q_name <- paste0("q", sprintf("%02d", round(q * 100))) + out[[q_name]] <- stats::quantile(value, probs = q, na.rm = TRUE, names = FALSE) + } + } + + out + }, by = final_group_vars] # return(summary_result) # Return a chainable object if (is.null(stats)) { + setorder(summed_data, date, disease_state, instance) # Return summary object for chaining return(MetaRVMSummary$new(summed_data, self$config, type = "instances")) } else { + setorderv(summary_result, c("date", group_by, "disease_state")) # Return summary object for chaining return(MetaRVMSummary$new(summary_result, self$config, type = "summary")) } @@ -785,18 +804,7 @@ MetaRVMSummary <- R6::R6Class( plot = function(ci_level = 0.95, theme = theme_minimal(), title = NULL) { columns <- names(self$data) - - # Check if we have the required columns for plotting - if (!("median_value" %in% columns)) { - stop("Plot method requires 'median_value' column. Please call summarize() first with stats = c('median', 'quantile')") - } - - # Check for quantile columns - quantile_cols <- grep("^q[0-9]", columns, value = TRUE) - if (length(quantile_cols) < 2) { - stop("Plot method requires quantile columns for confidence bands. Please call summarize() with stats = c('median', 'quantile')") - } - + # Detect grouping variables dynamically available_categories <- self$config$get_category_names() group_vars <- if (length(available_categories) > 0) { @@ -805,12 +813,53 @@ MetaRVMSummary <- R6::R6Class( character(0) } - if (length(group_vars) == 0) { - stop("Data must be grouped by at least one demographic variable") + # Instance trajectory plotting mode + if (identical(self$type, "instances")) { + color_var <- if ("disease_state" %in% columns) "disease_state" else "instance" + + p <- ggplot(self$data, aes(x = date, y = value, + color = rlang::.data[[color_var]], + group = interaction(rlang::.data[["instance"]], rlang::.data[[color_var]], drop = TRUE))) + + geom_line(alpha = 0.45, linewidth = 0.5) + + if (length(group_vars) == 1) { + p <- p + facet_wrap(stats::as.formula(paste("~", group_vars[1])), scales = "free_y") + } else if (length(group_vars) >= 2) { + p <- p + facet_grid(stats::as.formula(paste(group_vars[1], "~", group_vars[2])), scales = "free_y") + } + + p <- p + + labs( + title = title %||% "Simulation Trajectories", + x = "Date", + y = "Value", + color = tools::toTitleCase(gsub("_", " ", color_var)) + ) + + theme + + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + + return(p) + } + + # Check if we have the required columns for plotting + if (!("median_value" %in% columns)) { + stop("Plot method requires 'median_value' column. Please call summarize() first with stats = c('median', 'quantile')") } + # Check for quantile columns + quantile_cols <- grep("^q[0-9]+$", columns, value = TRUE) + quantile_levels <- suppressWarnings(as.numeric(sub("^q", "", quantile_cols))) + valid_q_idx <- which(!is.na(quantile_levels)) + if (length(valid_q_idx) < 2) { + stop("Plot method requires quantile columns for confidence bands. Please call summarize() with stats = c('median', 'quantile')") + } + quantile_cols <- quantile_cols[valid_q_idx][order(quantile_levels[valid_q_idx])] + # Create faceting strategy based on number of grouping variables - if (length(group_vars) == 1) { + if (length(group_vars) == 0) { + facet_formula <- NULL + color_var <- "disease_state" + } else if (length(group_vars) == 1) { # Single grouping variable: facet by that variable, color by disease_state facet_formula <- as.formula(paste("~", group_vars[1])) color_var <- "disease_state" @@ -827,16 +876,15 @@ MetaRVMSummary <- R6::R6Class( color_var <- group_vars[3] } - # Use first and last quantile columns for confidence bands + # Use lowest and highest quantile columns for confidence bands ci_lower_col <- quantile_cols[1] ci_upper_col <- quantile_cols[length(quantile_cols)] # Create the plot - p <- ggplot(self$data, aes(x = date, y = median_value, color = get(color_var))) + + p <- ggplot(self$data, aes(x = date, y = median_value, color = rlang::.data[[color_var]])) + geom_line(linewidth = 1) + - geom_ribbon(aes(ymin = get(ci_lower_col), ymax = get(ci_upper_col), - fill = get(color_var)), alpha = 0.2, color = NA) + - facet_grid(facet_formula, scales = "free_y") + + geom_ribbon(aes(ymin = rlang::.data[[ci_lower_col]], ymax = rlang::.data[[ci_upper_col]], + fill = rlang::.data[[color_var]]), alpha = 0.2, color = NA) + labs( title = title %||% paste0("Median Outcomes with ", ci_level*100, "% Empirical Quantiles"), x = "Date", @@ -846,6 +894,10 @@ MetaRVMSummary <- R6::R6Class( ) + theme + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + + if (!is.null(facet_formula)) { + p <- p + facet_grid(facet_formula, scales = "free_y") + } return(p) } diff --git a/R/metaODIN.R b/R/metaODIN.R index 6cdc634..38e439b 100644 --- a/R/metaODIN.R +++ b/R/metaODIN.R @@ -348,18 +348,19 @@ meta_sim <- function(N_pop, ts, ## ================================================= ## Draws from binomial distributions for numbers changing between ## compartments: - n_SE_eff[, ] <- if(S[i] <= 0) 0 else (if(stoch == 1) rbinom(S_eff_prod[j, i], p_SE[i]) else S_eff_prod[j, i] * p_SE[i]) - n_SE[] <- sum(n_SE_eff[i, ]) # rowSums + # n_SE_eff[, ] <- if(S[i] <= 0) 0 else (if(stoch == 1) rbinom(S_eff_prod[j, i], p_SE[i]) else S_eff_prod[j, i] * p_SE[i]) + n_SE_eff[, ] <- if(S[i] <= 0) 0 else (if(stoch == 1) rbinom(S_eff_prod[i, j], p_SE[j]) else S_eff_prod[i, j] * p_SE[j]) + n_SE[] <- sum(n_SE_eff[i, ]) # sum over destinations for each origin i n_EI[] <- if(E[i] == 0) 0 else (if(stoch == 1) rbinom(E[i], p_EIpresymp[i]) else E[i] * p_EIpresymp[i]) - n_EIpresymp[] <- n_EI[i] * (1 - pea[i]) + n_EIpresymp[] <- if(stoch == 1) rbinom(n_EI[i], 1 - pea[i]) else n_EI[i] * (1 - pea[i]) n_EIasymp[] <- n_EI[i] - n_EIpresymp[i] n_preIsymp[] <- if(I_presymp[i] == 0) 0 else (if(stoch == 1) rbinom(I_presymp[i], p_preIsymp[i]) else I_presymp[i] * p_preIsymp[i]) n_IasympR[] <- if(I_asymp[i] == 0) 0 else (if(stoch == 1) rbinom(I_asymp[i], p_IasympR[i]) else I_asymp[i] * p_IasympR[i]) n_IsympRH[] <- if(I_symp[i] == 0) 0 else (if(stoch == 1) rbinom(I_symp[i], p_IsympRH[i]) else I_symp[i] * p_IsympRH[i]) - n_IsympH[] <- n_IsympRH[i] * (1 - psr[i]) + n_IsympH[] <- if(stoch == 1) rbinom(n_IsympRH[i], 1 - psr[i]) else n_IsympRH[i] * (1 - psr[i]) n_IsympR[] <- n_IsympRH[i] - n_IsympH[i] n_HRD[] <- if(H[i] == 0) 0 else (if(stoch == 1) rbinom(H[i], p_HRD[i]) else H[i] * p_HRD[i]) - n_HR[] <- n_HRD[i] * phr[i] + n_HR[] <- if(stoch == 1) rbinom(n_HRD[i], phr[i]) else n_HRD[i] * phr[i] n_HD[] <- n_HRD[i] - n_HR[i] n_RS[] <- if(R[i] == 0) 0 else (if(stoch == 1) rbinom(R[i], p_RS[i]) else R[i] * p_RS[i]) n_VE_eff[, ] <- if(stoch == 1) rbinom(V_eff_prod[j, i], p_VE[i]) else V_eff_prod[j, i] * p_VE[i] @@ -387,6 +388,7 @@ meta_sim <- function(N_pop, ts, output(p_SE) <- TRUE output(p_VE) <- TRUE output(p_HRD) <- TRUE + output(p_RS) <- TRUE output(I_eff) <- TRUE output(n_SE) <- TRUE output(n_SV) <- TRUE @@ -402,6 +404,10 @@ meta_sim <- function(N_pop, ts, output(n_HR) <- TRUE output(n_HD) <- TRUE output(n_IasympR) <- TRUE + output(n_RS) <- TRUE + + output(n_SE_eff) <- TRUE + output(S_eff_prod) <- TRUE ## ================================================= diff --git a/man/MetaRVMResults.Rd b/man/MetaRVMResults.Rd index 61b5086..4ec0510 100644 --- a/man/MetaRVMResults.Rd +++ b/man/MetaRVMResults.Rd @@ -58,6 +58,7 @@ Arindam Fadikar \itemize{ \item \href{#method-MetaRVMResults-new}{\code{MetaRVMResults$new()}} \item \href{#method-MetaRVMResults-print}{\code{MetaRVMResults$print()}} +\item \href{#method-MetaRVMResults-plot}{\code{MetaRVMResults$plot()}} \item \href{#method-MetaRVMResults-subset_data}{\code{MetaRVMResults$subset_data()}} \item \href{#method-MetaRVMResults-summarize}{\code{MetaRVMResults$summarize()}} \item \href{#method-MetaRVMResults-clone}{\code{MetaRVMResults$clone()}} @@ -108,6 +109,55 @@ Self (invisible) } } \if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-MetaRVMResults-plot}{}}} +\subsection{Method \code{plot()}}{ +Plot simulation trajectories directly from results +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{MetaRVMResults$plot( + group_by = character(0), + disease_states = NULL, + date_range = NULL, + instances = NULL, + stats = NULL, + quantiles = c(0.25, 0.75), + exclude_p_columns = TRUE, + ci_level = 0.95, + theme = theme_minimal(), + title = NULL +)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{group_by}}{Vector of demographic category names to group/facet by} + +\item{\code{disease_states}}{Optional disease states to include} + +\item{\code{date_range}}{Optional date range for filtering} + +\item{\code{instances}}{Optional instance IDs to include} + +\item{\code{stats}}{Statistics for summary plotting. If NULL, plots raw trajectories.} + +\item{\code{quantiles}}{Quantiles for uncertainty bands when summary plotting} + +\item{\code{exclude_p_columns}}{Logical, whether to exclude p_ columns (default: TRUE)} + +\item{\code{ci_level}}{Confidence level label used in summary plot title} + +\item{\code{theme}}{ggplot2 theme function (default: theme_minimal())} + +\item{\code{title}}{Optional custom title} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +ggplot object +} +} +\if{html}{\out{
}} \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-MetaRVMResults-subset_data}{}}} \subsection{Method \code{subset_data()}}{ diff --git a/vignettes/checkpointing.Rmd b/vignettes/checkpointing.Rmd index beea3b9..cbaa7e3 100644 --- a/vignettes/checkpointing.Rmd +++ b/vignettes/checkpointing.Rmd @@ -32,7 +32,7 @@ The checkpointing system uses three main configuration options in the `simulatio ## Example -For this example, the checkpoint directory is set to a temporary directory, and a temporary YAML configuration file is created from a template. +For the purpose of this example, we will set the checkpointed directory to be a temporary directory, and will create a temporary yaml configuration file based on a template. ### Set Up @@ -48,7 +48,7 @@ example_config <- system.file("extdata", "example_config_checkpoint.yaml", # Create a temporary directory for checkpoints checkpoint_dir <- tempdir() cat("Checkpoint directory:", checkpoint_dir, "\n") -#> Checkpoint directory: /tmp/RtmpwQdGR9 +#> Checkpoint directory: /tmp/RtmpdtpnKS ``` ### Create a Configuration with Checkpointing @@ -63,9 +63,6 @@ yml_tmp <- data.table::copy(yml) # Assign checkpoint directory yml_tmp$simulation_config$checkpoint_dir <- checkpoint_dir -yml_tmp$simulation_config$nrep <- 1 -yml_tmp$simulation_config$simulation_mode <- "deterministic" -yml_tmp$simulation_config$random_seed <- 42 # Create a temporary config temp_config <- tempfile(tmpdir = dirname(example_config), fileext = ".yaml") @@ -125,7 +122,7 @@ print(head(results$results, 10)) ### Resume from a Checkpoint -A new configuration is created to restore from the checkpoint. The checkpoint file name must be provided, along with a start date that should be the day after the end date in the previous run. +We'll create a new configuration that restores from the checkpoint. We need to provide the checkpoint file name, and a start date which should be the next day of the end date in the previous run. ``` r @@ -142,11 +139,8 @@ if (length(checkpoint_files_full) > 0) { new_start_date <- "12/30/2024" yml_tmp <- data.table::copy(yml) yml_tmp$simulation_config$start_date <- new_start_date - yml_tmp$simulation_config$nrep <- 1 - yml_tmp$simulation_config$simulation_mode <- "deterministic" - yml_tmp$simulation_config$random_seed <- 42 yml_tmp$simulation_config$restore_from <- checkpoint_to_restore - yml_tmp$population_data$initialization <- "population_init_n24_only_mapping.csv" # ensure that the original initial conditions are not reapplied during restore + yml_tmp$population_data$initialization <- "population_init_n24_only_mapping.csv" # ensure that we don't want to reinitialize the population with the original initial conditions temp_config_resume <- tempfile(tmpdir = dirname(example_config), fileext = ".yaml") yaml::write_yaml(yml_tmp, temp_config_resume) @@ -163,7 +157,7 @@ if (length(checkpoint_files_full) > 0) { } else { cat("No checkpoint files found to resume from.\n") } -#> /tmp/RtmpwQdGR9/chk_2024-12-29_1.Rda +#> /tmp/RtmpdtpnKS/chk_2024-12-29_1.Rda #> #> Number of instances: 1 #> Date range: 2024-12-30 to 2025-03-29 @@ -184,7 +178,7 @@ ggplot(rbind(run1_hosp_sum, run2_hosp_sum), aes(date, total)) + ```
-plot of chunk plot +plot of chunk plot

plot of chunk plot

@@ -203,3 +197,5 @@ Where: For more information on configuring MetaRVM simulations, see the [YAML Configuration](yaml-configuration.html) vignette. + + From 75bb8a4acac8276c5461b9df443387a41cd75c66 Mon Sep 17 00:00:00 2001 From: arindam fadikar Date: Thu, 26 Mar 2026 12:27:16 -0400 Subject: [PATCH 3/5] Increment version number to 2.1.0 --- DESCRIPTION | 2 +- NEWS.md | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6de3a16..c8210df 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: MetaRVM Title: Meta-Population Compartmental Model for Respiratory Virus Diseases -Version: 2.0.0 +Version: 2.1.0 Authors@R: c( person("Arindam", "Fadikar", , "afadikar@anl.gov", role = c("aut", "cre", "cph"), comment = c(ORCID = "0000-0001-7396-0350")), person("Charles", "Macal", role = "ctb"), diff --git a/NEWS.md b/NEWS.md index 9c23cdc..f2eb67c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,5 @@ +# MetaRVM 2.1.0 + # MetaRVM 2.0.0 * Subpopulation categories are now user-defined from From fc4bbc0e1097dc775fc299d3ef53cc42dbc163a6 Mon Sep 17 00:00:00 2001 From: arindam fadikar Date: Thu, 26 Mar 2026 16:52:01 -0400 Subject: [PATCH 4/5] preparing for CRAN submission --- NEWS.md | 28 ++ R/classes.R | 10 +- R/metaODIN.R | 45 +++- README.Rmd | 8 + README.md | 8 + docs/404.html | 2 +- docs/LICENSE-text.html | 2 +- docs/LICENSE.html | 2 +- docs/articles/checkpointing.html | 19 +- docs/articles/getting-started.html | 48 ++-- docs/articles/index.html | 2 +- docs/articles/plot-1.png | Bin 17511 -> 17346 bytes docs/articles/running-a-simulation.html | 248 ++++++++++-------- .../figure-html/unnamed-chunk-12-1.png | Bin 89806 -> 89765 bytes .../figure-html/unnamed-chunk-17-1.png | Bin 0 -> 98651 bytes docs/articles/yaml-configuration.html | 140 +++++++--- docs/authors.html | 2 +- docs/index.html | 8 +- docs/news/index.html | 21 +- docs/pkgdown.yml | 2 +- docs/reference/MetaRVM-package.html | 2 +- docs/reference/MetaRVMCheck.html | 2 +- docs/reference/MetaRVMConfig.html | 2 +- docs/reference/MetaRVMResults-1.png | Bin 160545 -> 159681 bytes docs/reference/MetaRVMResults.html | 94 ++++++- docs/reference/MetaRVMSummary-1.png | Bin 200147 -> 199597 bytes docs/reference/MetaRVMSummary-2.png | Bin 204691 -> 204376 bytes docs/reference/MetaRVMSummary.html | 2 +- docs/reference/draw_sample.html | 2 +- docs/reference/format_metarvm_output.html | 2 +- docs/reference/grapes-or-or-grapes.html | 2 +- docs/reference/index.html | 2 +- docs/reference/metaRVM-1.png | Bin 181239 -> 179640 bytes docs/reference/metaRVM.html | 19 +- docs/reference/meta_sim.html | 2 +- docs/reference/parse_config.html | 32 ++- docs/reference/process_vac_data.html | 2 +- docs/search.json | 2 +- inst/extdata/example_config_stochastic.yaml | 2 +- man/figures/README-unnamed-chunk-3-1.png | Bin 23660 -> 23211 bytes tests/testthat/test-metaRVM.R | 49 ++++ vignettes/checkpointing.Rmd | 13 +- vignettes/plot-1.png | Bin 17511 -> 17346 bytes vignettes/running-a-simulation.Rmd | 17 +- 44 files changed, 580 insertions(+), 263 deletions(-) create mode 100644 docs/articles/running-a-simulation_files/figure-html/unnamed-chunk-17-1.png diff --git a/NEWS.md b/NEWS.md index f2eb67c..c0d2aea 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,33 @@ # MetaRVM 2.1.0 +* Added `simulation_mode` (`deterministic`/`stochastic`) and `nrep` support in + configuration parsing and run metadata. + +* Improved stochastic reproducibility: provided `random_seed` is honored, and a + seed is generated/stored when missing. + +* Updated stochastic transition splitting to use binomial draws for integer + compartment flows (`n_EIpresymp`, `n_IsympH`, `n_HR`). + +* Improved metapopulation infection-flow indexing consistency for destination- + specific force of infection in `n_SE_eff` and `n_VE_eff`. + +* Added bounded vaccination and transition safeguards to reduce overdraw risks + (for example, capping applied vaccination by available susceptible count and + guarding zero effective population in force-of-infection calculations). + +* Added integer row allocation for stochastic mixing flows (`S` and `V`) using + floor-plus-residual adjustment so row totals are conserved. + +* Added regression tests for stochastic seed behavior, `nsim * nrep` instance + counts, and row-allocation conservation checks for `S`/`V` mixing. + +* Enhanced results helper APIs: + - `MetaRVMResults$plot()` now provides direct trajectory plotting before + summarization. + - `MetaRVMSummary$plot()` now supports instance-trajectory plotting. + - `summarize()` uses a more efficient single-pass grouped aggregation. + # MetaRVM 2.0.0 * Subpopulation categories are now user-defined from diff --git a/R/classes.R b/R/classes.R index b7b6deb..2a9504a 100644 --- a/R/classes.R +++ b/R/classes.R @@ -818,8 +818,8 @@ MetaRVMSummary <- R6::R6Class( color_var <- if ("disease_state" %in% columns) "disease_state" else "instance" p <- ggplot(self$data, aes(x = date, y = value, - color = rlang::.data[[color_var]], - group = interaction(rlang::.data[["instance"]], rlang::.data[[color_var]], drop = TRUE))) + + color = .data[[color_var]], + group = interaction(.data[["instance"]], .data[[color_var]], drop = TRUE))) + geom_line(alpha = 0.45, linewidth = 0.5) if (length(group_vars) == 1) { @@ -881,10 +881,10 @@ MetaRVMSummary <- R6::R6Class( ci_upper_col <- quantile_cols[length(quantile_cols)] # Create the plot - p <- ggplot(self$data, aes(x = date, y = median_value, color = rlang::.data[[color_var]])) + + p <- ggplot(self$data, aes(x = date, y = median_value, color = .data[[color_var]])) + geom_line(linewidth = 1) + - geom_ribbon(aes(ymin = rlang::.data[[ci_lower_col]], ymax = rlang::.data[[ci_upper_col]], - fill = rlang::.data[[color_var]]), alpha = 0.2, color = NA) + + geom_ribbon(aes(ymin = .data[[ci_lower_col]], ymax = .data[[ci_upper_col]], + fill = .data[[color_var]]), alpha = 0.2, color = NA) + labs( title = title %||% paste0("Median Outcomes with ", ci_level*100, "% Empirical Quantiles"), x = "Date", diff --git a/R/metaODIN.R b/R/metaODIN.R index 38e439b..dc835d6 100644 --- a/R/metaODIN.R +++ b/R/metaODIN.R @@ -308,7 +308,7 @@ meta_sim <- function(N_pop, ts, dim(vac) <- user() # n_SV_eff[] <- n_SV[i] * vac_eff[i] - n_SV_eff[] <- n_SV[i] * 1 # all vaccinated people are moved to V + n_SV_eff[] <- if(n_SV[i] < 0) 0 else (if(n_SV[i] > S[i]) S[i] else n_SV[i]) # bounded by available S ## ================================================= # time/day specific mobility matrix @@ -334,22 +334,33 @@ meta_sim <- function(N_pop, ts, # first remove vaccinated people from S S_eff_prod[, ] <- m[i, j] * (S[i] - n_SV_eff[i]) - V_eff_prod[, ] <- m[i, j] * V[i] + S_src_int[] <- if((S[i] - n_SV_eff[i]) <= 0) 0 else floor(S[i] - n_SV_eff[i]) + S_target[, ] <- m[i, j] * S_src_int[i] + S_base[, ] <- floor(S_target[i, j]) + S_resid[] <- S_src_int[i] - sum(S_base[i, ]) + S_alloc[, ] <- if(j == N_pop) S_base[i, j] + S_resid[i] else S_base[i, j] + + V_src_int[] <- if(V[i] <= 0) 0 else floor(V[i]) + V_target[, ] <- m[i, j] * V_src_int[i] + V_base[, ] <- floor(V_target[i, j]) + V_resid[] <- V_src_int[i] - sum(V_base[i, ]) + V_alloc[, ] <- if(j == N_pop) V_base[i, j] + V_resid[i] else V_base[i, j] I_eff_prod[, ] <- m[i, j] * I_all[i] I_eff[] <- sum(I_eff_prod[, i]) # colSums ## ================================================= ## Force of infection - lambda_i[] <- beta_i[i] * I_eff[i] / P_eff[i] - lambda_v[] <- beta_i[i] * (1 - vac_eff[i]) * I_eff[i] / P_eff[i] + lambda_i[] <- if(P_eff[i] <= 0) 0 else beta_i[i] * I_eff[i] / P_eff[i] + lambda_v[] <- if(P_eff[i] <= 0) 0 else beta_i[i] * (1 - vac_eff[i]) * I_eff[i] / P_eff[i] ## ================================================= ## Draws from binomial distributions for numbers changing between ## compartments: - # n_SE_eff[, ] <- if(S[i] <= 0) 0 else (if(stoch == 1) rbinom(S_eff_prod[j, i], p_SE[i]) else S_eff_prod[j, i] * p_SE[i]) - n_SE_eff[, ] <- if(S[i] <= 0) 0 else (if(stoch == 1) rbinom(S_eff_prod[i, j], p_SE[j]) else S_eff_prod[i, j] * p_SE[j]) + n_SE_eff[, ] <- if((S[i] - n_SV_eff[i]) <= 0) 0 else ( + if(stoch == 1) rbinom(S_alloc[i, j], p_SE[j]) else S_eff_prod[i, j] * p_SE[j] + ) n_SE[] <- sum(n_SE_eff[i, ]) # sum over destinations for each origin i n_EI[] <- if(E[i] == 0) 0 else (if(stoch == 1) rbinom(E[i], p_EIpresymp[i]) else E[i] * p_EIpresymp[i]) n_EIpresymp[] <- if(stoch == 1) rbinom(n_EI[i], 1 - pea[i]) else n_EI[i] * (1 - pea[i]) @@ -363,9 +374,13 @@ meta_sim <- function(N_pop, ts, n_HR[] <- if(stoch == 1) rbinom(n_HRD[i], phr[i]) else n_HRD[i] * phr[i] n_HD[] <- n_HRD[i] - n_HR[i] n_RS[] <- if(R[i] == 0) 0 else (if(stoch == 1) rbinom(R[i], p_RS[i]) else R[i] * p_RS[i]) - n_VE_eff[, ] <- if(stoch == 1) rbinom(V_eff_prod[j, i], p_VE[i]) else V_eff_prod[j, i] * p_VE[i] + n_VE_eff[, ] <- if(V[i] <= 0) 0 else ( + if(stoch == 1) rbinom(V_alloc[i, j], p_VE[j]) else V_eff_prod[i, j] * p_VE[j] + ) n_VE[] <- sum(n_VE_eff[i, ]) # rowSums - n_VS[] <- if(stoch == 1) rbinom(V[i] - n_VE[i], p_VS[i]) else (V[i] - n_VE[i]) * p_VS[i] + n_VS[] <- if((V[i] - n_VE[i]) <= 0) 0 else ( + if(stoch == 1) rbinom(V[i] - n_VE[i], p_VS[i]) else (V[i] - n_VE[i]) * p_VS[i] + ) ## ================================================= ## Initial states: @@ -408,6 +423,10 @@ meta_sim <- function(N_pop, ts, output(n_SE_eff) <- TRUE output(S_eff_prod) <- TRUE + output(S_src_int) <- TRUE + output(S_alloc) <- TRUE + output(V_src_int) <- TRUE + output(V_alloc) <- TRUE ## ================================================= @@ -506,6 +525,16 @@ meta_sim <- function(N_pop, ts, dim(n_SV_eff) <- N_pop dim(n_SE_eff) <- c(N_pop, N_pop) dim(n_VE_eff) <- c(N_pop, N_pop) + dim(S_src_int) <- N_pop + dim(S_target) <- c(N_pop, N_pop) + dim(S_base) <- c(N_pop, N_pop) + dim(S_resid) <- N_pop + dim(S_alloc) <- c(N_pop, N_pop) + dim(V_src_int) <- N_pop + dim(V_target) <- c(N_pop, N_pop) + dim(V_base) <- c(N_pop, N_pop) + dim(V_resid) <- N_pop + dim(V_alloc) <- c(N_pop, N_pop) dim(lambda_i) <- N_pop dim(lambda_v) <- N_pop diff --git a/README.Rmd b/README.Rmd index 9a6df47..ea1fe7f 100644 --- a/README.Rmd +++ b/README.Rmd @@ -50,6 +50,14 @@ The current CRAN release is **2.0.0**. Install it with: install.packages("MetaRVM") ``` +To install the development branch for **MetaRVM 2.1.0** (branch: +`feature/simulation-mode-config`) directly from GitHub: + +``` r +# install.packages("remotes") +remotes::install_github("RESUME-Epi/MetaRVM@feature/simulation-mode-config") +``` + ### Running a simulation ```{r} diff --git a/README.md b/README.md index 666b59f..11d4d71 100644 --- a/README.md +++ b/README.md @@ -56,6 +56,14 @@ The current CRAN release is **2.0.0**. Install it with: install.packages("MetaRVM") ``` +To install the development branch for **MetaRVM 2.1.0** (branch: +`feature/simulation-mode-config`) directly from GitHub: + +``` r +# install.packages("remotes") +remotes::install_github("RESUME-Epi/MetaRVM@feature/simulation-mode-config") +``` + ### Running a simulation ``` r diff --git a/docs/404.html b/docs/404.html index 4d7532f..99bbfd1 100644 --- a/docs/404.html +++ b/docs/404.html @@ -28,7 +28,7 @@ MetaRVM - 2.0.0 + 2.1.0