Merge branch 'feat/rtss_epi_variable_cov' of https://github.com/mrc-i…

…de/malariasimulation into feat/rtss_epi_variable_cov
mrc-ide · Aug 2, 2022 · 1ba697e · 1ba697e
2 parents aca812b + 3fbba61
commit 1ba697e
Show file tree

Hide file tree

Showing 29 changed files with 767 additions and 277 deletions.
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
@@ -35,7 +35,7 @@ jobs:
         with:
           r-version: ${{ matrix.config.r }}
 
-      - uses: r-lib/actions/setup-pandoc@v1
+      - uses: r-lib/actions/setup-pandoc@v2
 
       - name: Query dependencies
         run: |

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -10,27 +10,27 @@ Authors@R: c(
   ),
   person(
     given = "Peter",
-    family = "Windskill",
+    family = "Winskill",
     role = c('aut'),
-    email = 'p.windskill@ic.ac.uk'
+    email = 'p.winskill@ic.ac.uk'
   ),
   person(
     given = "Hillary",
     family = "Topazian",
     role = c('aut'),
-    email = 'h.topazian@ic.ac.uk'
+    email = 'h.topazian@imperial.ac.uk'
   ),
   person(
     given = "Joseph",
     family = "Challenger",
     role = c('aut'),
-    email = 'j.challenger@ic.ac.uk'
+    email = 'j.challenger@imperial.ac.uk'
   ),
   person(
     given = "Richard",
     family = "Fitzjohn",
     role = c('aut'),
-    email = 'r.fitzjohn@ic.ac.uk'
+    email = 'r.fitzjohn@imperial.ac.uk'
   ),
   person(
     given = "Imperial College of Science, Technology and Medicine",

diff --git a/NAMESPACE b/NAMESPACE
@@ -12,6 +12,7 @@ export(get_parameters)
 export(parameterise_mosquito_equilibrium)
 export(parameterise_total_M)
 export(peak_season_offset)
+export(run_metapop_simulation)
 export(run_simulation)
 export(run_simulation_with_repetitions)
 export(set_bednets)

diff --git a/R/biting_process.R b/R/biting_process.R
@@ -8,8 +8,13 @@
 #' @param variables a list of all of the model variables
 #' @param events a list of all of the model events
 #' @param parameters model pararmeters
-#' @param lagged_infectivity a LaggedValue class with historical sums of infectivity
+#' @param lagged_infectivity a list of LaggedValue objects with historical sums
+#' of infectivity, one for every metapopulation
 #' @param lagged_eir a LaggedValue class with historical EIRs
+#' @param mixing a vector of mixing coefficients for the lagged_infectivity
+#' values (default: 1)
+#' @param mixing_index an index for this population's position in the
+#' lagged_infectivity list (default: 1)
 #' @noRd
 create_biting_process <- function(
   renderer,
@@ -19,7 +24,9 @@ create_biting_process <- function(
   events,
   parameters,
   lagged_infectivity,
-  lagged_eir
+  lagged_eir,
+  mixing = 1,
+  mixing_index = 1
   ) {
   function(timestep) {
     # Calculate combined EIR
@@ -35,7 +42,9 @@ create_biting_process <- function(
       parameters,
       timestep,
       lagged_infectivity,
-      lagged_eir
+      lagged_eir,
+      mixing,
+      mixing_index
     )
 
     simulate_infection(
@@ -61,7 +70,9 @@ simulate_bites <- function(
   parameters,
   timestep,
   lagged_infectivity,
-  lagged_eir
+  lagged_eir,
+  mixing = 1,
+  mixing_index = 1
   ) {
   bitten_humans <- individual::Bitset$new(parameters$human_population)
 
@@ -117,8 +128,17 @@ simulate_bites <- function(
       n_infectious <- calculate_infectious_compartmental(solver_states)
     }
 
-    lagged_eir[[s_i]]$save(n_infectious * a, timestep)
-    species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
+    # store the current population's EIR for later
+    lagged_eir[[mixing_index]][[s_i]]$save(n_infectious * a, timestep)
+
+    # calculated the EIR for this timestep after mixing
+    species_eir <- sum(
+      vnapply(
+        lagged_eir,
+        function(l) l[[s_i]]$get(timestep - parameters$de)
+      ) * mixing
+    )
+
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     expected_bites <- species_eir * mean(psi)
@@ -131,9 +151,15 @@ simulate_bites <- function(
       }
     }
 
-    infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)
-    lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
-    foim <- calculate_foim(a, infectivity)
+    infectivity <- vnapply(
+      lagged_infectivity,
+      function(l) l$get(timestep - parameters$delay_gam)
+    )
+    lagged_infectivity[[mixing_index]]$save(
+      sum(human_infectivity * .pi),
+      timestep
+    )
+    foim <- calculate_foim(a, infectivity, mixing)
     renderer$render(paste0('FOIM_', species_name), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
     renderer$render(paste0('mu_', species_name), mu, timestep)
@@ -276,8 +302,10 @@ unique_biting_rate <- function(age, parameters) {
 #' @title Calculate the force of infection towards mosquitoes
 #'
 #' @param a human blood meal rate
-#' @param infectivity_sum a sum of infectivity weighted by relative biting rate
+#' @param infectivity_sum a vector of sums of infectivity weighted by relative
+#' biting rate for each population
+#' @param mixing a vector of mixing coefficients for each population
 #' @noRd
-calculate_foim <- function(a, infectivity_sum) {
-  a * infectivity_sum
+calculate_foim <- function(a, infectivity_sum, mixing) {
+  a * sum(infectivity_sum * mixing)
 }
diff --git a/R/disease_progression.R b/R/disease_progression.R
@@ -7,31 +7,50 @@
 #' @param infectivity the handle for the infectivity variable
 #' @param new_infectivity the new infectivity of the progressed individuals
 #' @noRd
-create_infection_update_listener <- function(
- state,
- to_state,
- infectivity,
- new_infectivity
- ) {
-  function(timestep, to_move) {
-    state$queue_update(to_state, to_move)
-    infectivity$queue_update(new_infectivity, to_move)
-  }
+update_infection <- function(
+  state,
+  to_state,
+  infectivity,
+  new_infectivity,
+  to_move
+  ) {
+  state$queue_update(to_state, to_move)
+  infectivity$queue_update(new_infectivity, to_move)
 }
 
-create_progression_process <- function(event, state, from_state, rate) {
+create_progression_process <- function(
+  state,
+  from_state,
+  to_state,
+  rate,
+  infectivity,
+  new_infectivity
+  ) {
   function(timestep) {
-    event$schedule(state$get_index_of(from_state)$sample(1/rate), 0)
+    to_move <- state$get_index_of(from_state)$sample(1/rate)
+    update_infection(
+      state,
+      to_state,
+      infectivity,
+      new_infectivity,
+      to_move
+    )
   }
 }
 
-create_rate_listener <- function(from_state, to_state, renderer) {
-  renderer$set_default(paste0('rate_', from_state, '_', to_state), 0)
-  function(timestep, target) {
-    renderer$render(
-      paste0('rate_', from_state, '_', to_state),
-      target$size(),
-      timestep
+create_asymptomatic_progression_process <- function(
+  state,
+  rate,
+  variables,
+  parameters
+  ) {
+  function(timestep) {
+    to_move <- state$get_index_of('D')$sample(1/rate)
+    update_to_asymptomatic_infection(
+      variables,
+      parameters,
+      timestep,
+      to_move
     )
   }
 }
@@ -43,22 +62,25 @@ create_rate_listener <- function(from_state, to_state, renderer) {
 #' @param variables the available human variables
 #' @param parameters model parameters
 #' @noRd
-create_asymptomatic_update_listener <- function(variables, parameters) {
-  function(timestep, to_move) {
-    if (to_move$size() > 0) {
-      variables$state$queue_update('A', to_move)
-      new_infectivity <- asymptomatic_infectivity(
-        get_age(
-          variables$birth$get_values(to_move),
-          timestep
-        ),
-        variables$id$get_values(to_move),
-        parameters
-      )
-      variables$infectivity$queue_update(
-        new_infectivity,
-        to_move
-      )
-    }
+update_to_asymptomatic_infection <- function(
+  variables,
+  parameters,
+  timestep,
+  to_move
+  ) {
+  if (to_move$size() > 0) {
+    variables$state$queue_update('A', to_move)
+    new_infectivity <- asymptomatic_infectivity(
+      get_age(
+        variables$birth$get_values(to_move),
+        timestep
+      ),
+      variables$id$get_values(to_move),
+      parameters
+    )
+    variables$infectivity$queue_update(
+      new_infectivity,
+      to_move
+    )
   }
 }
diff --git a/R/events.R b/R/events.R
@@ -1,14 +1,5 @@
 create_events <- function(parameters) {
   events <- list(
-    # Disease progression events
-    asymptomatic_progression = individual::TargetedEvent$new(parameters$human_population),
-    subpatent_progression = individual::TargetedEvent$new(parameters$human_population),
-    recovery = individual::TargetedEvent$new(parameters$human_population),
-
-    # Human infection events
-    clinical_infection = individual::TargetedEvent$new(parameters$human_population),
-    asymptomatic_infection = individual::TargetedEvent$new(parameters$human_population),
-
     # MDA events
     mda_administer = individual::Event$new(),
     smc_administer = individual::Event$new(),
@@ -100,65 +91,6 @@ attach_event_listeners <- function(
   correlations,
   renderer
   ) {
-
-  # =============
-  # State updates
-  # =============
-  # When infection events fire, update the corresponding states and infectivity
-  # variables
-
-  # Infection events
-  events$clinical_infection$add_listener(
-    create_infection_update_listener(
-      variables$state,
-      'D',
-      variables$infectivity,
-      parameters$cd
-    )
-  )
-
-  events$asymptomatic_progression$add_listener(
-    create_asymptomatic_update_listener(
-      variables,
-      parameters
-    )
-  )
-  events$asymptomatic_progression$add_listener(
-    create_rate_listener('D', 'A', renderer)
-  )
-  events$asymptomatic_infection$add_listener(
-    create_asymptomatic_update_listener(
-      variables,
-      parameters
-    )
-  )
-
-  # Recovery events
-  events$subpatent_progression$add_listener(
-    create_infection_update_listener(
-      variables$state,
-      'U',
-      variables$infectivity,
-      parameters$cu
-    )
-  )
-
-  events$subpatent_progression$add_listener(
-    create_rate_listener('A', 'U', renderer)
-  )
-
-  events$recovery$add_listener(
-    create_infection_update_listener(
-      variables$state,
-      'S',
-      variables$infectivity,
-      0
-    )
-  )
-  events$recovery$add_listener(
-    create_rate_listener('U', 'S', renderer)
-  )
-
   # ===========
   # Progression
   # ===========
@@ -226,8 +158,8 @@ attach_event_listeners <- function(
       parameters$mda_drug,
       parameters$mda_timesteps,
       parameters$mda_coverages,
-      parameters$mda_min_age,
-      parameters$mda_max_age,
+      parameters$mda_min_ages,
+      parameters$mda_max_ages,
       correlations,
       'mda',
       parameters,
@@ -242,8 +174,8 @@ attach_event_listeners <- function(
       parameters$smc_drug,
       parameters$smc_timesteps,
       parameters$smc_coverages,
-      parameters$smc_min_age,
-      parameters$smc_max_age,
+      parameters$smc_min_ages,
+      parameters$smc_max_ages,
       correlations,
       'smc',
       parameters,