denim vs deSolve

Comparison between deSolve and denim

Using deSolve, we can replicate Erlang distribution and exponential distribution for testing

library(denim)
library(deSolve)

# --- 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

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

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

Execution time comparison

denim takes approximately 10.08 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),]

# ---- 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))