Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update unit testing #38

Open
wants to merge 16 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
143 changes: 83 additions & 60 deletions inst/tinytest/test-lfmcmc.R
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))
2 changes: 1 addition & 1 deletion inst/tinytest/test-multiple.R
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ saver <- make_saver(
run_multiple(
model_seircon,
ndays=days,
nsim=nsims,
nsims=nsims,
seed=1972,
saver=saver,
nthreads = 2L
Expand Down
38 changes: 38 additions & 0 deletions inst/tinytest/test-seir.R
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)
96 changes: 46 additions & 50 deletions inst/tinytest/test-sir.R
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"])
Loading
Loading