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)))
# mean(python_runs)
# 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.4081547 0.4145248 0.4115663 0.4107938 0.4103262 0.4193149 0.4170411
## [8] 0.4241791 0.4211700 0.4172759 0.4152267 0.4141960 0.4105721 0.4113030
## [15] 0.4091239 0.4162629 0.4119470 0.4095368 0.4096603 0.4062281 0.4068928
## [22] 0.4074152 0.4069958 0.4032819 0.4009092 0.4018691 0.4033310 0.4028161
## [29] 0.4023368 0.4104071 0.4114213 0.4117439 0.4253800 0.4138129 0.4080679
## [36] 0.4227841 0.4102170 0.4109051 0.4219439 0.4190531 0.4230320 0.4170420
## [43] 0.4159870 0.4140022 0.4180233 0.4103088 0.4142401 0.4103670 0.4150410
## [50] 0.4198253Median run time for uSEIR approach, Cython implementation: 0.4114938 seconds
deSolve
Model in R
Code for running SEIR in deSolve
## Warning: package 'deSolve' was built under R version 4.3.1
parameters <- c(scale_I = 4, shape_I=2,
scale_R = 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) ), {
gamma_rate_I = 1/scale_I
gamma_rate_R = 1/scale_R
tr = scale_R*shape_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] 3.33ms 3.17ms 3.36ms 3.28ms 5.43ms 3.18ms 3.13ms 3.08ms 3.09ms
## [10] 3.08ms 3.14ms 3.22ms 3.16ms 4.71ms 3.12ms 3.05ms 3.03ms 3.06ms
## [19] 3.05ms 3.06ms 3.07ms 3.08ms 4.57ms 3.1ms 3.1ms 3.07ms 3.06ms
## [28] 3.06ms 3.08ms 3.1ms 3.11ms 4.55ms 3.07ms 3.06ms 3.05ms 3.08ms
## [37] 3.1ms 3.15ms 3.11ms 3.1ms 23.66ms 3.21ms 3.13ms 3.14ms 3.2ms
## [46] 3.14ms 3.13ms 3.19ms 3.21ms 4.99msMedian run time for deSolve, with model defined in R: 0.0031036 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] 164µs 200µs 140µs 174µs 137µs 134µs 156µs 127µs 151µs 127µs 132µs 156µs
## [13] 130µs 149µs 125µs 125µs 147µs 119µs 141µs 120µs 122µs 142µs 120µs 142µs
## [25] 120µs 121µs 144µs 120µs 143µs 123µs 123µs 146µs 120µs 143µs 125µs 120µs
## [37] 144µs 120µs 148µs 122µs 122µs 143µs 119µs 146µs 119µs 119µs 145µs 128µs
## [49] 143µs 122µsMedian run time for deSolve, with model defined in C: 1.30626^{-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] 366µs 195µs 218µs 179µs 177µs 204µs 172µs 190µs 170µs 174µs 191µs 173µs
## [13] 191µs 170µs 171µs 188µs 171µs 189µs 168µs 172µs 194µs 169µs 190µs 168µs
## [25] 168µs 190µs 165µs 196µs 165µs 168µs 190µs 174µs 189µs 167µs 168µs 188µs
## [37] 171µs 189µs 171µs 168µs 199µs 166µs 197µs 169µs 162µs 183µs 168µs 183µs
## [49] 163µs 162µsMedian run time for deSolve, with model defined in Fortran: 1.734505^{-4} seconds
denim
Code for running SEIR in denim
library(denim)
denim_model <- list(
"S -> E" = "(R0/tr) * timeStepDur * S * (I/N)", # formulate according that of uSEIR method
"E -> I" = d_gamma(scale = 4, shape = 2),
"I -> R" = d_gamma(scale = 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, timeStepDur = 0.01)
# ---- Get runtimes ----
denim_runs <- bench::mark(
sim(transitions = denim_model,
initialValues = initialValues,
parameters = parameters,
simulationDuration = sim_duration, timeStep = 0.01),
iterations = total_runs
)
# ---- Get output ----
denim_out <- sim(transitions = denim_model,
initialValues = initialValues,
parameters = parameters,
simulationDuration = sim_duration, timeStep = 0.01)Run time for denim implementation
denim_runs$time## [[1]]
## [1] 819ms 810ms 821ms 820ms 818ms 827ms 818ms 820ms 821ms 812ms 808ms 815ms
## [13] 805ms 814ms 819ms 819ms 815ms 817ms 795ms 818ms 818ms 812ms 821ms 821ms
## [25] 817ms 803ms 817ms 817ms 816ms 818ms 808ms 819ms 826ms 806ms 815ms 808ms
## [37] 805ms 809ms 813ms 823ms 827ms 811ms 847ms 820ms 826ms 834ms 825ms 819ms
## [49] 819ms 830msMedian run time for denim: 0.8180356 seconds
Compare run time
The following plot shows run time for 50 runs (with horizontal line showing median run time) of each approach.
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.