To assess denim’s performance, we benchmark it against other approaches and packages.
The benchmarking process was conducted on a Macbook with M2 Pro chip, 16 GBs of RAM and 10 cores (6 performance and 4 efficiency).
Benchmark settings
All approaches will simulate the following SEIR model, with the same simulation duration of 180
Each approach will be run 50 times
total_runs <- 50L # number of runs
sim_duration <- 180 # duration of simulationuSEIR
Simulate model using uSEIR approach (Hernández et al. 2021)
Source code: https://github.com/jjgomezcadenas/useirn/blob/master/nb/uSEIR.ipynb
load useir implementation in pure Python
## Warning: package 'reticulate' was built under R version 4.3.3
# use_python("/opt/anaconda3/envs/bnn/bin/python", required = TRUE)
use_condaenv(condaenv='bnn', required = TRUE)
matplotlib <- import("matplotlib")
matplotlib$use("Agg", force = TRUE)
py_run_file("../supplements/useir_python.py")Python implementation
Code for running uSEIR in pure Python
import time
import concurrent.futures
import pickle
import os
from statistics import mean
python_runs = []
def get_python_runtime(n):
start = time.time()
df = solve_uSeir(ti_shape = 2,
ti_scale = 4,
tr_shape = 2,
tr_scale = 3,
R0 = 3.5)
end = time.time()
return end - start
# load cached result if available instead of rerun due to long run time
cached_python_runs = "../supplements/python_runs.pkl"
if os.path.exists(cached_python_runs):
# If the file exists, load the Python list from the file
with open(cached_python_runs, 'rb') as f:
python_runs = pickle.load(f)
else:
print("no cache found")
# multithread instead for quicker result
with concurrent.futures.ProcessPoolExecutor(max_workers=8) as executor:
python_runs = list(executor.map(get_python_runtime, range(r.total_runs)))
# Save to cache
with open(cached_python_runs, 'wb') as f:
pickle.dump(python_runs, f)
# plot_useir((df,), ('G',), T = 'uSEIR', figsize=(14,8))
# print(f'python solve_seir call: dr = {end-start}')Run time for uSEIR approach, Python implementation (in seconds)
py$python_runs## [1] 54.79370 54.23037 54.61863 54.13476 54.73439 54.35928 54.48391 54.56821
## [9] 53.26567 53.35437 53.44551 53.31425 53.51092 53.19878 53.92493 53.84889
## [17] 54.05378 53.66922 53.65536 53.54520 53.73345 53.88690 54.23988 54.19460
## [25] 53.64543 53.55217 53.66061 53.67430 53.59397 53.65503 54.32232 53.98363
## [33] 53.53084 53.79021 53.85689 53.45313 53.61043 53.53294 54.00833 54.11180
## [41] 53.94760 53.61952 53.80676 53.74364 53.68474 53.63701 53.59615 53.81793
## [49] 44.74661 44.77127Median run time for uSEIR approach, Python implementation: 53.7090955 seconds
Cython implementation
Code for running uSEIR in Cython (C backend)
# import precompiled cython module
import sys
sys.path.insert(0, "../supplements")
import useir
import time
import pyarrow as pa
cython_runs = []
# --- Get runtime ----
for i in range(r.total_runs):
start = time.time()
df = useir.csolve_uSeir(dist = "gamma",
ti_shape = 2,
ti_scale = 4,
tr_shape = 2,
tr_scale = 3,
R0 = 3.5
)
end = time.time()
cython_runs = cython_runs + [end - start]## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334
## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334# ---- Get output for uSEIR -----
df = useir.csolve_uSeir(dist = "gamma",
ti_shape = 2,
ti_scale = 4,
tr_shape = 2,
tr_scale = 3,
R0 = 3.5
)## Function compute_pde with sampling = Fine, time epsilon = 0.01
## statistical distribution = gamma, ti = 8.0, tr = 6.0
## number of exposed compartments = 4702, infected compartments = 3526
## len pdE = 4702, max(pdE) =0.0009196976461058881 len pdi = 3526 max(pdI) =0.0012262625368251734
## prob = 0.5833333333333334, pn = 0.005833333333333334Run time for uSEIR approach, Cython implementation (in seconds)
py$cython_runs## [1] 0.4006190 0.4024279 0.4006343 0.4025679 0.3997817 0.4068038 0.4022102
## [8] 0.4043179 0.4043519 0.4014809 0.4099109 0.4061849 0.3975039 0.3988829
## [15] 0.4015999 0.4003770 0.4012938 0.4009390 0.4056711 0.4066510 0.4015381
## [22] 0.3969321 0.4054148 0.4033730 0.4022348 0.4026082 0.4015408 0.4071460
## [29] 0.4019132 0.4040194 0.4105749 0.4004490 0.4023907 0.4013791 0.4040160
## [36] 0.4061792 0.4076910 0.4065399 0.4004521 0.3999770 0.4053557 0.4042711
## [43] 0.4062769 0.4112499 0.4052560 0.4052248 0.4040771 0.4084358 0.4083548
## [50] 0.4056518Median run time for uSEIR approach, Cython implementation: 0.4036945 seconds
deSolve
Model in R
Code for running SEIR in deSolve
## Warning: package 'deSolve' was built under R version 4.3.1
parameters <- c(gamma_rate_I = 1/4, shape_I=2,
gamma_rate_R = 1/3, shape_R = 2,
R0 = 3.5, N = 1e6)
initialValues <- c(S = 999999, E1 = 1,
E2 = 0, E = 0, I1=0,
I2=0, I=0, R=0
)
# --- Transition def for deSolve
transition_func <- function(t, state, param){
with(as.list( c(state, param) ), {
tr = shape_R*(1/gamma_rate_R)
dS = - (R0/tr) * S * I/N
# apply linear chain trick
dE1 = (R0/tr) * S * I/N - gamma_rate_I*E1
dE2 = gamma_rate_I*E1 - gamma_rate_I*E2
dE = dE1 + dE2
dI1 = gamma_rate_I*E2 - gamma_rate_R*I1
dI2 = gamma_rate_R*I1 - gamma_rate_R*I2
dI = dI1 + dI2
dR = gamma_rate_R*I2
list(c(dS, dE1, dE2, dE, dI1, dI2, dI, dR))
})
}
times <- seq(0, sim_duration, 1)
# ------ Compute run time ------
desolve_runs <- bench::mark(
ode(y = initialValues, times = times, parms = parameters, func = transition_func),
iterations = total_runs
)
ode_mod <- ode(y = initialValues, times = times, parms = parameters, func = transition_func)
ode_mod <- as.data.frame(ode_mod)Run time for deSolve implementation
desolve_runs$time## [[1]]
## [1] 5.57ms 3.16ms 3.14ms 3.1ms 3.1ms 4.72ms 3.1ms 3.05ms 2.99ms
## [10] 3.04ms 28.32ms 3.19ms 3.09ms 3.1ms 3.1ms 3.06ms 5.43ms 3.1ms
## [19] 3.1ms 3.07ms 3.02ms 4.44ms 3.15ms 3.01ms 2.99ms 2.98ms 3ms
## [28] 4.58ms 3.25ms 3.22ms 3.03ms 3.03ms 3.27ms 5.96ms 3.48ms 3.43ms
## [37] 3.34ms 3.27ms 5.38ms 3.14ms 3.21ms 3.05ms 3.07ms 3.06ms 4.81ms
## [46] 3.1ms 3.06ms 3.04ms 3.03ms 4.33msMedian run time for deSolve, with model defined in R: 0.0030976 seconds
Model in C
Code for running model defined in C
# compile model
# system("R CMD SHLIB supplements/desolve_mod/benchmark_mod.c")
# compiled file on Windows will have .dll extension instead of .so
dyn.load("../supplements/desolve_mod/benchmark_mod.so")
initialValues <- c(S = 999999, E1 = 1,
E2 = 0, E = 0, I1=0,
I2=0, I=0, R=0
)
parameters <- c(R0 = 3.5, scale_I = 4, shape_I=2,
scale_R = 3, shape_R = 2, N = 1e6)
deSolve_c_runs <- bench::mark(
# run model defined in C
ode(initialValues, times, func = "derivs", parms = parameters,
dllname = "benchmark_mod", initfunc = "initmod"),
iterations = total_runs
)
dyn.unload("../supplements/desolve_mod/benchmark_mod.so")Run time for deSolve with model defined in C
deSolve_c_runs$time## [[1]]
## [1] 200µs 244µs 142µs 159µs 118µs 117µs 147µs 112µs 140µs 110µs 109µs 134µs
## [13] 109µs 135µs 115µs 115µs 131µs 112µs 132µs 112µs 115µs 148µs 114µs 139µs
## [25] 124µs 121µs 136µs 115µs 141µs 114µs 115µs 131µs 118µs 136µs 113µs 114µs
## [37] 135µs 114µs 134µs 117µs 115µs 141µs 118µs 135µs 115µs 122µs 141µs 114µs
## [49] 141µs 117µsMedian run time for deSolve, with model defined in C: 1.1916651^{-4} seconds
Model in Fortran
Code for running model defined in C
# compile model in fortran
# system("R CMD SHLIB supplements/desolve_mod/benchmark_mod_fortran.f")
dyn.load("../supplements/desolve_mod/benchmark_mod_fortran.so")
initialValues <- c(S = 999999, E1 = 1,
E2 = 0, E = 0, I1=0,
I2=0, I=0, R=0
)
parameters <- c(R0 = 3.5, scale_I = 4, shape_I=2,
scale_R = 3, shape_R = 2, N = 1e6)
deSolve_fortran_runs <- bench::mark(
# run model defined in C
ode(initialValues, times, func = "derivs", parms = parameters,
dllname = "benchmark_mod_fortran", initfunc = "initmod"),
iterations = total_runs
)
dyn.unload("../supplements/desolve_mod/benchmark_mod_fortran.so")
deSolve_fortran_runs$time## [[1]]
## [1] 368µs 197µs 222µs 184µs 177µs 197µs 171µs 193µs 170µs 177µs 192µs 168µs
## [13] 191µs 169µs 171µs 189µs 167µs 191µs 171µs 171µs 189µs 167µs 199µs 176µs
## [25] 170µs 189µs 170µs 190µs 167µs 169µs 188µs 168µs 186µs 169µs 167µs 186µs
## [37] 167µs 183µs 161µs 163µs 184µs 160µs 187µs 162µs 167µs 186µs 161µs 188µs
## [49] 160µs 163µsMedian run time for deSolve, with model defined in Fortran: 1.7345051^{-4} seconds
odin
We also implement uSEIR model with denim’s algorithm using
odin package for comparison
Code for running SEIR in odin
# ---- Install packages -----
# install.packages(
# "odin2",
# repos = c("https://mrc-ide.r-universe.dev", "https://cloud.r-project.org"))
# install.packages(
# "dust2",
# repos = c("https://mrc-ide.r-universe.dev", "https://cloud.r-project.org"))
library(odin2)
odin_mod <- odin2::odin(
{
# ----- Define algo to update compartments here ---------
update(S) <- S - dt * (R0/tr) * S * sum(I)/N
# --- E compartment ------
update(E[1]) <- dt * (R0/tr) * S * sum(I)/N
# starting from 2: to simulate individuals staying in E for another timestep
update(E[2:e_maxtime]) <- E[i-1]*(1-e_transprob[i-1])
# compute total population from E -> I
dim(E_to_I) <- e_maxtime
E_to_I[1:e_maxtime] <- e_transprob[i]*E[i]
sum_E_to_I <- sum(E_to_I)
# --- I compartment ------
update(I[1]) <- sum_E_to_I
update(I[2:i_maxtime]) <- I[i-1]*(1-i_transprob[i-1])
# compute total population from I -> R
dim(I_to_R) <- i_maxtime
I_to_R[1:i_maxtime] <- i_transprob[i]*I[i]
sum_I_to_R <- sum(I_to_R)
# --- R compartment ------
update(R) <- R + sum_I_to_R
# initialize population from input
initial(S) <- S_init
initial(E[]) <- E_init[i]
dim(E) <- e_maxtime
initial(I[]) <- I_init[i]
dim(I) <- i_maxtime
initial(R) <- R_init
# ----- Inputs -------
R0 <- parameter()
tr <- parameter()
# transition prob of E
e_transprob<- parameter()
e_maxtime <- parameter()
dim(e_transprob) <- e_maxtime
# transition prob of I
i_transprob <- parameter()
i_maxtime <- parameter()
dim(i_transprob) <- i_maxtime
# initial populations
S_init <- user()
E_init <- user()
dim(E_init) <- e_maxtime
I_init <- user()
dim(I_init) <- i_maxtime
R_init <- user()
N <- parameter(1000)
}
)## Warning in odin2::odin({: Found 4 compatibility issues
## Replace calls to 'user()' with 'parameter()'
## ✖ S_init <- user()
## ✔ S_init <- parameter()
## ✖ E_init <- user()
## ✔ E_init <- parameter()
## ✖ I_init <- user()
## ✔ I_init <- parameter()
## ✖ R_init <- user()
## ✔ R_init <- parameter()## ✔ Wrote 'DESCRIPTION'## ✔ Wrote 'NAMESPACE'## ✔ Wrote 'R/dust.R'## ✔ Wrote 'src/dust.cpp'## ✔ Wrote 'src/Makevars'## ℹ 12 functions decorated with [[cpp11::register]]## ✔ generated file cpp11.R## ✔ generated file cpp11.cpp## ℹ Re-compiling odin.system79e020ca## ── R CMD INSTALL ───────────────────────────────────────────────────────────────
## * installing *source* package ‘odin.system79e020ca’ ...
## ** using staged installation
## ** libs
## using C++ compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.3)’
## using SDK: ‘NA’
## clang++ -arch arm64 -std=gnu++17 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -O2 -I'/Users/anhptq/Library/R/arm64/4.3/library/cpp11/include' -I'/Users/anhptq/Library/R/arm64/4.3/library/dust2/include' -I'/Users/anhptq/Library/R/arm64/4.3/library/monty/include' -I/opt/R/arm64/include -DHAVE_INLINE -fPIC -isystem /Library/Developer/CommandLineTools/SDKs/MacOSX15.4.sdk/usr/include/c++/v1 -O2 -c cpp11.cpp -o cpp11.o
## clang++ -arch arm64 -std=gnu++17 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -O2 -I'/Users/anhptq/Library/R/arm64/4.3/library/cpp11/include' -I'/Users/anhptq/Library/R/arm64/4.3/library/dust2/include' -I'/Users/anhptq/Library/R/arm64/4.3/library/monty/include' -I/opt/R/arm64/include -DHAVE_INLINE -fPIC -isystem /Library/Developer/CommandLineTools/SDKs/MacOSX15.4.sdk/usr/include/c++/v1 -O2 -c dust.cpp -o dust.o
## clang++ -arch arm64 -std=gnu++17 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/Library/Frameworks/R.framework/Resources/lib -L/Library/Developer/CommandLineTools/SDKs/MacOSX15.4.sdk/usr/lib -o odin.system79e020ca.so cpp11.o dust.o -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation
## ld: warning: -single_module is obsolete
## ld: warning: -multiply_defined is obsolete
## installing to /private/var/folders/rf/dwxhm19j1ws1mfmsvj9yfp140000gr/T/RtmpOyjPQe/devtools_install_2fe61bd96486/00LOCK-dust_2fe62c73c9dd/00new/odin.system79e020ca/libs
## ** checking absolute paths in shared objects and dynamic libraries
## * DONE (odin.system79e020ca)## ℹ Loading odin.system79e020ca
compute_transprob <- function(dist_func,..., timestep=0.05, error_tolerance=0.001){
maxtime <- timestep
prev_prob <- 0
transprob <- numeric()
cumulative_dist <- numeric()
prob_dist <- numeric()
while(TRUE){
# get current cumulative prob and check whether it is sufficiently close to 1
temp_prob <- ifelse(
dist_func(maxtime, ...) < (1 - error_tolerance),
dist_func(maxtime, ...),
1);
cumulative_dist <- c(cumulative_dist, temp_prob)
# get f(t)
curr_prob <- temp_prob - prev_prob
prob_dist <- c(prob_dist, curr_prob)
# compute transprob
curr_transprob <- curr_prob/(1-prev_prob)
transprob <- c(transprob, curr_transprob)
prev_prob <- temp_prob
maxtime <- maxtime + timestep
if(temp_prob == 1){
break
}
}
data.frame(
prob_dist = prob_dist,
cumulative_dist = cumulative_dist,
transprob = transprob
)
}Run model for bench mark. Note that the process of computing the transition probability is also included as part of the benchmark for a fair comparison with denim.
timeStep <- 0.01
errorTolerance <- 0.001
# ---- Get runtimes ----
odin_runs <- bench::mark(
{
# ---- Compute transprob -----
e_transprob <- compute_transprob(pgamma, rate=1/4, shape=2,
timestep = timeStep, error_tolerance = errorTolerance)$transprob
i_transprob <- compute_transprob(pgamma, rate=1/3, shape=2,
timestep = timeStep, error_tolerance = errorTolerance)$transprob
# ---- Run model and plot -----
# initialize params
odin_pars <- list(
R0 = 3.5,
tr = 3*2, # compute mean recovery time, for gamma it's scale*shape
N = 1e6,
e_transprob = e_transprob,
e_maxtime = length(e_transprob),
i_transprob = i_transprob,
i_maxtime = length(i_transprob),
S_init = 999999,
E_init = array( c(1, rep(0, length(e_transprob) - 1) ) ),
I_init = array( rep(0, length(i_transprob)) ),
R_init = 0
)
# run model
t_seq <- seq(0, sim_duration, 0.25)
odin_seir <- dust2::dust_system_create(odin_mod, odin_pars, dt = timeStep)
dust2::dust_system_set_state_initial(odin_seir)
out <- dust2::dust_system_simulate(odin_seir, t_seq)
out <- dust2::dust_unpack_state(odin_seir, out)
},
iterations = total_runs
)## Warning: Some expressions had a GC in every iteration; so filtering is
## disabled.
odin_out <- data.frame(
t = t_seq,
S = out$S,
E = colSums(out$E),
I = colSums(out$I),
R = out$R
)
odin_runs$time## [[1]]
## [1] 428ms 421ms 315ms 411ms 326ms 333ms 430ms 331ms 409ms 325ms 316ms 324ms
## [13] 392ms 314ms 317ms 393ms 315ms 319ms 317ms 319ms 395ms 314ms 394ms 312ms
## [25] 334ms 319ms 384ms 318ms 388ms 319ms 329ms 315ms 386ms 315ms 320ms 384ms
## [37] 397ms 320ms 319ms 327ms 390ms 314ms 318ms 391ms 316ms 329ms 395ms 314ms
## [49] 324ms 316msMedian run time for odin: 0.3246013 seconds
denim
Parametric
Code for running SEIR in denim
timeStep <- 0.01
errorTolerance <- 0.001
library(denim)
denim_model <- denim_dsl({
S -> E = (R0/tr) * timeStep * S * (I/N) # formulate according that of uSEIR method
E -> I = d_gamma(rate = 1/4, shape = 2)
I -> R = d_gamma(rate = 1/3, shape = 2)
})
initialValues <- c(S = 999999, E = 1, I= 0, R= 0)
parameters <- c(R0 = 3.5,
tr = 3*2, # compute mean recovery time, for gamma it's scale*shape
N = 1e6)
# ---- Get runtimes ----
denim_runs <- bench::mark(
sim(
transitions = denim_model,
initialValues = initialValues,
parameters = parameters,
simulationDuration = sim_duration,
timeStep = timeStep,
errorTolerance = errorTolerance
),
iterations = total_runs
)
# ---- Get output ----
denim_out <- sim(transitions = denim_model,
initialValues = initialValues,
parameters = parameters,
simulationDuration = sim_duration, timeStep = timeStep)Run time for denim implementation
denim_runs$time## [[1]]
## [1] 791ms 802ms 792ms 788ms 775ms 776ms 775ms 776ms 774ms 775ms 776ms 780ms
## [13] 783ms 792ms 797ms 788ms 790ms 797ms 790ms 788ms 784ms 770ms 768ms 771ms
## [25] 770ms 768ms 769ms 770ms 773ms 788ms 783ms 808ms 797ms 791ms 784ms 799ms
## [37] 793ms 768ms 770ms 773ms 770ms 768ms 771ms 770ms 780ms 792ms 784ms 790ms
## [49] 793ms 791msMedian run time for denim: 0.7830995 seconds
Nonparametric
We can also define the SEIR model in denim using function
nonparametric() where the dwell-time distribution will be
pre-computed using the helper function from Section @ref(sec-odin)
Code for running SEIR in denim
timeStep <- 0.01
errorTolerance <- 0.001
denim_nonparametric_model <- denim_dsl({
S -> E = (R0/tr) * timeStep * S * (I/N) # formulate according that of uSEIR method
E -> I = nonparametric(ei_dist) #ei_dist is considered a model parameter
I -> R = nonparametric(ir_dist) #ir_dist is also a model parameter
})
initialValues2 <- c(S = 999999, E = 1, I= 0, R= 0)
ei_dist <- compute_transprob(pgamma, rate = 1/4, shape = 2,
timestep = timeStep, error_tolerance = errorTolerance)$prob_dist
ir_dist <- compute_transprob(pgamma, rate = 1/3, shape = 2,
timestep = timeStep, error_tolerance = errorTolerance)$prob_dist
parameters2 <- list(R0 = 3.5,
tr = 3*2, # compute mean recovery time, for gamma it's scale*shape
N = 1e6,
ei_dist = ei_dist,
ir_dist = ir_dist)
# ---- Get runtimes ----
denim_nonparametric_runs <- bench::mark(
sim(transitions = denim_nonparametric_model,
initialValues = initialValues2,
parameters = parameters2,
simulationDuration = sim_duration, timeStep = timeStep),
iterations = total_runs
)## Warning: Some expressions had a GC in every iteration; so filtering is
## disabled.
# ---- Get output ----
denim_nonparametric_out <- sim(transitions = denim_nonparametric_model,
initialValues = initialValues2,
parameters = parameters2,
simulationDuration = sim_duration, timeStep = timeStep)Run time for denim, using nonparametric()
with pre-computed distribution
denim_nonparametric_runs$time## [[1]]
## [1] 1.25s 1.24s 1.25s 1.22s 1.22s 1.25s 1.16s 1.24s 1.24s 1.26s 1.24s 1.25s
## [13] 1.24s 1.21s 1.13s 1.22s 1.22s 1.22s 1.25s 1.26s 1.24s 1.24s 1.14s 1.23s
## [25] 1.24s 1.24s 1.24s 1.22s 1.24s 1.23s 1.18s 1.25s 1.25s 1.24s 1.25s 1.27s
## [37] 1.22s 1.26s 1.17s 1.23s 1.24s 1.24s 1.25s 1.25s 1.23s 1.24s 1.15s 1.22s
## [49] 1.21s 1.22sMedian run time for denim using nonparametric() with
pre-computed distribution: 1.235576 seconds. This longer run time
compared to parametric approach is due to the overhead of interfacing
large vectors (ei_dist and ir_dist in this
example) between R and C++.For this reason, it is recommended to only
use nonparametric() when the observed distribution cannot
be adequately represented by one of the available parametric
transitions.
Compare run time
The following plot shows run time for 50 runs (with horizontal line showing median run time) of each approach.
Run time scaling in denim
It is worth noting that runtime for denim is also dependent on
duration of time step (timeStep parameter for
sim).
The following plot demonstrates how run time changes as value for
timeStep changes, using the same model for benchmarking. The values for
timeStep being evaluated are
[0.01, 0.02, 0.05, 0.1, 0.25, 0.5, 1].
Hernández, Pilar, Carlos Pena, Alberto Ramos, and Juan José
Gómez-Cadenas. 2021. “A New Formulation of Compartmental Epidemic
Modelling for Arbitrary Distributions of Incubation and Removal
Times.” Edited by Eric Forgoston. PLOS ONE 16 (2):
e0244107. https://doi.org/10.1371/journal.pone.0244107.