Implement R only models

[pratikunterwegs]

This PR is a work in progress towards implementing R only versions of the library models.

1. Fixes #79 by adding epidemic_default_r(),
2. Fixes #80 by adding epidemic_vacamole_r(),
3. Fixes #84 by adding a simple vignette with speed benchmarks on the default and epidemic models.

[TimTaylor]

I've had a quick gander at the default R version. I suspect it can be sped up considerably by:

• Taking a closure based approach where constants are pulled from the input outside of the ODE solver function. Otherwise there is a lot of overhead accrued via [[.
• Precalculating vaccination and intervention effects. With vaccinations this should be straightforward but a little more involved with the interventions.
• Removing internal error checking, especially when functions are being called multiple times, e.g. in intervention_on_cm() (that check shouldn't be needed if input's are already sanitized?)

[pratikunterwegs]

Thanks, could you give me an example of what you mean by the first two points?

[TimTaylor]

On point one something similar to:

user_facing <- function(params, ...) {

    ### missing stuff ###

    p1 <- params[[1L]]
    p2 <- params[[2L]]
    solver <- function() {
        dx <- #something something p1 p2 dy
        dy <- # something something p1 p2 dx
        list(dx, dy)  
    }

    ### missing stuff ###

    deSolve::lsoda(init, times, solver, list())
}

[TimTaylor]

on point (2) - IIANM - the number of resultant cumulative intervention contact matrices can be precalculated. They don't need calculating on every timestep but a time can be mapped uniquely to one of these. This would increase memory usage and be a little trickier but I think it could speed up both the R and C++.

[pratikunterwegs]

Here's a reprex on point 1 for the default model. It's got maybe one more step than what you suggest so that I can make use of the pre-written argument preparation function. The intervention on contact matrix function has had input checking removed. Is this what you were thinking of?

devtools::load_all(".")
#> ℹ Loading epidemics
library(deSolve)

# prepare arguments
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "United Kingdom",
  age.limits = c(0, 20, 40),
  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

# Prepare some initial objects
uk_population <- population(
  name = "UK population",
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = matrix(
    c(0.9999, 0.0001, 0, 0, 0),
    nrow = nrow(contact_matrix), ncol = 5L,
    byrow = TRUE
  )
)

# Prepare epi parameters
pandemic <- infection(
  r0 = 3,
  preinfectious_period = 3,
  infectious_period = 7
)

# create a multi intervention
multi_intervention <- c(
  intervention(
    time_begin = 50, time_end = 100,
    contact_reduction = matrix(
      0.2, nrow(contact_matrix), 1
    )
  ),
  intervention(
    time_begin = 70, time_end = 90,
    contact_reduction = matrix(
      0.3, nrow(contact_matrix), 1
    )
  )
)

# create a vaccination regime
vax_regime <- vaccination(
  time_begin = matrix(20, nrow(contact_matrix), 1),
  time_end = matrix(100, nrow(contact_matrix), 1),
  nu = matrix(0.01, nrow(contact_matrix), 1)
)

# Define a closure that generates the SIR model function
create_model <- function(beta, gamma, alpha, contact_matrix, n_age_groups,
                         npi_time_begin, npi_time_end, npi_cr,
                         vax_time_begin, vax_time_end, vax_nu) {
  model <- function(t, y, params) {
    # create a matrix
    y <- matrix(y, nrow = n_age_groups, ncol = 5L, byrow = FALSE)

    # scale the contact matrix if within the intervention period
    contact_matrix_ <- intervention_on_cm(
      t = t,
      cm = contact_matrix,
      time_begin = npi_time_begin,
      time_end = npi_time_end,
      cr = npi_cr
    )

    # modify the vaccination rate depending on the regime
    # the number of doses is already checked before passing
    current_nu <- vax_nu * ((vax_time_begin < t) & (vax_time_end > t))

    # calculate transitions
    sToE <- (beta * y[, 1] * contact_matrix_ %*% y[, 3])
    eToI <- alpha * y[, 2]
    iToR <- gamma * y[, 3]
    sToV <- current_nu * y[, 1]

    # define compartmental changes
    dS <- -sToE - sToV
    dE <- sToE - eToI
    dI <- eToI - iToR
    dR <- iToR
    dV <- sToV

    # return a list
    list(c(dS, dE, dI, dR, dV))
  }

  return(model)
}

model_args = list(
  population = uk_population, infection = pandemic,
  intervention = multi_intervention, vaccination = vax_regime,
  time_end = 100, increment = 1
)

model_args = .check_args_epidemic_default(model_args)
model_args = .prepare_args_epidemic_default(model_args)

# Create the SIR model closure
modelfn <- create_model(
  beta = model_args$beta,
  alpha = model_args$alpha,
  gamma = model_args$gamma,
  contact_matrix = model_args$contact_matrix, 
  n_age_groups = 3, npi_time_begin = model_args$npi_time_begin,
  npi_time_end = model_args$npi_time_end, npi_cr = model_args$npi_cr,
  vax_time_begin = model_args$vax_time_begin,
  vax_time_end = model_args$vax_time_end,
  vax_nu = model_args$vax_nu
)

# Initial conditions
initial <- model_args$initial_state

# Time points for simulation
times <- seq(0, 100, by = 1)

# Solve the ODE using the ode function
solution <- deSolve::lsoda(initial, times, modelfn)

microbenchmark::microbenchmark(
  deSolve::lsoda(initial, times, modelfn)
)
#> Warning in microbenchmark::microbenchmark(deSolve::lsoda(initial, times, : less
#> accurate nanosecond times to avoid potential integer overflows
#> Unit: milliseconds
#>                                     expr      min       lq     mean   median
#>  deSolve::lsoda(initial, times, modelfn) 18.65984 19.09831 20.83849 19.40196
#>        uq      max neval
#>  21.49933 105.3003   100

^{Created on 2023-08-16 by the reprex package (v2.0.1)}

[TimTaylor]

Any timings pre this change?

[TimTaylor]

From self timings, closure approach has only a small difference (comparable to the total runtime of the rcpp version).

[pratikunterwegs]

From self timings, closure approach has only a small difference (comparable to the total runtime of the rcpp version).

Okay thanks - can still move to implementing this if there's interest.

on point (2) - IIANM - the number of resultant cumulative intervention contact matrices can be precalculated. They don't need calculating on every timestep but a time can be mapped uniquely to one of these. This would increase memory usage and be a little trickier but I think it could speed up both the R and C++.

This could also be sorted by using events (in both R and C++). This forms issue #82, and I was intending to make a separate PR as that is likely to be quite a big change (see related issue #83).

[TimTaylor]

This could also be sorted by using events (in both R and C++). This forms issue https://github.com/epiverse-trace/epidemics/issues/82, and I was intending to make a separate PR as that is likely to be quite a big change (see related issue https://github.com/epiverse-trace/epidemics/issues/83).

Cool - I suspect either way that it's implemented it's the fiddly pre-calculation that will make most difference timewise but may not be worth the effort.

[pratikunterwegs]

This could also be sorted by using events (in both R and C++). This forms issue #82, and I was intending to make a separate PR as that is likely to be quite a big change (see related issue #83).

Cool - I suspect either way that it's implemented it's the fiddly pre-calculation that will make most difference timewise but may not be worth the effort.

I'm mistaken - Boost does not offer event handling, per this review.

Working off @rozeggo's feedback on Slack, I'll keep the R models as they are for now. I can explore event handling for the R models using deSolve::events, but the documentation suggests that implementing complex interventions and vaccinations with age-specific effects would not be easy.

If there's agreement on this, I will merge this PR - any feedback for related features or documentation etc. to add to this PR is welcome.

[pratikunterwegs]

Merging this for now to keep development moving, as #90 is in the pipeline already.