From df56a3c919d99bdc80e831a01294ec09ca702cb3 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 16 Jun 2022 15:42:33 +0100
Subject: [PATCH 01/12] Implement metapopulation model:

 * e2e, integration and unit tests for metapopulation modelling
 * run_metapop_simulation function for setting up and integrating
   multiple models
 * update calculate_foim for multiple infectivity values
---
 R/biting_process.R                          | 39 +++++++++---
 R/model.R                                   | 70 ++++++++++++++++++++-
 R/processes.R                               | 40 ++++++++----
 R/variables.R                               |  5 ++
 tests/testthat/test-biting-integration.R    | 45 ++++++++++++-
 tests/testthat/test-biting.R                | 11 ++++
 tests/testthat/test-infection-integration.R |  3 -
 tests/testthat/test-simulation-e2e.R        | 25 ++++++++
 8 files changed, 211 insertions(+), 27 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index 97d801f0..291f982f 100644
--- 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)
 
@@ -131,9 +142,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 +293,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/model.R b/R/model.R
index 6c88e551..73de57b1 100644
--- a/R/model.R
+++ b/R/model.R
@@ -111,7 +111,8 @@ run_simulation <- function(
       parameters,
       vector_models,
       solvers,
-      correlations
+      correlations,
+      list(create_lagged_infectivity(variables, parameters))
     ),
     variables = variables,
     events = unlist(events),
@@ -120,6 +121,73 @@ run_simulation <- function(
   renderer$to_dataframe()
 }
 
+#' @title Run a metapopulation model
+#'
+#' @param timesteps the number of timesteps to run the simulation for (in days)
+#' @param parameters a named list of parameters to use
+#' @param correlations correlation parameters
+#' @param mixing matrix of mixing coefficients for infectivity
+#' @return dataframe of results
+#' @export
+run_metapop_simulation <- function(
+  timesteps,
+  parameters,
+  correlations = NULL,
+  mixing
+  ) {
+  random_seed(ceiling(runif(1) * .Machine$integer.max))
+  if (is.null(correlations)) {
+    correlations <- lapply(parameters, get_correlation_parameters)
+  }
+  variables <- lapply(parameters, create_variables)
+  events <- lapply(parameters, create_events)
+  renderer <- lapply(parameters, function(.) individual::Render$new(timesteps))
+  for (i in seq_along(parameters)) {
+    initialise_events(events[[i]], variables[[i]], parameters[[i]])
+    attach_event_listeners(
+      events[[i]],
+      variables[[i]],
+      parameters[[i]],
+      correlations[[i]],
+      renderer[[i]]
+    )
+  }
+  vector_models <- lapply(parameters, parameterise_mosquito_models)
+  solvers <- lapply(
+    seq_along(parameters),
+    function(i) parameterise_solvers(vector_models[[i]], parameters[[i]])
+  )
+  lagged_infectivity <- lapply(
+    seq_along(parameters),
+    function(i) create_lagged_infectivity(variables[[i]], parameters[[i]])
+  )
+  processes <- lapply(
+    seq_along(parameters),
+    function(i) {
+      create_processes(
+        renderer[[i]],
+        variables[[i]],
+        events[[i]],
+        parameters[[i]],
+        vector_models[[i]],
+        solvers[[i]],
+        correlations[[i]],
+        lagged_infectivity,
+        mixing[i,],
+        i
+      )
+    }
+  )
+  individual::simulation_loop(
+    processes = unlist(processes),
+    variables = unlist(variables),
+    events = unlist(events),
+    timesteps = timesteps
+  )
+  
+  lapply(renderer, function(r) r$to_dataframe())
+}
+
 #' @title Run the simulation with repetitions
 #'
 #' @param timesteps the number of timesteps to run the simulation for
diff --git a/R/processes.R b/R/processes.R
index f5b4b6ee..a22a8ee3 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -10,6 +10,12 @@
 #' @param models a list of vector models, one for each species
 #' @param solvers a list of ode solvers, one for each species
 #' @param correlations the intervention correlations object
+#' @param lagged_infectivity a list of LaggedValue objects for each population
+#' in the simulation
+#' @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_processes <- function(
   renderer,
@@ -18,7 +24,10 @@ create_processes <- function(
   parameters,
   models,
   solvers,
-  correlations
+  correlations,
+  lagged_infectivity,
+  mixing = 1,
+  mixing_index = 1
   ) {
   # ========
   # Immunity
@@ -70,15 +79,6 @@ create_processes <- function(
     }
   )
 
-  age <- get_age(variables$birth$get_values(), 0)
-  psi <- unique_biting_rate(age, parameters)
-  .pi <- human_pi(psi, variables$zeta$get_values())
-  init_infectivity <- sum(.pi * variables$infectivity$get_values())
-  lagged_infectivity <- LaggedValue$new(
-    max_lag = parameters$delay_gam + 2,
-    default = init_infectivity
-  )
-
   processes <- c(
     processes,
     create_biting_process(
@@ -89,7 +89,9 @@ create_processes <- function(
       events,
       parameters,
       lagged_infectivity,
-      lagged_eir
+      lagged_eir,
+      mixing,
+      mixing_index
     ),
     create_mortality_process(variables, events, renderer, parameters),
     create_progression_process(
@@ -230,3 +232,19 @@ create_exponential_decay_process <- function(variable, rate) {
   decay_rate <- exp(-1/rate)
   function(timestep) variable$queue_update(variable$get_values() * decay_rate)
 }
+
+#' @title Create and initialise lagged_infectivity object
+#'
+#' @param variables model variables for initialisation
+#' @param parameters model parameters 
+#' @noRd
+create_lagged_infectivity <- function(variables, parameters) {
+  age <- get_age(variables$birth$get_values(), 0)
+  psi <- unique_biting_rate(age, parameters)
+  .pi <- human_pi(psi, variables$zeta$get_values())
+  init_infectivity <- sum(.pi * variables$infectivity$get_values())
+  LaggedValue$new(
+    max_lag = parameters$delay_gam + 2,
+    default = init_infectivity
+  )
+}
diff --git a/R/variables.R b/R/variables.R
index 09d4f56e..95acb018 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -307,6 +307,11 @@ create_variables <- function(parameters) {
   variables
 }
 
+
+create_export_variable <- function(metapop_params) {
+  individual::DoubleVariable$new(rep(0, length(metapop_params$x)))
+}
+
 # =========
 # Utilities
 # =========
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 2b82a1ba..30541c7e 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -45,7 +45,9 @@ test_that('biting_process integrates mosquito effects and human infection', {
     parameters,
     timestep,
     lagged_foim,
-    lagged_eir
+    lagged_eir,
+    1,
+    1
   )
 
   mockery::expect_args(
@@ -98,7 +100,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
   mockery::stub(simulate_bites, 'aquatic_mosquito_model_update', eqs_update)
   models <- parameterise_mosquito_models(parameters)
   solvers <- parameterise_solvers(models, parameters)
-  lagged_foim <- LaggedValue$new(12.5, .001)
+  lagged_foim <- list(LaggedValue$new(12.5, .001))
   lagged_eir <- list(LaggedValue$new(12, 10))
   bitten <- simulate_bites(
     renderer,
@@ -142,3 +144,42 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     c(.5, .5, .5, .5)
   )
 })
+
+
+test_that('simulate_bites works with mixed populations', {
+  population <- 4
+  timestep <- 5
+  renderer <- individual::Render$new(5)
+  parameters <- get_parameters(
+    list(human_population = population)
+  )
+  events <- create_events(parameters)
+  variables <- create_variables(parameters)
+
+  mock_foim <- mockery::mock(1)
+  mock_a <- mockery::mock(.3)
+
+  mockery::stub(simulate_bites, 'calculate_foim', mock_foim)
+  mockery::stub(simulate_bites, '.human_blood_meal_rate', mock_a)
+  models <- parameterise_mosquito_models(parameters)
+  solvers <- parameterise_solvers(models, parameters)
+  lagged_foim <- list(LaggedValue$new(12.5, .001), LaggedValue$new(12.5, .01))
+  lagged_eir <- list(LaggedValue$new(12, 10))
+  age <- c(20, 24, 5, 39) * 365
+  bitten <- simulate_bites(
+    renderer,
+    solvers,
+    models,
+    variables,
+    events,
+    age,
+    parameters,
+    timestep,
+    lagged_foim,
+    lagged_eir,
+    c(0.2, 0.8),
+    2
+  )
+
+  mockery::expect_args(mock_foim, 1, .3, c(.001, .01), c(.2, .8))
+})
diff --git a/tests/testthat/test-biting.R b/tests/testthat/test-biting.R
index e6c19fb9..91a60be4 100644
--- a/tests/testthat/test-biting.R
+++ b/tests/testthat/test-biting.R
@@ -10,3 +10,14 @@ test_that('gonotrophic_cycle cannot be negative', {
   vparams$blood_meal_rates <- 5
   expect_error(set_species(params, list(vparams), 1))
 })
+
+test_that('calculate_foim mixes correctly', {
+  expect_equal( # check 1 case
+    calculate_foim(a=0.5, infectivity=6, mixing=1),
+    3
+  )
+  expect_equal( # check 2 cases
+    calculate_foim(a=0.5, infectivity=c(6, 4), mixing=c(.5, .5)),
+    2.5
+  )
+})
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 819c5a69..eea1e4b8 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -16,9 +16,6 @@ test_that('simulate_infection integrates different types of infection and schedu
     state = list(get_index_of = mockery::mock(asymptomatics))
   )
 
-  total_eir <- 5
-  eir <- rep(total_eir / population, population)
-
   bitten <- individual::Bitset$new(population)$insert(c(1, 3, 5, 7))
   boost_immunity_mock <- mockery::mock()
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
diff --git a/tests/testthat/test-simulation-e2e.R b/tests/testthat/test-simulation-e2e.R
index f5f1b8ab..0e7f7a6c 100644
--- a/tests/testthat/test-simulation-e2e.R
+++ b/tests/testthat/test-simulation-e2e.R
@@ -2,3 +2,28 @@ test_that('Simulation runs for a few timesteps', {
   sim <- run_simulation(100)
   expect_equal(nrow(sim), 100)
 })
+
+test_that('run_metapop_simulation fails with incorrect mixing matrix', {
+  population <- 4
+  timestep <- 5
+  parameters <- get_parameters(list(human_population = population))
+  paramlist <- list(parameters, parameters)
+  # incorrect params
+  mixing <- matrix(c(1, 1), nrow = 1, ncol = 2)
+  expect_error(
+    run_metapop_simulation(timesteps, parameters, NULL, mixing)
+  )
+})
+
+test_that('run_metapop_simulation integrates two models correctly', {
+  population <- 4
+  timesteps <- 5
+  parameters <- get_parameters(list(human_population = population))
+  parametersets <- list(parameters, parameters)
+  mixing <- diag(1, nrow = 2, ncol = 2)
+
+  outputs <- run_metapop_simulation(timesteps, parametersets, NULL, mixing)
+  expect_equal(length(outputs), 2)
+  expect_equal(nrow(outputs[[1]]), 5)
+  expect_equal(nrow(outputs[[2]]), 5)
+})

From 5dd0d033df642db2f95e88e5238f9012111be35c Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 16 Jun 2022 15:58:05 +0100
Subject: [PATCH 02/12] Add documentation and input validation for metapop

---
 NAMESPACE                     |  1 +
 R/model.R                     | 16 ++++++++++++----
 man/run_metapop_simulation.Rd | 25 +++++++++++++++++++++++++
 man/set_equilibrium.Rd        |  3 ++-
 4 files changed, 40 insertions(+), 5 deletions(-)
 create mode 100644 man/run_metapop_simulation.Rd

diff --git a/NAMESPACE b/NAMESPACE
index e8f26fc0..c0d7331a 100644
--- 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/model.R b/R/model.R
index 73de57b1..17064790 100644
--- a/R/model.R
+++ b/R/model.R
@@ -124,10 +124,12 @@ run_simulation <- function(
 #' @title Run a metapopulation model
 #'
 #' @param timesteps the number of timesteps to run the simulation for (in days)
-#' @param parameters a named list of parameters to use
-#' @param correlations correlation parameters
-#' @param mixing matrix of mixing coefficients for infectivity
-#' @return dataframe of results
+#' @param parameters a list of model parameter lists for each population
+#' @param correlations a list of correlation parameters for each population
+#' (default: NULL)
+#' @param mixing matrix of mixing coefficients for infectivity towards
+#' mosquitoes
+#' @return a list of dataframe of results
 #' @export
 run_metapop_simulation <- function(
   timesteps,
@@ -136,6 +138,12 @@ run_metapop_simulation <- function(
   mixing
   ) {
   random_seed(ceiling(runif(1) * .Machine$integer.max))
+  if (nrow(mixing) != ncol(mixing)) {
+    stop('mixing matrix must be square')
+  }
+  if (nrow(mixing) != length(parameters)) {
+    stop('mixing matrix rows must match length of parameters')
+  }
   if (is.null(correlations)) {
     correlations <- lapply(parameters, get_correlation_parameters)
   }
diff --git a/man/run_metapop_simulation.Rd b/man/run_metapop_simulation.Rd
new file mode 100644
index 00000000..83b83968
--- /dev/null
+++ b/man/run_metapop_simulation.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/model.R
+\name{run_metapop_simulation}
+\alias{run_metapop_simulation}
+\title{Run a metapopulation model}
+\usage{
+run_metapop_simulation(timesteps, parameters, correlations = NULL, mixing)
+}
+\arguments{
+\item{timesteps}{the number of timesteps to run the simulation for (in days)}
+
+\item{parameters}{a list of model parameter lists for each population}
+
+\item{correlations}{a list of correlation parameters for each population
+(default: NULL)}
+
+\item{mixing}{matrix of mixing coefficients for infectivity towards
+mosquitoes}
+}
+\value{
+a list of dataframe of results
+}
+\description{
+Run a metapopulation model
+}
diff --git a/man/set_equilibrium.Rd b/man/set_equilibrium.Rd
index f41905bd..f1c83c54 100644
--- a/man/set_equilibrium.Rd
+++ b/man/set_equilibrium.Rd
@@ -9,7 +9,8 @@ set_equilibrium(parameters, init_EIR, eq_params = NULL)
 \arguments{
 \item{parameters}{model parameters to update}
 
-\item{init_EIR}{the desired initial EIR}
+\item{init_EIR}{the desired initial EIR (infectious bites per day over the entire human
+population)}
 
 \item{eq_params}{parameters from the malariaEquilibrium package, if null.
 The default malariaEquilibrium parameters will be used}

From bd479fad821050c47b70c7682e83d142249129ba Mon Sep 17 00:00:00 2001
From: Hillary Topazian <80535782+htopazian@users.noreply.github.com>
Date: Mon, 20 Jun 2022 11:31:22 +0100
Subject: [PATCH 03/12] Update DESCRIPTION

---
 DESCRIPTION | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 49cc359f..d7385c26 100644
--- 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@imperial.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",

From e9193c523be2614e5cf03f0e5defdbecccb9f40f Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Tue, 28 Jun 2022 13:42:06 +0100
Subject: [PATCH 04/12] Fix disease progression bug:

 * remove overcomplicated events from disease progression
 * update regression tests
---
 R/disease_progression.R                     | 94 +++++++++++++--------
 R/events.R                                  | 68 ---------------
 R/human_infection.R                         | 29 ++++---
 R/mortality_processes.R                     |  6 --
 R/processes.R                               | 26 +++---
 tests/testthat/test-infection-integration.R | 52 +++++-------
 tests/testthat/test-progression.R           | 13 +--
 7 files changed, 122 insertions(+), 166 deletions(-)

diff --git a/R/disease_progression.R b/R/disease_progression.R
index e7d8a654..05367252 100644
--- 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
index e8e5704c..d3425097 100644
--- 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
   # ===========
diff --git a/R/human_infection.R b/R/human_infection.R
index bd76ac97..a7474382 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -73,7 +73,6 @@ simulate_infection <- function(
   treated <- calculate_treated(
     variables,
     clinical_infections,
-    events$recovery,
     parameters,
     timestep,
     renderer
@@ -83,11 +82,12 @@ simulate_infection <- function(
   renderer$render('n_infections', infected_humans$size(), timestep)
 
   schedule_infections(
-    events,
+    variables,
     clinical_infections,
     treated,
     infected_humans,
-    parameters
+    parameters,
+    timestep
   )
 }
 
@@ -255,7 +255,6 @@ update_severe_disease <- function(
 #' Sample treated humans from the clinically infected
 #' @param variables a list of all of the model variables
 #' @param clinical_infections a bitset of clinically infected humans
-#' @param recovery the recovery event
 #' @param parameters model parameters
 #' @param timestep the current timestep
 #' @param renderer simulation renderer
@@ -263,7 +262,6 @@ update_severe_disease <- function(
 calculate_treated <- function(
   variables,
   clinical_infections,
-  recovery,
   parameters,
   timestep,
   renderer
@@ -319,12 +317,12 @@ calculate_treated <- function(
 #' @param parameters model parameters
 #' @noRd
 schedule_infections <- function(
-  events,
+  variables,
   clinical_infections,
   treated,
   infections,
   parameters,
-  asymptomatics
+  timestep
   ) {
   included <- treated$not(TRUE)
 
@@ -334,13 +332,22 @@ schedule_infections <- function(
   )
 
   if(to_infect$size() > 0) {
-    infection_times <- log_uniform(to_infect$size(), parameters$de)
-    events$clinical_infection$schedule(to_infect, 0)
+    update_infection(
+      variables$state,
+      'D',
+      variables$infectivity,
+      parameters$cd,
+      to_infect
+    )
   }
 
   if(to_infect_asym$size() > 0) {
-    infection_times <- log_uniform(to_infect_asym$size(), parameters$de)
-    events$asymptomatic_infection$schedule(to_infect_asym, 0)
+    update_to_asymptomatic_infection(
+      variables,
+      parameters,
+      timestep,
+      to_infect_asym
+    )
   }
 }
 
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index 4abc2053..5884c928 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -78,12 +78,6 @@ reset_target <- function(variables, events, target, state, timestep) {
   if (target$size() > 0) {
     # clear events
     to_clear <- c(
-      'asymptomatic_progression',
-      'subpatent_progression',
-      'recovery',
-      'clinical_infection',
-      'asymptomatic_infection',
-      'detection',
       'throw_away_net',
       'rtss_mass_doses',
       'rtss_mass_booster',
diff --git a/R/processes.R b/R/processes.R
index a22a8ee3..d8d52e17 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -94,29 +94,35 @@ create_processes <- function(
       mixing_index
     ),
     create_mortality_process(variables, events, renderer, parameters),
-    create_progression_process(
-      events$asymptomatic_progression,
+    create_asymptomatic_progression_process(
       variables$state,
-      'D',
-      parameters$dd
+      parameters$dd,
+      variables,
+      parameters
     ),
     create_progression_process(
-      events$subpatent_progression,
       variables$state,
       'A',
-      parameters$da
+      'U',
+      parameters$da,
+      variables$infectivity,
+      parameters$cu
     ),
     create_progression_process(
-      events$recovery,
       variables$state,
       'U',
-      parameters$du
+      'S',
+      parameters$du,
+      variables$infectivity,
+      0
     ),
     create_progression_process(
-      events$recovery,
       variables$state,
       'Tr',
-      parameters$dt
+      'S',
+      parameters$dt,
+      variables$infectivity,
+      0
     )
   )
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index eea1e4b8..9ff16b3a 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -91,7 +91,6 @@ test_that('simulate_infection integrates different types of infection and schedu
     1,
     variables,
     clinical,
-    events$recovery,
     parameters,
     timestep,
     renderer
@@ -100,11 +99,12 @@ test_that('simulate_infection integrates different types of infection and schedu
   mockery::expect_args(
     schedule_mock,
     1,
-    events,
+    variables,
     clinical,
     treated,
     infected,
-    parameters
+    parameters,
+    timestep
   )
 })
 
@@ -265,7 +265,6 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   calculate_treated(
     variables,
     clinical_infections,
-    events$recovery,
     parameters,
     timestep,
     mock_render(timestep)
@@ -297,44 +296,39 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
 })
 
 test_that('schedule_infections correctly schedules new infections', {
-  parameters <- get_parameters()
-
-  clinical_mock <- mockery::mock()
-  asym_mock <- mockery::mock()
-
-  events <- list(
-    clinical_infection = list(
-      schedule = clinical_mock
-    ),
-    asymptomatic_infection = list(
-      schedule = asym_mock
-    ),
-    subpatent_infection = list(clear_schedule = mockery::mock()),
-    recovery = list(clear_schedule = mockery::mock())
-  )
+  parameters <- get_parameters(list(human_population = 20))
+  variables <- create_variables(parameters)
 
   infections <- individual::Bitset$new(20)$insert(1:20)
   clinical_infections <- individual::Bitset$new(20)$insert(5:15)
   treated <- individual::Bitset$new(20)$insert(7:12)
 
+  infection_mock <- mockery::mock()
+  asymp_mock <- mockery::mock()
+
+  mockery::stub(schedule_infections, 'update_infection', infection_mock)
+  mockery::stub(schedule_infections, 'update_to_asymptomatic_infection', asymp_mock)
+
   schedule_infections(
-    events,
+    variables,
     clinical_infections,
     treated,
     infections,
-    parameters
+    parameters,
+    42 
   )
 
-  expect_bitset_schedule(
-    clinical_mock,
-    c(5, 6, 13, 14, 15),
-    0
+  actual_infected <- mockery::mock_args(infection_mock)[[1]][[5]]$to_vector()
+  actual_asymp_infected <- mockery::mock_args(asymp_mock)[[1]][[4]]$to_vector()
+
+  expect_equal(
+    actual_infected,
+    c(5, 6, 13, 14, 15)
   )
 
-  expect_bitset_schedule(
-    asym_mock,
-    c(1, 2, 3, 4, 16, 17, 18, 19, 20),
-    0,
+  expect_equal(
+    actual_asymp_infected,
+    c(1, 2, 3, 4, 16, 17, 18, 19, 20)
   )
 })
 
diff --git a/tests/testthat/test-progression.R b/tests/testthat/test-progression.R
index d4ee3459..a8b7415c 100644
--- a/tests/testthat/test-progression.R
+++ b/tests/testthat/test-progression.R
@@ -2,16 +2,17 @@
 test_that('Asymptomatic infection listener updates infectivity correctly', {
   parameters <- get_parameters(list(human_population = 4))
   variables <- create_variables(parameters)
-  listener <- create_asymptomatic_update_listener(variables, parameters)
-
-  mockery::stub(listener, 'asymptomatic_infectivity', mockery::mock(c(.2, .5)))
 
   update_mock <- mockery::mock()
   variables$infectivity <- list(queue_update = update_mock)
 
-  to_move <- individual::Bitset$new(4)
-  to_move$insert(c(2, 4))
-  listener(1, to_move)
+  to_move <- individual::Bitset$new(4)$insert(c(2, 4))
+  mockery::stub(
+    update_to_asymptomatic_infection,
+    'asymptomatic_infectivity',
+    mockery::mock(c(.2, .5))
+  )
+  update_to_asymptomatic_infection(variables, parameters, 1, to_move)
 
   expect_bitset_update(
     update_mock,

From 5c5edcc1e7f0ced08685be51010146d99517d05f Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Wed, 13 Jul 2022 17:45:04 +0100
Subject: [PATCH 05/12] Fix events bug with metapopulation model:

 * Add diagnostic rendering and documentation for net_usage
 * Make run_metapop_simulation forceAndCall when initialising and attaching event
   listeners
---
 .github/workflows/R-CMD-check.yaml |  2 +-
 R/model.R                          | 15 +++++++++++++--
 R/processes.R                      |  3 ++-
 R/vector_control.R                 | 10 ++++++++++
 4 files changed, 26 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index 4b6a405e..b0d19a8d 100644
--- 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/R/model.R b/R/model.R
index 17064790..9c9862c0 100644
--- a/R/model.R
+++ b/R/model.R
@@ -71,6 +71,7 @@
 #' subpatent
 #'  * rate_U_S: rate that humans transition from subpatent to
 #' susceptible
+#'  * net_usage: the proportion of the population protected by a bed net
 #'  * mosquito_deaths: number of adult female mosquitoes who die this timestep
 #'
 #' @param timesteps the number of timesteps to run the simulation for (in days)
@@ -151,8 +152,18 @@ run_metapop_simulation <- function(
   events <- lapply(parameters, create_events)
   renderer <- lapply(parameters, function(.) individual::Render$new(timesteps))
   for (i in seq_along(parameters)) {
-    initialise_events(events[[i]], variables[[i]], parameters[[i]])
-    attach_event_listeners(
+    # NOTE: forceAndCall is necessary here to make sure i refers to the current
+    # iteration
+    forceAndCall(
+      3,
+      initialise_events,
+      events[[i]],
+      variables[[i]],
+      parameters[[i]]
+    )
+    forceAndCall(
+      5,
+      attach_event_listeners,
       events[[i]],
       variables[[i]],
       parameters[[i]],
diff --git a/R/processes.R b/R/processes.R
index d8d52e17..2f55b914 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -208,7 +208,8 @@ create_processes <- function(
         events$throw_away_net,
         parameters,
         correlations
-      )
+      ),
+      net_usage_renderer(variables$net_time, renderer)
     )
   }
 
diff --git a/R/vector_control.R b/R/vector_control.R
index f197a7d1..0537d746 100644
--- a/R/vector_control.R
+++ b/R/vector_control.R
@@ -185,3 +185,13 @@ bednet_decay <- function(t, gamma) {
 spraying_decay <- function(t, theta, gamma) {
   1 / (1 + exp(-(theta + gamma * t)))
 }
+
+net_usage_renderer <- function(net_time, renderer) {
+  function(t) {
+    renderer$render(
+      'net_usage',
+      net_time$get_index_of(-1)$not(TRUE)$size(),
+      t
+    )
+  }
+}

From 6ff7feb887c8bba668e62d69d5319ed887517c6b Mon Sep 17 00:00:00 2001
From: Peter Winskill <p.winskill@imperial.ac.uk>
Date: Thu, 21 Jul 2022 10:14:44 +0100
Subject: [PATCH 06/12] Name typo correction

---
 DESCRIPTION | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 49cc359f..02dffb8d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -10,9 +10,9 @@ 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",
@@ -66,7 +66,7 @@ Suggests:
   ggplot2,
   covr,
   mgcv 
-RoxygenNote: 7.1.2
+RoxygenNote: 7.2.1
 Roxygen: list(markdown = TRUE)
 LinkingTo:
   Rcpp,

From eb4fcec4328ca1e4287caf67b81cb3239bb933db Mon Sep 17 00:00:00 2001
From: Peter Winskill <p.winskill@imperial.ac.uk>
Date: Thu, 21 Jul 2022 11:04:12 +0100
Subject: [PATCH 07/12] Bump version

---
 DESCRIPTION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 02dffb8d..1fb1c563 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: malariasimulation
 Title: An individual based model for malaria
-Version: 1.3.0
+Version: 1.4.0
 Authors@R: c(
   person(
     given = "Giovanni",

From 7a842f54c6060c6b5af127d6403b24caaa70972a Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 21 Jul 2022 12:00:15 +0100
Subject: [PATCH 08/12] Mix EIR betweeen metapopulations

---
 R/biting_process.R                       | 13 ++++++-
 R/model.R                                |  8 +++-
 R/processes.R                            | 49 +++++++++++++++---------
 man/run_simulation.Rd                    |  1 +
 tests/testthat/test-biting-integration.R |  4 +-
 5 files changed, 51 insertions(+), 24 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index 291f982f..9dbb2d41 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -128,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)
diff --git a/R/model.R b/R/model.R
index 9c9862c0..31b84d3e 100644
--- a/R/model.R
+++ b/R/model.R
@@ -71,7 +71,7 @@
 #' subpatent
 #'  * rate_U_S: rate that humans transition from subpatent to
 #' susceptible
-#'  * net_usage: the proportion of the population protected by a bed net
+#'  * net_usage: the number people protected by a bed net
 #'  * mosquito_deaths: number of adult female mosquitoes who die this timestep
 #'
 #' @param timesteps the number of timesteps to run the simulation for (in days)
@@ -113,6 +113,7 @@ run_simulation <- function(
       vector_models,
       solvers,
       correlations,
+      list(create_lagged_eir(variables, solvers, parameters)),
       list(create_lagged_infectivity(variables, parameters))
     ),
     variables = variables,
@@ -176,6 +177,10 @@ run_metapop_simulation <- function(
     seq_along(parameters),
     function(i) parameterise_solvers(vector_models[[i]], parameters[[i]])
   )
+  lagged_eir <- lapply(
+    seq_along(parameters),
+    function(i) create_lagged_eir(variables[[i]], solvers[[i]], parameters[[i]])
+  )
   lagged_infectivity <- lapply(
     seq_along(parameters),
     function(i) create_lagged_infectivity(variables[[i]], parameters[[i]])
@@ -191,6 +196,7 @@ run_metapop_simulation <- function(
         vector_models[[i]],
         solvers[[i]],
         correlations[[i]],
+        lagged_eir,
         lagged_infectivity,
         mixing[i,],
         i
diff --git a/R/processes.R b/R/processes.R
index 2f55b914..beedff6c 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -10,12 +10,14 @@
 #' @param models a list of vector models, one for each species
 #' @param solvers a list of ode solvers, one for each species
 #' @param correlations the intervention correlations object
-#' @param lagged_infectivity a list of LaggedValue objects for each population
+#' @param lagged_eir a list of list of LaggedValue objects for EIR for each
+#' population and species in the simulation
+#' @param lagged_infectivity a list of LaggedValue objects for FOIM for each population
 #' in the simulation
-#' @param mixing a vector of mixing coefficients for the lagged_infectivity
+#' @param mixing a vector of mixing coefficients for the lagged transmission
 #' values (default: 1)
 #' @param mixing_index an index for this population's position in the
-#' lagged_infectivity list (default: 1)
+#' lagged transmission lists (default: 1)
 #' @noRd
 create_processes <- function(
   renderer,
@@ -25,6 +27,7 @@ create_processes <- function(
   models,
   solvers,
   correlations,
+  lagged_eir,
   lagged_infectivity,
   mixing = 1,
   mixing_index = 1
@@ -63,22 +66,6 @@ create_processes <- function(
   # schedule infections for humans and set last_boosted_*
   # move mosquitoes into incubating state
   # kill mosquitoes caught in vector control
-  lagged_eir <- lapply(
-    seq_along(parameters$species),
-    function(species) {
-      LaggedValue$new(
-        max_lag = parameters$de + 1,
-        default = calculate_eir(
-          species,
-          solvers,
-          variables,
-          parameters,
-          0
-        )
-      )
-    }
-  )
-
   processes <- c(
     processes,
     create_biting_process(
@@ -255,3 +242,27 @@ create_lagged_infectivity <- function(variables, parameters) {
     default = init_infectivity
   )
 }
+
+#' @title Create and initialise a list of lagged_eir objects per species
+#'
+#' @param variables model variables for initialisation
+#' @param solvers model differential equation solvers
+#' @param parameters model parameters 
+#' @noRd
+create_lagged_eir <- function(variables, solvers, parameters) {
+  lapply(
+    seq_along(parameters$species),
+    function(species) {
+      LaggedValue$new(
+        max_lag = parameters$de + 1,
+        default = calculate_eir(
+          species,
+          solvers,
+          variables,
+          parameters,
+          0
+        )
+      )
+    }
+  )
+}
diff --git a/man/run_simulation.Rd b/man/run_simulation.Rd
index 8c371f8e..7e4f795c 100644
--- a/man/run_simulation.Rd
+++ b/man/run_simulation.Rd
@@ -88,6 +88,7 @@ asymptomatic
 subpatent
 \item rate_U_S: rate that humans transition from subpatent to
 susceptible
+\item net_usage: the proportion of the population protected by a bed net
 \item mosquito_deaths: number of adult female mosquitoes who die this timestep
 }
 }
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 30541c7e..26e2da6f 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -101,7 +101,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
   models <- parameterise_mosquito_models(parameters)
   solvers <- parameterise_solvers(models, parameters)
   lagged_foim <- list(LaggedValue$new(12.5, .001))
-  lagged_eir <- list(LaggedValue$new(12, 10))
+  lagged_eir <- list(list(LaggedValue$new(12, 10)))
   bitten <- simulate_bites(
     renderer,
     solvers,
@@ -164,7 +164,7 @@ test_that('simulate_bites works with mixed populations', {
   models <- parameterise_mosquito_models(parameters)
   solvers <- parameterise_solvers(models, parameters)
   lagged_foim <- list(LaggedValue$new(12.5, .001), LaggedValue$new(12.5, .01))
-  lagged_eir <- list(LaggedValue$new(12, 10))
+  lagged_eir <- list(list(LaggedValue$new(12, 10)), list(LaggedValue$new(12, 10)))
   age <- c(20, 24, 5, 39) * 365
   bitten <- simulate_bites(
     renderer,

From 788a2c70b2b3950d5757893f5d87f20ffe821d7e Mon Sep 17 00:00:00 2001
From: Peter Winskill <p.winskill@imperial.ac.uk>
Date: Thu, 21 Jul 2022 17:27:03 +0100
Subject: [PATCH 09/12] Time-varying MDA age range

---
 R/events.R                                  |   8 +-
 R/mda_parameters.R                          |  50 +++++--
 R/mda_processes.R                           |  10 +-
 R/parameters.R                              |   8 +-
 man/CorrelationParameters.Rd                |  36 ++---
 man/set_equilibrium.Rd                      |   3 +-
 man/set_mda.Rd                              |   8 +-
 man/set_smc.Rd                              |   8 +-
 tests/testthat/test-mda.R                   | 144 ++++++++++++++++++--
 tests/testthat/test-mortality-integration.R |   4 +-
 vignettes/MDA.Rmd                           |  10 +-
 11 files changed, 214 insertions(+), 75 deletions(-)

diff --git a/R/events.R b/R/events.R
index e8e5704c..dbcdc8a3 100644
--- a/R/events.R
+++ b/R/events.R
@@ -226,8 +226,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 +242,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,
diff --git a/R/mda_parameters.R b/R/mda_parameters.R
index 166f8a1c..287ab61b 100644
--- a/R/mda_parameters.R
+++ b/R/mda_parameters.R
@@ -2,10 +2,10 @@
 #' @param parameters a list of parameters to modify
 #' @param drug the index of the drug to administer
 #' @param timesteps a vector of timesteps for each round of mda
-#' @param coverages the proportion of the target population who recieve each
+#' @param coverages a vector of the proportion of the target population who receive each
 #' round
-#' @param min_age the minimum age of the target population exclusive (in timesteps)
-#' @param max_age the maximum age of the target population exclusive (in timesteps)
+#' @param min_ages a vector of minimum ages of the target population for each round exclusive (in timesteps)
+#' @param max_ages a vector of maximum ages of the target population for each round exclusive (in timesteps)
 #' drug
 #' @export
 set_mda <- function(
@@ -13,15 +13,26 @@ set_mda <- function(
   drug,
   timesteps,
   coverages,
-  min_age,
-  max_age
+  min_ages,
+  max_ages
   ) {
+  
+  if(length(coverages) != length(timesteps)){
+    stop("coverages and timesteps do no align")
+  }
+  if(length(min_ages) != length(timesteps)){
+    stop("minimum ages and timesteps do no align")
+  }
+  if(length(max_ages) != length(timesteps)){
+    stop("maximum ages and timesteps do no align")
+  }
+  
   parameters$mda <- TRUE
   parameters$mda_drug <- drug
   parameters$mda_timesteps <- timesteps
   parameters$mda_coverages <- coverages
-  parameters$mda_min_age <- min_age
-  parameters$mda_max_age <- max_age
+  parameters$mda_min_ages <- min_ages
+  parameters$mda_max_ages <- max_ages
   parameters
 }
 
@@ -29,10 +40,10 @@ set_mda <- function(
 #' @param parameters a list of parameters to modify
 #' @param drug the index of the drug to administer
 #' @param timesteps a vector of timesteps for each round of smc
-#' @param coverages the proportion of the target population who recieve each
+#' @param coverages a vector of the proportion of the target population who receive each
 #' round
-#' @param min_age the minimum age of the target population exclusive (in timesteps)
-#' @param max_age the maximum age of the target population exclusive (in timesteps)
+#' @param min_ages a vector of minimum ages of the target population for each round exclusive (in timesteps)
+#' @param max_ages a vector of maximum ages of the target population for each round exclusive (in timesteps)
 #' drug
 #' @export
 set_smc <- function(
@@ -40,14 +51,25 @@ set_smc <- function(
   drug,
   timesteps,
   coverages,
-  min_age,
-  max_age
+  min_ages,
+  max_ages
   ) {
+  
+  if(length(coverages) != length(timesteps)){
+    stop("coverages and timesteps do no align")
+  }
+  if(length(min_ages) != length(timesteps)){
+    stop("minimum ages and timesteps do no align")
+  }
+  if(length(max_ages) != length(timesteps)){
+    stop("maximum ages and timesteps do no align")
+  }
+  
   parameters$smc <- TRUE
   parameters$smc_drug <- drug
   parameters$smc_timesteps <- timesteps
   parameters$smc_coverages <- coverages
-  parameters$smc_min_age <- min_age
-  parameters$smc_max_age <- max_age
+  parameters$smc_min_ages <- min_ages
+  parameters$smc_max_ages <- max_ages
   parameters
 }
diff --git a/R/mda_processes.R b/R/mda_processes.R
index 3833a152..348e224b 100644
--- a/R/mda_processes.R
+++ b/R/mda_processes.R
@@ -4,8 +4,8 @@
 #' @param drug the drug to administer
 #' @param timesteps timesteps for each round
 #' @param coverages the coverage for each round
-#' @param min_age minimum age for the target population
-#' @param max_age maximum age for the target population
+#' @param min_ages minimum age for the target population for each round
+#' @param max_ages maximum age for the target population for each round
 #' @param correlations correlation parameters
 #' @param int_name the name of this intervention (either 'smc' or 'mda')
 #' @param parameters the model parameters
@@ -18,8 +18,8 @@ create_mda_listeners <- function(
   drug,
   timesteps,
   coverages,
-  min_age,
-  max_age,
+  min_ages,
+  max_ages,
   correlations,
   int_name,
   parameters,
@@ -30,7 +30,7 @@ create_mda_listeners <- function(
     coverage <- coverages[[time_index]]
     age <- get_age(variables$birth$get_values(), timestep)
 
-    in_age <- which((age > min_age) & (age < max_age))
+    in_age <- which((age > min_ages[[time_index]]) & (age < max_ages[[time_index]]))
     target <- in_age[sample_intervention(in_age, int_name, coverage, correlations)]
 
     successful_treatments <- bernoulli(
diff --git a/R/parameters.R b/R/parameters.R
index 44b4a512..309b0265 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -364,14 +364,14 @@ get_parameters <- function(overrides = list()) {
     mda_drug = 0,
     mda_timesteps = NULL,
     mda_coverages = NULL,
-    mda_min_age = -1,
-    mda_max_age = -1,
+    mda_min_ages = -1,
+    mda_max_ages = -1,
     smc = FALSE,
     smc_drug = 0,
     smc_timesteps = NULL,
     smc_coverages = NULL,
-    smc_min_age = -1,
-    smc_max_age = -1,
+    smc_min_ages = -1,
+    smc_max_ages = -1,
     # tbv
     tbv = FALSE,
     tbv_mt = 35,
diff --git a/man/CorrelationParameters.Rd b/man/CorrelationParameters.Rd
index b1d22578..c2e6ada7 100644
--- a/man/CorrelationParameters.Rd
+++ b/man/CorrelationParameters.Rd
@@ -14,17 +14,17 @@ Describes an event in the simulation
 \section{Methods}{
 \subsection{Public methods}{
 \itemize{
-\item \href{#method-new}{\code{CorrelationParameters$new()}}
-\item \href{#method-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}}
-\item \href{#method-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}}
-\item \href{#method-sigma}{\code{CorrelationParameters$sigma()}}
-\item \href{#method-mvnorm}{\code{CorrelationParameters$mvnorm()}}
-\item \href{#method-clone}{\code{CorrelationParameters$clone()}}
+\item \href{#method-CorrelationParameters-new}{\code{CorrelationParameters$new()}}
+\item \href{#method-CorrelationParameters-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}}
+\item \href{#method-CorrelationParameters-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}}
+\item \href{#method-CorrelationParameters-sigma}{\code{CorrelationParameters$sigma()}}
+\item \href{#method-CorrelationParameters-mvnorm}{\code{CorrelationParameters$mvnorm()}}
+\item \href{#method-CorrelationParameters-clone}{\code{CorrelationParameters$clone()}}
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-new"></a>}}
-\if{latex}{\out{\hypertarget{method-new}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-new"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-new}{}}}
 \subsection{Method \code{new()}}{
 initialise correlation parameters
 \subsection{Usage}{
@@ -40,8 +40,8 @@ initialise correlation parameters
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-inter_round_rho"></a>}}
-\if{latex}{\out{\hypertarget{method-inter_round_rho}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-inter_round_rho"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_round_rho}{}}}
 \subsection{Method \code{inter_round_rho()}}{
 Add rho between rounds
 \subsection{Usage}{
@@ -60,8 +60,8 @@ the intervention}
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-inter_intervention_rho"></a>}}
-\if{latex}{\out{\hypertarget{method-inter_intervention_rho}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-inter_intervention_rho"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_intervention_rho}{}}}
 \subsection{Method \code{inter_intervention_rho()}}{
 Add rho between interventions
 \subsection{Usage}{
@@ -83,8 +83,8 @@ the intervention}
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-sigma"></a>}}
-\if{latex}{\out{\hypertarget{method-sigma}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-sigma"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-sigma}{}}}
 \subsection{Method \code{sigma()}}{
 Standard deviation of each intervention between rounds
 \subsection{Usage}{
@@ -93,8 +93,8 @@ Standard deviation of each intervention between rounds
 
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-mvnorm"></a>}}
-\if{latex}{\out{\hypertarget{method-mvnorm}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-mvnorm"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-mvnorm}{}}}
 \subsection{Method \code{mvnorm()}}{
 multivariate norm draws for these parameters
 \subsection{Usage}{
@@ -103,8 +103,8 @@ multivariate norm draws for these parameters
 
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-clone"></a>}}
-\if{latex}{\out{\hypertarget{method-clone}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-clone"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-clone}{}}}
 \subsection{Method \code{clone()}}{
 The objects of this class are cloneable with this method.
 \subsection{Usage}{
diff --git a/man/set_equilibrium.Rd b/man/set_equilibrium.Rd
index f41905bd..f1c83c54 100644
--- a/man/set_equilibrium.Rd
+++ b/man/set_equilibrium.Rd
@@ -9,7 +9,8 @@ set_equilibrium(parameters, init_EIR, eq_params = NULL)
 \arguments{
 \item{parameters}{model parameters to update}
 
-\item{init_EIR}{the desired initial EIR}
+\item{init_EIR}{the desired initial EIR (infectious bites per day over the entire human
+population)}
 
 \item{eq_params}{parameters from the malariaEquilibrium package, if null.
 The default malariaEquilibrium parameters will be used}
diff --git a/man/set_mda.Rd b/man/set_mda.Rd
index 68db2958..0b6f81c0 100644
--- a/man/set_mda.Rd
+++ b/man/set_mda.Rd
@@ -4,7 +4,7 @@
 \alias{set_mda}
 \title{Parameterise a Mass Drug Administration}
 \usage{
-set_mda(parameters, drug, timesteps, coverages, min_age, max_age)
+set_mda(parameters, drug, timesteps, coverages, min_ages, max_ages)
 }
 \arguments{
 \item{parameters}{a list of parameters to modify}
@@ -13,12 +13,12 @@ set_mda(parameters, drug, timesteps, coverages, min_age, max_age)
 
 \item{timesteps}{a vector of timesteps for each round of mda}
 
-\item{coverages}{the proportion of the target population who recieve each
+\item{coverages}{a vector of the proportion of the target population who receive each
 round}
 
-\item{min_age}{the minimum age of the target population exclusive (in timesteps)}
+\item{min_ages}{a vector of minimum ages of the target population for each round exclusive (in timesteps)}
 
-\item{max_age}{the maximum age of the target population exclusive (in timesteps)
+\item{max_ages}{a vector of maximum ages of the target population for each round exclusive (in timesteps)
 drug}
 }
 \description{
diff --git a/man/set_smc.Rd b/man/set_smc.Rd
index adf66701..95004a51 100644
--- a/man/set_smc.Rd
+++ b/man/set_smc.Rd
@@ -4,7 +4,7 @@
 \alias{set_smc}
 \title{Parameterise a Seasonal Malaria Chemoprevention}
 \usage{
-set_smc(parameters, drug, timesteps, coverages, min_age, max_age)
+set_smc(parameters, drug, timesteps, coverages, min_ages, max_ages)
 }
 \arguments{
 \item{parameters}{a list of parameters to modify}
@@ -13,12 +13,12 @@ set_smc(parameters, drug, timesteps, coverages, min_age, max_age)
 
 \item{timesteps}{a vector of timesteps for each round of smc}
 
-\item{coverages}{the proportion of the target population who recieve each
+\item{coverages}{a vector of the proportion of the target population who receive each
 round}
 
-\item{min_age}{the minimum age of the target population exclusive (in timesteps)}
+\item{min_ages}{a vector of minimum ages of the target population for each round exclusive (in timesteps)}
 
-\item{max_age}{the maximum age of the target population exclusive (in timesteps)
+\item{max_ages}{a vector of maximum ages of the target population for each round exclusive (in timesteps)
 drug}
 }
 \description{
diff --git a/tests/testthat/test-mda.R b/tests/testthat/test-mda.R
index d9d65ada..034501ec 100644
--- a/tests/testthat/test-mda.R
+++ b/tests/testthat/test-mda.R
@@ -1,3 +1,36 @@
+test_that('MDA inputs align', {
+  parameters <- get_parameters()
+  expect_error(
+    parameters <- set_mda(
+      parameters,
+      drug = 1,
+      timesteps = c(50, 150),
+      coverages= 1,
+      min_ages = c(5* 365, 5 * 365),
+      max_ages = c(10 * 365, 10 * 365)
+    ), "coverages and timesteps do no align")
+  
+  expect_error(
+    parameters <- set_mda(
+      parameters,
+      drug = 1,
+      timesteps = c(50, 150),
+      coverages = c(1, 1),
+      min_ages = 5 * 365,
+      max_ages = c(10 * 365, 10 * 365)
+    ), "minimum ages and timesteps do no align")
+  
+  expect_error(
+    parameters <- set_mda(
+      parameters,
+      drug = 1,
+      timesteps = c(50, 150),
+      coverages= c(1, 1),
+      min_ages = c(5* 365, 5 * 365),
+      max_ages = 10 * 365
+    ), "maximum ages and timesteps do no align")
+ 
+})
 
 test_that('MDA moves the diseased and non-diseased population correctly', {
   timestep <- 50
@@ -9,11 +42,11 @@ test_that('MDA moves the diseased and non-diseased population correctly', {
     drug = 1,
     timesteps = c(50, 150),
     coverages= c(1, 1),
-    min_age = 5 * 365,
-    max_age = 10 * 365
+    min_ages = c(5 * 365, 0),
+    max_ages = c(10 * 365, 100 * 365)
   )
   events <- create_events(parameters)
-
+  
   variables <- list(
     state = mock_category(
       c('D', 'S', 'A', 'U', 'Tr'),
@@ -24,28 +57,28 @@ test_that('MDA moves the diseased and non-diseased population correctly', {
     drug_time = mock_double(c(1, 2, 3, 4)),
     drug = mock_double(c(1, 2, 1, 2))
   )
-
+  
   events$mda_administer <- mock_event(events$mda_administer)
-
+  
   listener <- create_mda_listeners(
     variables,
     events$mda_administer,
     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,
     get_correlation_parameters(parameters),
     'mda',
     parameters,
     renderer
   )
-
+  
   mockery::stub(listener, 'bernoulli', mockery::mock(c(TRUE, TRUE)))
   mock_correlation <- mockery::mock(c(TRUE, TRUE))
   mockery::stub(listener, 'sample_intervention', mock_correlation)
   listener(timestep)
-
+  
   expect_equal(
     mockery::mock_args(mock_correlation)[[1]][[1]],
     c(3, 4)
@@ -55,35 +88,118 @@ test_that('MDA moves the diseased and non-diseased population correctly', {
     'Tr',
     3
   )
-
+  
   expect_bitset_update(
     variables$state$queue_update_mock(),
     'S',
     4,
     call = 2
   )
-
+  
   expect_bitset_update(
     variables$infectivity$queue_update_mock(),
     c(.3, .4) * SP_AQ_params[[2]],
     c(3, 4)
   )
-
+  
   expect_bitset_update(
     variables$drug$queue_update_mock(),
     1,
     c(3, 4)
   )
-
+  
   expect_bitset_update(
     variables$drug_time$queue_update_mock(),
     50,
     c(3, 4)
   )
-
+  
   mockery::expect_args(
     events$mda_administer$schedule,
     1,
     100
   )
 })
+
+test_that('MDA moves the diseased and non-diseased population correctly - second round, varying age range', {
+  timestep <- 150
+  renderer <- individual::Render$new(timestep)
+  parameters <- get_parameters(list(human_population = 4))
+  parameters <- set_drugs(parameters, list(SP_AQ_params))
+  parameters <- set_mda(
+    parameters,
+    drug = 1,
+    timesteps = c(50, 150),
+    coverages= c(1, 1),
+    min_ages = c(5 * 365, 0),
+    max_ages = c(10 * 365, 100 * 365)
+  )
+  events <- create_events(parameters)
+  
+  variables <- list(
+    state = mock_category(
+      c('D', 'S', 'A', 'U', 'Tr'),
+      c('D', 'S', 'A', 'U')
+    ),
+    birth = mock_double(-365 * c(2, 20, 5, 7)),
+    infectivity = mock_double(c(.1, .2, .3, .4)),
+    drug_time = mock_double(c(1, 2, 3, 4)),
+    drug = mock_double(c(1, 2, 1, 2))
+  )
+  
+  events$mda_administer <- mock_event(events$mda_administer)
+  
+  listener <- create_mda_listeners(
+    variables,
+    events$mda_administer,
+    parameters$mda_drug,
+    parameters$mda_timesteps,
+    parameters$mda_coverages,
+    parameters$mda_min_ages,
+    parameters$mda_max_ages,
+    get_correlation_parameters(parameters),
+    'mda',
+    parameters,
+    renderer
+  )
+  
+  mockery::stub(listener, 'bernoulli', mockery::mock(c(TRUE, TRUE, TRUE, TRUE)))
+  mock_correlation <- mockery::mock(c(TRUE, TRUE, TRUE, TRUE))
+  mockery::stub(listener, 'sample_intervention', mock_correlation)
+  listener(timestep)
+  
+  expect_equal(
+    mockery::mock_args(mock_correlation)[[1]][[1]],
+    c(1, 2, 3, 4)
+  )
+  expect_bitset_update(
+    variables$state$queue_update_mock(),
+    'Tr',
+    c(1, 3)
+  )
+  
+  expect_bitset_update(
+    variables$state$queue_update_mock(),
+    'S',
+    c(2, 4),
+    call = 2
+  )
+  
+  expect_bitset_update(
+    variables$infectivity$queue_update_mock(),
+    c(.1, .2, .3, .4) * SP_AQ_params[[2]],
+    c(1, 2, 3, 4)
+  )
+  
+  expect_bitset_update(
+    variables$drug$queue_update_mock(),
+    1,
+    c(1, 2, 3, 4)
+  )
+  
+  expect_bitset_update(
+    variables$drug_time$queue_update_mock(),
+    150,
+    c(1, 2, 3, 4)
+  )
+})
diff --git a/tests/testthat/test-mortality-integration.R b/tests/testthat/test-mortality-integration.R
index b4ce6e04..138f3d19 100644
--- a/tests/testthat/test-mortality-integration.R
+++ b/tests/testthat/test-mortality-integration.R
@@ -7,8 +7,8 @@ test_that('mortality_process resets humans correctly', {
     drug = 1,
     timesteps = c(50, 150),
     coverages= c(1, 1),
-    min_age = 5 * 365,
-    max_age = 10 * 365
+    min_ages = c(5 * 365, 5 * 365),
+    max_ages = c(10 * 365, 10 * 365)
   )
   events <- create_events(parameters)
   variables <- create_variables(parameters)
diff --git a/vignettes/MDA.Rmd b/vignettes/MDA.Rmd
index 1b1bba09..7ac3f20c 100644
--- a/vignettes/MDA.Rmd
+++ b/vignettes/MDA.Rmd
@@ -109,8 +109,8 @@ mdaparams <- set_mda(
   drug = 1,
   timesteps = mda_events$timestep,
   coverages = rep(.8, 2),
-  min_age = 0,
-  max_age = 200 * 365
+  min_ages = rep(0, length(mda_events$timestep)),
+  max_ages = rep(200 * 365, length(mda_events$timestep))
 )
 
 output <- run_simulation(sim_length, mdaparams)
@@ -126,7 +126,7 @@ add_mda_lines(plot_prevalence(output), mda_events)
 
 ## SMC
 
-This is a dose of SP-AQ to 90% of 2 - 11 year olds once a year a month before the peak season for mosquitos.
+This is a dose of SP-AQ to 90% of 2 - 11 year olds once a year a month before the peak season for mosquitoes.
 
 ```{r}
 smcparams <- simparams
@@ -142,8 +142,8 @@ smcparams <- set_smc(
   drug = 1,
   timesteps = smc_events$timestep,
   coverages = rep(.9, 2),
-  min_age = 2 * 365 - 1,
-  max_age = 11 * 365
+  min_ages = rep(2 * 365 - 1, length(smc_events$timestep)),
+  max_ages = rep(11 * 365, length(smc_events$timestep))
 )
 
 output <- run_simulation(sim_length, smcparams)

From 6e08e74b828694d7e2cc858644d82e8b39274f6a Mon Sep 17 00:00:00 2001
From: Hillary Topazian <80535782+htopazian@users.noreply.github.com>
Date: Fri, 22 Jul 2022 09:27:03 +0100
Subject: [PATCH 10/12] mixing matrix row & col tests

- rows and cols must sum to 1
- sums rounded to one decimal place
---
 R/model.R | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/R/model.R b/R/model.R
index 31b84d3e..86262e49 100644
--- a/R/model.R
+++ b/R/model.R
@@ -130,7 +130,7 @@ run_simulation <- function(
 #' @param correlations a list of correlation parameters for each population
 #' (default: NULL)
 #' @param mixing matrix of mixing coefficients for infectivity towards
-#' mosquitoes
+#' mosquitoes. Each element must be between 0 and 1 and all rows and columns must sum to 1.
 #' @return a list of dataframe of results
 #' @export
 run_metapop_simulation <- function(
@@ -146,6 +146,12 @@ run_metapop_simulation <- function(
   if (nrow(mixing) != length(parameters)) {
     stop('mixing matrix rows must match length of parameters')
   }
+  if (!all(round(rowSums(mixing), 1) == 1)) {
+    stop('all mixing matrix rows must sum to 1')
+  }
+  if (!all(round(colSums(mixing), 1) == 1)) {
+    stop('all mixing matrix columns must sum to 1')
+  }
   if (is.null(correlations)) {
     correlations <- lapply(parameters, get_correlation_parameters)
   }

From 9860a419686c526defec6c809299cebbcc8ca142 Mon Sep 17 00:00:00 2001
From: Hillary Topazian <80535782+htopazian@users.noreply.github.com>
Date: Fri, 22 Jul 2022 12:38:44 +0100
Subject: [PATCH 11/12] metapopulation vignette

---
 vignettes/Metapopulation.Rmd | 122 +++++++++++++++++++++++++++++++++++
 1 file changed, 122 insertions(+)
 create mode 100644 vignettes/Metapopulation.Rmd

diff --git a/vignettes/Metapopulation.Rmd b/vignettes/Metapopulation.Rmd
new file mode 100644
index 00000000..bcfff5de
--- /dev/null
+++ b/vignettes/Metapopulation.Rmd
@@ -0,0 +1,122 @@
+---
+title: "Metapopulation"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Treatment}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+suppressPackageStartupMessages(library(ggplot2))
+library(malariasimulation)
+library(malariaEquilibrium)
+suppressPackageStartupMessages(library(dplyr))
+library(purrr)
+library(data.table)
+```
+
+# Parameterisation
+
+The metapopulation model runs the individual-based model for multiple parameterized units (countries, regions, admins, etc.) simultaneously. The inputted mixing matrix allows for the transmission of one unit to affect the transmission of another unit. 'Mixing' in the model occurs through the variables `foim_m` and `EIR`. 
+
+Here we will set up a case study of three distinct units and compare output with various transmission mixing patterns. 
+
+```{r}
+# set variables
+year <- 365
+human_population <- 5000
+sim_length <- 5 * year
+EIR_vector <- c(1, 2, 5)
+
+# get parameters
+ms_parameterize <- function(x){ # index of EIR
+  
+  params <- get_parameters(list(human_population = human_population,
+                                model_seasonality = FALSE,
+                                individual_mosquitoes = FALSE))
+  
+  # setting treatment
+  params <- set_drugs(params, list(AL_params))
+  params <- set_clinical_treatment(params, drug = 1, timesteps = 1, coverages = 0.40)
+
+  params <- set_equilibrium(params, init_EIR = EIR_vector[x])
+  
+  return(params)
+  
+}
+
+# creating a list of three parameter lists
+paramslist <- lapply(seq(1, length(EIR_vector), 1), ms_parameterize)
+
+```
+
+# Modelling 
+
+Our parameters for the three distinct units are stored in the object `paramslist`. Next we will run the metapopulation model with these parameters. We will plug in three mixing matrices - A) isolated, B) semi-mixed, C) perfectly mixed. 
+
+```{r, warning = F}
+# isolated
+mix_1 <- diag(length(EIR_vector))
+
+# semi-mixed
+mix_2 <- matrix(c(0.8, 0.1, 0.1, 
+                  0.1, 0.8, 0.1,
+                  0.1, 0.1, 0.8),
+               nrow = 3, ncol = 3)
+
+# perfectly-mixed
+mix_3 <- matrix(rep(.33, 9), nrow = 3, ncol = 3)
+
+# run model
+set.seed(123)
+
+metapop_loop <- function(mixing, mixname){ # mixing matrix and matrix name
+
+  output <- run_metapop_simulation(timesteps = sim_length,
+                                   parameters = paramslist,
+                                   correlations = NULL,
+                                   mixing = mixing)
+  
+  # convert to dataframe and label EIR and mixing matrix type
+  output <- data.table::rbindlist(output, idcol = 'EIR') %>%
+  mutate(EIR = c(sort(rep(EIR_vector, sim_length))),
+         mix = mixname)
+  
+  return(output)
+  
+}
+
+output <- map2_dfr(list(mix_1, mix_2, mix_3),
+                   c('isolated', 'semi-mixed', 'perfectly-mixed'),
+                   metapop_loop)
+
+# get mean PfPR 2-10 by year
+output <- output %>% 
+  mutate(prev2to10 = p_detect_730_3650 / n_730_3650) %>%
+  mutate(year = ceiling(timestep/year)) %>%
+  group_by(mix, EIR, year) %>%
+  summarize(prev2to10 = mean(prev2to10, na.rm = T),
+         mix = factor(mix, levels = c('isolated', 'semi-mixed', 'perfectly-mixed')))
+
+```
+
+Now let's visualize the results of mixing on PfPR2-10:
+
+```{r}
+# plot
+ggplot(data = output) + 
+  geom_line(aes(x = year, y = prev2to10, color = factor(EIR))) +
+  facet_wrap(~ mix) +
+  labs(x = 'time (years)',
+       y = 'PfPR 2-10 (month)',
+       color = 'EIR') +
+  theme_classic()
+```

From 01789fd4c4b07c551e63b5370fa7e0476ae3bbe3 Mon Sep 17 00:00:00 2001
From: Hillary Topazian <80535782+htopazian@users.noreply.github.com>
Date: Fri, 22 Jul 2022 13:57:41 +0100
Subject: [PATCH 12/12] taking out vignette dependencies

---
 vignettes/Metapopulation.Rmd | 38 +++++++++++++++++++-----------------
 1 file changed, 20 insertions(+), 18 deletions(-)

diff --git a/vignettes/Metapopulation.Rmd b/vignettes/Metapopulation.Rmd
index bcfff5de..c3674a90 100644
--- a/vignettes/Metapopulation.Rmd
+++ b/vignettes/Metapopulation.Rmd
@@ -2,7 +2,7 @@
 title: "Metapopulation"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{Treatment}
+  %\VignetteIndexEntry{Metapopulation}
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
 ---
@@ -18,9 +18,6 @@ knitr::opts_chunk$set(
 suppressPackageStartupMessages(library(ggplot2))
 library(malariasimulation)
 library(malariaEquilibrium)
-suppressPackageStartupMessages(library(dplyr))
-library(purrr)
-library(data.table)
 ```
 
 # Parameterisation
@@ -32,7 +29,7 @@ Here we will set up a case study of three distinct units and compare output with
 ```{r}
 # set variables
 year <- 365
-human_population <- 5000
+human_population <- 1000
 sim_length <- 5 * year
 EIR_vector <- c(1, 2, 5)
 
@@ -78,7 +75,7 @@ mix_3 <- matrix(rep(.33, 9), nrow = 3, ncol = 3)
 # run model
 set.seed(123)
 
-metapop_loop <- function(mixing, mixname){ # mixing matrix and matrix name
+metapop_loop <- function(mixing, mixnam){ # mixing matrix 
 
   output <- run_metapop_simulation(timesteps = sim_length,
                                    parameters = paramslist,
@@ -86,25 +83,29 @@ metapop_loop <- function(mixing, mixname){ # mixing matrix and matrix name
                                    mixing = mixing)
   
   # convert to dataframe and label EIR and mixing matrix type
-  output <- data.table::rbindlist(output, idcol = 'EIR') %>%
-  mutate(EIR = c(sort(rep(EIR_vector, sim_length))),
-         mix = mixname)
+  output <-  do.call('rbind', output)
+  output$EIR <- c(sort(rep(EIR_vector, sim_length)))
   
   return(output)
   
 }
 
-output <- map2_dfr(list(mix_1, mix_2, mix_3),
-                   c('isolated', 'semi-mixed', 'perfectly-mixed'),
-                   metapop_loop)
+output1 <- metapop_loop(mix_1)
+output1$mix <- 'isolated'
+  
+output2 <- metapop_loop(mix_2)
+output2$mix <- 'semi-mixed'
+  
+output3 <- metapop_loop(mix_3)
+output3$mix <- 'perfectly-mixed'
+  
+output <- rbind(output1, output2, output3)
 
 # get mean PfPR 2-10 by year
-output <- output %>% 
-  mutate(prev2to10 = p_detect_730_3650 / n_730_3650) %>%
-  mutate(year = ceiling(timestep/year)) %>%
-  group_by(mix, EIR, year) %>%
-  summarize(prev2to10 = mean(prev2to10, na.rm = T),
-         mix = factor(mix, levels = c('isolated', 'semi-mixed', 'perfectly-mixed')))
+output$prev2to10 = output$p_detect_730_3650 / output$n_730_3650
+output$year = ceiling(output$timestep / 365)
+output$mix = factor(output$mix, levels = c('isolated', 'semi-mixed', 'perfectly-mixed'))
+output <- aggregate(prev2to10 ~ mix + EIR + year, data = output, FUN = mean)
 
 ```
 
@@ -115,6 +116,7 @@ Now let's visualize the results of mixing on PfPR2-10:
 ggplot(data = output) + 
   geom_line(aes(x = year, y = prev2to10, color = factor(EIR))) +
   facet_wrap(~ mix) +
+  scale_y_continuous(limits = c(0, 0.35)) + 
   labs(x = 'time (years)',
        y = 'PfPR 2-10 (month)',
        color = 'EIR') +