Closed pratikunterwegs closed 10 months ago
Posting an initial benchmarking of the deafult model in {epidemics} started by @bahadzie in https://github.com/bahadzie/epidemics/blob/benchmarks/vignettes/COMPARE.Rmd.
Using both the {bench} and {microbenchmark} packages the Cpp version is quicker, with the difference larger when using {bench} (not sure why this is).
This benchmarking is on a Mac with an M2 chip.
r
pak::pak("bahadzie/epidemics@benchmarks")
#> ℹ Loading metadata database
#> ℹ Loading metadata database
#> ✔ Loading metadata database ... done
#> ✔ Loading metadata database ... done
#>
#>
#>
#> ℹ No downloads are needed
#> ℹ No downloads are needed
#> ✔ 1 pkg + 10 deps: kept 8 [5s]
#> ✔ 1 pkg + 10 deps: kept 8 [5s]
library(epidemics)
# Create stratified populations based on supplied age limits
pop <- function(age_limits = c(0, 20, 40)) {
# load contact and population data from socialmixr::polymod
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
polymod,
countries = "United Kingdom",
age.limits = age_limits,
symmetric = TRUE
)
# prepare contact matrix
contact_matrix <- t(contact_data$matrix)
contact_matrix
# prepare the demography vector
demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)
# initial conditions: one in every 1 million is infected
initial_i <- 1e-6
initial_condition <- c(
S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
)
ncol <- length(initial_condition)
nrow <- length(rownames(contact_matrix))
initial_conditions <- matrix(rep(initial_condition, nrow), ncol = ncol)
colnames(initial_conditions) <- names(initial_condition)
rownames(initial_conditions) <- rownames(contact_matrix)
population(
name = "UK",
contact_matrix = contact_matrix,
demography_vector = demography_vector,
initial_conditions = initial_conditions
)
}
uk_pop <- pop(c(0, 20, 65)) # ages beyond 80?
#> Using POLYMOD social contact data. To cite this in a publication, use the 'cite' function
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option
# simulate a pandemic, with an R0,
# an infectious period, and an pre-infectious period
pandemic_influenza <- infection(
r0 = 1.5,
preinfectious_period = 3,
infectious_period = 7
)
# an intervention to close schools
close_schools <- intervention(
time_begin = 200,
time_end = 260,
contact_reduction = round(runif(length(uk_pop$demography_vector)), 2)
)
# a vaccination to close schools
vax <- vaccination(
time_begin = 100,
time_end = 400,
nu = runif(length(uk_pop$demography_vector)) * 1e-4
)
# benchmarking using {bench}
bnch <- bench::mark(
min_iterations = 100, max_iterations = 100,
"epidemic_cpp" = epidemic(
model = "default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
"epidemic_r" = r_epidemic(
model = "r_default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
check = FALSE # setting false for now
)
summary(bnch)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 epidemic_cpp 6.73ms 7.33ms 89.0 1.56MB 2.75
#> 2 epidemic_r 26.64ms 28.86ms 34.4 2.69MB 6.07
# benchmarking using {microbenchmark}
mbench <- microbenchmark::microbenchmark(
"epidemic_cpp" = epidemic(
model = "default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
"epidemic_r" = r_epidemic(
model = "r_default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
times = 100
)
#> Warning in microbenchmark::microbenchmark(epidemic_cpp = epidemic(model =
#> "default", : less accurate nanosecond times to avoid potential integer overflows
summary(mbench)
#> expr min lq mean median uq max neval
#> 1 epidemic_cpp 14.78911 18.92238 19.52856 19.78615 20.32046 24.46855 100
#> 2 epidemic_r 24.46983 28.73629 31.10285 30.25005 31.57761 122.99360 100
Created on 2023-06-05 by the reprex package (v2.0.1)
Benchmarking the default model with 1000 iterations on an M2 chip; the Rcpp implementation is roughly 1.5x faster.
pak::pak("bahadzie/epidemics@benchmarks")
#> ℹ Loading metadata database
#> ℹ Loading metadata database
#> ✔ Loading metadata database ... done
#> ✔ Loading metadata database ... done
#>
#>
#>
#> ℹ No downloads are needed
#> ℹ No downloads are needed
#> ✔ 1 pkg + 10 deps: kept 8 [5.6s]
#> ✔ 1 pkg + 10 deps: kept 8 [5.6s]
library(epidemics)
# Create stratified populations based on supplied age limits
pop <- function(age_limits = c(0, 20, 40)) {
# load contact and population data from socialmixr::polymod
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
polymod,
countries = "United Kingdom",
age.limits = age_limits,
symmetric = TRUE
)
# prepare contact matrix
contact_matrix <- t(contact_data$matrix)
contact_matrix
# prepare the demography vector
demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)
# initial conditions: one in every 1 million is infected
initial_i <- 1e-6
initial_condition <- c(
S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
)
ncol <- length(initial_condition)
nrow <- length(rownames(contact_matrix))
initial_conditions <- matrix(rep(initial_condition, nrow), ncol = ncol)
colnames(initial_conditions) <- names(initial_condition)
rownames(initial_conditions) <- rownames(contact_matrix)
population(
name = "UK",
contact_matrix = contact_matrix,
demography_vector = demography_vector,
initial_conditions = initial_conditions
)
}
uk_pop <- pop(c(0, 20, 65)) # ages beyond 80?
#> Using POLYMOD social contact data. To cite this in a publication, use the 'cite' function
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option
# simulate a pandemic, with an R0,
# an infectious period, and an pre-infectious period
pandemic_influenza <- infection(
r0 = 1.5,
preinfectious_period = 3,
infectious_period = 7
)
# an intervention to close schools
close_schools <- intervention(
time_begin = 200,
time_end = 260,
contact_reduction = round(runif(length(uk_pop$demography_vector)), 2)
)
# a vaccination to close schools
vax <- vaccination(
time_begin = 100,
time_end = 400,
nu = runif(length(uk_pop$demography_vector)) * 1e-4
)
# benchmarking using {bench}
bnch <- bench::mark(
min_iterations = 1000, max_iterations = 1000,
"epidemic_cpp" = epidemic(
model = "default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
"epidemic_r" = r_epidemic(
model = "r_default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
check = FALSE # setting false for now
)
summary(bnch)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 epidemic_cpp 6.6ms 22.5ms 47.0 1.56MB 1.35
#> 2 epidemic_r 15.8ms 28.5ms 35.8 2.69MB 6.02
# benchmarking using {microbenchmark}
mbench <- microbenchmark::microbenchmark(
"epidemic_cpp" = epidemic(
model = "default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
"epidemic_r" = r_epidemic(
model = "r_default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
times = 1000
)
#> Warning in microbenchmark::microbenchmark(epidemic_cpp = epidemic(model =
#> "default", : less accurate nanosecond times to avoid potential integer overflows
summary(mbench)
#> expr min lq mean median uq max neval
#> 1 epidemic_cpp 7.24142 18.40187 19.33215 19.78393 20.29320 116.95586 1000
#> 2 epidemic_r 19.26077 28.41257 29.92554 29.64601 31.18464 57.17938 1000
Created on 2023-06-05 by the reprex package (v2.0.1)
Using the reprex from https://github.com/epiverse-trace/epidemics/issues/62#issuecomment-1576399861, with 500 iterations:
pak::pak("bahadzie/epidemics@benchmarks")
#> ℹ Loading metadata database
#> ✔ Loading metadata database ... done
#>
#>
#> ℹ No downloads are needed
#> ✔ 1 pkg + 10 deps: kept 10 [3.5s]
library(epidemics)
# Create stratified populations based on supplied age limits
pop <- function(age_limits = c(0, 20, 40)) {
# load contact and population data from socialmixr::polymod
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
polymod,
countries = "United Kingdom",
age.limits = age_limits,
symmetric = TRUE
)
# prepare contact matrix
contact_matrix <- t(contact_data$matrix)
contact_matrix
# prepare the demography vector
demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)
# initial conditions: one in every 1 million is infected
initial_i <- 1e-6
initial_condition <- c(
S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
)
ncol <- length(initial_condition)
nrow <- length(rownames(contact_matrix))
initial_conditions <- matrix(rep(initial_condition, nrow), ncol = ncol)
colnames(initial_conditions) <- names(initial_condition)
rownames(initial_conditions) <- rownames(contact_matrix)
population(
name = "UK",
contact_matrix = contact_matrix,
demography_vector = demography_vector,
initial_conditions = initial_conditions
)
}
uk_pop <- pop(c(0, 20, 65)) # ages beyond 80?
#> Using POLYMOD social contact data. To cite this in a publication, use the 'cite' function
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option
# simulate a pandemic, with an R0,
# an infectious period, and an pre-infectious period
pandemic_influenza <- infection(
r0 = 1.5,
preinfectious_period = 3,
infectious_period = 7
)
# an intervention to close schools
close_schools <- intervention(
time_begin = 200,
time_end = 260,
contact_reduction = round(runif(length(uk_pop$demography_vector)), 2)
)
# a vaccination to close schools
vax <- vaccination(
time_begin = 100,
time_end = 400,
nu = runif(length(uk_pop$demography_vector)) * 1e-4
)
# benchmarking using {bench}
bnch <- bench::mark(
min_iterations = 500,
"epidemic_cpp" = epidemic(
model = "default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
"epidemic_r" = r_epidemic(
model = "r_default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
check = FALSE # setting false for now
)
summary(bnch)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 epidemic_cpp 4.21ms 5.25ms 182. 1.56MB 10.8
#> 2 epidemic_r 22.1ms 25.34ms 39.0 2.71MB 8.68
# benchmarking using {microbenchmark}
mbench <- microbenchmark::microbenchmark(
"epidemic_cpp" = epidemic(
model = "default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
"epidemic_r" = r_epidemic(
model = "r_default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
times = 500
)
summary(mbench)
#> expr min lq mean median uq max neval
#> 1 epidemic_cpp 4.254519 4.873327 5.502002 5.227099 5.762491 16.5548 500
#> 2 epidemic_r 21.874709 24.381577 27.292796 25.971410 28.353433 173.4696 500
Created on 2023-06-05 with reprex v2.0.2
Checking with options(matprod = "internal")
in case %*%
on MacOS is the issue - does not seem to be, as the results are fairly similar to earlier runs on MacOS.
pak::pak("bahadzie/epidemics@benchmarks")
#>
#> ✔ Updated metadata database: 4.86 MB in 4 files.
#> ✔ Updated metadata database: 4.86 MB in 4 files.
#>
#> ℹ Updating metadata database
#> ℹ Updating metadata database
#> ✔ Updating metadata database ... done
#> ✔ Updating metadata database ... done
#>
#>
#>
#> → Will update 1 package.
#> → Will update 1 package.
#> → Will download 1 package with unknown size.
#> → Will download 1 package with unknown size.
#> + epidemics 0.0.0.9000 → 0.0.0.9000 👷🏿♂️🔧 ⬇ (GitHub: 0df5c96)
#> + epidemics 0.0.0.9000 → 0.0.0.9000 👷🏿♂️🔧 ⬇ (GitHub: 0df5c96)
#> ℹ Getting 1 pkg with unknown size
#> ℹ Getting 1 pkg with unknown size
#> ✔ Cached copy of epidemics 0.0.0.9000 (source) is the latest build
#> ✔ Cached copy of epidemics 0.0.0.9000 (source) is the latest build
#> ℹ Packaging epidemics 0.0.0.9000
#> ℹ Packaging epidemics 0.0.0.9000
#> ✔ Packaged epidemics 0.0.0.9000 (1.8s)
#> ✔ Packaged epidemics 0.0.0.9000 (1.8s)
#> ℹ Building epidemics 0.0.0.9000
#> ℹ Building epidemics 0.0.0.9000
#> ✔ Built epidemics 0.0.0.9000 (10s)
#> ✔ Built epidemics 0.0.0.9000 (10s)
#> ✔ Installed epidemics 0.0.0.9000 (github::bahadzie/epidemics@0df5c96) (61ms)
#> ✔ Installed epidemics 0.0.0.9000 (github::bahadzie/epidemics@0df5c96) (61ms)
#> ✔ 1 pkg + 11 deps: kept 7, upd 1 [27.1s]
#> ✔ 1 pkg + 11 deps: kept 7, upd 1 [27.1s]
library(epidemics)
options(matprod = "internal")
# Create stratified populations based on supplied age limits
pop <- function(age_limits = c(0, 20, 40)) {
# load contact and population data from socialmixr::polymod
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
polymod,
countries = "United Kingdom",
age.limits = age_limits,
symmetric = TRUE
)
# prepare contact matrix
contact_matrix <- t(contact_data$matrix)
contact_matrix
# prepare the demography vector
demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)
# initial conditions: one in every 1 million is infected
initial_i <- 1e-6
initial_condition <- c(
S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
)
ncol <- length(initial_condition)
nrow <- length(rownames(contact_matrix))
initial_conditions <- matrix(rep(initial_condition, nrow), ncol = ncol)
colnames(initial_conditions) <- names(initial_condition)
rownames(initial_conditions) <- rownames(contact_matrix)
population(
name = "UK",
contact_matrix = contact_matrix,
demography_vector = demography_vector,
initial_conditions = initial_conditions
)
}
uk_pop <- pop(c(0, 20, 65)) # ages beyond 80?
#> Using POLYMOD social contact data. To cite this in a publication, use the 'cite' function
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option
# simulate a pandemic, with an R0,
# an infectious period, and an pre-infectious period
pandemic_influenza <- infection(
r0 = 1.5,
preinfectious_period = 3,
infectious_period = 7
)
# an intervention to close schools
close_schools <- intervention(
time_begin = 200,
time_end = 260,
contact_reduction = round(runif(length(uk_pop$demography_vector)), 2)
)
# a vaccination to close schools
vax <- vaccination(
time_begin = 100,
time_end = 400,
nu = runif(length(uk_pop$demography_vector)) * 1e-4
)
# benchmarking using {bench}
bnch <- bench::mark(
min_iterations = 1000, max_iterations = 1000,
"epidemic_cpp" = epidemic(
model = "default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
"epidemic_r" = r_epidemic(
model = "r_default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
check = FALSE # setting false for now
)
summary(bnch)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 epidemic_cpp 6.69ms 23.6ms 44.7 1.56MB 1.33
#> 2 epidemic_r 18.93ms 30.4ms 33.6 2.69MB 5.66
# benchmarking using {microbenchmark}
mbench <- microbenchmark::microbenchmark(
"epidemic_cpp" = epidemic(
model = "default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
"epidemic_r" = r_epidemic(
model = "r_default",
population = uk_pop,
infection = pandemic_influenza,
intervention = close_schools,
vaccination = vax,
time_end = 600, increment = 1.0
),
times = 1000
)
#> Warning in microbenchmark::microbenchmark(epidemic_cpp = epidemic(model =
#> "default", : less accurate nanosecond times to avoid potential integer overflows
summary(mbench)
#> expr min lq mean median uq max neval
#> 1 epidemic_cpp 10.76369 19.99843 20.58345 20.82152 21.21734 118.59910 1000
#> 2 epidemic_r 20.97130 31.24284 32.58153 32.33664 33.76649 50.74873 1000
Created on 2023-06-12 by the reprex package (v2.0.1)
Posting a benchmark of the Vacamole model after some speed-ups implemented in #86. Now achieving an approx 2x speedup on an M2 Mac. Similar speed-up expected for the default model which has received the same edits.
library(epidemics)
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
polymod,
countries = "United Kingdom",
age.limits = c(0, 60),
symmetric = TRUE
)
#> Using POLYMOD social contact data. To cite this in a publication, use the 'cite' function
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option
contact_matrix <- t(contact_data$matrix)
demography_vector <- contact_data$demography$population
# // 0| 1| 2|3| 4|5| 6|7| 8|9|10
# // S|V1|V2|E|EV|I|IV|H|HV|D|R
# make initial conditions - order is important
initial_conditions <- c(
S = 1 - 1e-6,
V1 = 0, V2 = 0,
E = 0, EV = 0,
I = 1e-6, IV = 0,
H = 0, HV = 0, D = 0, R = 0
)
initial_conditions <- rbind(
initial_conditions,
initial_conditions
)
# create a population
uk_population <- population(
name = "UK population",
contact_matrix = matrix(1, 2, 2),
demography_vector = 67e6 * c(0.4, 0.6),
initial_conditions = initial_conditions
)
# prepare a two dose vaccination regime for three age groups
double_vaccination <- vaccination(
name = "double_vaccination",
nu = matrix(
1e-3,
nrow = 2, ncol = 2
),
time_begin = matrix(
00,
nrow = 2, ncol = 2
),
time_end = matrix(
200,
nrow = 2, ncol = 2
)
)
# make infection class for Vacamole model
# note extra arguments passed as ...
infect <- infection(
name = "covid", r0 = 5, infectious_period = 10,
preinfectious_period = 5,
eta = 1 / 1000, omega = 1 / 1000,
susc_reduction_vax = 0.5,
hosp_reduction_vax = 0.7,
mort_reduction_vax = 0.9
)
microbenchmark::microbenchmark(
"vacamole_new_cpp" = epidemic_vacamole_cpp(
population = uk_population,
infection = infect,
vaccination = double_vaccination,
time_end = 500, increment = 1
)
)
#> Warning in microbenchmark::microbenchmark(vacamole_new_cpp =
#> epidemic_vacamole_cpp(population = uk_population, : less accurate nanosecond
#> times to avoid potential integer overflows
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> vacamole_new_cpp 7.920052 8.636712 8.831481 8.753028 8.850362 13.34374 100
bench::mark(
"vacamole_new_cpp" = epidemic_vacamole_cpp(
population = uk_population,
infection = infect,
vaccination = double_vaccination,
time_end = 500, increment = 1
)
)
#> # A tibble: 1 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 vacamole_new_cpp 4.67ms 8.04ms 142. 894KB 4.12
Created on 2023-08-09 by the reprex package (v2.0.1)
Closing in favour of moving this to a dedicated Discussion item in the repo.
This is an open Discussion issue to post benchmarks of different implementations of {epidemics}.