-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* Add tests to test-sir.R * Update tests in test-sir.R * Add queuing timing test to test-sir.R * Clean up test-sir.R * Add add tests to test-lfmcmc.R * Add tests for factory methods and missing parameters to test-lfmcmc.R * Create test-sirconn.R and populate with tests * Create test-sird.R and populate with tests * Lower acceptance tolerance for expected check in test-sirconn.R * Set tolerance to 0.1 in test-sirconn.R * Create test-sis.R and populate with tests * Add checks to test-lfmcmc, -sir, -sirconn, -sird, -sis for unexpected inheritance * Create test-seir.R and populate with tests * Add check for warning from queuing_on() in test-sirconn.R * Crate test-tool.R and populate with tests * Clean up test-tool.R
- Loading branch information
1 parent
8b817d1
commit ad6aca2
Showing
8 changed files
with
356 additions
and
111 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,96 +1,119 @@ | ||
# Create Model to use in LFMCMC simulation | ||
model_seed <- 122 | ||
# Test just this file: tinytest::run_test_file("inst/tinytest/test-lfmcmc.R") | ||
|
||
model_sir <- ModelSIR( | ||
name = "COVID-19", | ||
prevalence = .1, | ||
transmission_rate = .9, | ||
recovery_rate = .3 | ||
) | ||
|
||
agents_smallworld( | ||
model_sir, | ||
n = 1000, | ||
k = 5, | ||
d = FALSE, | ||
p = 0.01 | ||
) | ||
# Set model parameters | ||
model_seed <- 122 | ||
|
||
# Create and run SIR Model for LFMCMC simulation ------------------------------- | ||
model_sir <- ModelSIR(name = "COVID-19", prevalence = .1, | ||
transmission_rate = .9, recovery_rate = .3) | ||
agents_smallworld(model_sir, n = 1000, k = 5, d = FALSE, p = 0.01) | ||
verbose_off(model_sir) | ||
run(model_sir, ndays = 50, seed = model_seed) | ||
|
||
run( | ||
model_sir, | ||
ndays = 50, | ||
seed = model_seed | ||
) | ||
# Check bad init of LFMCMC model ----------------------------------------------- | ||
expect_error(lfmcmc_bad <- LFMCMC(), 'argument "model" is missing') | ||
expect_error(lfmcmc_bad <- LFMCMC(c("not_a_model")), "model should be of class 'epiworld_model'") | ||
|
||
# Create LFMCMC model ---------------------------------------------------------- | ||
expect_silent(lfmcmc_model <- LFMCMC(model_sir)) | ||
|
||
# Check initialization | ||
expect_inherits(lfmcmc_model, "epiworld_lfmcmc") | ||
expect_length(class(lfmcmc_model), 1) | ||
|
||
# Setup LFMCMC | ||
## Extract the observed data from the model | ||
# Extract observed data from the model | ||
obs_data <- unname(as.integer(get_today_total(model_sir))) | ||
|
||
## Define the LFMCMC simulation function | ||
simfun <- function(params) { | ||
expect_silent(set_observed_data(lfmcmc_model, obs_data)) | ||
|
||
# Define LFMCMC functions | ||
simfun <- function(params) { | ||
set_param(model_sir, "Recovery rate", params[1]) | ||
set_param(model_sir, "Transmission rate", params[2]) | ||
|
||
run( | ||
model_sir, | ||
ndays = 50 | ||
) | ||
|
||
run(model_sir, ndays = 50) | ||
res <- unname(as.integer(get_today_total(model_sir))) | ||
return(res) | ||
} | ||
|
||
## Define the LFMCMC summary function | ||
sumfun <- function(dat) { | ||
return(dat) | ||
} | ||
sumfun <- function(dat) { return(dat) } | ||
|
||
## Define the LFMCMC proposal function | ||
## - Based on proposal_fun_normal from lfmcmc-meat.hpp | ||
propfun <- function(params_prev) { | ||
res <- params_prev + rnorm(length(params_prev), ) | ||
return(res) | ||
} | ||
|
||
## Define the LFMCMC kernel function | ||
## - Based on kernel_fun_uniform from lfmcmc-meat.hpp | ||
kernelfun <- function(stats_now, stats_obs, epsilon) { | ||
|
||
ans <- sum(mapply(function(v1, v2) (v1 - v2)^2, | ||
stats_obs, | ||
stats_now)) | ||
|
||
ans <- sum(mapply(function(v1, v2) (v1 - v2)^2, stats_obs, stats_now)) | ||
return(ifelse(sqrt(ans) < epsilon, 1.0, 0.0)) | ||
} | ||
|
||
## Create the LFMCMC model | ||
lfmcmc_model <- LFMCMC(model_sir) |> | ||
set_simulation_fun(simfun) |> | ||
set_summary_fun(sumfun) |> | ||
set_proposal_fun(propfun) |> | ||
set_kernel_fun(kernelfun) |> | ||
set_observed_data(obs_data) | ||
# Check adding functions to LFMCMC | ||
expect_silent(set_simulation_fun(lfmcmc_model, simfun)) | ||
expect_silent(set_summary_fun(lfmcmc_model, sumfun)) | ||
expect_silent(set_proposal_fun(lfmcmc_model, propfun)) | ||
expect_silent(set_kernel_fun(lfmcmc_model, kernelfun)) | ||
|
||
# Run LFMCMC simulation | ||
## Set initial parameters | ||
# Run LFMCMC simulation -------------------------------------------------------- | ||
# Initial parameters | ||
par0 <- as.double(c(0.1, 0.5)) | ||
n_samp <- 2000 | ||
epsil <- as.double(1.0) | ||
|
||
## Run the LFMCMC simulation | ||
run_lfmcmc( | ||
expect_silent(run_lfmcmc( | ||
lfmcmc = lfmcmc_model, | ||
params_init_ = par0, | ||
n_samples_ = n_samp, | ||
epsilon_ = epsil, | ||
seed = model_seed | ||
)) | ||
|
||
expect_silent(set_stats_names(lfmcmc_model, get_states(model_sir))) | ||
expect_silent(set_par_names(lfmcmc_model, c("Immune recovery", "Infectiousness"))) | ||
|
||
expect_stdout(print(lfmcmc_model)) | ||
|
||
# Check LFMCMC using factory functions ----------------------------------------- | ||
expect_silent(use_proposal_norm_reflective(lfmcmc_model)) | ||
expect_silent(use_kernel_fun_gaussian(lfmcmc_model)) | ||
|
||
expect_silent(run_lfmcmc( | ||
lfmcmc = lfmcmc_model, | ||
params_init_ = par0, | ||
n_samples_ = n_samp, | ||
epsilon_ = epsil, | ||
seed = model_seed | ||
) | ||
)) | ||
|
||
# Check running LFMCMC with missing parameters --------------------------------- | ||
expect_silent(run_lfmcmc( | ||
lfmcmc = lfmcmc_model, | ||
params_init_ = par0, | ||
n_samples_ = n_samp, | ||
epsilon_ = epsil | ||
)) | ||
|
||
expect_error(run_lfmcmc( | ||
lfmcmc = lfmcmc_model, | ||
params_init_ = par0, | ||
n_samples_ = n_samp | ||
)) | ||
|
||
# Print the results | ||
set_stats_names(lfmcmc_model, get_states(model_sir)) | ||
set_par_names(lfmcmc_model, c("Immune recovery", "Infectiousness")) | ||
expect_error(run_lfmcmc( | ||
lfmcmc = lfmcmc_model, | ||
params_init_ = par0, | ||
epsilon_ = epsil | ||
)) | ||
|
||
expect_error(run_lfmcmc( | ||
lfmcmc = lfmcmc_model, | ||
n_samples_ = n_samp, | ||
epsilon_ = epsil | ||
)) | ||
|
||
expect_error(run_lfmcmc( | ||
params_init_ = par0, | ||
n_samples_ = n_samp, | ||
epsilon_ = epsil | ||
)) | ||
|
||
print(lfmcmc_model) | ||
expect_error(run_lfmcmc(lfmcmc = lfmcmc_model)) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,38 @@ | ||
# Test just this file: tinytest::run_test_file("inst/tinytest/test-seir.R") | ||
|
||
# Create small world population SEIR Model ------------------------------------- | ||
expect_silent(seir_0 <- ModelSEIR( | ||
name = "A Virus", | ||
prevalence = .01, | ||
transmission_rate = .5, | ||
incubation_days = 7.0, | ||
recovery_rate = .1 | ||
)) | ||
|
||
# Check model initialization | ||
expect_inherits(seir_0, "epiworld_seir") | ||
expect_inherits(seir_0, "epiworld_model") | ||
expect_length(class(seir_0), 2) | ||
expect_silent(agents_smallworld( | ||
seir_0, | ||
n = 10000, | ||
k = 5, | ||
d = FALSE, | ||
p = .01 | ||
)) | ||
|
||
|
||
# Check model run without queuing ---------------------------------------------- | ||
expect_silent(verbose_off(seir_0)) | ||
expect_silent(queuing_off(seir_0)) | ||
expect_silent(initial_states(seir_0, c(.3, .5))) | ||
expect_warning(expect_error(plot(seir_0))) # Plot fails before model is run | ||
expect_silent(run(seir_0, ndays = 100, seed = 1231)) | ||
expect_silent(plot(seir_0)) # Plot succeeds after model is run | ||
|
||
hist_0 <- get_hist_total(seir_0) | ||
|
||
expect_equal(hist_0[1,3], 4950) | ||
expect_equal(hist_0[2,3], 70) | ||
expect_equal(hist_0[3,3], 30) | ||
expect_equal(hist_0[4,3], 4950) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,65 +1,61 @@ | ||
# Adding a Small world population without queuing ------------------------------ | ||
sir_0 <- ModelSIR( | ||
name = "COVID-19", | ||
prevalence = .01, | ||
transmission_rate = .9, | ||
recovery_rate = .3 | ||
) | ||
|
||
agents_smallworld( | ||
sir_0, | ||
n = 50000, | ||
k = 5, | ||
d = FALSE, | ||
p = .01 | ||
) | ||
|
||
# Initializing | ||
queuing_off(sir_0) | ||
|
||
# Running and printing | ||
run(sir_0, ndays = 50, seed = 1912) | ||
|
||
tmat_0 <- get_transition_probability(sir_0) | ||
# Test just this file: tinytest::run_test_file("inst/tinytest/test-sir.R") | ||
|
||
# Function to test transition probability matrix ------------------------ | ||
test_tmat_matches_expected <- function(tmat) { | ||
tmat_expected <- structure( | ||
c( | ||
0.961843, 0, 0, | ||
0.03815696, 0.69985167, 0, | ||
0, 0.3001483, 1 | ||
), | ||
dim = c(3L, 3L), | ||
dimnames = list( | ||
c("Susceptible", "Infected", "Recovered"), | ||
c("Susceptible", "Infected", "Recovered") | ||
) | ||
) | ||
|
||
expect_equal(tmat, tmat_expected, tolerance = 0.0000001) | ||
} | ||
|
||
# Creating a SIR model with queuing -------------------------------------------- | ||
sir_1 <- ModelSIR( | ||
# Create small world population SIR Model -------------------------------------- | ||
expect_silent(sir_0 <- ModelSIR( | ||
name = "COVID-19", | ||
prevalence = .01, | ||
transmission_rate = .9, | ||
recovery_rate = .3 | ||
) | ||
)) | ||
|
||
# Adding a Small world population | ||
agents_smallworld( | ||
sir_1, | ||
# Check model initialization | ||
expect_inherits(sir_0, "epiworld_sir") | ||
expect_inherits(sir_0, "epiworld_model") | ||
expect_length(class(sir_0), 2) | ||
expect_silent(agents_smallworld( | ||
sir_0, | ||
n = 50000, | ||
k = 5, | ||
d = FALSE, | ||
p = .01 | ||
) | ||
|
||
# Running and printing | ||
run(sir_1, ndays = 50, seed = 1912) | ||
)) | ||
|
||
tmat_1 <- get_transition_probability(sir_1) | ||
# Check model run with queuing ------------------------------------------------- | ||
expect_silent(verbose_off(sir_0)) | ||
expect_warning(expect_error(plot(sir_0))) # Plot fails before model is run | ||
expect_silent(run(sir_0, ndays = 50, seed = 1912)) | ||
expect_silent(plot(sir_0)) # Plot succeeds after model is run | ||
|
||
# Expected | ||
tmat_expected <- structure( | ||
c( | ||
0.962139427661896, 0, 0, | ||
0.0378605611622334, 0.70049238204956, | ||
0, 0, 0.299507558345795, 1 | ||
), | ||
dim = c(3L, 3L), | ||
dimnames = list( | ||
c("Susceptible", "Infected", "Recovered"), | ||
c("Susceptible", "Infected", "Recovered") | ||
) | ||
) | ||
tmat_queuing <- get_transition_probability(sir_0) | ||
test_tmat_matches_expected(tmat_queuing) | ||
|
||
expect_equivalent(tmat_0, tmat_1) | ||
expect_equivalent(tmat_0, tmat_expected) | ||
expect_equivalent(tmat_1, tmat_expected) | ||
# Check model run without queuing ---------------------------------------------- | ||
expect_silent(queuing_off(sir_0)) | ||
run(sir_0, ndays = 50, seed = 1912) | ||
|
||
tmat_noqueuing <- get_transition_probability(sir_0) | ||
expect_identical(tmat_noqueuing, tmat_queuing) | ||
|
||
# Check queuing is faster ------------------------------------------------------ | ||
runtime_noqueuing <- system.time(run(sir_0, ndays = 50, seed = 1912)) | ||
queuing_on(sir_0) | ||
runtime_queuing <- system.time(run(sir_0, ndays = 50, seed = 1912)) | ||
expect_true(runtime_queuing["elapsed"] < runtime_noqueuing["elapsed"]) |
Oops, something went wrong.