Merge pull request #264 from mrc-ide/bug/min_wait

Fix bug with pev min_wait

NEWS.md

@@ -2,6 +2,8 @@
 
  * Fix bug in competing hazards between mass and EPI vaccines. Where individuals
    can be enrolled onto both strategies if applied on the same timestep.
+ * Fix bug with min_wait. Min wait was working off of the final primary dose. It
+   now works of of the first dose.
 
 # malariasimulation 1.6.0
 

R/human_infection.R

@@ -129,9 +129,10 @@ calculate_infections <- function(
 
   # calculate vaccine efficacy
   vaccine_efficacy <- rep(0, length(source_vector))
-  vaccine_times <- variables$pev_timestep$get_values(source_vector)
-  vaccinated <- vaccine_times > -1
+  vaccine_times <- variables$last_eff_pev_timestep$get_values(source_vector)
   pev_profile <- variables$pev_profile$get_values(source_vector)
+  # get vector of individuals who have received their 3rd dose
+  vaccinated <- vaccine_times > -1
   pev_profile <- pev_profile[vaccinated]
   if (length(vaccinated) > 0) {
     antibodies <- calculate_pev_antibodies(

R/mortality_processes.R

@@ -106,7 +106,8 @@ reset_target <- function(variables, events, target, state, timestep) {
     variables$drug_time$queue_update(-1, target)
 
     # vaccination
-    variables$pev_timestep$queue_update(-1, target)
+    variables$last_pev_timestep$queue_update(-1, target)
+    variables$last_eff_pev_timestep$queue_update(-1, target)
     variables$pev_profile$queue_update(-1, target)
     variables$tbv_vaccinated$queue_update(-1, target)
 

R/pev.R

@@ -40,7 +40,7 @@ create_epi_pev_process <- function(
     if (parameters$pev_epi_min_wait == 0) {
       target <- to_vaccinate$to_vector()
     } else {
-      not_recently_vaccinated <- variables$pev_timestep$get_index_of(
+      not_recently_vaccinated <- variables$last_pev_timestep$get_index_of(
         a = max(timestep - parameters$pev_epi_min_wait, 0),
         b = timestep
       )$not(TRUE)
@@ -56,6 +56,9 @@ create_epi_pev_process <- function(
       )
     ]
 
+    # Update the latest vaccination time
+    variables$last_pev_timestep$queue_update(timestep, target)
+
     schedule_vaccination(
       target,
       events,
@@ -91,7 +94,7 @@ create_mass_pev_listener <- function(
     if (parameters$mass_pev_min_wait == 0) {
       target <- in_age_group
     } else {
-      not_recently_vaccinated <- variables$pev_timestep$get_index_of(
+      not_recently_vaccinated <- variables$last_pev_timestep$get_index_of(
         a = max(timestep - parameters$mass_pev_min_wait, 0),
         b = timestep
       )$not(TRUE)
@@ -116,6 +119,11 @@ create_mass_pev_listener <- function(
         correlations
       )
     ]
+
+    # Update the latest vaccination time
+    variables$last_pev_timestep$queue_update(timestep, target)
+
+    # Schedule future doses
     schedule_vaccination(
       target,
       events,
@@ -162,7 +170,7 @@ schedule_vaccination <- function(
 create_pev_efficacy_listener <- function(variables, pev_profile_index) {
   function(timestep, target) {
     if (target$size() > 0) {
-      variables$pev_timestep$queue_update(timestep, target)
+      variables$last_eff_pev_timestep$queue_update(timestep, target)
       variables$pev_profile$queue_update(pev_profile_index, target)
     }
   }
@@ -185,7 +193,8 @@ create_pev_booster_listener <- function(
   force(coverage)
   function(timestep, target) {
     target <- sample_bitset(target, coverage)
-    variables$pev_timestep$queue_update(timestep, target)
+    variables$last_pev_timestep$queue_update(timestep, target)
+    variables$last_eff_pev_timestep$queue_update(timestep, target)
     variables$pev_profile$queue_update(pev_profile_index, target)
     renderer$render(render_name, target$size(), timestep)
 
@@ -240,6 +249,12 @@ attach_pev_dose_listeners <- function(
     dose_events[[d]]$add_listener(
       create_dosage_renderer(renderer, strategy, d)
     )
+    # update last vaccination on every primary dose
+    dose_events[[d]]$add_listener(
+      function(t, target) {
+        variables$last_pev_timestep$queue_update(t, target)
+      }
+    )
     if (d == length(dose_events)) {
       dose_events[[d]]$add_listener(
         create_pev_efficacy_listener(

R/variables.R

@@ -1,5 +1,4 @@
-#' @title Define model variables
-#' @description
+#' @title Define model variables #' @description
 #' create_variables creates the human and mosquito variables for
 #' the model. Variables are used to track real data for each individual over
 #' time, they are read and updated by processes
@@ -18,10 +17,13 @@
 #' * ID - Acquired immunity to detectability
 #' * zeta - Heterogeneity of human individuals
 #' * zeta_group - Discretised heterogeneity of human individuals
-#' * pev_timestep - The timestep of the last pev vaccination (-1 if there
-#' haven't been any)
-#' * pev_profile - The index of the profile of the last administered pev vaccine 
-#' (-1 if there haven't been any)
+#' * last_pev_timestep - The timestep of the last pev vaccination (-1 if there
+#' * last_eff_pev_timestep - The timestep of the last efficacious pev
+#' vaccination, including final primary dose and booster doses (-1 if there have not been any)
+#' * pev_profile - The index of the efficacy profile of any pev vaccinations.
+#' Not set until the final primary dose.
+#' This is only set on the final primary dose and subsequent booster doses
+#' (-1 otherwise)
 #' * tbv_vaccinated - The timstep of the last tbv vaccination (-1 if there
 #' haven't been any
 #' * net_time - The timestep when a net was last put up (-1 if never)
@@ -188,7 +190,8 @@ create_variables <- function(parameters) {
   drug <- individual::IntegerVariable$new(rep(0, size))
   drug_time <- individual::IntegerVariable$new(rep(-1, size))
 
-  pev_timestep <- individual::IntegerVariable$new(rep(-1, size))
+  last_pev_timestep <- individual::IntegerVariable$new(rep(-1, size))
+  last_eff_pev_timestep <- individual::IntegerVariable$new(rep(-1, size))
   pev_profile <- individual::IntegerVariable$new(rep(-1, size))
 
   tbv_vaccinated <- individual::DoubleVariable$new(rep(-1, size))
@@ -215,7 +218,8 @@ create_variables <- function(parameters) {
     infectivity = infectivity,
     drug = drug,
     drug_time = drug_time,
-    pev_timestep = pev_timestep,
+    last_pev_timestep = last_pev_timestep,
+    last_eff_pev_timestep = last_eff_pev_timestep,
     pev_profile = pev_profile,
     tbv_vaccinated = tbv_vaccinated,
     net_time = net_time,

tests/testthat/test-infection-integration.R

@@ -133,7 +133,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     ),
     drug = individual::DoubleVariable$new(c(1, 2, 0, 0)),
     drug_time = individual::DoubleVariable$new(c(20, 30, -1, -1)),
-    pev_timestep = individual::DoubleVariable$new(c(-1, 10, 40, -1)),
+    last_eff_pev_timestep = individual::DoubleVariable$new(c(-1, 10, 40, -1)),
     pev_profile = individual::IntegerVariable$new(c(-1, 1, 2, -1)),
     ib = individual::DoubleVariable$new(c(.2, .3, .5, .9))
   )
@@ -366,7 +366,7 @@ test_that('prophylaxis is considered for medicated humans', {
     ),
     drug = individual::DoubleVariable$new(c(0, 2, 1, 0)),
     drug_time = individual::DoubleVariable$new(c(-1, 49, 40, -1)),
-    pev_timestep = individual::DoubleVariable$new(c(-1, -1, -1, -1)),
+    last_eff_pev_timestep = individual::DoubleVariable$new(c(-1, -1, -1, -1)),
     pev_profile = individual::IntegerVariable$new(c(-1, -1, -1, -1)),
     ib = individual::DoubleVariable$new(c(.2, .3, .5, .9))
   )

tests/testthat/test-pev-epi.R

@@ -86,49 +86,48 @@ test_that('pev epi targets correct age and respects min_wait', {
   variables$birth <- individual::IntegerVariable$new(
     -c(18, 18, 2.9, 18, 18) * 365 + timestep
   )
-  variables$pev_timestep <- mock_integer(
+  variables$last_pev_timestep <- mock_integer(
     c(50, -1, -1, 4*365, -1)
   )
+  variables$pev_profile <- mock_integer(
+    c(1, -1, -1, 1, -1)
+  )
 
-  events$pev_epi_doses <- lapply(events$pev_epi_doses, mock_event)
-
+  correlations <- get_correlation_parameters(parameters)
   process <- create_epi_pev_process(
     variables,
     events,
     parameters,
-    get_correlation_parameters(parameters),
+    correlations,
     parameters$pev_epi_coverages,
     parameters$pev_epi_timesteps
   )
 
+  sample_mock <- mockery::mock(c(TRUE, TRUE, FALSE))
   mockery::stub(
     process,
     'sample_intervention',
-    mockery::mock(c(TRUE, TRUE, FALSE))
+    sample_mock
   )
 
   process(timestep)
 
   mockery::expect_args(
-    events$pev_epi_doses[[1]]$schedule,
+    sample_mock,
     1,
-    c(1, 2),
-    parameters$pev_doses[[1]]
+    c(1, 2, 5),
+    'pev',
+    .8,
+    correlations
   )
 
   mockery::expect_args(
-    events$pev_epi_doses[[2]]$schedule,
+    variables$last_pev_timestep$queue_update_mock(),
     1,
-    c(1, 2),
-    parameters$pev_doses[[2]]
+    timestep,
+    c(1, 2)
   )
 
-  mockery::expect_args(
-    events$pev_epi_doses[[3]]$schedule,
-    1,
-    c(1, 2),
-    parameters$pev_doses[[3]]
-  )
 })
 
 test_that('EPI ignores individuals scheduled for mass vaccination', {

tests/testthat/test-pev.R

@@ -97,7 +97,7 @@ test_that('Infection considers pev efficacy', {
   variables$birth <- individual::IntegerVariable$new(
     -c(8, 2.9, 3.2, 18.4) * 365 - 100
   )
-  variables$pev_timestep <- individual::IntegerVariable$new(
+  variables$last_eff_pev_timestep <- individual::IntegerVariable$new(
     c(-1, -1, 50, 50 + 30)
   )
   variables$pev_profile <- individual::IntegerVariable$new(
@@ -299,7 +299,10 @@ test_that('Mass boosters update profile params and reschedule correctly', {
   variables$birth <- individual::IntegerVariable$new(
     -c(2.9, 3.2, 18.4) * 365 + 100
   )
-  variables$pev_timestep <- mock_double(
+  variables$last_pev_timestep <- mock_double(
+    c(50, 50, 50)
+  )
+  variables$last_eff_pev_timestep <- mock_double(
     c(50, 50, 50)
   )
   variables$pev_profile <- mock_integer(
@@ -326,7 +329,13 @@ test_that('Mass boosters update profile params and reschedule correctly', {
   listener(timestep, individual::Bitset$new(3)$insert(c(1, 2, 3)))
 
   expect_bitset_update(
-    variables$pev_timestep$queue_update_mock(),
+    variables$last_pev_timestep$queue_update_mock(),
+    timestep,
+    c(1, 2, 3)
+  )
+
+  expect_bitset_update(
+    variables$last_eff_pev_timestep$queue_update_mock(),
     timestep,
     c(1, 2, 3)
   )
@@ -368,7 +377,10 @@ test_that('Mass booster coverages sample subpopulations correctly', {
   variables$birth <- individual::IntegerVariable$new(
     -c(2.9, 3.2, 18.4) * 365 + 100
   )
-  variables$pev_timestep <- mock_double(
+  variables$last_pev_timestep <- mock_double(
+    c(50, 50, 50)
+  )
+  variables$last_eff_pev_timestep <- mock_double(
     c(50, 50, 50)
   )
   variables$pev_profile <- mock_integer(
@@ -395,7 +407,13 @@ test_that('Mass booster coverages sample subpopulations correctly', {
   mockery::expect_args(sample_mock, 1, target, .9)
 
   expect_bitset_update(
-    variables$pev_timestep$queue_update_mock(),
+    variables$last_pev_timestep$queue_update_mock(),
+    timestep,
+    c(2, 3)
+  )
+
+  expect_bitset_update(
+    variables$last_eff_pev_timestep$queue_update_mock(),
     timestep,
     c(2, 3)
   )
@@ -413,6 +431,68 @@ test_that('Mass booster coverages sample subpopulations correctly', {
   )
 })
 
+test_that('mass pev targets correct age and respects min_wait', {
+  timestep <- 5*365 
+  parameters <- get_parameters(list(human_population = 5))
+  parameters <- set_mass_pev(
+    parameters,
+    profile = rtss_profile,
+    timesteps = c(4, 5) * 365,
+    coverages = c(0.8, 0.8),
+    min_ages = 0,
+    max_ages = 19 * 365,
+    min_wait = 2*365,
+    booster_timestep = c(1, 6) * 30,
+    booster_coverage = c(.9, .8),
+    booster_profile = list(rtss_booster_profile, rtss_booster_profile)
+  )
+  events <- create_events(parameters)
+  variables <- create_variables(parameters)
+  variables$birth <- individual::IntegerVariable$new(
+    -c(18, 18, 30, 18, 18) * 365 + timestep
+  )
+  variables$last_pev_timestep <- mock_integer(
+    c(50, -1, -1, 4*365, -1)
+  )
+
+  variables$pev_profile <- mock_integer(
+    c(1, -1, -1, 1, -1)
+  )
+
+  correlations <- get_correlation_parameters(parameters)
+  listener <- create_mass_pev_listener(
+    variables,
+    events,
+    parameters,
+    get_correlation_parameters(parameters)
+  )
+
+  sample_mock <- mockery::mock(c(TRUE, TRUE, FALSE))
+  mockery::stub(
+    listener,
+    'sample_intervention',
+    sample_mock
+  )
+
+  listener(timestep)
+
+  mockery::expect_args(
+    sample_mock,
+    1,
+    c(1, 2, 5),
+    'pev',
+    .8,
+    correlations
+  )
+
+  mockery::expect_args(
+    variables$last_pev_timestep$queue_update_mock(),
+    1,
+    timestep,
+    c(1, 2)
+  )
+})
+
 test_that('Mass efficacy listener works correctly', {
   timestep <- 50
   parameters <- get_parameters()
@@ -430,15 +510,15 @@ test_that('Mass efficacy listener works correctly', {
   )
 
   variables <- create_variables(parameters)
-  variables$pev_timestep <- mock_integer(c(-1, -1, -1))
+  variables$last_eff_pev_timestep <- mock_integer(c(-1, -1, -1))
   variables$pev_profile <- mock_integer(c(-1, -1, -1))
   listener <- create_pev_efficacy_listener(variables, 1)
 
   listener(timestep, individual::Bitset$new(3)$insert(c(1, 2, 3)))
 
   # vaccinated time
   expect_bitset_update(
-    variables$pev_timestep$queue_update_mock(),
+    variables$last_eff_pev_timestep$queue_update_mock(),
     timestep,
     c(1, 2, 3)
   )