diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index b4ddd76..24aef03 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ -Version: 2.0.0 -Date: 2026-03-19 19:14:10 UTC -SHA: 4795478e0caef4ec96b0827d2c20fa4275bcf2f3 +Version: 2.1.0 +Date: 2026-03-26 20:54:41 UTC +SHA: fc4bbc0e1097dc775fc299d3ef53cc42dbc163a6 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..c0d2aea 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +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 9667639..2a9504a 100644 --- a/R/classes.R +++ b/R/classes.R @@ -342,11 +342,57 @@ 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) }, + #' @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) @@ -435,6 +481,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), @@ -510,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")) } @@ -766,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) { @@ -786,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 = .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) { + 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" @@ -808,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 = .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 = .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", @@ -827,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..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,37 +334,53 @@ 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[] <- sum(n_SE_eff[i, ]) # rowSums + 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[] <- 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] + 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: @@ -387,6 +403,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 +419,14 @@ 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 + output(S_src_int) <- TRUE + output(S_alloc) <- TRUE + output(V_src_int) <- TRUE + output(V_alloc) <- TRUE ## ================================================= @@ -500,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/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..4469882 100644 --- a/README.Rmd +++ b/README.Rmd @@ -18,6 +18,10 @@ knitr::opts_chunk$set( ``` +[![R-CMD-check](https://github.com/RESUME-Epi/MetaRVM/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/RESUME-Epi/MetaRVM/actions/workflows/R-CMD-check.yaml) +[![Codecov test coverage](https://codecov.io/gh/RESUME-Epi/MetaRVM/branch/main/graph/badge.svg)](https://app.codecov.io/gh/RESUME-Epi/MetaRVM) +[![Lifecycle: maturing](https://img.shields.io/badge/lifecycle-maturing-blue.svg)](https://lifecycle.r-lib.org/articles/stages.html) +[![pkgdown site](https://img.shields.io/badge/docs-pkgdown-blue.svg)](https://RESUME-Epi.github.io/MetaRVM/) [![CRAN status](https://www.r-pkg.org/badges/version/MetaRVM)](https://CRAN.R-project.org/package=MetaRVM) [![CRAN downloads](https://cranlogs.r-pkg.org/badges/MetaRVM)](https://cran.r-project.org/package=MetaRVM) @@ -88,6 +92,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 +130,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 +223,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..eb4728a 100644 --- a/README.md +++ b/README.md @@ -5,6 +5,13 @@ +[![R-CMD-check](https://github.com/RESUME-Epi/MetaRVM/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/RESUME-Epi/MetaRVM/actions/workflows/R-CMD-check.yaml) +[![Codecov test +coverage](https://codecov.io/gh/RESUME-Epi/MetaRVM/branch/main/graph/badge.svg)](https://app.codecov.io/gh/RESUME-Epi/MetaRVM) +[![Lifecycle: +maturing](https://img.shields.io/badge/lifecycle-maturing-blue.svg)](https://lifecycle.r-lib.org/articles/stages.html) +[![pkgdown +site](https://img.shields.io/badge/docs-pkgdown-blue.svg)](https://RESUME-Epi.github.io/MetaRVM/) [![CRAN status](https://www.r-pkg.org/badges/version/MetaRVM)](https://CRAN.R-project.org/package=MetaRVM) [![CRAN @@ -95,6 +102,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 +160,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/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