From cbfe929fb6c286d1491d5ae6f712bc65c416e461 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 26 Aug 2021 11:31:49 +0100
Subject: [PATCH 01/23] Fixes to stop EIR being < 0:

* Validate gonotrophic cycle
* Parameterise foraging time per species
---
 R/biting_process.R           |  4 ++--
 R/mosquito_biology.R         |  4 ++--
 R/vector_parameters.R        | 16 ++++++++++++++++
 tests/testthat/test-biting.R | 12 ++++++++++++
 4 files changed, 32 insertions(+), 4 deletions(-)
 create mode 100644 tests/testthat/test-biting.R

diff --git a/R/biting_process.R b/R/biting_process.R
index 15f65e2c..7da01dc0 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -221,7 +221,7 @@ calculate_infectious_individual <- function(
 }
 
 calculate_infectious_compartmental <- function(solver_states) {
-  solver_states[[ADULT_ODE_INDICES['Im']]]
+  max(solver_states[[ADULT_ODE_INDICES['Im']]], 0)
 }
 
 intervention_coefficient <- function(p_bitten) {
@@ -234,7 +234,7 @@ human_pi <- function(zeta, psi) {
 
 blood_meal_rate <- function(v, z, parameters) {
   gonotrophic_cycle <- get_gonotrophic_cycle(v, parameters)
-  interrupted_foraging_time <- parameters$foraging_time / (1 - z)
+  interrupted_foraging_time <- parameters$foraging_time[[v]] / (1 - z)
   1 / (interrupted_foraging_time + gonotrophic_cycle)
 }
 
diff --git a/R/mosquito_biology.R b/R/mosquito_biology.R
index d30f3714..0b1076cc 100644
--- a/R/mosquito_biology.R
+++ b/R/mosquito_biology.R
@@ -148,7 +148,7 @@ peak_season_offset <- function(parameters) {
 #' @noRd
 death_rate <- function(f, W, Z, species, parameters) {
   mum <- parameters$mum[[species]]
-  p1_0 <- exp(-mum * parameters$foraging_time)
+  p1_0 <- exp(-mum * parameters$foraging_time[[species]])
   gonotrophic_cycle <- get_gonotrophic_cycle(species, parameters)
   p2 <- exp(-mum * gonotrophic_cycle)
   p1 <- p1_0 * W / (1 - Z * p1_0)
@@ -157,7 +157,7 @@ death_rate <- function(f, W, Z, species, parameters) {
 
 get_gonotrophic_cycle <- function(v, parameters) {
   f <- parameters$blood_meal_rates[[v]]
-  gonotrophic_cycle <- 1 / f - parameters$foraging_time
+  gonotrophic_cycle <- 1 / f - parameters$foraging_time[[v]]
 }
 
 #' @title Update the individual mosquito model after biting
diff --git a/R/vector_parameters.R b/R/vector_parameters.R
index 888477e1..8ff25c54 100644
--- a/R/vector_parameters.R
+++ b/R/vector_parameters.R
@@ -2,6 +2,7 @@
 #' @details Default parameters:
 #' species: "gamb"
 #' blood_meal_rates: 0.3333333
+#' foraging_time: .69
 #' Q0: 0.92
 #' phi_bednets: 0.85
 #' phi_indoors: 0.90
@@ -13,6 +14,7 @@
 gamb_params <- list(
   species = 'gamb',
   blood_meal_rates = 1/3,
+  foraging_time = .69,
   Q0 = .92,
   phi_bednets = .85,
   phi_indoors = .90,
@@ -23,6 +25,7 @@ gamb_params <- list(
 #' @details Default parameters:
 #' species: "arab"
 #' blood_meal_rates: 0.3333333
+#' foraging_time: .69
 #' Q0: 0.71
 #' phi_bednets: 0.8
 #' phi_indoors: 0.86
@@ -34,6 +37,7 @@ gamb_params <- list(
 arab_params <- list(
   species = 'arab',
   blood_meal_rates = 1/3,
+  foraging_time = .69,
   Q0 = .71,
   phi_bednets = .8,
   phi_indoors = .86,
@@ -44,6 +48,7 @@ arab_params <- list(
 #' @details Default parameters:
 #' species: "fun"
 #' blood_meal_rates: 0.3333333
+#' foraging_time: .69
 #' Q0: 0.94
 #' phi_bednets: 0.78
 #' phi_indoors: 0.87
@@ -55,6 +60,7 @@ arab_params <- list(
 fun_params <- list(
   species = 'fun',
   blood_meal_rates = 1/3,
+  foraging_time = .69,
   Q0 = .94,
   phi_bednets = .78,
   phi_indoors = .87,
@@ -77,6 +83,7 @@ set_species <- function(parameters, species, proportions) {
   keys <- c(
     'species',
     'blood_meal_rates',
+    'foraging_time',
     'Q0',
     'phi_bednets',
     'phi_indoors',
@@ -86,6 +93,15 @@ set_species <- function(parameters, species, proportions) {
     parameters[[key]] <- rep(NA, length(species))
   }
   for (v in seq_along(species)) {
+    if (species[[v]]$foraging_time > 1 / species[[v]]$blood_meal_rates) {
+      stop(
+        sprintf(
+          "blood meal time (%f) must be >= foraging time (%f)",
+          1 / species[[v]]$blood_meal_rates,
+          species[[v]]$foraging_time
+        )
+      )
+    }
     for (key in keys) {
       parameters[[key]][[v]] <- species[[v]][[key]]
     }
diff --git a/tests/testthat/test-biting.R b/tests/testthat/test-biting.R
new file mode 100644
index 00000000..e6c19fb9
--- /dev/null
+++ b/tests/testthat/test-biting.R
@@ -0,0 +1,12 @@
+test_that('compartmental always gives positive infectious', {
+  solver_states <- rep(0, length(ODE_INDICES) + length(ADULT_ODE_INDICES))
+  solver_states[[ADULT_ODE_INDICES['Im']]] <- -1e-10
+  expect_gte(calculate_infectious_compartmental(solver_states), 0)
+})
+
+test_that('gonotrophic_cycle cannot be negative', {
+  params <- get_parameters()
+  vparams <- gamb_params
+  vparams$blood_meal_rates <- 5
+  expect_error(set_species(params, list(vparams), 1))
+})

From f6f2b5fea2b4001c191f29d9258b93534db009c9 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 9 Sep 2021 11:31:45 +0100
Subject: [PATCH 02/23] Add treated to n_inc_clinical

---
 R/events.R              | 21 ----------
 R/human_infection.R     | 24 +++++++----
 R/model.R               |  2 +-
 R/mortality_processes.R |  1 -
 R/render.R              | 89 +++++++++++++++++++++++------------------
 5 files changed, 68 insertions(+), 69 deletions(-)

diff --git a/R/events.R b/R/events.R
index 66905042..908a00e0 100644
--- a/R/events.R
+++ b/R/events.R
@@ -9,9 +9,6 @@ create_events <- function(parameters) {
     clinical_infection = individual::TargetedEvent$new(parameters$human_population),
     asymptomatic_infection = individual::TargetedEvent$new(parameters$human_population),
 
-    # whether the infection is detected
-    detection = individual::TargetedEvent$new(parameters$human_population), 
-
     # MDA events
     mda_administer = individual::Event$new(),
     smc_administer = individual::Event$new(),
@@ -165,24 +162,6 @@ attach_event_listeners <- function(
   # Progression
   # ===========
   # When infection events fire, schedule the next stages of infection
-
-  events$clinical_infection$add_listener(
-    create_clinical_incidence_renderer(
-      variables$birth,
-      parameters,
-      renderer
-    )
-  )
-
-  events$detection$add_listener(
-    create_incidence_renderer(
-      variables$birth,
-      variables$is_severe,
-      parameters,
-      renderer
-    )
-  )
-
   if (parameters$individual_mosquitoes) {
     events$mosquito_infection$add_listener(
       individual::update_category_listener(variables$mosquito_state, 'Im')
diff --git a/R/human_infection.R b/R/human_infection.R
index 625a96ad..a41a966e 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -36,6 +36,15 @@ simulate_infection <- function(
     timestep
   )
 
+  incidence_renderer(
+    variables$birth,
+    variables$is_severe,
+    parameters,
+    renderer,
+    timestep,
+    infected_humans
+  )
+
   if (infected_humans$size() > 0) {
     boost_immunity(
       variables$ica,
@@ -59,6 +68,14 @@ simulate_infection <- function(
     parameters
   )
 
+  clinical_incidence_renderer(
+    variables$birth,
+    parameters,
+    renderer,
+    timestep,
+    clinical_infections
+  )
+
   if (parameters$severe_enabled) {
     update_severe_disease(
       timestep,
@@ -69,12 +86,10 @@ simulate_infection <- function(
     )
   }
 
-
   treated <- calculate_treated(
     variables,
     clinical_infections,
     events$recovery,
-    events$detection,
     parameters,
     timestep,
     renderer
@@ -224,7 +239,6 @@ calculate_treated <- function(
   variables,
   clinical_infections,
   recovery,
-  detection,
   parameters,
   timestep,
   renderer
@@ -266,7 +280,6 @@ calculate_treated <- function(
       timestep,
       treated_index
     )
-    detection$schedule(treated_index, 0)
   }
   treated_index
 }
@@ -295,17 +308,14 @@ schedule_infections <- function(
     included
   )
 
-  # change to symptomatic
   if(to_infect$size() > 0) {
     infection_times <- log_uniform(to_infect$size(), parameters$de)
     events$clinical_infection$schedule(to_infect, 0)
-    events$detection$schedule(to_infect, 0)
   }
 
   if(to_infect_asym$size() > 0) {
     infection_times <- log_uniform(to_infect_asym$size(), parameters$de)
     events$asymptomatic_infection$schedule(to_infect_asym, 0)
-    events$detection$schedule(to_infect_asym, 0)
   }
 }
 
diff --git a/R/model.R b/R/model.R
index ce8c0d86..0d540d9d 100644
--- a/R/model.R
+++ b/R/model.R
@@ -32,7 +32,7 @@
 #'  * n: number of humans between an inclusive age range at this timestep. This
 #' defaults to n_730_3650. Other age ranges can be set with
 #' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
-#' * n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This
+#'  * n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This
 #' defaults to n_detect_730_3650. Other age ranges can be set with
 #' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
 #'  * n_severe: number of humans with a severe infection detectable by microscopy 
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index 7c07b911..2a350a94 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -72,7 +72,6 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
         'recovery',
         'clinical_infection',
         'asymptomatic_infection',
-        'detection',
         'throw_away_net',
         'rtss_mass_doses',
         'rtss_mass_booster',
diff --git a/R/render.R b/R/render.R
index 4fd2cb15..ab30f579 100644
--- a/R/render.R
+++ b/R/render.R
@@ -81,37 +81,44 @@ create_prevelance_renderer <- function(
 #' @param is_severe variable for if individual has severe malaria
 #' @param parameters model parameters
 #' @param renderer model renderer
+#' @param timestep current target
+#' @param target newly infected population
 #' 
 #' @noRd
-create_incidence_renderer <- function(birth, is_severe, parameters, renderer) {
-  function(timestep, target) {
-    severe <- is_severe$get_index_of('yes')$and(target)
-    for (i in seq_along(parameters$incidence_rendering_min_ages)) {
-      lower <- parameters$incidence_rendering_min_ages[[i]]
-      upper <- parameters$incidence_rendering_max_ages[[i]]
-      in_age <- in_age_range(birth, timestep, lower, upper)
-      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
-      renderer$render(
-        paste0('n_inc_', lower, '_', upper),
-        in_age$and(target)$size(),
-        timestep
-      )
-    }
-    for (i in seq_along(parameters$severe_incidence_rendering_min_ages)) {
-      lower <- parameters$severe_incidence_rendering_min_ages[[i]]
-      upper <- parameters$severe_incidence_rendering_max_ages[[i]]
-      in_age <- in_age_range(birth, timestep, lower, upper)
-      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
-      renderer$render(
-        paste0('n_inc_severe_', lower, '_', upper),
-        in_age$and(target)$and(severe)$size(),
-        timestep
-      )
-    }
+incidence_renderer <- function(
+  birth,
+  is_severe,
+  parameters,
+  renderer,
+  timestep,
+  target
+  ) {
+  severe <- is_severe$get_index_of('yes')$and(target)
+  for (i in seq_along(parameters$incidence_rendering_min_ages)) {
+    lower <- parameters$incidence_rendering_min_ages[[i]]
+    upper <- parameters$incidence_rendering_max_ages[[i]]
+    in_age <- in_age_range(birth, timestep, lower, upper)
+    renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+    renderer$render(
+      paste0('n_inc_', lower, '_', upper),
+      in_age$and(target)$size(),
+      timestep
+    )
+  }
+  for (i in seq_along(parameters$severe_incidence_rendering_min_ages)) {
+    lower <- parameters$severe_incidence_rendering_min_ages[[i]]
+    upper <- parameters$severe_incidence_rendering_max_ages[[i]]
+    in_age <- in_age_range(birth, timestep, lower, upper)
+    renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+    renderer$render(
+      paste0('n_inc_severe_', lower, '_', upper),
+      in_age$and(target)$and(severe)$size(),
+      timestep
+    )
   }
 }
 
-#' @title Render incidence statistics
+#' @title Render clinical incidence statistics
 #' 
 #' @description renders clinical incidence (new for this timestep)
 #' 
@@ -120,19 +127,23 @@ create_incidence_renderer <- function(birth, is_severe, parameters, renderer) {
 #' @param renderer model renderer
 #' 
 #' @noRd
-create_clinical_incidence_renderer <- function(birth, parameters, renderer) {
-  function(timestep, target) {
-    for (i in seq_along(parameters$clinical_incidence_rendering_min_ages)) {
-      lower <- parameters$clinical_incidence_rendering_min_ages[[i]]
-      upper <- parameters$clinical_incidence_rendering_max_ages[[i]]
-      in_age <- in_age_range(birth, timestep, lower, upper)
-      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
-      renderer$render(
-        paste0('n_inc_clinical_', lower, '_', upper),
-        in_age$and(target)$size(),
-        timestep
-      )
-    }
+clinical_incidence_renderer <- function(
+  birth,
+  parameters,
+  renderer,
+  target,
+  timestep
+  ) {
+  for (i in seq_along(parameters$clinical_incidence_rendering_min_ages)) {
+    lower <- parameters$clinical_incidence_rendering_min_ages[[i]]
+    upper <- parameters$clinical_incidence_rendering_max_ages[[i]]
+    in_age <- in_age_range(birth, timestep, lower, upper)
+    renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+    renderer$render(
+      paste0('n_inc_clinical_', lower, '_', upper),
+      in_age$and(target)$size(),
+      timestep
+    )
   }
 }
 

From 8fc5eb30607eb71d07368b7a80cba4099b8e7ac9 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 9 Sep 2021 13:19:46 +0100
Subject: [PATCH 03/23] Add tests for new incidence rendering

---
 R/human_infection.R                         |  8 +++---
 R/render.R                                  |  4 +--
 tests/testthat/test-infection-integration.R | 27 ++---------------
 tests/testthat/test-render.R                | 32 +++++++++++++++------
 4 files changed, 31 insertions(+), 40 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index a41a966e..f5910991 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -41,8 +41,8 @@ simulate_infection <- function(
     variables$is_severe,
     parameters,
     renderer,
-    timestep,
-    infected_humans
+    infected_humans,
+    timestep
   )
 
   if (infected_humans$size() > 0) {
@@ -72,8 +72,8 @@ simulate_infection <- function(
     variables$birth,
     parameters,
     renderer,
-    timestep,
-    clinical_infections
+    clinical_infections,
+    timestep
   )
 
   if (parameters$severe_enabled) {
diff --git a/R/render.R b/R/render.R
index ab30f579..0757aa47 100644
--- a/R/render.R
+++ b/R/render.R
@@ -90,8 +90,8 @@ incidence_renderer <- function(
   is_severe,
   parameters,
   renderer,
-  timestep,
-  target
+  target,
+  timestep
   ) {
   severe <- is_severe$get_index_of('yes')$and(target)
   for (i in seq_along(parameters$incidence_rendering_min_ages)) {
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index f55e484c..df9f6231 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -38,6 +38,8 @@ test_that('simulate_infection integrates different types of infection and schedu
   mockery::stub(simulate_infection, 'update_severe_disease', severe_infection_mock)
   mockery::stub(simulate_infection, 'calculate_treated', treated_mock)
   mockery::stub(simulate_infection, 'schedule_infections', schedule_mock)
+  mockery::stub(simulate_infection, 'incidence_renderer', mockery::mock())
+  mockery::stub(simulate_infection, 'clinical_incidence_renderer', mockery::mock())
   simulate_infection(
     variables,
     events,
@@ -91,7 +93,6 @@ test_that('simulate_infection integrates different types of infection and schedu
     variables,
     clinical,
     events$recovery,
-    events$detection,
     parameters,
     timestep,
     renderer
@@ -245,9 +246,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   )
 
   recovery_mock <- mockery::mock()
-  detection_mock <- mockery::mock()
   mockery::stub(calculate_treated, 'recovery$schedule', recovery_mock)
-  mockery::stub(calculate_treated, 'detection$schedule', detection_mock)
 
   seek_treatment <- individual::Bitset$new(4)$insert(c(1, 2, 4))
   mockery::stub(
@@ -267,7 +266,6 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     variables,
     clinical_infections,
     events$recovery,
-    events$detection,
     parameters,
     timestep,
     mock_render(timestep)
@@ -296,12 +294,6 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   )
   expect_bitset_update(variables$drug$queue_update, c(2, 1), c(1, 4))
   expect_bitset_update(variables$drug_time$queue_update, 5, c(1, 4))
-
-  expect_bitset_schedule(
-    detection_mock,
-    c(1, 4),
-    0
-  )
 })
 
 test_that('schedule_infections correctly schedules new infections', {
@@ -309,10 +301,8 @@ test_that('schedule_infections correctly schedules new infections', {
 
   clinical_mock <- mockery::mock()
   asym_mock <- mockery::mock()
-  detection_mock <- mockery::mock()
 
   events <- list(
-    detection = list(schedule = detection_mock),
     clinical_infection = list(
       schedule = clinical_mock
     ),
@@ -341,24 +331,11 @@ test_that('schedule_infections correctly schedules new infections', {
     0
   )
 
-  expect_bitset_schedule(
-    detection_mock,
-    c(5, 6, 13, 14, 15),
-    0
-  )
-
   expect_bitset_schedule(
     asym_mock,
     c(1, 2, 3, 4, 16, 17, 18, 19, 20),
     0,
   )
-
-  expect_bitset_schedule(
-    detection_mock,
-    c(1, 2, 3, 4, 16, 17, 18, 19, 20),
-    0,
-    call = 2
-  )
 })
 
 test_that('prophylaxis is considered for medicated humans', {
diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R
index ea980230..94ed32f7 100644
--- a/tests/testthat/test-render.R
+++ b/tests/testthat/test-render.R
@@ -14,7 +14,6 @@ test_that('that default rendering works', {
     c('no', 'no', 'yes', 'no')
   )
 
-
   renderer <- mock_render(1)
   process <- create_prevelance_renderer(
     state,
@@ -163,9 +162,13 @@ test_that('that clinical incidence rendering works', {
   )
 
   renderer <- mock_render(1)
-  process <- create_clinical_incidence_renderer(birth, parameters, renderer)
-
-  process(timestep, individual::Bitset$new(4)$insert(c(1, 2, 4)))
+  clinical_incidence_renderer(
+    birth,
+    parameters,
+    renderer,
+    individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    timestep
+  )
 
   mockery::expect_args( 
     renderer$render_mock(),
@@ -221,9 +224,14 @@ test_that('that incidence rendering works', {
   )
 
   renderer <- mock_render(1)
-  process <- create_incidence_renderer(birth, is_severe, parameters, renderer)
-
-  process(timestep, individual::Bitset$new(4)$insert(c(1, 2, 4)))
+  process <- incidence_renderer(
+    birth,
+    is_severe,
+    parameters,
+    renderer,
+    individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    timestep
+  )
 
   mockery::expect_args(
     renderer$render_mock(),
@@ -278,9 +286,15 @@ test_that('that severe incidence rendering works', {
   )
 
   renderer <- mock_render(1)
-  process <- create_incidence_renderer(birth, is_severe, parameters, renderer)
+  process <- incidence_renderer(
+    birth,
+    is_severe,
+    parameters,
+    renderer,
+    individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    timestep
+  )
 
-  process(timestep, individual::Bitset$new(4)$insert(c(1, 2, 4)))
 
   mockery::expect_args(
     renderer$render_mock(),

From eacf1e43807de1d1a7faf30135734d2e8a5b45e1 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 9 Sep 2021 13:22:16 +0100
Subject: [PATCH 04/23] Bump version and add news

---
 DESCRIPTION | 2 +-
 NEWS.md     | 4 ++++
 2 files changed, 5 insertions(+), 1 deletion(-)
 create mode 100644 NEWS.md

diff --git a/DESCRIPTION b/DESCRIPTION
index ebf727d7..07be5235 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: malariasimulation
 Title: An individual based model for malaria
-Version: 1.1.0
+Version: 1.2.0
 Authors@R: c(
   person(
     given = "Giovanni",
diff --git a/NEWS.md b/NEWS.md
new file mode 100644
index 00000000..726b4e27
--- /dev/null
+++ b/NEWS.md
@@ -0,0 +1,4 @@
+# malariasimulation 1.2.0
+
+  * Added a `NEWS.md` file to track changes to the package.
+  * n_inc_clinical includes treated humans

From 10f277d0f55bf9a300eda04c8b5ef384e186bde7 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 9 Sep 2021 13:32:53 +0100
Subject: [PATCH 05/23] Update documentation

---
 man/arab_params.Rd    |  3 ++-
 man/fun_params.Rd     |  3 ++-
 man/gamb_params.Rd    |  3 ++-
 man/get_parameters.Rd |  5 ++---
 man/run_simulation.Rd | 18 ++++++++++--------
 man/set_mass_rtss.Rd  |  4 +++-
 man/set_rtss_epi.Rd   |  6 ++++--
 7 files changed, 25 insertions(+), 17 deletions(-)

diff --git a/man/arab_params.Rd b/man/arab_params.Rd
index b5655b25..144bbf70 100644
--- a/man/arab_params.Rd
+++ b/man/arab_params.Rd
@@ -5,7 +5,7 @@
 \alias{arab_params}
 \title{Preset parameters for the An. arabiensis vector}
 \format{
-An object of class \code{list} of length 6.
+An object of class \code{list} of length 7.
 }
 \usage{
 arab_params
@@ -17,6 +17,7 @@ Preset parameters for the An. arabiensis vector
 Default parameters:
 species: "arab"
 blood_meal_rates: 0.3333333
+foraging_time: .69
 Q0: 0.71
 phi_bednets: 0.8
 phi_indoors: 0.86
diff --git a/man/fun_params.Rd b/man/fun_params.Rd
index b5d0eaff..6063bfbc 100644
--- a/man/fun_params.Rd
+++ b/man/fun_params.Rd
@@ -5,7 +5,7 @@
 \alias{fun_params}
 \title{Preset parameters for the An. funestus vector}
 \format{
-An object of class \code{list} of length 6.
+An object of class \code{list} of length 7.
 }
 \usage{
 fun_params
@@ -17,6 +17,7 @@ Preset parameters for the An. funestus vector
 Default parameters:
 species: "fun"
 blood_meal_rates: 0.3333333
+foraging_time: .69
 Q0: 0.94
 phi_bednets: 0.78
 phi_indoors: 0.87
diff --git a/man/gamb_params.Rd b/man/gamb_params.Rd
index 434c46fa..fef4b29f 100644
--- a/man/gamb_params.Rd
+++ b/man/gamb_params.Rd
@@ -5,7 +5,7 @@
 \alias{gamb_params}
 \title{Preset parameters for the An. gambiae s.s vector}
 \format{
-An object of class \code{list} of length 6.
+An object of class \code{list} of length 7.
 }
 \usage{
 gamb_params
@@ -17,6 +17,7 @@ Preset parameters for the An. gambiae s.s vector
 Default parameters:
 species: "gamb"
 blood_meal_rates: 0.3333333
+foraging_time: .69
 Q0: 0.92
 phi_bednets: 0.85
 phi_indoors: 0.90
diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd
index 0f6f8d11..5fea735c 100644
--- a/man/get_parameters.Rd
+++ b/man/get_parameters.Rd
@@ -233,8 +233,8 @@ outputs; default = 730
 \item incidence_rendering_min_ages - the minimum ages for incidence
 outputs (includes asymptomatic microscopy +); default = turned off
 \item incidence_rendering_max_ages - the corresponding max ages; default = turned off
-clinical_incidence_rendering_min_ages - the minimum ages for clinical incidence outputs (symptomatic); default = 0
-clinical_incidence_rendering_max_ages - the corresponding max ages; default = 1825
+\item clinical_incidence_rendering_min_ages - the minimum ages for clinical incidence outputs (symptomatic); default = 0
+\item clinical_incidence_rendering_max_ages - the corresponding max ages; default = 1825
 \item severe_prevalence_rendering_min_ages - the minimum ages for severe
 prevalence outputs; default = turned off
 \item severe_prevalence_rendering_max_ages - the corresponding max ages; default = turned off
@@ -253,7 +253,6 @@ individually or compartmentally; default = TRUE
 \item r_tol - the relative tolerance for the ode solver; default = 1e-4
 \item a_tol - the absolute tolerance for the ode solver; default = 1e-4
 \item ode_max_steps - the max number of steps for the solver; default = 1e6
-\item ode_debug - whether to print debug output for ode parameter updates; default = FALSE
 \item enable_heterogeneity - boolean whether to include heterogeneity in biting
 rates; default = TRUE
 }}
diff --git a/man/run_simulation.Rd b/man/run_simulation.Rd
index d25a6526..f991b373 100644
--- a/man/run_simulation.Rd
+++ b/man/run_simulation.Rd
@@ -46,20 +46,22 @@ population)
 \item icm_mean: the mean maternal immunity to clinical infection over the population of humans
 \item ib_mean: the mean blood immunity to all infection over the population of humans
 \item id_mean: the mean immunity from detection through microscopy over the population of humans
-\item pv: prevalence for humans between an inclusive age range (in timesteps). This
-defaults to pv_730_3650. Other prevalence columns can be set with
+\item n: number of humans between an inclusive age range at this timestep. This
+defaults to n_730_3650. Other age ranges can be set with
 prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
-\item pv_severe: prevalence for severe malaria in humans between an inclusive age range (in timesteps).
-These columns can be set with the
+\item n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This
+defaults to n_detect_730_3650. Other age ranges can be set with
+prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
+\item n_severe: number of humans with a severe infection detectable by microscopy
+between an inclusive age range at this timestep. Age ranges can be set with
 severe_prevalence_rendering_min_ages and severe_prevalence_rendering_max_ages parameters.
-\item inc: incidence for humans between an inclusive age range (in timesteps).
+\item n_inc: number of new infections for humans between an inclusive age range at this timestep.
 incidence columns can be set with
 incidence_rendering_min_ages and incidence_rendering_max_ages parameters.
-\item clin_inc: clinical incidence for humans between an inclusive age range (in timesteps)
+\item n_inc_clinical: number of new clinical infections for humans between an inclusive age range at this timestep.
 clinical incidence columns can be set with
 clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters.
-\item inc_severe: severe incidence for humans between an inclusive age range (in
-timesteps).
+\item n_inc_severe: number of new severe infections for humans between an inclusive age range at this timestep.
 severe incidence columns can be set with
 severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
 \item severe_deaths: number of deaths due to severe malaria. severe_enabled must be
diff --git a/man/set_mass_rtss.Rd b/man/set_mass_rtss.Rd
index 8622e912..aa1a479d 100644
--- a/man/set_mass_rtss.Rd
+++ b/man/set_mass_rtss.Rd
@@ -26,7 +26,9 @@ set_mass_rtss(
 
 \item{max_ages}{for the target population, inclusive (in timesteps)}
 
-\item{min_wait}{the minimum acceptable time since the last vaccination (in timesteps);}
+\item{min_wait}{the minimum acceptable time since the last vaccination (in timesteps);
+When using both set_mass_rtss and set_rtss_epi, this represents the minimum
+time between an individual being vaccinated under one scheme and vaccinated under another.}
 
 \item{boosters}{the timesteps (following the initial vaccination) at which booster vaccinations are administered}
 
diff --git a/man/set_rtss_epi.Rd b/man/set_rtss_epi.Rd
index 77b82309..3d759ef0 100644
--- a/man/set_rtss_epi.Rd
+++ b/man/set_rtss_epi.Rd
@@ -28,8 +28,10 @@ set_rtss_epi(
 \item{age}{for the target population, (in timesteps)}
 
 \item{min_wait}{the minimum acceptable time since the last vaccination (in
-timesteps); This includes seasonal boosters and vaccinations between EPI and mass
-strategies}
+timesteps); When seasonal_boosters = TRUE, this represents the minimum time
+between an individual receiving the final dose and the first booster. When using
+both set_mass_rtss and set_rtss_epi, this represents the minimum time between
+an individual being vaccinated under one scheme and vaccinated under another.}
 
 \item{boosters}{the timesteps (following the final dose) at which booster vaccinations are administered}
 

From de5546d53caa1c8c9f4234a7951abf1e79a39c46 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 16 Sep 2021 11:23:27 +0100
Subject: [PATCH 06/23] Disaggregate vector species outputs:

  * Disaggregate EIR and state counts by mosquito species
  * Remove redundant Total_M and and global EIR outputs
---
 NEWS.md               |  6 ++++--
 R/biting_process.R    |  7 ++++---
 R/compartmental.R     | 21 ++++++++++-----------
 R/model.R             | 25 +++++++++++--------------
 R/processes.R         | 10 +++-------
 R/render.R            | 25 +++++++++++++------------
 man/run_simulation.Rd | 25 +++++++++++--------------
 7 files changed, 56 insertions(+), 63 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 726b4e27..aed84d5a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,6 @@
-# malariasimulation 1.2.0
+# malariasimulation 1.2.0 (wip)
 
-  * Added a `NEWS.md` file to track changes to the package.
+  * added a `news.md` file to track changes to the package.
   * n_inc_clinical includes treated humans
+  * disaggregate EIR and state counts by mosquito species
+  * remove redundant Total_M and and global EIR outputs
diff --git a/R/biting_process.R b/R/biting_process.R
index 7da01dc0..87537a41 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -91,6 +91,7 @@ simulate_bites <- function(
   EIR <- 0
 
   for (s_i in seq_along(parameters$species)) {
+    species_name <- parameters$species[[s_i]]
     solver_states <- solver_get_states(solvers[[s_i]])
     p_bitten <- prob_bitten(timestep, variables, s_i, parameters)
     Q0 <- parameters$Q0[[s_i]]
@@ -118,6 +119,7 @@ simulate_bites <- function(
 
     lagged_eir[[s_i]]$save(n_infectious * a, timestep)
     species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
+    renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     n_bites <- rpois(1, species_eir * mean(psi))
     if (n_bites > 0) {
@@ -129,9 +131,9 @@ simulate_bites <- function(
     infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)
     lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
     foim <- calculate_foim(a, infectivity)
-    renderer$render(paste0('FOIM_', s_i), foim, timestep)
+    renderer$render(paste0('FOIM_', species_name), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
-    renderer$render(paste0('mu_', s_i), mu, timestep)
+    renderer$render(paste0('mu_', species_name), mu, timestep)
 
     if (parameters$individual_mosquitoes) {
       # update the ODE with stats for ovoposition calculations
@@ -167,7 +169,6 @@ simulate_bites <- function(
     }
   }
 
-  renderer$render('EIR', EIR, timestep)
   renderer$render('n_bitten', bitten_humans$size(), timestep)
   bitten_humans
 }
diff --git a/R/compartmental.R b/R/compartmental.R
index 79d52d76..1e712b49 100644
--- a/R/compartmental.R
+++ b/R/compartmental.R
@@ -77,7 +77,7 @@ parameterise_solvers <- function(models, parameters) {
   )
 }
 
-create_ode_rendering_process <- function(renderer, solvers, parameters) {
+create_compartmental_rendering_process <- function(renderer, solvers, parameters) {
   if (parameters$individual_mosquitoes) {
     indices <- ODE_INDICES
   } else {
@@ -86,18 +86,17 @@ create_ode_rendering_process <- function(renderer, solvers, parameters) {
 
   function(timestep) {
     counts <- rep(0, length(indices))
-    for (i in seq_along(solvers)) {
-      row <- solver_get_states(solvers[[i]])
+    for (s_i in seq_along(solvers)) {
+      row <- solver_get_states(solvers[[s_i]])
+      for (i in seq_along(indices)) {
+        renderer$render(
+          paste0(names(indices)[[i]], '_', parameters$species[[s_i]], '_count'),
+          row[[i]],
+          timestep
+        )
+      }
       counts <- counts + row
     }
-
-    for (i in seq_along(indices)) {
-      renderer$render(
-        paste0(names(indices)[[i]], '_count'),
-        counts[[i]],
-        timestep
-      )
-    }
   }
 }
 
diff --git a/R/model.R b/R/model.R
index 0d540d9d..6f3405e8 100644
--- a/R/model.R
+++ b/R/model.R
@@ -8,14 +8,10 @@
 #'
 #'  * timestep: the timestep for the row
 #'  * infectivity: the infectivity from humans towards mosquitoes
-#'  * lambda: the effective biting rate on humans (per timestep) (per
-#'species). This defaults to lambda_All
-#'  * normal_lambda: the effective biting rate on adult humans (per
-#' timestep) (per species). This defaults to normal_lambda_All
 #'  * FOIM: the force of infection towards mosquitoes (per species)
 #'  * mu: the death rate of adult mosquitoes (per species)
-#'  * EIR: the Entomological Inoculation Rate (per timestep, over the whole
-#'  population)
+#'  * EIR: the Entomological Inoculation Rate (per timestep, per species, over 
+#'  the whole population)
 #'  * n_bitten: number of humans bitten by an infectious mosquito
 #'  * n_treated: number of humans treated for clinical or severe malaria this timestep
 #'  * n_infections: number of humans who get an asymptomatic, clinical or severe malaria this timestep
@@ -49,14 +45,15 @@
 #' severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
 #'  * severe_deaths: number of deaths due to severe malaria. severe_enabled must be
 #' set to TRUE
-#'  * E_count: number of mosquitoes in the early larval stage
-#'  * L_count: number of mosquitoes in the late larval stage
-#'  * P_count: number of mosquitoes in the pupal stage
-#'  * Sm_count: number of adult female mosquitoes who are Susceptible
-#'  * Pm_count: number of adult female mosquitoes who are incubating
-#'  * Im_count: number of adult female mosquitoes who are infectious
-#'  * total_M: number of adult female mosquitoes. Variables with suffixes total_M_# refer to the number of adult 
-#' female mosquitoes from the different indicies of the set_species() vector
+#'  * E_count: number of mosquitoes in the early larval stage (per species)
+#'  * L_count: number of mosquitoes in the late larval stage (per species)
+#'  * P_count: number of mosquitoes in the pupal stage (per species)
+#'  * Sm_count: number of adult female mosquitoes who are Susceptible (per
+#'  species)
+#'  * Pm_count: number of adult female mosquitoes who are incubating (per
+#'  species)
+#'  * Im_count: number of adult female mosquitoes who are infectious (per
+#'  species)
 #'  * rate_D_A: rate that humans transition from clinical disease to
 #' asymptomatic
 #'  * rate_A_U: rate that humans transition from asymptomatic to
diff --git a/R/processes.R b/R/processes.R
index 6a7621a8..713ced0f 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -164,20 +164,16 @@ create_processes <- function(
       parameters,
       renderer
     ),
-    create_ode_rendering_process(renderer, solvers, parameters)
+    create_compartmental_rendering_process(renderer, solvers, parameters)
   )
 
   if (parameters$individual_mosquitoes) {
     processes <- c(
       processes,
-      individual::categorical_count_renderer_process(
-        renderer,
-        variables$mosquito_state,
-        c('Sm', 'Pm', 'Im')
-      ),
-      create_total_M_renderer_individual(
+      create_vector_count_renderer_individual(
         variables$mosquito_state,
         variables$species,
+        variables$mosquito_state,
         renderer,
         parameters
       )
diff --git a/R/render.R b/R/render.R
index 0757aa47..a957c6e5 100644
--- a/R/render.R
+++ b/R/render.R
@@ -163,36 +163,37 @@ create_variable_mean_renderer_process <- function(
   }
 }
 
-create_total_M_renderer_individual <- function(
+create_vector_count_renderer_individual <- function(
   mosquito_state,
   species,
+  state,
   renderer,
   parameters
   ) {
   function(timestep) {
     adult <- mosquito_state$get_index_of('NonExistent')$not()
     for (i in seq_along(parameters$species)) {
-      renderer$render(
-        paste0('total_M_', i),
-        species$get_index_of(parameters$species[[i]])$and(adult)$size(),
-        timestep
-      )
+      species_name <- parameters$species[[i]]
+      species_index <- species$get_index_of(species_name)
+      for (s in state$get_categories()) {
+        renderer$render(
+          paste0(s, '_', species_name, '_count'),
+          state$get_index_of(s)$and(species_index)$size(),
+          timestep
+        )
+      }
     }
-
-    renderer$render('total_M', adult$size(), timestep)
   }
 }
 
-create_total_M_renderer_compartmental <- function(renderer, solvers) {
+create_total_M_renderer_compartmental <- function(renderer, solvers, parameters) {
   function(timestep) {
     total_M <- 0
     for (i in seq_along(solvers)) {
       row <- solver_get_states(solvers[[i]])
       species_M <- sum(row[ADULT_ODE_INDICES])
       total_M <- total_M + species_M
-      renderer$render(paste0('total_M_', i), species_M, timestep)
+      renderer$render(paste0('total_M_', parameters$species[[i]]), species_M, timestep)
     }
-
-    renderer$render('total_M', total_M, timestep)
   }
 }
diff --git a/man/run_simulation.Rd b/man/run_simulation.Rd
index f991b373..e4bd832e 100644
--- a/man/run_simulation.Rd
+++ b/man/run_simulation.Rd
@@ -25,14 +25,10 @@ The resulting dataframe contains the following columns:
 \itemize{
 \item timestep: the timestep for the row
 \item infectivity: the infectivity from humans towards mosquitoes
-\item lambda: the effective biting rate on humans (per timestep) (per
-species). This defaults to lambda_All
-\item normal_lambda: the effective biting rate on adult humans (per
-timestep) (per species). This defaults to normal_lambda_All
 \item FOIM: the force of infection towards mosquitoes (per species)
 \item mu: the death rate of adult mosquitoes (per species)
-\item EIR: the Entomological Inoculation Rate (per timestep, over the whole
-population)
+\item EIR: the Entomological Inoculation Rate (per timestep, per species, over
+the whole population)
 \item n_bitten: number of humans bitten by an infectious mosquito
 \item n_treated: number of humans treated for clinical or severe malaria this timestep
 \item n_infections: number of humans who get an asymptomatic, clinical or severe malaria this timestep
@@ -66,14 +62,15 @@ severe incidence columns can be set with
 severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
 \item severe_deaths: number of deaths due to severe malaria. severe_enabled must be
 set to TRUE
-\item E_count: number of mosquitoes in the early larval stage
-\item L_count: number of mosquitoes in the late larval stage
-\item P_count: number of mosquitoes in the pupal stage
-\item Sm_count: number of adult female mosquitoes who are Susceptible
-\item Pm_count: number of adult female mosquitoes who are incubating
-\item Im_count: number of adult female mosquitoes who are infectious
-\item total_M: number of adult female mosquitoes. Variables with suffixes total_M_# refer to the number of adult
-female mosquitoes from the different indicies of the set_species() vector
+\item E_count: number of mosquitoes in the early larval stage (per species)
+\item L_count: number of mosquitoes in the late larval stage (per species)
+\item P_count: number of mosquitoes in the pupal stage (per species)
+\item Sm_count: number of adult female mosquitoes who are Susceptible (per
+species)
+\item Pm_count: number of adult female mosquitoes who are incubating (per
+species)
+\item Im_count: number of adult female mosquitoes who are infectious (per
+species)
 \item rate_D_A: rate that humans transition from clinical disease to
 asymptomatic
 \item rate_A_U: rate that humans transition from asymptomatic to

From 26235b71b4d8adcc133364a5eaf42a091d43f0cb Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 23 Sep 2021 09:06:04 +0100
Subject: [PATCH 07/23] Do not bite if EIR is 0

---
 R/biting_process.R | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index 87537a41..fbe58572 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -121,11 +121,14 @@ simulate_bites <- function(
     species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
-    n_bites <- rpois(1, species_eir * mean(psi))
-    if (n_bites > 0) {
-      bitten_humans$insert(
-        fast_weighted_sample(n_bites, lambda)
-      )
+    expected_bites <- species_eir * mean(psi)
+    if (expected_bites > 0) {
+      n_bites <- rpois(1, expected_bites)
+      if (n_bites > 0) {
+        bitten_humans$insert(
+          fast_weighted_sample(n_bites, lambda)
+        )
+      }
     }
 
     infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)

From 17311c0e8ba729136aef835c0cee9489fa30d493 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 23 Sep 2021 09:45:15 +0100
Subject: [PATCH 08/23] Fix compartmental rendering bug

---
 R/processes.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/processes.R b/R/processes.R
index 713ced0f..58b5a3b9 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -183,7 +183,8 @@ create_processes <- function(
       processes,
       create_total_M_renderer_compartmental(
         renderer,
-        solvers
+        solvers,
+        parameters
       )
     )
   }

From 212f0e4e974c7e6476ac51521db0128e37af3ca9 Mon Sep 17 00:00:00 2001
From: Hillary Topazian <80535782+htopazian@users.noreply.github.com>
Date: Thu, 23 Sep 2021 12:28:50 +0100
Subject: [PATCH 09/23] prophylaxis spelling correction

---
 DESCRIPTION    | 2 +-
 R/parameters.R | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 07be5235..a2029297 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -40,7 +40,7 @@ Suggests:
   cowplot,
   ggplot2,
   covr
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.2
 Roxygen: list(markdown = TRUE)
 LinkingTo:
   Rcpp,
diff --git a/R/parameters.R b/R/parameters.R
index d3ed8d6a..847e2b5f 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -346,8 +346,8 @@ get_parameters <- function(overrides = list()) {
     # treatment
     drug_efficacy          = numeric(0),
     drug_rel_c             = numeric(0),
-    drug_prophilaxis_shape = numeric(0),
-    drug_prophilaxis_scale = numeric(0),
+    drug_prophylaxis_shape = numeric(0),
+    drug_prophylaxis_scale = numeric(0),
     clinical_treatment_drugs     = list(),
     clinical_treatment_timesteps = list(),
     clinical_treatment_coverages = list(),

From c486ffff027e7aff255f03135c0d04766c412736 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Tue, 12 Oct 2021 13:43:39 +0100
Subject: [PATCH 10/23] Specify inplace bitset.not operations

---
 R/biting_process.R      | 4 ++--
 R/human_infection.R     | 4 ++--
 R/mda_processes.R       | 2 +-
 R/mortality_processes.R | 2 +-
 R/render.R              | 2 +-
 R/rtss.R                | 4 ++--
 6 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index 87537a41..07b7d93b 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -85,7 +85,7 @@ simulate_bites <- function(
   if (parameters$individual_mosquitoes) {
     infectious_index <- variables$mosquito_state$get_index_of('Im')
     susceptible_index <- variables$mosquito_state$get_index_of('Sm')
-    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not()
+    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not(TRUE)
   }
 
   EIR <- 0
@@ -192,7 +192,7 @@ effective_biting_rates <- function(a, .pi, p_bitten) {
 
 calculate_infectious <- function(species, solvers, variables, parameters) {
   if (parameters$individual_mosquitoes) {
-    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not()
+    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not(TRUE)
     return(
       calculate_infectious_individual(
         species,
diff --git a/R/human_infection.R b/R/human_infection.R
index f5910991..93ef8744 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -301,10 +301,10 @@ schedule_infections <- function(
   parameters,
   asymptomatics
   ) {
-  included <- treated$not()
+  included <- treated$not(TRUE)
 
   to_infect <- clinical_infections$and(included)
-  to_infect_asym <- clinical_infections$not()$and(infections)$and(
+  to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(
     included
   )
 
diff --git a/R/mda_processes.R b/R/mda_processes.R
index 002df03c..3833a152 100644
--- a/R/mda_processes.R
+++ b/R/mda_processes.R
@@ -50,7 +50,7 @@ create_mda_listeners <- function(
       }
 
       # Move everyone else
-      other <- to_move$copy()$and(diseased$not())
+      other <- to_move$copy()$and(diseased$not(TRUE))
       if (other$size() > 0) {
         variables$state$queue_update('S', other)
       }
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index 2a350a94..af1bd9e0 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -16,7 +16,7 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
     ) / 365)
 
     died <- individual::Bitset$new(parameters$human_population)
-    died <- sample_bitset(died$not(), 1 / parameters$average_age)
+    died <- sample_bitset(died$not(TRUE), 1 / parameters$average_age)
 
     renderer$render('natural_deaths', died$size(), timestep)
 
diff --git a/R/render.R b/R/render.R
index a957c6e5..5d0f5b0f 100644
--- a/R/render.R
+++ b/R/render.R
@@ -171,7 +171,7 @@ create_vector_count_renderer_individual <- function(
   parameters
   ) {
   function(timestep) {
-    adult <- mosquito_state$get_index_of('NonExistent')$not()
+    adult <- mosquito_state$get_index_of('NonExistent')$not(TRUE)
     for (i in seq_along(parameters$species)) {
       species_name <- parameters$species[[i]]
       species_index <- species$get_index_of(species_name)
diff --git a/R/rtss.R b/R/rtss.R
index 638deb7a..6d603777 100644
--- a/R/rtss.R
+++ b/R/rtss.R
@@ -28,7 +28,7 @@ create_rtss_epi_process <- function(
       not_recently_vaccinated <- variables$rtss_vaccinated$get_index_of(
         a = max(timestep - parameters$rtss_epi_min_wait, 0),
         b = timestep
-      )$not()
+      )$not(TRUE)
       target <- to_vaccinate$and(not_recently_vaccinated)$to_vector()
     }
 
@@ -78,7 +78,7 @@ create_rtss_mass_listener <- function(
       not_recently_vaccinated <- variables$rtss_vaccinated$get_index_of(
         a = max(timestep - parameters$rtss_mass_min_wait, 0),
         b = timestep
-      )$not()
+      )$not(TRUE)
       target <- in_age_group$and(not_recently_vaccinated)$to_vector()
     }
     

From 7c1867d2ed1faafd11e153ed57b991ca5f476b97 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Tue, 12 Oct 2021 16:17:45 +0100
Subject: [PATCH 11/23] Fix regressions and peg individual version:

 * mock classes copy lazy bitsets
 * individual >= 0.1.6
 * prob_bitten uses inplace not
---
 DESCRIPTION                         |  5 +++--
 R/vector_control.R                  |  3 +--
 tests/testthat/helper-integration.R | 18 +++++++++++++++++-
 tests/testthat/test-mda.R           |  1 -
 4 files changed, 21 insertions(+), 6 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 07be5235..664a9a4d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -19,9 +19,10 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 Remotes:
-  mrc-ide/malariaEquilibrium
+  mrc-ide/malariaEquilibrium,
+  mrc-ide/individual
 Imports:
-  individual,
+  individual (>= 0.1.6),
   malariaEquilibrium,
   Rcpp,
   statmod,
diff --git a/R/vector_control.R b/R/vector_control.R
index dc39dd1c..f197a7d1 100644
--- a/R/vector_control.R
+++ b/R/vector_control.R
@@ -39,8 +39,7 @@ prob_bitten <- function(
 
   if (parameters$spraying) {
     phi_indoors <- parameters$phi_indoors[[species]]
-    unprotected <- variables$spray_time$get_index_of(set=-1)
-    protected <- unprotected$not()
+    protected <- variables$spray_time$get_index_of(set=-1)$not(TRUE)
     spray_time <- variables$spray_time$get_values(protected)
     matches <- match(spray_time, parameters$spraying_timesteps)
     ls_theta <- parameters$spraying_ls_theta[matches, species]
diff --git a/tests/testthat/helper-integration.R b/tests/testthat/helper-integration.R
index 8a3d6f4e..b69b7dcb 100644
--- a/tests/testthat/helper-integration.R
+++ b/tests/testthat/helper-integration.R
@@ -23,7 +23,23 @@ mock_method <- function(class, method, mock) {
     paste0('Mock', class$classname),
     inherit = class
   )
-  MockClass$set('public', method, function(...) mock(...))
+  MockClass$set(
+    'public',
+    method,
+    function(...) {
+      args <- lapply(
+        list(...),
+        function(arg) {
+          # copy a any bitsets because these will be lazily evaluated
+          if (inherits(arg, 'Bitset')) {
+            return(arg$copy())
+          }
+          arg
+        }
+      )
+      do.call(mock, args)
+    }
+  )
   MockClass$set('public', paste0(method, '_mock'), function() mock)
   MockClass
 }
diff --git a/tests/testthat/test-mda.R b/tests/testthat/test-mda.R
index 689cad3c..d9d65ada 100644
--- a/tests/testthat/test-mda.R
+++ b/tests/testthat/test-mda.R
@@ -50,7 +50,6 @@ test_that('MDA moves the diseased and non-diseased population correctly', {
     mockery::mock_args(mock_correlation)[[1]][[1]],
     c(3, 4)
   )
-
   expect_bitset_update(
     variables$state$queue_update_mock(),
     'Tr',

From df6536a16b5a76fc2a5745dff2818566f431c859 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Tue, 28 Sep 2021 11:35:47 +0100
Subject: [PATCH 12/23] Quick and dirty patch for resetting severe cases

---
 R/events.R              | 10 ++++++++++
 R/mortality_processes.R |  2 +-
 2 files changed, 11 insertions(+), 1 deletion(-)

diff --git a/R/events.R b/R/events.R
index 908a00e0..6afa03b8 100644
--- a/R/events.R
+++ b/R/events.R
@@ -132,6 +132,11 @@ attach_event_listeners <- function(
       parameters
     )
   )
+  events$asymptomatic_progression$add_listener(
+    function(timestep, target) {
+      variables$is_severe$queue_update('no', target)
+    }
+  )
 
   # Recovery events
   events$subpatent_progression$add_listener(
@@ -157,6 +162,11 @@ attach_event_listeners <- function(
   events$recovery$add_listener(
     create_rate_listener('U', 'S', renderer)
   )
+  events$recovery$add_listener(
+    function(timestep, target) {
+      variables$is_severe$queue_update('no', target)
+    }
+  )
 
   # ===========
   # Progression
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index af1bd9e0..2bcb1253 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -122,5 +122,5 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
 died_from_severe <- function(severe, diseased, v, treated, fvt) {
   at_risk <- severe$copy()$and(diseased)
   treated_at_risk <- severe$copy()$and(treated)
-  sample_bitset(at_risk, v)$or(sample_bitset(treated_at_risk, fvt))
+  sample_bitset(at_risk$or(sample_bitset(treated_at_risk, fvt)), v)
 }

From 516a684ddcb2a1a7e0e528a10c40d586008af253 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Tue, 28 Sep 2021 17:55:52 +0100
Subject: [PATCH 13/23] Generalised incidence rendering

---
 R/human_infection.R | 31 +++++++++++++++-------
 R/render.R          | 63 +++++++++------------------------------------
 2 files changed, 34 insertions(+), 60 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index 93ef8744..24d62b17 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -38,10 +38,11 @@ simulate_infection <- function(
 
   incidence_renderer(
     variables$birth,
-    variables$is_severe,
-    parameters,
     renderer,
     infected_humans,
+    'n_inc_',
+    parameters$incidence_rendering_min_ages,
+    parameters$incidence_rendering_max_ages,
     timestep
   )
 
@@ -68,11 +69,13 @@ simulate_infection <- function(
     parameters
   )
 
-  clinical_incidence_renderer(
+  incidence_renderer(
     variables$birth,
-    parameters,
     renderer,
     clinical_infections,
+    'n_inc_clinical_',
+    parameters$clinical_incidence_rendering_min_ages,
+    parameters$clinical_incidence_rendering_max_ages,
     timestep
   )
 
@@ -82,7 +85,8 @@ simulate_infection <- function(
       clinical_infections,
       variables,
       infected_humans,
-      parameters
+      parameters,
+      renderer
     )
   }
 
@@ -193,13 +197,15 @@ calculate_clinical_infections <- function(variables, infections, parameters) {
 #' @param variables a list of all of the model variables
 #' @param infections indices of all infected humans (for immunity boosting)
 #' @param parameters model parameters
+#' @param renderer model outputs
 #' @noRd
 update_severe_disease <- function(
   timestep,
   clinical_infections,
   variables,
   infections,
-  parameters
+  parameters,
+  renderer
   ) {
   if (clinical_infections$size() > 0) {
     age <- get_age(variables$birth$get_values(clinical_infections), timestep)
@@ -211,9 +217,16 @@ update_severe_disease <- function(
       parameters
     )
     develop_severe <- bernoulli_multi_p(theta)
-    variables$is_severe$queue_update(
-      'yes',
-      bitset_at(clinical_infections, develop_severe)
+    severe_infections <- bitset_at(clinical_infections, develop_severe)
+    variables$is_severe$queue_update('yes', severe_infections)
+    incidence_renderer(
+      variables$birth,
+      renderer,
+      severe_infections,
+      'n_inc_severe_',
+      parameters$severe_incidence_rendering_min_ages,
+      parameters$severe_incidence_rendering_max_ages,
+      timestep
     )
     boost_immunity(
       variables$iva,
diff --git a/R/render.R b/R/render.R
index 5d0f5b0f..47790255 100644
--- a/R/render.R
+++ b/R/render.R
@@ -78,69 +78,30 @@ create_prevelance_renderer <- function(
 #' detected by microscopy and with severe malaria
 #' 
 #' @param birth variable for birth of the individual
-#' @param is_severe variable for if individual has severe malaria
-#' @param parameters model parameters
-#' @param renderer model renderer
+#' @param renderer object for model outputs
+#' @param target incidence population
+#' @param prefix for model outputs
+#' @param lowers age bounds
+#' @param uppers age bounds
 #' @param timestep current target
-#' @param target newly infected population
 #' 
 #' @noRd
 incidence_renderer <- function(
   birth,
-  is_severe,
-  parameters,
-  renderer,
-  target,
-  timestep
-  ) {
-  severe <- is_severe$get_index_of('yes')$and(target)
-  for (i in seq_along(parameters$incidence_rendering_min_ages)) {
-    lower <- parameters$incidence_rendering_min_ages[[i]]
-    upper <- parameters$incidence_rendering_max_ages[[i]]
-    in_age <- in_age_range(birth, timestep, lower, upper)
-    renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
-    renderer$render(
-      paste0('n_inc_', lower, '_', upper),
-      in_age$and(target)$size(),
-      timestep
-    )
-  }
-  for (i in seq_along(parameters$severe_incidence_rendering_min_ages)) {
-    lower <- parameters$severe_incidence_rendering_min_ages[[i]]
-    upper <- parameters$severe_incidence_rendering_max_ages[[i]]
-    in_age <- in_age_range(birth, timestep, lower, upper)
-    renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
-    renderer$render(
-      paste0('n_inc_severe_', lower, '_', upper),
-      in_age$and(target)$and(severe)$size(),
-      timestep
-    )
-  }
-}
-
-#' @title Render clinical incidence statistics
-#' 
-#' @description renders clinical incidence (new for this timestep)
-#' 
-#' @param birth variable for birth of the individual
-#' @param parameters model parameters
-#' @param renderer model renderer
-#' 
-#' @noRd
-clinical_incidence_renderer <- function(
-  birth,
-  parameters,
   renderer,
   target,
+  prefix,
+  lowers,
+  uppers,
   timestep
   ) {
-  for (i in seq_along(parameters$clinical_incidence_rendering_min_ages)) {
-    lower <- parameters$clinical_incidence_rendering_min_ages[[i]]
-    upper <- parameters$clinical_incidence_rendering_max_ages[[i]]
+  for (i in seq_along(lowers)) {
+    lower <- lowers[[i]]
+    upper <- uppers[[i]]
     in_age <- in_age_range(birth, timestep, lower, upper)
     renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
     renderer$render(
-      paste0('n_inc_clinical_', lower, '_', upper),
+      paste0(prefix, lower, '_', upper),
       in_age$and(target)$size(),
       timestep
     )

From cb43c1f223a78ddbfe74eb9e36e2951c0bbaa811 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Thu, 7 Oct 2021 11:21:09 +0100
Subject: [PATCH 14/23] Fix tests and docs for new incidence

---
 R/render.R                                  |   1 -
 tests/testthat/test-infection-integration.R |   3 +-
 tests/testthat/test-render.R                | 139 +-------------------
 3 files changed, 6 insertions(+), 137 deletions(-)

diff --git a/R/render.R b/R/render.R
index 47790255..6796b6ee 100644
--- a/R/render.R
+++ b/R/render.R
@@ -75,7 +75,6 @@ create_prevelance_renderer <- function(
 #' @title Render incidence statistics
 #' 
 #' @description renders incidence (new for this timestep) for indivduals
-#' detected by microscopy and with severe malaria
 #' 
 #' @param birth variable for birth of the individual
 #' @param renderer object for model outputs
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index df9f6231..161ec357 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -84,7 +84,8 @@ test_that('simulate_infection integrates different types of infection and schedu
     clinical,
     variables,
     infected,
-    parameters
+    parameters,
+    renderer
   )
 
   mockery::expect_args(
diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R
index 94ed32f7..7545cba1 100644
--- a/tests/testthat/test-render.R
+++ b/tests/testthat/test-render.R
@@ -150,23 +150,18 @@ test_that('that severe rendering works', {
 test_that('that clinical incidence rendering works', {
   timestep <- 0
   year <- 365
-  parameters <- get_parameters(list(
-    clinical_incidence_rendering_min_ages = c(0, 2) * year,
-    clinical_incidence_rendering_max_ages = c(5, 10) * year,
-    prevalence_rendering_min_ages = NULL,
-    prevalence_rendering_max_ages = NULL
-  ))
-
   birth <- individual::IntegerVariable$new(
     -c(2, 5, 10, 11) * year
   )
 
   renderer <- mock_render(1)
-  clinical_incidence_renderer(
+  incidence_renderer(
     birth,
-    parameters,
     renderer,
     individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    'n_inc_clinical_',
+    c(0, 2) * year,
+    c(5, 10) * year,
     timestep
   )
 
@@ -202,129 +197,3 @@ test_that('that clinical incidence rendering works', {
     timestep
   )
 })
-
-
-test_that('that incidence rendering works', {
-  timestep <- 0
-  year <- 365
-  parameters <- get_parameters(list(
-    incidence_rendering_min_ages = c(0, 2) * year,
-    incidence_rendering_max_ages = c(5, 10) * year,
-    prevalence_rendering_min_ages = NULL,
-    prevalence_rendering_max_ages = NULL
-  ))
-
-  birth <- individual::IntegerVariable$new(
-    -c(2, 5, 10, 11) * 365
-  )
-
-  is_severe <- individual::CategoricalVariable$new(
-    c('yes', 'no'),
-    c('no', 'yes', 'no', 'no')
-  )
-
-  renderer <- mock_render(1)
-  process <- incidence_renderer(
-    birth,
-    is_severe,
-    parameters,
-    renderer,
-    individual::Bitset$new(4)$insert(c(1, 2, 4)),
-    timestep
-  )
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    1,
-    'n_0_1825',
-    2,
-    timestep
-  )
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    2,
-    'n_inc_0_1825',
-    2,
-    timestep
-  )
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    3,
-    'n_730_3650',
-    3,
-    timestep
-  )
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    4,
-    'n_inc_730_3650',
-    2,
-    timestep
-  )
-})
-
-test_that('that severe incidence rendering works', {
-  timestep <- 0
-  year <- 365
-  parameters <- get_parameters(list(
-    severe_incidence_rendering_min_ages = c(0, 2) * year,
-    severe_incidence_rendering_max_ages = c(5, 10) * year,
-    prevalence_rendering_min_ages = NULL,
-    prevalence_rendering_max_ages = NULL
-  ))
-
-  birth <- individual::IntegerVariable$new(
-    -c(2, 6, 10, 11) * 365
-  )
-
-  is_severe <- individual::CategoricalVariable$new(
-    c('yes', 'no'),
-    c('no', 'yes', 'no', 'no')
-  )
-
-  renderer <- mock_render(1)
-  process <- incidence_renderer(
-    birth,
-    is_severe,
-    parameters,
-    renderer,
-    individual::Bitset$new(4)$insert(c(1, 2, 4)),
-    timestep
-  )
-
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    1,
-    'n_0_1825',
-    1,
-    timestep
-  )
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    2,
-    'n_inc_severe_0_1825',
-    0,
-    timestep
-  )
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    3,
-    'n_730_3650',
-    3,
-    timestep
-  )
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    4,
-    'n_inc_severe_730_3650',
-    1,
-    timestep
-  )
-})

From 09f4dfa055dd9942c74f934a7fd828f737338fee Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Mon, 11 Oct 2021 16:32:52 +0100
Subject: [PATCH 15/23] Add optimised versions of p_* outputs

---
 R/RcppExports.R                             |  4 ++
 R/human_infection.R                         | 69 ++++++++++++++-------
 R/model.R                                   | 17 ++++-
 R/render.R                                  | 43 ++++++++-----
 R/utils.R                                   |  9 +++
 man/run_simulation.Rd                       | 17 ++++-
 src/RcppExports.cpp                         | 13 ++++
 src/test-helper.cpp                         | 12 ----
 src/test-helper.h                           | 15 -----
 src/test-mock.h                             | 38 ------------
 src/utils.cpp                               | 17 +++++
 tests/testthat/test-infection-integration.R | 14 ++++-
 tests/testthat/test-render.R                | 39 ++++++++++--
 tests/testthat/test-rtss.R                  |  1 +
 14 files changed, 193 insertions(+), 115 deletions(-)
 delete mode 100644 src/test-helper.cpp
 delete mode 100644 src/test-helper.h
 delete mode 100644 src/test-mock.h

diff --git a/R/RcppExports.R b/R/RcppExports.R
index 249955b9..ccf76a71 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -65,6 +65,10 @@ bernoulli_multi_p_cpp <- function(p) {
     .Call(`_malariasimulation_bernoulli_multi_p_cpp`, p)
 }
 
+bitset_index_cpp <- function(a, b) {
+    .Call(`_malariasimulation_bitset_index_cpp`, a, b)
+}
+
 fast_weighted_sample <- function(size, probs) {
     .Call(`_malariasimulation_fast_weighted_sample`, size, probs)
 }
diff --git a/R/human_infection.R b/R/human_infection.R
index 24d62b17..c0b246c3 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -33,16 +33,7 @@ simulate_infection <- function(
     variables,
     bitten_humans,
     parameters,
-    timestep
-  )
-
-  incidence_renderer(
-    variables$birth,
     renderer,
-    infected_humans,
-    'n_inc_',
-    parameters$incidence_rendering_min_ages,
-    parameters$incidence_rendering_max_ages,
     timestep
   )
 
@@ -66,16 +57,8 @@ simulate_infection <- function(
   clinical_infections <- calculate_clinical_infections(
     variables,
     infected_humans,
-    parameters
-  )
-
-  incidence_renderer(
-    variables$birth,
+    parameters,
     renderer,
-    clinical_infections,
-    'n_inc_clinical_',
-    parameters$clinical_incidence_rendering_min_ages,
-    parameters$clinical_incidence_rendering_max_ages,
     timestep
   )
 
@@ -117,12 +100,14 @@ simulate_infection <- function(
 #' @param variables a list of all of the model variables
 #' @param bitten_humans bitset of bitten humans
 #' @param parameters model parameters
+#' @param renderer model render object
 #' @param timestep current timestep
 #' @noRd
 calculate_infections <- function(
   variables,
   bitten_humans,
   parameters,
+  renderer,
   timestep
   ) {
   source_humans <- variables$state$get_index_of(
@@ -169,10 +154,22 @@ calculate_infections <- function(
     )
   }
 
-  bitset_at(
+  prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+  infected <- bitset_at(source_humans, bernoulli_multi_p(prob))
+
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    infected,
     source_humans,
-    bernoulli_multi_p(b * (1 - prophylaxis) * (1 - vaccine_efficacy))
+    prob,
+    'inc_',
+    parameters$incidence_rendering_min_ages,
+    parameters$incidence_rendering_max_ages,
+    timestep
   )
+
+  infected
 }
 
 #' @title Calculate clinical infections
@@ -181,12 +178,32 @@ calculate_infections <- function(
 #' @param variables a list of all of the model variables
 #' @param infections bitset of infected humans
 #' @param parameters model parameters
+#' @param renderer model render
+#' @param timestep current timestep
 #' @noRd
-calculate_clinical_infections <- function(variables, infections, parameters) {
+calculate_clinical_infections <- function(
+  variables,
+  infections,
+  parameters,
+  renderer,
+  timestep
+  ) {
   ica <- variables$ica$get_values(infections)
   icm <- variables$icm$get_values(infections)
   phi <- clinical_immunity(ica, icm, parameters)
-  bitset_at(infections, bernoulli_multi_p(phi))
+  clinical_infections <- bitset_at(infections, bernoulli_multi_p(phi))
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    clinical_infections,
+    infections,
+    phi,
+    'inc_clinical_',
+    parameters$clinical_incidence_rendering_min_ages,
+    parameters$clinical_incidence_rendering_max_ages,
+    timestep
+  )
+  clinical_infections
 }
 
 #' @title Calculate severe infections
@@ -219,11 +236,17 @@ update_severe_disease <- function(
     develop_severe <- bernoulli_multi_p(theta)
     severe_infections <- bitset_at(clinical_infections, develop_severe)
     variables$is_severe$queue_update('yes', severe_infections)
+    variables$is_severe$queue_update(
+      'no',
+      severe_infections$copy()$not(TRUE)$and(infections)
+    ) # Prevent any unwanted severe reinfections
     incidence_renderer(
       variables$birth,
       renderer,
       severe_infections,
-      'n_inc_severe_',
+      clinical_infections,
+      theta,
+      'inc_severe_',
       parameters$severe_incidence_rendering_min_ages,
       parameters$severe_incidence_rendering_max_ages,
       timestep
diff --git a/R/model.R b/R/model.R
index 6f3405e8..2b636991 100644
--- a/R/model.R
+++ b/R/model.R
@@ -31,18 +31,31 @@
 #'  * n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This
 #' defaults to n_detect_730_3650. Other age ranges can be set with
 #' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
-#'  * n_severe: number of humans with a severe infection detectable by microscopy 
-#'  between an inclusive age range at this timestep. Age ranges can be set with
+#'  * p_detect: the sum of probabilities of detection by microscopy between an
+#' inclusive age range at this timestep. This
+#' defaults to p_detect_730_3650. Other age ranges can be set with
+#' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
+#'  * n_severe: number of humans with a severe infection between an inclusive
+#' age range at this timestep. Age ranges can be set with
 #' severe_prevalence_rendering_min_ages and severe_prevalence_rendering_max_ages parameters.
 #'  * n_inc: number of new infections for humans between an inclusive age range at this timestep.
 #' incidence columns can be set with
 #' incidence_rendering_min_ages and incidence_rendering_max_ages parameters.
+#'  * p_inc: sum of probabilities of infection for humans between an inclusive age range at this timestep.
+#' incidence columns can be set with
+#' incidence_rendering_min_ages and incidence_rendering_max_ages parameters.
 #'  * n_inc_clinical: number of new clinical infections for humans between an inclusive age range at this timestep. 
 #' clinical incidence columns can be set with
 #' clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters.
+#'  * p_inc_clinical: sub of probabilities of clinical infection for humans between an inclusive age range at this timestep. 
+#' clinical incidence columns can be set with
+#' clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters.
 #'  * n_inc_severe: number of new severe infections for humans between an inclusive age range at this timestep.
 #' severe incidence columns can be set with
 #' severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
+#'  * p_inc_severe: the sum of probabilities of severe infection for humans between an inclusive age range at this timestep.
+#' severe incidence columns can be set with
+#' severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
 #'  * severe_deaths: number of deaths due to severe malaria. severe_enabled must be
 #' set to TRUE
 #'  * E_count: number of mosquitoes in the early larval stage (per species)
diff --git a/R/render.R b/R/render.R
index 6796b6ee..90755891 100644
--- a/R/render.R
+++ b/R/render.R
@@ -24,19 +24,15 @@ create_prevelance_renderer <- function(
   renderer
   ) {
   function(timestep) {
-    age <- get_age(birth$get_values(), timestep)
     asymptomatic <- state$get_index_of('A')
-    asymptomatic_vector <- asymptomatic$to_vector()
-    asymptomatic_detected <- bernoulli_multi_p(
-      probability_of_detection(
-        age[asymptomatic_vector],
-        immunity$get_values(asymptomatic),
-        parameters
-      )
-    )
-    detected <- state$get_index_of(c('Tr', 'D'))$or(
-      bitset_at(asymptomatic, asymptomatic_detected)
+    prob <- probability_of_detection(
+      get_age(birth$get_values(asymptomatic), timestep),
+      immunity$get_values(asymptomatic),
+      parameters
     )
+    asymptomatic_detected <- bitset_at(asymptomatic, bernoulli_multi_p(prob))
+    clinically_detected <- state$get_index_of(c('Tr', 'D'))
+    detected <- clinically_detected$copy()$or(asymptomatic_detected)
 
     severe <- is_severe$get_index_of('yes')$and(detected)
     for (i in seq_along(parameters$prevalence_rendering_min_ages)) {
@@ -50,7 +46,14 @@ create_prevelance_renderer <- function(
       ) 
       renderer$render(
         paste0('n_detect_', lower, '_', upper),
-        in_age$and(detected)$size(),
+        in_age$copy()$and(detected)$size(),
+        timestep
+      )
+      renderer$render(
+        paste0('p_detect_', lower, '_', upper),
+        in_age$copy()$and(clinically_detected)$size() + sum(
+          prob[bitset_index(asymptomatic, in_age)]
+        ),
         timestep
       )
     }
@@ -62,7 +65,7 @@ create_prevelance_renderer <- function(
         paste0('n_', lower, '_', upper),
         in_age$size(),
         timestep
-      ) 
+      )
       renderer$render(
         paste0('n_severe_', lower, '_', upper),
         in_age$and(severe)$size(),
@@ -79,6 +82,8 @@ create_prevelance_renderer <- function(
 #' @param birth variable for birth of the individual
 #' @param renderer object for model outputs
 #' @param target incidence population
+#' @param source_pop the population which is sampled for infection
+#' @param prob probability of infection
 #' @param prefix for model outputs
 #' @param lowers age bounds
 #' @param uppers age bounds
@@ -89,6 +94,8 @@ incidence_renderer <- function(
   birth,
   renderer,
   target,
+  source_pop,
+  prob,
   prefix,
   lowers,
   uppers,
@@ -100,8 +107,14 @@ incidence_renderer <- function(
     in_age <- in_age_range(birth, timestep, lower, upper)
     renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
     renderer$render(
-      paste0(prefix, lower, '_', upper),
-      in_age$and(target)$size(),
+      paste0('n_', prefix, lower, '_', upper),
+      in_age$copy()$and(target)$size(),
+      timestep
+    )
+
+    renderer$render(
+      paste0('p_', prefix, lower, '_', upper),
+      sum(prob[bitset_index(source_pop, in_age)]),
       timestep
     )
   }
diff --git a/R/utils.R b/R/utils.R
index f9d046f9..b946c7ec 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -15,6 +15,15 @@ bitset_at <- function(b, i) {
 
 bernoulli_multi_p <- function(p) bernoulli_multi_p_cpp(p)
 
+
+#' @title find the indices of a where it intersects with b
+#' @description synonymous with \code{which(a$to_vector() %in%
+#' b$to_vector())} but faster
+#' @param a the bitset to index
+#' @param b the bitset to check
+#' @noRd
+bitset_index <- function(a, b) bitset_index_cpp(a$.bitset, b$.bitset)
+
 #' @importFrom stats runif
 log_uniform <- function(size, rate) -rate * log(runif(size))
 
diff --git a/man/run_simulation.Rd b/man/run_simulation.Rd
index e4bd832e..5dd1c8ec 100644
--- a/man/run_simulation.Rd
+++ b/man/run_simulation.Rd
@@ -48,18 +48,31 @@ prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
 \item n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This
 defaults to n_detect_730_3650. Other age ranges can be set with
 prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
-\item n_severe: number of humans with a severe infection detectable by microscopy
-between an inclusive age range at this timestep. Age ranges can be set with
+\item p_detect: the sum of probabilities of detection by microscopy between an
+inclusive age range at this timestep. This
+defaults to p_detect_730_3650. Other age ranges can be set with
+prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
+\item n_severe: number of humans with a severe infection between an inclusive
+age range at this timestep. Age ranges can be set with
 severe_prevalence_rendering_min_ages and severe_prevalence_rendering_max_ages parameters.
 \item n_inc: number of new infections for humans between an inclusive age range at this timestep.
 incidence columns can be set with
 incidence_rendering_min_ages and incidence_rendering_max_ages parameters.
+\item p_inc: sum of probabilities of infection for humans between an inclusive age range at this timestep.
+incidence columns can be set with
+incidence_rendering_min_ages and incidence_rendering_max_ages parameters.
 \item n_inc_clinical: number of new clinical infections for humans between an inclusive age range at this timestep.
 clinical incidence columns can be set with
 clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters.
+\item p_inc_clinical: sub of probabilities of clinical infection for humans between an inclusive age range at this timestep.
+clinical incidence columns can be set with
+clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters.
 \item n_inc_severe: number of new severe infections for humans between an inclusive age range at this timestep.
 severe incidence columns can be set with
 severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
+\item p_inc_severe: the sum of probabilities of severe infection for humans between an inclusive age range at this timestep.
+severe incidence columns can be set with
+severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
 \item severe_deaths: number of deaths due to severe malaria. severe_enabled must be
 set to TRUE
 \item E_count: number of mosquitoes in the early larval stage (per species)
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 0d778d51..a73a67bf 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -232,6 +232,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// bitset_index_cpp
+std::vector<size_t> bitset_index_cpp(Rcpp::XPtr<individual_index_t> a, Rcpp::XPtr<individual_index_t> b);
+RcppExport SEXP _malariasimulation_bitset_index_cpp(SEXP aSEXP, SEXP bSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Rcpp::XPtr<individual_index_t> >::type a(aSEXP);
+    Rcpp::traits::input_parameter< Rcpp::XPtr<individual_index_t> >::type b(bSEXP);
+    rcpp_result_gen = Rcpp::wrap(bitset_index_cpp(a, b));
+    return rcpp_result_gen;
+END_RCPP
+}
 // fast_weighted_sample
 Rcpp::IntegerVector fast_weighted_sample(size_t size, std::vector<double> probs);
 RcppExport SEXP _malariasimulation_fast_weighted_sample(SEXP sizeSEXP, SEXP probsSEXP) {
@@ -263,6 +275,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_malariasimulation_solver_step", (DL_FUNC) &_malariasimulation_solver_step, 1},
     {"_malariasimulation_random_seed", (DL_FUNC) &_malariasimulation_random_seed, 1},
     {"_malariasimulation_bernoulli_multi_p_cpp", (DL_FUNC) &_malariasimulation_bernoulli_multi_p_cpp, 1},
+    {"_malariasimulation_bitset_index_cpp", (DL_FUNC) &_malariasimulation_bitset_index_cpp, 2},
     {"_malariasimulation_fast_weighted_sample", (DL_FUNC) &_malariasimulation_fast_weighted_sample, 2},
     {"run_testthat_tests", (DL_FUNC) &run_testthat_tests, 0},
     {NULL, NULL, 0}
diff --git a/src/test-helper.cpp b/src/test-helper.cpp
deleted file mode 100644
index 898ae2bb..00000000
--- a/src/test-helper.cpp
+++ /dev/null
@@ -1,12 +0,0 @@
-/*
- * test-helper.cpp
- *
- *  Created on: 20 Aug 2020
- *      Author: gc1610
- */
-
-#include "test-helper.h"
-
-individual_index_t individual_index(size_t size, std::vector<size_t> i) {
-    return individual_index_t(size, i.cbegin(), i.cend());
-}
diff --git a/src/test-helper.h b/src/test-helper.h
deleted file mode 100644
index b1065d8a..00000000
--- a/src/test-helper.h
+++ /dev/null
@@ -1,15 +0,0 @@
-/*
- * test-helper.h
- *
- *  Created on: 20 Aug 2020
- *      Author: gc1610
- */
-
-#ifndef SRC_TEST_HELPER_H_
-#define SRC_TEST_HELPER_H_
-
-#include <individual.h>
-
-individual_index_t individual_index(size_t size, std::vector<size_t> i);
-
-#endif /* SRC_TEST_HELPER_H_ */
diff --git a/src/test-mock.h b/src/test-mock.h
deleted file mode 100644
index dcc5ffad..00000000
--- a/src/test-mock.h
+++ /dev/null
@@ -1,38 +0,0 @@
-// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
-/*
- * test-emergence.h
- *
- *  Created on: 5 Aug 2020
- *      Author: gc1610
- */
-
-#ifndef SRC_TEST_MOCK_H_
-#define SRC_TEST_MOCK_H_
-
-#include "solver.h"
-#include <testthat/trompeloeil.hpp>
-#include <individual.h>
-#include "Random.h"
-
-struct MockCategory : public CategoricalVariable {
-    MockCategory(
-        const std::vector<std::string> categories,
-        const std::vector<std::string> values)
-        : CategoricalVariable(categories, values) {}
-    MAKE_MOCK2(queue_update, void(const std::string, const individual_index_t&), override);
-};
-
-class MockRandom : public RandomInterface {
-public:
-    MAKE_MOCK2(bernoulli, std::vector<size_t>(size_t, double), override);
-    MAKE_MOCK3(sample, std::vector<size_t>(size_t, size_t, bool), override);
-};
-
-
-integration_function_t mock_integration = [](
-    const state_t&,
-    state_t&,
-    double t
-) {};
-
-#endif /* SRC_TEST_MOCK_H_ */
diff --git a/src/utils.cpp b/src/utils.cpp
index 585a75b7..838a9530 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -1,6 +1,7 @@
 
 #include <Rcpp.h>
 #include "Random.h"
+#include <individual.h>
 
 //[[Rcpp::export]]
 void random_seed(size_t seed) {
@@ -16,6 +17,22 @@ std::vector<size_t> bernoulli_multi_p_cpp(const std::vector<double> p) {
     return values;
 }
 
+//[[Rcpp::export]]
+std::vector<size_t> bitset_index_cpp(
+    Rcpp::XPtr<individual_index_t> a,
+    Rcpp::XPtr<individual_index_t> b
+    ) {
+    auto values = std::vector<size_t>();
+    auto i = 1u;
+    for (const auto& v : *a) {
+        if (b->find(v) != b->cend()) {
+            values.push_back(i);
+        }
+        ++i;
+    }
+    return values;
+}
+
 //[[Rcpp::export(rng = false)]]
 Rcpp::IntegerVector fast_weighted_sample(
     size_t size,
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 161ec357..2dbd2ad5 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -66,6 +66,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     variables,
     bitten,
     parameters,
+    renderer,
     timestep
   )
 
@@ -74,7 +75,9 @@ test_that('simulate_infection integrates different types of infection and schedu
     1,
     variables,
     infected,
-    parameters
+    parameters,
+    renderer,
+    timestep
   )
 
   mockery::expect_args(
@@ -149,6 +152,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     variables,
     bitten_humans, 
     parameters,
+    mock_render(timestep),
     timestep
   )
 
@@ -361,7 +365,13 @@ test_that('prophylaxis is considered for medicated humans', {
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
 
-  calculate_infections(variables, bitten, parameters, timestep)
+  calculate_infections(
+    variables,
+    bitten,
+    parameters,
+    mock_render(timestep),
+    timestep
+  )
 
   expect_equal(
     mockery::mock_args(m)[[1]][[1]],
diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R
index 7545cba1..a1e68f43 100644
--- a/tests/testthat/test-render.R
+++ b/tests/testthat/test-render.R
@@ -24,6 +24,7 @@ test_that('that default rendering works', {
     renderer
   )
 
+  mockery::stub(process, 'probability_of_detection', mockery::mock(.5))
   mockery::stub(process, 'bernoulli_multi_p', mockery::mock(1))
   process(timestep)
 
@@ -43,6 +44,14 @@ test_that('that default rendering works', {
     timestep
   )
 
+  mockery::expect_args(
+    renderer$render_mock(),
+    3,
+    'p_detect_730_3650',
+    1.5,
+    timestep
+  )
+
 })
 
 test_that('that default rendering works when no one is in the age range', {
@@ -159,13 +168,15 @@ test_that('that clinical incidence rendering works', {
     birth,
     renderer,
     individual::Bitset$new(4)$insert(c(1, 2, 4)),
-    'n_inc_clinical_',
+    individual::Bitset$new(4)$insert(seq(4)),
+    c(.1, .2, .3, .4),
+    'inc_clinical_',
     c(0, 2) * year,
     c(5, 10) * year,
     timestep
   )
 
-  mockery::expect_args( 
+  mockery::expect_args(
     renderer$render_mock(),
     1,
     'n_0_1825',
@@ -173,7 +184,7 @@ test_that('that clinical incidence rendering works', {
     timestep
   )
 
-  mockery::expect_args( 
+  mockery::expect_args(
     renderer$render_mock(),
     2,
     'n_inc_clinical_0_1825',
@@ -181,19 +192,35 @@ test_that('that clinical incidence rendering works', {
     timestep
   )
 
-  mockery::expect_args( 
+  mockery::expect_args(
     renderer$render_mock(),
     3,
+    'p_inc_clinical_0_1825',
+    .3,
+    timestep
+  )
+
+  mockery::expect_args(
+    renderer$render_mock(),
+    4,
     'n_730_3650',
     3,
     timestep
   )
 
-  mockery::expect_args( 
+  mockery::expect_args(
     renderer$render_mock(),
-    4,
+    5,
     'n_inc_clinical_730_3650',
     2,
     timestep
   )
+
+  mockery::expect_args(
+    renderer$render_mock(),
+    6,
+    'p_inc_clinical_730_3650',
+    .6,
+    timestep
+  )
 })
diff --git a/tests/testthat/test-rtss.R b/tests/testthat/test-rtss.R
index 060a8a21..3ee0619a 100644
--- a/tests/testthat/test-rtss.R
+++ b/tests/testthat/test-rtss.R
@@ -90,6 +90,7 @@ test_that('Infection considers vaccine efficacy', {
     variables,
     bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
     parameters,
+    mock_render(timestep),
     timestep
   )
 

From 24d9df0cb94126681d450ab12fd5f75be5c8a778 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Tue, 12 Oct 2021 11:08:30 +0100
Subject: [PATCH 16/23] Remove redundant listener

---
 R/events.R | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/R/events.R b/R/events.R
index 6afa03b8..941de05d 100644
--- a/R/events.R
+++ b/R/events.R
@@ -162,11 +162,6 @@ attach_event_listeners <- function(
   events$recovery$add_listener(
     create_rate_listener('U', 'S', renderer)
   )
-  events$recovery$add_listener(
-    function(timestep, target) {
-      variables$is_severe$queue_update('no', target)
-    }
-  )
 
   # ===========
   # Progression

From f4650800f1a3628c51cad5a06feaab9b958ff7f3 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Tue, 12 Oct 2021 15:48:03 +0100
Subject: [PATCH 17/23] Add vignette on variation

---
 R/processes.R           |  3 +-
 vignettes/Vaccines.Rmd  |  1 -
 vignettes/Variation.Rmd | 91 +++++++++++++++++++++++++++++++++++++++++
 3 files changed, 93 insertions(+), 2 deletions(-)
 create mode 100644 vignettes/Variation.Rmd

diff --git a/R/processes.R b/R/processes.R
index 713ced0f..58b5a3b9 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -183,7 +183,8 @@ create_processes <- function(
       processes,
       create_total_M_renderer_compartmental(
         renderer,
-        solvers
+        solvers,
+        parameters
       )
     )
   }
diff --git a/vignettes/Vaccines.Rmd b/vignettes/Vaccines.Rmd
index 994af5c0..fb434ed0 100644
--- a/vignettes/Vaccines.Rmd
+++ b/vignettes/Vaccines.Rmd
@@ -17,7 +17,6 @@ knitr::opts_chunk$set(
 ```{r setup}
 suppressPackageStartupMessages(library(ggplot2))
 library(malariasimulation)
-library(malariaEquilibrium)
 ```
 
 # Parameterisation
diff --git a/vignettes/Variation.Rmd b/vignettes/Variation.Rmd
new file mode 100644
index 00000000..e8f39746
--- /dev/null
+++ b/vignettes/Variation.Rmd
@@ -0,0 +1,91 @@
+---
+title: "Variation"
+output: html_document
+vignette: >
+  %\VignetteIndexEntry{Variation}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include=FALSE}
+library(malariasimulation)
+library(cowplot)
+library(ggplot2)
+```
+
+## Variation in outputs
+
+Malariasimulation is a stochastic model, so there will be be variation in model outputs. For example, we could look at detected malaria cases over a year...
+
+```{r}
+# A small population
+p <- get_parameters(list(
+  human_population = 1000,
+  individual_mosquitoes = FALSE
+))
+p <- set_equilibrium(p, 2)
+small_o <- run_simulation(365, p)
+p <- get_parameters(list(
+  human_population = 10000,
+  individual_mosquitoes = FALSE
+))
+p <- set_equilibrium(p, 2)
+big_o <- run_simulation(365, p)
+
+plot_grid(
+  ggplot(small_o) + geom_line(aes(timestep, n_detect_730_3650)) + ggtitle('n_detect at n = 1,000'),
+  ggplot(small_o) + geom_line(aes(timestep, p_detect_730_3650)) + ggtitle('p_detect at n = 1,000'),
+  ggplot(big_o) + geom_line(aes(timestep, n_detect_730_3650)) + ggtitle('n_detect at n = 10,000'),
+  ggplot(big_o) + geom_line(aes(timestep, p_detect_730_3650)) + ggtitle('p_detect at n = 10,000')
+)
+```
+
+The n_detect output shows the result of sampling individuals who would be detected by microscopy.
+
+While the p_detect output shows the sum of probabilities of detection in the population.
+
+Notice that the p_detect output is slightly smoother. That's because it forgoes the sampling step. At low population sizes, p_detect will be smoother than the n_detect counterpart.
+
+## Estimating variation 
+
+We can estimate the variation in the number of detectable cases by repeating the simulation several times...
+
+```{r}
+outputs <- run_simulation_with_repetitions(
+  timesteps = 365,
+  repetitions = 10,
+  overrides = p,
+  parallel=TRUE
+)
+
+df <- aggregate(
+  outputs$n_detect_730_3650,
+  by=list(outputs$timestep),
+  FUN = function(x) {
+    c(
+      median = median(x),
+      lowq = unname(quantile(x, probs = .25)),
+      highq = unname(quantile(x, probs = .75)),
+      mean = mean(x),
+      lowci = mean(x) + 1.96*sd(x),
+      highci = mean(x) - 1.96*sd(x)
+    )
+  }
+)
+
+df <- data.frame(cbind(t = df$Group.1, df$x))
+
+plot_grid(
+  ggplot(df)
+  + geom_line(aes(x = t, y = median, group = 1))
+  + geom_ribbon(aes(x = t, ymin = lowq, ymax = highq), alpha = .2)
+  + xlab('timestep') + ylab('n_detect') + ggtitle('IQR spread'),
+  ggplot(df)
+  + geom_line(aes(x = t, y = mean, group = 1))
+  + geom_ribbon(aes(x = t, ymin = lowci, ymax = highci), alpha = .2)
+  + xlab('timestep') + ylab('n_detect') + ggtitle('95% confidence interval')
+)
+
+```
+
+These are useful techniques for comparing the expected variance from model outputs to outcomes from trials with lower population sizes.
\ No newline at end of file

From fb91ed40937d56d7f0646fd5cc35f3136e8bbf71 Mon Sep 17 00:00:00 2001
From: Hillary Topazian <80535782+htopazian@users.noreply.github.com>
Date: Tue, 19 Oct 2021 10:57:58 +0100
Subject: [PATCH 18/23] adding vignette for matching prev to EIR

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

diff --git a/vignettes/EIRprevmatch.Rmd b/vignettes/EIRprevmatch.Rmd
new file mode 100644
index 00000000..b8b111f6
--- /dev/null
+++ b/vignettes/EIRprevmatch.Rmd
@@ -0,0 +1,161 @@
+---
+title: "Matching PfPR2-10 to EIR"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{EIRprevmatch}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup, message=F, warning=F}
+library(tidyverse)
+library(malariasimulation)
+library(mgcv)
+library(fuzzyjoin)
+```
+
+Begin by parameterizing the model you wish to run. Be sure to include any seasonality factors and interventions you want to include at baseline. 
+
+```{r parameters}
+year <- 365
+human_population <- 1000
+
+params <- get_parameters(list(
+  human_population = human_population,
+  average_age = 8453.323,     # to match flat_demog
+  model_seasonality = TRUE,   # assign seasonality
+  g0 = 0.284596,                              
+  g = c(-0.317878, -0.0017527, 0.116455),     
+  h = c(-0.331361, 0.293128, -0.0617547),     
+  prevalence_rendering_min_ages = 2 * year,    # set prev age range
+  prevalence_rendering_max_ages = 10 * year,
+  fvt = 0,
+  v = 0,
+  severe_enabled = T,
+  individual_mosquitoes = FALSE))
+
+  # set species / drugs / treatment parameters
+params <- set_species(params, species = list(arab_params, fun_params, gamb_params), 
+                      proportions = c(0.25, 0.25, 0.5))
+params <- set_drugs(params, list(AL_params))
+params <- set_clinical_treatment(params, 1, c(1), c(0.45))
+
+```
+
+Now run `malariasimulation` over a range of EIRs. The goal is to run enough points to generate a curve to which you can match PfPR2-10 to EIR effectively.
+
+```{r model, message=F, warning=F}
+# loop over malariasimulation runs
+init_EIR <- c(0.01, 0.1, seq(1,9,1), seq(10, 50, by=5)) # set EIR values
+
+# run model
+outputs <- lapply(
+  init_EIR,
+  function(init) {
+    p_i <- set_equilibrium(params, init)
+    run_simulation(5 * year, p_i) # sim time = 5 years
+  }
+)
+
+# output EIR values
+EIR <- lapply(
+  outputs,
+  function(output) {
+    mean(
+      rowSums(
+        output[
+          output$timestep %in% seq(4 * 365, 5 * 365), # just use data from the last year
+          grepl('EIR_', names(output))
+        ] / human_population * year
+      )
+    )
+  }
+)
+
+# output prev 2-10 values
+prev <- lapply(
+  outputs,
+  function(output) {
+    mean(
+      output[
+        output$timestep %in% seq(4 * 365, 5 * 365),
+        'n_detect_730_3650'
+      ] / output[
+        output$timestep %in% seq(4 * 365, 5 * 365),
+        'n_730_3650'
+      ]
+    )
+  }
+)
+
+# create dataframe of initial EIR, output EIR, and prev 2-10 results
+EIR_prev <- cbind(init_EIR, unlist(EIR), unlist(prev)) %>% 
+  as_tibble() %>% rename(output_EIR = V2, prev = V3)
+
+```
+
+Now plot your results! Code is included to compare the results of matching PfPR2-10 to EIR based on `malariaEquilibrium` (blue line) versus matching based on parameterized `malariasimulation` runs (red line). Notice that the generated points do not form a smooth curve. We ran `malariasimulation` using a population of just 1,000. Increasing the population to 10,000 or even 100,000 will generate more accurate estimates, but will take longer to run. 
+
+```{r plot, warning=F}
+# calculate EIR / prev 2-10 relationship from malariEquilibrium
+eir <- seq(from = 0.1, to = 50, by=.5)
+eq_params <- malariaEquilibrium::load_parameter_set("Jamie_parameters.rds")
+
+prev <- vapply( # calculate prevalence between 2:10 for a bunch of EIRs
+  eir,
+  function(eir) {
+    eq <- malariaEquilibrium::human_equilibrium(
+      eir,
+      ft=0,
+      p=eq_params,
+      age=0:100
+    )
+    sum(eq$states[2:10, 'pos_M']) / sum(eq$states[2:10, 'prop'])
+  },
+  numeric(1)
+)
+
+prevmatch <- as_tibble(cbind(eir,prev))
+
+colors <- c("malariaEquilibrium" = "blue", "malariasimulation" = "red")
+
+p <- ggplot(data=EIR_prev) + 
+  geom_point(aes(x=init_EIR, y=prev), color='black') + 
+  stat_smooth(aes(x=init_EIR, y=prev, color='malariasimulation'), method = 'gam', n=1000) + 
+  geom_line(data=prevmatch, aes(x=eir, y=prev, color = 'malariaEquilibrium')) + 
+  labs(x='initial EIR', y=expression(paste(italic(Pf),"PR"[2-10])), color='model') + 
+  scale_color_manual(values = colors) + 
+  theme_classic()
+
+p
+
+```
+
+Now extract values from the fitted line to generate EIR estimates at your desired PfPR2-10 values.
+
+```{r table, warning=F}
+# extract points from stat_smooth
+p2 <- ggplot_build(p)
+p2 <- p2$data[[2]]
+
+p2 <- p2 %>% as_tibble() %>% select(x,y) %>% rename(init_EIR=x, pred=y) 
+
+# Pre-intervention baseline PfPR2-10 starting at values 10, 25, 35, 55 
+PfPR <- as_tibble_col(c(.10, .25, .35, .55), column_name="pfpr")
+
+
+# match via stat_smooth predictions
+match <- PfPR %>%
+  fuzzyjoin::difference_left_join(p2, by=c("pfpr"="pred"), 
+                                  max_dist=1, distance_col="dist") %>%
+  group_by(pfpr) %>% slice_min(dist) %>% rename(PfPR = pfpr, EIR = init_EIR, predicted_PfPR = pred, difference = dist)
+
+match
+```

From 3f3b0b40d4e3d8b4d0a12c3e919c0c7b8fd24d07 Mon Sep 17 00:00:00 2001
From: Hillary Topazian <80535782+htopazian@users.noreply.github.com>
Date: Wed, 20 Oct 2021 09:16:38 +0100
Subject: [PATCH 19/23] Apply suggestions from code review

re: Giovanni suggestions to simplify / eliminate dplyr

Co-authored-by: Giovanni <giovanni.charles@gmail.com>
---
 vignettes/EIRprevmatch.Rmd | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/vignettes/EIRprevmatch.Rmd b/vignettes/EIRprevmatch.Rmd
index b8b111f6..5aa4ecfe 100644
--- a/vignettes/EIRprevmatch.Rmd
+++ b/vignettes/EIRprevmatch.Rmd
@@ -96,8 +96,7 @@ prev <- lapply(
 )
 
 # create dataframe of initial EIR, output EIR, and prev 2-10 results
-EIR_prev <- cbind(init_EIR, unlist(EIR), unlist(prev)) %>% 
-  as_tibble() %>% rename(output_EIR = V2, prev = V3)
+EIR_prev <- cbind.data.frame(init_EIR, output_EIR = unlist(EIR), prev = unlist(prev))
 
 ```
 
@@ -122,11 +121,11 @@ prev <- vapply( # calculate prevalence between 2:10 for a bunch of EIRs
   numeric(1)
 )
 
-prevmatch <- as_tibble(cbind(eir,prev))
+prevmatch <- cbind.data.frame(eir, prev)
 
 colors <- c("malariaEquilibrium" = "blue", "malariasimulation" = "red")
 
-p <- ggplot(data=EIR_prev) + 
+ggplot(data=EIR_prev) + 
   geom_point(aes(x=init_EIR, y=prev), color='black') + 
   stat_smooth(aes(x=init_EIR, y=prev, color='malariasimulation'), method = 'gam', n=1000) + 
   geom_line(data=prevmatch, aes(x=eir, y=prev, color = 'malariaEquilibrium')) + 
@@ -134,8 +133,6 @@ p <- ggplot(data=EIR_prev) +
   scale_color_manual(values = colors) + 
   theme_classic()
 
-p
-
 ```
 
 Now extract values from the fitted line to generate EIR estimates at your desired PfPR2-10 values.
@@ -145,7 +142,7 @@ Now extract values from the fitted line to generate EIR estimates at your desire
 p2 <- ggplot_build(p)
 p2 <- p2$data[[2]]
 
-p2 <- p2 %>% as_tibble() %>% select(x,y) %>% rename(init_EIR=x, pred=y) 
+p2 <- data.frame(init_EIR = p2$x, pred = p2$y)
 
 # Pre-intervention baseline PfPR2-10 starting at values 10, 25, 35, 55 
 PfPR <- as_tibble_col(c(.10, .25, .35, .55), column_name="pfpr")

From 5528e6e4d662f857a4763c9749182773951d880a Mon Sep 17 00:00:00 2001
From: Hillary Topazian <80535782+htopazian@users.noreply.github.com>
Date: Wed, 20 Oct 2021 11:43:17 +0100
Subject: [PATCH 20/23] eliminating dependencies

---
 vignettes/EIRprevmatch.Rmd | 72 ++++++++++++++++++++++++--------------
 1 file changed, 45 insertions(+), 27 deletions(-)

diff --git a/vignettes/EIRprevmatch.Rmd b/vignettes/EIRprevmatch.Rmd
index 5aa4ecfe..a30b0f0d 100644
--- a/vignettes/EIRprevmatch.Rmd
+++ b/vignettes/EIRprevmatch.Rmd
@@ -15,17 +15,15 @@ knitr::opts_chunk$set(
 ```
 
 ```{r setup, message=F, warning=F}
-library(tidyverse)
 library(malariasimulation)
 library(mgcv)
-library(fuzzyjoin)
 ```
 
 Begin by parameterizing the model you wish to run. Be sure to include any seasonality factors and interventions you want to include at baseline. 
 
 ```{r parameters}
 year <- 365
-human_population <- 1000
+human_population <- 5000
 
 params <- get_parameters(list(
   human_population = human_population,
@@ -100,10 +98,10 @@ EIR_prev <- cbind.data.frame(init_EIR, output_EIR = unlist(EIR), prev = unlist(p
 
 ```
 
-Now plot your results! Code is included to compare the results of matching PfPR2-10 to EIR based on `malariaEquilibrium` (blue line) versus matching based on parameterized `malariasimulation` runs (red line). Notice that the generated points do not form a smooth curve. We ran `malariasimulation` using a population of just 1,000. Increasing the population to 10,000 or even 100,000 will generate more accurate estimates, but will take longer to run. 
+Now plot your results! Code is included to compare the results of matching PfPR2-10 to EIR based on `malariaEquilibrium` (blue line) versus matching based on parameterized `malariasimulation` runs (red line). Notice that the generated points do not form a smooth curve. We ran `malariasimulation` using a population of just 5,000. Increasing the population to 10,000 or even 100,000 will generate more accurate estimates, but will take longer to run. 
 
 ```{r plot, warning=F}
-# calculate EIR / prev 2-10 relationship from malariEquilibrium
+# calculate EIR / prev 2-10 relationship from malariaEquilibrium
 eir <- seq(from = 0.1, to = 50, by=.5)
 eq_params <- malariaEquilibrium::load_parameter_set("Jamie_parameters.rds")
 
@@ -123,36 +121,56 @@ prev <- vapply( # calculate prevalence between 2:10 for a bunch of EIRs
 
 prevmatch <- cbind.data.frame(eir, prev)
 
-colors <- c("malariaEquilibrium" = "blue", "malariasimulation" = "red")
-
-ggplot(data=EIR_prev) + 
-  geom_point(aes(x=init_EIR, y=prev), color='black') + 
-  stat_smooth(aes(x=init_EIR, y=prev, color='malariasimulation'), method = 'gam', n=1000) + 
-  geom_line(data=prevmatch, aes(x=eir, y=prev, color = 'malariaEquilibrium')) + 
-  labs(x='initial EIR', y=expression(paste(italic(Pf),"PR"[2-10])), color='model') + 
-  scale_color_manual(values = colors) + 
-  theme_classic()
+# calculate best fit line through malariasimulation data
+fit <- predict(gam(prev~s(init_EIR, k=-1), data=EIR_prev), 
+               newdata = data.frame(init_EIR=c(0,seq(0.1,50,0.1))),  type="response")
+fit <- cbind(fit, data.frame(init_EIR=c(0,seq(0.1,50,0.1))))
+
+# plot
+plot(x=1, type = "n", 
+     frame = F, 
+     xlab = "initial EIR", 
+     ylab = expression(paste(italic(Pf),"PR"[2-10])),
+     xlim = c(0,50), 
+     ylim = c(0,.7))
+
+points(x=EIR_prev$init_EIR, 
+       y=EIR_prev$prev,
+       pch = 19,
+       col = 'black')
+
+lines(x=fit$init_EIR, 
+      y=fit$fit, 
+      col = "red", 
+      type = "l", 
+      lty = 1)
+
+lines(x=prevmatch$eir, 
+      y=prevmatch$prev, 
+      col = "blue", 
+      type = "l", 
+      lty = 1)
+
+legend("bottomright", legend=c("malariasimulation", "malariaEquilibrium"), col=c("red", "blue"), lty = c(1,1), cex=0.8)
 
 ```
 
 Now extract values from the fitted line to generate EIR estimates at your desired PfPR2-10 values.
 
 ```{r table, warning=F}
-# extract points from stat_smooth
-p2 <- ggplot_build(p)
-p2 <- p2$data[[2]]
-
-p2 <- data.frame(init_EIR = p2$x, pred = p2$y)
-
 # Pre-intervention baseline PfPR2-10 starting at values 10, 25, 35, 55 
-PfPR <- as_tibble_col(c(.10, .25, .35, .55), column_name="pfpr")
-
+PfPR <- c(.10, .25, .35, .45)
 
 # match via stat_smooth predictions
-match <- PfPR %>%
-  fuzzyjoin::difference_left_join(p2, by=c("pfpr"="pred"), 
-                                  max_dist=1, distance_col="dist") %>%
-  group_by(pfpr) %>% slice_min(dist) %>% rename(PfPR = pfpr, EIR = init_EIR, predicted_PfPR = pred, difference = dist)
+match <- function(x){
+  
+  m <- which.min(abs(fit$fit-x)) # index of closest prev match
+  fit[m,2]                       # print match
+  
+}
+
+eir <- unlist(lapply(PfPR, match))
+
+cbind.data.frame(PfPR, eir)
 
-match
 ```

From 0a3de425f13d922aaaeafbb2d7552a218d194783 Mon Sep 17 00:00:00 2001
From: Hillary Topazian <80535782+htopazian@users.noreply.github.com>
Date: Wed, 20 Oct 2021 13:41:38 +0100
Subject: [PATCH 21/23] adding mgcv to description

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

diff --git a/DESCRIPTION b/DESCRIPTION
index 36e7bd70..1e65b852 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -40,7 +40,8 @@ Suggests:
   DiagrammeR,
   cowplot,
   ggplot2,
-  covr
+  covr,
+  mgcv 
 RoxygenNote: 7.1.2
 Roxygen: list(markdown = TRUE)
 LinkingTo:

From 7a2cc46d24cb59fb67baa382955c0cc6b8d082e0 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Fri, 29 Oct 2021 15:36:50 +0100
Subject: [PATCH 22/23] Faster vignettes

---
 .Rbuildignore              |  1 -
 vignettes/EIRprevmatch.Rmd |  8 ++++----
 vignettes/Treatment.Rmd    | 13 ++++---------
 3 files changed, 8 insertions(+), 14 deletions(-)

diff --git a/.Rbuildignore b/.Rbuildignore
index 14fb78e2..d969fb5a 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -9,4 +9,3 @@ docker
 LICENSE.md
 codecov.yml
 .github
-vignettes/*
diff --git a/vignettes/EIRprevmatch.Rmd b/vignettes/EIRprevmatch.Rmd
index a30b0f0d..db608a14 100644
--- a/vignettes/EIRprevmatch.Rmd
+++ b/vignettes/EIRprevmatch.Rmd
@@ -49,9 +49,9 @@ params <- set_clinical_treatment(params, 1, c(1), c(0.45))
 
 Now run `malariasimulation` over a range of EIRs. The goal is to run enough points to generate a curve to which you can match PfPR2-10 to EIR effectively.
 
-```{r model, message=F, warning=F}
+```{r model}
 # loop over malariasimulation runs
-init_EIR <- c(0.01, 0.1, seq(1,9,1), seq(10, 50, by=5)) # set EIR values
+init_EIR <- c(0.01, 0.1, 1, 5, 10, 25, 50) # set EIR values
 
 # run model
 outputs <- lapply(
@@ -100,7 +100,7 @@ EIR_prev <- cbind.data.frame(init_EIR, output_EIR = unlist(EIR), prev = unlist(p
 
 Now plot your results! Code is included to compare the results of matching PfPR2-10 to EIR based on `malariaEquilibrium` (blue line) versus matching based on parameterized `malariasimulation` runs (red line). Notice that the generated points do not form a smooth curve. We ran `malariasimulation` using a population of just 5,000. Increasing the population to 10,000 or even 100,000 will generate more accurate estimates, but will take longer to run. 
 
-```{r plot, warning=F}
+```{r plot}
 # calculate EIR / prev 2-10 relationship from malariaEquilibrium
 eir <- seq(from = 0.1, to = 50, by=.5)
 eq_params <- malariaEquilibrium::load_parameter_set("Jamie_parameters.rds")
@@ -122,7 +122,7 @@ prev <- vapply( # calculate prevalence between 2:10 for a bunch of EIRs
 prevmatch <- cbind.data.frame(eir, prev)
 
 # calculate best fit line through malariasimulation data
-fit <- predict(gam(prev~s(init_EIR, k=-1), data=EIR_prev), 
+fit <- predict(gam(prev~s(init_EIR, k=5), data=EIR_prev), 
                newdata = data.frame(init_EIR=c(0,seq(0.1,50,0.1))),  type="response")
 fit <- cbind(fit, data.frame(init_EIR=c(0,seq(0.1,50,0.1))))
 
diff --git a/vignettes/Treatment.Rmd b/vignettes/Treatment.Rmd
index d2be6bfd..3986c202 100644
--- a/vignettes/Treatment.Rmd
+++ b/vignettes/Treatment.Rmd
@@ -24,17 +24,11 @@ Set up some baseline parameters:
 
 ```{r}
 year <- 365
-sim_length <- 1 * year
+sim_length <- 2 * year
 human_population <- 1000
 
 simparams <- get_parameters(list(human_population = human_population))
 simparams <- set_drugs(simparams, list(AL_params, DHA_PQP_params))
-
-# NOTE: this equilibrium is calculated using a different prophylaxis model
-# it is just used as a starting point
-get_equilibrium <- function(EIR, ft) {
-  human_equilibrium(EIR = EIR, ft = ft, p = jamie_params, age = 0:99)
-}
 ```
 
 Run the model for some starting EIRs and fts:
@@ -45,11 +39,12 @@ starting_EIR <- 10
 fts <- c(0, .25, .5, .75, 1.)
 
 for (ft in fts) {
-  eq_params <- load_parameter_set()
   simparams <- set_clinical_treatment(simparams, 1, 1, ft / 2)
   simparams <- set_clinical_treatment(simparams, 2, 1, ft / 2)
-  simparams <- set_equilibrium(simparams, starting_EIR, eq_params)
+  simparams <- set_equilibrium(simparams, starting_EIR)
   output <- run_simulation(sim_length, simparams)
+  # cut out the warmup
+  output <- output[output$timestep > year,]
   outputs[[length(outputs) + 1]] <- output
 }
 ```

From 86c36b968bb4d1fe30d48fda95077da026277708 Mon Sep 17 00:00:00 2001
From: Giovanni Charles <giovanni.charles@gmail.com>
Date: Fri, 29 Oct 2021 15:39:21 +0100
Subject: [PATCH 23/23] Fix individual dependency and NEWS

---
 DESCRIPTION | 3 +--
 NEWS.md     | 3 ++-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 1e65b852..1affbf8b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -20,9 +20,8 @@ Encoding: UTF-8
 LazyData: true
 Remotes:
   mrc-ide/malariaEquilibrium,
-  mrc-ide/individual
 Imports:
-  individual (>= 0.1.6),
+  individual (>= 0.1.7),
   malariaEquilibrium,
   Rcpp,
   statmod,
diff --git a/NEWS.md b/NEWS.md
index aed84d5a..0a813cf8 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,7 @@
-# malariasimulation 1.2.0 (wip)
+# malariasimulation 1.2.0
 
   * added a `news.md` file to track changes to the package.
   * n_inc_clinical includes treated humans
   * disaggregate EIR and state counts by mosquito species
   * remove redundant Total_M and and global EIR outputs
+  * new vignette on how to calibrate EIR to prevalence