Merge pull request #334 from mrc-ide/progression_rates_mass_admin

MDA and PMC code now consistent with each other, correctly update individual variables and take into account antimalarial resistance.

R/disease_progression.R

@@ -1,34 +1,34 @@
-#' @title Calculate recovery rates
-#' @description Calculates recovery rates for each individual in the population
+#' @title Calculate disease progression rates
+#' @description Calculates disease progression rates for each individual in the population
 #' for storage in competing hazards object and subsequent resolution
 #'
 #' @param variables the available human variables
-#' @param recovery_outcome competing hazards object for recovery rates
+#' @param progression_outcome competing hazards object for disease progression rates
 #' @noRd
-create_recovery_rates_process <- function(
+create_progression_rates_process <- function(
   variables,
-  recovery_outcome
+  progression_outcome
 ) {
   function(timestep){
     target <- variables$state$get_index_of("S")$not()
-    recovery_outcome$set_rates(
+    progression_outcome$set_rates(
       target,
-      variables$recovery_rates$get_values(target))
+      variables$progression_rates$get_values(target))
   }
 }
 
 
-#' @title Disease progression outcomes (recovery)
+#' @title Disease progression outcomes
 #' @description Following resolution of competing hazards, update state and
 #' infectivity of sampled individuals
 #'
 #' @param timestep the current timestep
-#' @param target the sampled recovering individuals
+#' @param target the sampled progressing individuals
 #' @param variables the available human variables
 #' @param parameters model parameters
-#' @param renderer competing hazards object for recovery rates
+#' @param renderer competing hazards object for disease progression rates
 #' @noRd
-recovery_outcome_process <- function(
+progression_outcome_process <- function(
     timestep,
     target,
     variables,
@@ -48,7 +48,7 @@ recovery_outcome_process <- function(
     "U",
     variables$infectivity,
     parameters$cu,
-    variables$recovery_rates,
+    variables$progression_rates,
     1/parameters$du,
     variables$state$get_index_of("A")$and(target)
   )
@@ -58,7 +58,7 @@ recovery_outcome_process <- function(
     "S",
     variables$infectivity,
     0,
-    variables$recovery_rates,
+    variables$progression_rates,
     0,
     variables$state$get_index_of(c("U","Tr"))$and(target)
   )
@@ -73,19 +73,21 @@ recovery_outcome_process <- function(
 #' @param to_state the destination disease state
 #' @param infectivity the handle for the infectivity variable
 #' @param new_infectivity the new infectivity of the progressed individuals
+#' @param progression_rates the handle for the progression_rates variable
+#' @param new_progression the new disease progression rate of the progressed individuals
 #' @noRd
 update_infection <- function(
     state,
     to_state,
     infectivity,
     new_infectivity,
-    recovery_rates,
-    new_recovery_rate,
+    progression_rates,
+    new_progression_rate,
     to_move
 ) {
   state$queue_update(to_state, to_move)
   infectivity$queue_update(new_infectivity, to_move)
-  recovery_rates$queue_update(new_recovery_rate, to_move)
+  progression_rates$queue_update(new_progression_rate, to_move)
 }
 
 #' @title Modelling the progression to asymptomatic disease
@@ -115,7 +117,7 @@ update_to_asymptomatic_infection <- function(
       new_infectivity,
       to_move
     )
-    variables$recovery_rates$queue_update(
+    variables$progression_rates$queue_update(
       1/parameters$da,
       to_move
     )

R/human_infection.R

@@ -347,13 +347,73 @@ calculate_treated <- function(
     )
   ])
 
+  successfully_treated <- calculate_successful_treatments(
+    parameters,
+    seek_treatment,
+    drugs,
+    timestep,
+    renderer,
+    ""
+  )
+
+  if (successfully_treated$successfully_treated$size() > 0) {
+
+    if(parameters$antimalarial_resistance) {
+      dt_update_vector <- successfully_treated$dt_spc_combined
+    } else {
+      dt_update_vector <- parameters$dt
+    }
+
+    update_infection(
+      variables$state,
+      'Tr',
+      variables$infectivity,
+      parameters$cd * parameters$drug_rel_c[successfully_treated$drugs],
+      variables$progression_rates,
+      1/dt_update_vector,
+      successfully_treated$successfully_treated
+    )
+
+    variables$drug$queue_update(
+      successfully_treated$drugs,
+      successfully_treated$successfully_treated
+    )
+    variables$drug_time$queue_update(
+      timestep,
+      successfully_treated$successfully_treated
+    )
+  }
+
+  successfully_treated$successfully_treated
+
+}
+
+
+#' @title Calculate successfully treated humans
+#' @description
+#' Sample successful treatments based on drug efficacy and antimalarial resistance
+#' @param parameters model parameters
+#' @param target bitset of treated humans
+#' @param drugs drug index
+#' @param timestep the current timestep
+#' @param renderer simulation renderer
+#' @param int_name the intervention name to use for rendering, use "" for frontline treatment
+#' @noRd
+calculate_successful_treatments <- function(
+    parameters,
+    target,
+    drugs,
+    timestep,
+    renderer,
+    int_name){
+
   #+++ DRUG EFFICACY +++#
   #+++++++++++++++++++++#
   effectively_treated_index <- bernoulli_multi_p(parameters$drug_efficacy[drugs])
-  effectively_treated <- bitset_at(seek_treatment, effectively_treated_index)
+  effectively_treated <- bitset_at(target, effectively_treated_index)
   drugs <- drugs[effectively_treated_index]
-  n_drug_efficacy_failures <- n_treat - effectively_treated$size()
-  renderer$render('n_drug_efficacy_failures', n_drug_efficacy_failures, timestep)
+  n_drug_efficacy_failures <- target$size() - effectively_treated$size()
+  renderer$render(paste0('n_', int_name, 'drug_efficacy_failures'), n_drug_efficacy_failures, timestep)
 
   #+++ ANTIMALARIAL RESISTANCE +++#
   #+++++++++++++++++++++++++++++++#
@@ -370,7 +430,7 @@ calculate_treated <- function(
     successfully_treated_indices <- bernoulli_multi_p(p = 1 - early_treatment_failure_probability)
     successfully_treated <- bitset_at(effectively_treated, successfully_treated_indices)
     n_early_treatment_failure <- effectively_treated$size() - successfully_treated$size()
-    renderer$render('n_early_treatment_failure', n_early_treatment_failure, timestep)
+    renderer$render(paste0('n_', int_name, 'early_treatment_failure'), n_early_treatment_failure, timestep)
     drugs <- drugs[successfully_treated_indices]
     dt_slow_parasite_clearance <- resistance_parameters$dt_slow_parasite_clearance[successfully_treated_indices]
 
@@ -380,51 +440,30 @@ calculate_treated <- function(
       resistance_parameters$slow_parasite_clearance_probability[successfully_treated_indices]
     slow_parasite_clearance_indices <- bernoulli_multi_p(p = slow_parasite_clearance_probability)
     slow_parasite_clearance_individuals <- bitset_at(successfully_treated, slow_parasite_clearance_indices)
-    renderer$render('n_slow_parasite_clearance', slow_parasite_clearance_individuals$size(), timestep)
+    renderer$render(paste0('n_', int_name, 'slow_parasite_clearance'), slow_parasite_clearance_individuals$size(), timestep)
     non_slow_parasite_clearance_individuals <- successfully_treated$copy()$set_difference(slow_parasite_clearance_individuals)
-    renderer$render('n_successfully_treated', successfully_treated$size(), timestep)
+    renderer$render(paste0('n_', int_name, 'successfully_treated'), successfully_treated$size(), timestep)
     dt_slow_parasite_clearance <- dt_slow_parasite_clearance[slow_parasite_clearance_indices]
 
+    dt_spc_combined <- rep(parameters$dt, length(successfully_treated_indices))
+    dt_spc_combined[slow_parasite_clearance_indices] <- dt_slow_parasite_clearance
+
+    successfully_treated_list <- list(
+      drugs = drugs,
+      successfully_treated = successfully_treated,
+      dt_spc_combined = dt_spc_combined)
+
   } else {
 
     successfully_treated <- effectively_treated
-    renderer$render('n_successfully_treated', successfully_treated$size(), timestep)
+    renderer$render(paste0('n_', int_name, 'successfully_treated'), successfully_treated$size(), timestep)
+
+    successfully_treated_list <- list(
+      drugs = drugs,
+      successfully_treated = successfully_treated)
 
   }
-
-  if (successfully_treated$size() > 0) {
-    variables$state$queue_update("Tr", successfully_treated)
-    variables$infectivity$queue_update(
-      parameters$cd * parameters$drug_rel_c[drugs],
-      successfully_treated
-    )
-    variables$drug$queue_update(
-      drugs,
-      successfully_treated
-    )
-    variables$drug_time$queue_update(
-      timestep,
-      successfully_treated
-    )
-    if(parameters$antimalarial_resistance) {
-      variables$recovery_rates$queue_update(
-        1/parameters$dt,
-        non_slow_parasite_clearance_individuals
-      )
-      variables$recovery_rates$queue_update(
-        1/dt_slow_parasite_clearance,
-        slow_parasite_clearance_individuals
-      )
-    } else {
-      variables$recovery_rates$queue_update(
-        1/parameters$dt,
-        successfully_treated
-      )
-    }
-  }
-
-  successfully_treated
-
+  successfully_treated_list
 }
 
 #' @title Schedule infections
@@ -444,6 +483,7 @@ schedule_infections <- function(
     parameters,
     timestep
 ) {
+
   included <- treated$not(TRUE)
 
   to_infect <- clinical_infections$and(included)
@@ -457,7 +497,7 @@ schedule_infections <- function(
       'D',
       variables$infectivity,
       parameters$cd,
-      variables$recovery_rates,
+      variables$progression_rates,
       1/parameters$dd,
       to_infect
     )