Comparison between deSolve and denim
Using deSolve, we can replicate Erlang
distribution and exponential distribution for
testing
library(denim)
library(deSolve)
## Warning: package 'deSolve' was built under R version 4.3.1
# --- Transition def for denim
transitions <- list(
"S -> I" = d_exponential(0.2),
"I -> R" = d_gamma(3, 2)
)
parameters <- c(rate = 0.2, scale = 3, shape=2)
initialValues <- c(S = 999, I = 1, I1 = 1, I2=0, R=0)
# --- Transition def for deSolve
transition_func <- function(t, state, param){
with(as.list( c(state, param) ), {
gamma_rate = 1/scale
dS = -rate*S
# apply linear chain trick
dI1 = rate*S - gamma_rate*I1
dI2 = gamma_rate*I1 - gamma_rate*I2
dI = dI1 + dI2
dR = gamma_rate*I2
list(c(dS, dI, dI1, dI2, dR))
})
}
# --- Timestep definition
simulationDuration <- 20
timestep <- 0.001 # small timestep required for comparison
denim_vs_deSolve.R
Run simulation with denim
denim_start <- Sys.time()
mod <- sim(transitions = transitions, initialValues = initialValues, parameters, simulationDuration = simulationDuration, timeStep = timestep)
denim_end <- Sys.time()
# --- show output
head(mod[mod$Time %in% 1:simulationDuration,])
## Time S I R
## 1001 1 817.9120 179.0627 3.025308
## 2001 2 669.6497 310.8811 19.469173
## 3001 3 548.2628 398.4336 53.303539
## 4001 4 448.8796 448.0932 103.027204
## 5001 5 367.5116 467.6504 164.838060
## 6001 6 300.8930 464.7307 234.376244
denim_vs_deSolve.R
Run simulation with deSolve
times <- seq(0, simulationDuration, timestep)
desolve_start <- Sys.time()
ode_mod <- ode(y = initialValues, times = times, parms = parameters, func = transition_func)
desolve_end <- Sys.time()
# --- show output
ode_mod <- as.data.frame(ode_mod)
head(ode_mod[ode_mod$time %in% 1:simulationDuration, c("time", "S", "I", "R")])
## time S I R
## 1001 1 817.9120 179.0585 3.029466
## 2001 2 669.6497 310.8686 19.481654
## 3001 3 548.2628 398.4125 53.324630
## 4001 4 448.8796 448.0650 103.055392
## 5001 5 367.5116 467.6172 164.871207
## 6001 6 300.8930 464.6948 234.412204
denim_vs_deSolve.R
Execution time comparison
denim takes approximately 16.38 times as long as
deSolve to compute the result with the given specifications
.
This significant difference can be attributed to the difference in
approaches: deSolve solves a system of ODEs while
denim iterates through each timestep and updates the
population in each compartment
While the approach in denim allow more flexibility in
types of dwell time distributions, the computation time scales up as
timestep grows smaller (in O(n) time complexity).
Plot the result
# increase timestep before plotting
mod <- mod[mod$Time %in% seq(0, simulationDuration, 0.2),]
ode_mod <- ode_mod[ode_mod$time %in% seq(0, simulationDuration, 0.2),]
denim_vs_deSolve.R
# ---- Plot S compartment
plot(x = mod$Time, y = mod$S,xlab = "Time", ylab = "Count", main="S compartment",
col = "#4876ff", type="l", lwd=3)
lines(ode_mod$time, ode_mod$S, lwd=3, lty=3)
legend(x = 15, y = 900,legend=c("denim", "deSolve"), col = c("#4876ff", "black"), lty=c(1,3))
# ---- Plot I compartment
plot(x = mod$Time, y = mod$I, xlab = "Time", ylab = "Count", main="I compartment",
col = "#4876ff", type="l", lwd=2)
lines(ode_mod$time, ode_mod$I, lwd=3, lty=3)
legend(x = 15, y = 350,legend=c("denim", "deSolve"), col = c("#4876ff", "black"), lty=c(1,3))
# ---- Plot R compartment
plot(x = mod$Time, y = mod$R, xlab = "Time", ylab = "Count", main="R compartment",
col = "#4876ff", type="l", lwd=2)
lines(ode_mod$time, ode_mod$R, lwd=3, lty=3)
legend(x = 15, y = 300,legend=c("denim", "deSolve"), col = c("#4876ff", "black"), lty=c(1,3))
denim_vs_deSolve.R
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.