Fix vaccination leading to negative susceptibles

[pratikunterwegs]

This PR fixes #226. A reprex copied from #224 where this issue was first reported is attached.

library(epidemics)
library(socialmixr)
#> 
#> Attaching package: 'socialmixr'
#> The following object is masked from 'package:utils':
#> 
#>     cite
library(tidyverse)

preinfectious_period <- 3.0
infectious_period <- 7.0
basic_reproduction <- 1.46

infectiousness_rate <- 1.0/preinfectious_period
recovery_rate <- 1.0/infectious_period
transmission_rate <- basic_reproduction/infectious_period

polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "Italy",
  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 without age information. To change this behaviour, set the 'missing.participant.age' option
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option

# prepare contact matrix
contact_matrix <- t(contact_data$matrix)

# prepare the demography vector
demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)

# initial conditions: one in every 1 million is infected
initial_i <- 1e-6
initial_conditions <- c(
  S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
)

initial_conditions <- matrix(
  rep(initial_conditions, dim(contact_matrix)[1]),
  ncol = 5, byrow = TRUE
)

italy_population <- population(
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = initial_conditions
)

vaccinate <- vaccination(
  name = "vaccinate all",
  time_begin = matrix(200, nrow(contact_matrix)),
  time_end = matrix(200 + 150, nrow(contact_matrix)),
  nu = matrix(c(0.01, 0.0, 0.01))
)

intervention_vaccine <- model_default(
  # population
  population = italy_population,
  # rates
  transmission_rate = transmission_rate,
  infectiousness_rate = infectiousness_rate,
  recovery_rate = recovery_rate,
  # intervention
  vaccination = vaccinate,
  # time
  time_end = 600,increment = 1
)

intervention_vaccine %>% 
  filter(compartment == "infectious") %>% 
  ggplot(aes(x = time, y = value, colour = demography_group)) +
  geom_line()

^{Created on 2024-05-24 with reprex v2.0.2}

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 68696ea is merged into main:

• ✔️default_ode: 10.7ms -> 10.7ms [-0.86%, +0.79%]
• ✔️default_ode_interventions: 96.8ms -> 96.5ms [-1.44%, +0.88%]
• ✔️default_ode_param_vec: 825ms -> 822ms [-0.95%, +0.22%]
• ✔️default_ode_paramvec_intervs: 6.17s -> 6.18s [-0.06%, +0.64%]
• ✔️ebola: 552ms -> 550ms [-1.27%, +0.68%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 68696ea is merged into main:

• ✔️default_ode: 11.1ms -> 11.1ms [-0.82%, +2.02%]
• ✔️default_ode_interventions: 98.7ms -> 98.8ms [-1.06%, +1.12%]
• ✔️default_ode_param_vec: 832ms -> 836ms [-0.04%, +0.92%]
• ✔️default_ode_paramvec_intervs: 6.26s -> 6.27s [-0.36%, +0.67%]
• ✔️ebola: 560ms -> 560ms [-0.63%, +0.66%]
  Further explanation regarding interpretation and methodology can be found in the documentation.