biting_process renamed as create_biting_process, where functionality …

…now occurs in biting_process.R. Similarly, create_recovery_rates_process functionality is now in disease_progression.R stash_source_humans in competing hazard resolution is now redundant and has been removed from infection processes. recovery rates are now assigned to a variable which is updated as infections, recoveries (via update_infection) and deaths occur. This also takes into account recovery time of treated due to slow parasite clearance and replaced the dt variable. I've also noticed (I think) that Tr->S and U->S recoveries can be done in a single step, so have attempted to combine these at risk of making the code more complicated! Added a metapopulation incidence rendering wrapper function for initialising multiple renderers. Rewrote create_age_group_renderer for clarity following feedback. Some tests have been adjusted to reflect these changes.
mrc-ide · Jun 24, 2024 · 352fbae · 352fbae
1 parent 1b30562
commit 352fbae
Show file tree

Hide file tree

Showing 12 changed files with 181 additions and 187 deletions.
diff --git a/R/biting_process.R b/R/biting_process.R
@@ -18,7 +18,7 @@
 #' @param infection_outcome competing hazards object for infection rates
 #' @param timestep the current timestep
 #' @noRd
-biting_process <- function(
+create_biting_process <- function(
     renderer,
     solvers,
     models,
@@ -29,37 +29,36 @@ biting_process <- function(
     lagged_eir,
     mixing = 1,
     mixing_index = 1,
-    infection_outcome,
-    timestep
-) {
-
-  age <- get_age(variables$birth$get_values(), timestep)
-
-  bitten_humans <- simulate_bites(
-    renderer,
-    solvers,
-    models,
-    variables,
-    events,
-    age,
-    parameters,
-    timestep,
-    lagged_infectivity,
-    lagged_eir,
-    mixing,
-    mixing_index
-  )
-
-  simulate_infection(
-    variables,
-    events,
-    bitten_humans,
-    age,
-    parameters,
-    timestep,
-    renderer,
     infection_outcome
-  )
+) {
+  function(timestep){
+    age <- get_age(variables$birth$get_values(), timestep)
+    bitten_humans <- simulate_bites(
+      renderer,
+      solvers,
+      models,
+      variables,
+      events,
+      age,
+      parameters,
+      timestep,
+      lagged_infectivity,
+      lagged_eir,
+      mixing,
+      mixing_index
+    )
+
+    simulate_infection(
+      variables,
+      events,
+      bitten_humans,
+      age,
+      parameters,
+      timestep,
+      renderer,
+      infection_outcome
+    )
+  }
 }
 
 #' @importFrom stats rpois

diff --git a/R/competing_hazards.R b/R/competing_hazards.R
@@ -18,15 +18,11 @@ CompetingOutcome <- R6::R6Class(
     set_rates = function(rates){
       self$rates <- rates
     },
-    stash_source_humans = function(source_humans){
-      self$source_humans <- source_humans
-    },
     execute = function(t, target){
       private$targeted_process(t, target)
       self$rates <- rep(0, length(self$rates))
     },
-    rates = NULL,
-    source_humans = NULL
+    rates = NULL
   )
 )
 

diff --git a/R/disease_progression.R b/R/disease_progression.R
@@ -1,28 +1,20 @@
 #' @title Calculate recovery rates
 #' @description Calculates recovery rates for each individual in the population
-#' for storage in competing hazards object and future resolution
+#' for storage in competing hazards object and subsequent resolution
 #'
 #' @param variables the available human variables
-#' @param parameters model parameters
 #' @param recovery_outcome competing hazards object for recovery rates
 #' @noRd
-calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt_input){
-
-  # Get correct input for dt depending on antimalarial resistance
-  if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){
-    dt_v <- dt_input
-  } else if (isTRUE(parameters$antimalarial_resistance)) {
-    dt_v <- dt_input$get_values(variables$state$get_index_of("Tr"))
+create_recovery_rates_process <- function(
+  variables,
+  recovery_outcome
+) {
+  function(timestep){
+    recovery_outcome$set_rates(variables$recovery_rates$get_values())
   }
-
-  recovery_rates <- numeric(length = parameters$human_population)
-  recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd
-  recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da
-  recovery_rates[variables$state$get_index_of("U")$to_vector()] <- 1/parameters$du
-  recovery_rates[variables$state$get_index_of("Tr")$to_vector()] <- 1/dt_v
-  recovery_outcome$set_rates(recovery_rates)
 }
 
+
 #' @title Disease progression outcomes (recovery)
 #' @description Following resolution of competing hazards, update state and
 #' infectivity of sampled individuals
@@ -33,7 +25,7 @@ calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt
 #' @param parameters model parameters
 #' @param renderer competing hazards object for recovery rates
 #' @noRd
-recovery_process_resolved_hazard <- function(
+recovery_outcome_process <- function(
     timestep,
     target,
     variables,
@@ -53,25 +45,31 @@ recovery_process_resolved_hazard <- function(
     "U",
     variables$infectivity,
     parameters$cu,
+    variables$recovery_rates,
+    1/parameters$du,
     variables$state$get_index_of("A")$and(target)
   )
 
+  # Is there a reason we aren't doing this in one step...? These all recover to S...
   update_infection(
     variables$state,
     "S",
     variables$infectivity,
     0,
-    variables$state$get_index_of("U")$and(target)
-  )
-
-  update_infection(
-    variables$state,
-    "S",
-    variables$infectivity,
+    variables$recovery_rates,
     0,
-    variables$state$get_index_of("Tr")$and(target)
+    variables$state$get_index_of(c("U","Tr"))$and(target)
   )
 
+  # update_infection(
+  #   variables$state,
+  #   "S",
+  #   variables$infectivity,
+  #   0,
+  #   variables$recovery_rates,
+  #   0,
+  #   variables$state$get_index_of("Tr")$and(target)
+  # )
 }
 
 #' @title Update the state of an individual as infection events occur
@@ -84,14 +82,17 @@ recovery_process_resolved_hazard <- function(
 #' @param new_infectivity the new infectivity of the progressed individuals
 #' @noRd
 update_infection <- function(
-  state,
-  to_state,
-  infectivity,
-  new_infectivity,
-  to_move
-  ) {
+    state,
+    to_state,
+    infectivity,
+    new_infectivity,
+    recovery_rates,
+    new_recovery_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)
 }
 
 #' @title Modelling the progression to asymptomatic disease
@@ -102,11 +103,11 @@ update_infection <- function(
 #' @param parameters model parameters
 #' @noRd
 update_to_asymptomatic_infection <- function(
-  variables,
-  parameters,
-  timestep,
-  to_move
-  ) {
+    variables,
+    parameters,
+    timestep,
+    to_move
+) {
   if (to_move$size() > 0) {
     variables$state$queue_update('A', to_move)
     new_infectivity <- asymptomatic_infectivity(
@@ -121,5 +122,9 @@ update_to_asymptomatic_infection <- function(
       new_infectivity,
       to_move
     )
+    variables$recovery_rates$queue_update(
+      1/parameters$da,
+      to_move
+    )
   }
 }
diff --git a/R/human_infection.R b/R/human_infection.R
@@ -1,7 +1,7 @@
 #' @title Simulate malaria infection in humans
 #' @description
 #' This function ends with the assignment of rates of infection to the competing
-#' hazard resolution object.  Boosts immunity given infectious bites.
+#' hazard resolution object and boosts immunity given infectious bites.
 #' @param variables a list of all of the model variables
 #' @param events a list of all of the model events
 #' @param bitten_humans a bitset of bitten humans
@@ -60,7 +60,6 @@ calculate_infections <- function(
 ) {
   source_humans <- variables$state$get_index_of(
     c('S', 'A', 'U'))$and(bitten_humans)
-  infection_outcome$stash_source_humans(source_humans)
 
   b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
 
@@ -113,7 +112,7 @@ calculate_infections <- function(
   incidence_probability_renderer(
     variables$birth,
     renderer,
-    source_pop,
+    source_humans,
     prob,
     "inc_",
     parameters$incidence_min_ages,
@@ -122,7 +121,7 @@ calculate_infections <- function(
   )
 
   ## capture infection rates to resolve in competing hazards
-  infection_rates <- numeric(length = parameters$human_population)
+  infection_rates <- rep(0, length = parameters$human_population)
   infection_rates[source_vector] <- prob_to_rate(prob)
   infection_outcome$set_rates(infection_rates)
 }
@@ -133,16 +132,14 @@ calculate_infections <- function(
 #' and treated malaria; and resulting boosts in immunity
 #' @param timestep current timestep
 #' @param infected_humans bitset of infected humans
-#' @param source_humans bitset of humans with infection potential
 #' @param variables a list of all of the model variables
 #' @param renderer model render object
 #' @param parameters model parameters
 #' @param prob vector of population probabilities of infection
 #' @noRd
-infection_process_resolved_hazard <- function(
+infection_outcome_process <- function(
     timestep,
     infected_humans,
-    source_humans,
     variables,
     renderer,
     parameters,
@@ -410,14 +407,19 @@ calculate_treated <- function(
       successfully_treated
     )
     if(parameters$antimalarial_resistance) {
-      variables$dt$queue_update(
-        parameters$dt,
+      variables$recovery_rates$queue_update(
+        1/parameters$dt,
         non_slow_parasite_clearance_individuals
       )
-      variables$dt$queue_update(
-        dt_slow_parasite_clearance,
+      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
+      )
     }
   }
 
@@ -455,6 +457,8 @@ schedule_infections <- function(
       'D',
       variables$infectivity,
       parameters$cd,
+      variables$recovery_rates,
+      1/parameters$dd,
       to_infect
     )
   }
@@ -523,7 +527,7 @@ severe_immunity <- function(age, acquired_immunity, maternal_immunity, parameter
   parameters$theta0 * (parameters$theta1 + (1 - parameters$theta1) / (
     1 + fv * (
       (acquired_immunity + maternal_immunity) / parameters$iv0) ** parameters$kv
-  )
+    )
   )
 }
 

diff --git a/R/model.R b/R/model.R
@@ -220,6 +220,7 @@ run_metapop_simulation <- function(
   variables <- lapply(parameters, create_variables)
   events <- lapply(parameters, create_events)
   renderer <- lapply(parameters, function(.) individual::Render$new(timesteps))
+  populate_metapopulation_incidence_rendering_columns(renderer, parameters)
   for (i in seq_along(parameters)) {
     # NOTE: forceAndCall is necessary here to make sure i refers to the current
     # iteration

diff --git a/R/mortality_processes.R b/R/mortality_processes.R
@@ -113,12 +113,8 @@ reset_target <- function(variables, events, target, state, parameters, timestep)
 
     # onwards infectiousness
     variables$infectivity$queue_update(0, target)
+    variables$recovery_rates$queue_update(0, target)
 
-    # treated compartment residence time:
-    if(!is.null(variables$dt)) {
-      variables$dt$queue_update(parameters$dt, target)
-    }
-
     # zeta and zeta group and vector controls survive rebirth
   }
 }