From f6d8e3e437b7fc88cb5b04cec8c85c69fed2ba80 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Fri, 28 Jun 2024 16:40:21 +0100
Subject: [PATCH 01/14] p.v allows subpatent infections to occur, so we now
 calculate patent/subpatent infections, render these and schedule them.

---
 R/human_infection.R                         | 146 ++++++++++++++++----
 tests/testthat/test-infection-integration.R |   2 +
 tests/testthat/test-vivax.R                 |  65 +++++++++
 3 files changed, 184 insertions(+), 29 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index 113acb98..6540e7ff 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -180,6 +180,25 @@ infection_outcome_process <- function(
         timestep,
         parameters$ud
       )
+      
+      clinical_infections <- calculate_clinical_infections(
+        variables,
+        infected_humans,
+        parameters,
+        renderer,
+        timestep
+      )
+      
+      update_severe_disease(
+        timestep,
+        infected_humans,
+        variables,
+        parameters,
+        renderer
+      )
+      
+      patent_infections <- NULL
+      
     } else if (parameters$parasite == "vivax"){
       boost_immunity(
         variables$iaa,
@@ -188,27 +207,27 @@ infection_outcome_process <- function(
         timestep,
         parameters$ua
       )
+      
+      ## Only S and U infections need to be split using the patent infection function
+      patent_infections <- calculate_patent_infections(
+        variables,
+        variables$state$get_index_of(c("S","U"))$and(infected_humans),
+        parameters,
+        renderer,
+        timestep
+      )
+      
+      # Patent level infected S and U, and all A infections to get clinical infections
+      clinical_infections <- calculate_clinical_infections(
+        variables,
+        variables$state$get_index_of("A")$and(infected_humans)$or(patent_infections),
+        parameters,
+        renderer,
+        timestep
+      )
     }
   }
   
-  clinical_infections <- calculate_clinical_infections(
-    variables,
-    infected_humans,
-    parameters,
-    renderer,
-    timestep
-  )
-  
-  if(parameters$parasite == "falciparum"){
-    update_severe_disease(
-      timestep,
-      infected_humans,
-      variables,
-      parameters,
-      renderer
-    )
-  }
-  
   treated <- calculate_treated(
     variables,
     clinical_infections,
@@ -221,6 +240,7 @@ infection_outcome_process <- function(
   
   schedule_infections(
     variables,
+    patent_infections,
     clinical_infections,
     treated,
     infected_humans,
@@ -229,6 +249,53 @@ infection_outcome_process <- function(
   )
 }
 
+#' @title Calculate patent infections (p.v only)
+#' @description
+#' Sample patent infections from all infections
+#' @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_patent_infections <- function(
+    variables,
+    infections,
+    parameters,
+    renderer,
+    timestep
+) {
+  
+  iaa <- variables$iaa$get_values(infections)
+  iam <- variables$iam$get_values(infections)
+  
+  philm <- anti_parasite_immunity(
+    min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
+    k = parameters$klm, iaa = iaa, iam = iam)
+  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    patent_infections,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  incidence_probability_renderer(
+    variables$birth,
+    renderer,
+    infections,
+    philm,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  patent_infections
+}
+
 #' @title Calculate clinical infections
 #' @description
 #' Sample clinical infections from all infections
@@ -459,20 +526,18 @@ calculate_treated <- function(
 #' @noRd
 schedule_infections <- function(
     variables,
+    patent_infections,
     clinical_infections,
     treated,
     infections,
     parameters,
     timestep
 ) {
+  
   included <- treated$not(TRUE)
+  to_infect_clinical <- clinical_infections$and(included)
 
-  to_infect <- clinical_infections$and(included)
-  to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(
-    included
-  )
-
-  if(to_infect$size() > 0) {
+  if(to_infect_clinical$size() > 0) {
     update_infection(
       variables$state,
       'D',
@@ -480,12 +545,14 @@ schedule_infections <- function(
       parameters$cd,
       variables$recovery_rates,
       1/parameters$dd,
-      to_infect
+      to_infect_clinical
     )
   }
 
-  if(to_infect_asym$size() > 0) {
-    if(parameters$parasite == "falciparum"){
+  if(parameters$parasite == "falciparum"){
+    to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(included)
+    
+    if(to_infect_asym$size() > 0){
       # p.f has immunity-determined asymptomatic infectivity
       update_to_asymptomatic_infection(
         variables,
@@ -493,7 +560,12 @@ schedule_infections <- function(
         timestep,
         to_infect_asym
       )
-    } else if (parameters$parasite == "vivax"){
+    }
+  } else if (parameters$parasite == "vivax"){
+    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(patent_infections$not(FALSE))
+    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(patent_infections)$and(clinical_infections$not(FALSE))
+    
+    if(to_infect_asym$size() > 0){
       # p.v has constant asymptomatic infectivity
       update_infection(
         variables$state,
@@ -502,7 +574,23 @@ schedule_infections <- function(
         parameters$ca,
         variables$recovery_rates,
         1/parameters$da,
-        to_infect
+        to_infect_asym
+      )
+    }
+    if(to_infect_subpatent$size() > 0){
+      # p.v subpatent recovery rate is immunity dependent
+      update_infection(
+        variables$state,
+        'U',
+        variables$infectivity,
+        parameters$cu,
+        variables$recovery_rates,
+        1/anti_parasite_immunity(
+          parameters$dpcr_min, parameters$dpcr_max, parameters$apcr50, parameters$kpcr,
+          variables$iaa$get_values(to_infect_subpatent),
+          variables$iam$get_values(to_infect_subpatent)
+        ),
+        to_infect_subpatent
       )
     }
   }
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index d399a527..17d0e958 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -157,6 +157,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     schedule_mock,
     1,
     variables,
+    NULL,
     clinical,
     treated,
     infected,
@@ -659,6 +660,7 @@ test_that('schedule_infections correctly schedules new infections', {
   
   schedule_infections(
     variables,
+    NULL,
     clinical_infections,
     treated,
     infections,
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index db443891..289848ce 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -151,3 +151,68 @@ test_that('that vivax patent prevalence rendering works', {
   )
   
 })
+
+test_that('Test default vivax incidence rendering works', {
+  
+  timestep <- 0
+  year <- 365
+  birth <- individual::IntegerVariable$new(
+    -c(2, 5, 10, 11) * year
+  )
+  vivax_parameters <- get_parameters(
+    parasite = "vivax")
+  
+  renderer <- mock_render(1)
+  incidence_renderer(
+    birth,
+    renderer,
+    individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    'inc_patent_',
+    c(0, 2) * year,
+    c(5, 10) * year,
+    timestep
+  )
+  
+  incidence_probability_renderer(
+    birth,
+    renderer,
+    individual::Bitset$new(4)$insert(seq(4)),
+    c(.1, .2, .3, .4),
+    'inc_patent_',
+    c(0, 2) * year,
+    c(5, 10) * year,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_inc_patent_0_1825',
+    2,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    2,
+    'n_inc_patent_730_3650',
+    2,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    3,
+    'p_inc_patent_0_1825',
+    0.3,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    4,
+    'p_inc_patent_730_3650',
+    .6,
+    timestep
+  )
+})

From 6b577290cc213feeb5fda8f252e7c3d46d024588 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Fri, 28 Jun 2024 16:40:21 +0100
Subject: [PATCH 02/14] p.v allows subpatent infections to occur, so we now
 calculate patent/subpatent infections, render these and schedule them.

---
 R/human_infection.R                         | 144 ++++++++++++++++----
 tests/testthat/test-infection-integration.R |   2 +
 tests/testthat/test-vivax.R                 |  65 +++++++++
 3 files changed, 183 insertions(+), 28 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index d4450cf2..6540e7ff 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -180,6 +180,25 @@ infection_outcome_process <- function(
         timestep,
         parameters$ud
       )
+      
+      clinical_infections <- calculate_clinical_infections(
+        variables,
+        infected_humans,
+        parameters,
+        renderer,
+        timestep
+      )
+      
+      update_severe_disease(
+        timestep,
+        infected_humans,
+        variables,
+        parameters,
+        renderer
+      )
+      
+      patent_infections <- NULL
+      
     } else if (parameters$parasite == "vivax"){
       boost_immunity(
         variables$iaa,
@@ -188,27 +207,27 @@ infection_outcome_process <- function(
         timestep,
         parameters$ua
       )
+      
+      ## Only S and U infections need to be split using the patent infection function
+      patent_infections <- calculate_patent_infections(
+        variables,
+        variables$state$get_index_of(c("S","U"))$and(infected_humans),
+        parameters,
+        renderer,
+        timestep
+      )
+      
+      # Patent level infected S and U, and all A infections to get clinical infections
+      clinical_infections <- calculate_clinical_infections(
+        variables,
+        variables$state$get_index_of("A")$and(infected_humans)$or(patent_infections),
+        parameters,
+        renderer,
+        timestep
+      )
     }
   }
   
-  clinical_infections <- calculate_clinical_infections(
-    variables,
-    infected_humans,
-    parameters,
-    renderer,
-    timestep
-  )
-  
-  if(parameters$parasite == "falciparum"){
-    update_severe_disease(
-      timestep,
-      infected_humans,
-      variables,
-      parameters,
-      renderer
-    )
-  }
-  
   treated <- calculate_treated(
     variables,
     clinical_infections,
@@ -221,6 +240,7 @@ infection_outcome_process <- function(
   
   schedule_infections(
     variables,
+    patent_infections,
     clinical_infections,
     treated,
     infected_humans,
@@ -229,6 +249,53 @@ infection_outcome_process <- function(
   )
 }
 
+#' @title Calculate patent infections (p.v only)
+#' @description
+#' Sample patent infections from all infections
+#' @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_patent_infections <- function(
+    variables,
+    infections,
+    parameters,
+    renderer,
+    timestep
+) {
+  
+  iaa <- variables$iaa$get_values(infections)
+  iam <- variables$iam$get_values(infections)
+  
+  philm <- anti_parasite_immunity(
+    min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
+    k = parameters$klm, iaa = iaa, iam = iam)
+  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    patent_infections,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  incidence_probability_renderer(
+    variables$birth,
+    renderer,
+    infections,
+    philm,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  patent_infections
+}
+
 #' @title Calculate clinical infections
 #' @description
 #' Sample clinical infections from all infections
@@ -459,20 +526,18 @@ calculate_treated <- function(
 #' @noRd
 schedule_infections <- function(
     variables,
+    patent_infections,
     clinical_infections,
     treated,
     infections,
     parameters,
     timestep
 ) {
+  
   included <- treated$not(TRUE)
+  to_infect_clinical <- clinical_infections$and(included)
 
-  to_infect <- clinical_infections$and(included)
-  to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(
-    included
-  )
-
-  if(to_infect$size() > 0) {
+  if(to_infect_clinical$size() > 0) {
     update_infection(
       variables$state,
       'D',
@@ -480,12 +545,14 @@ schedule_infections <- function(
       parameters$cd,
       variables$recovery_rates,
       1/parameters$dd,
-      to_infect
+      to_infect_clinical
     )
   }
 
-  if(to_infect_asym$size() > 0) {
-    if(parameters$parasite == "falciparum"){
+  if(parameters$parasite == "falciparum"){
+    to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(included)
+    
+    if(to_infect_asym$size() > 0){
       # p.f has immunity-determined asymptomatic infectivity
       update_to_asymptomatic_infection(
         variables,
@@ -493,7 +560,12 @@ schedule_infections <- function(
         timestep,
         to_infect_asym
       )
-    } else if (parameters$parasite == "vivax"){
+    }
+  } else if (parameters$parasite == "vivax"){
+    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(patent_infections$not(FALSE))
+    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(patent_infections)$and(clinical_infections$not(FALSE))
+    
+    if(to_infect_asym$size() > 0){
       # p.v has constant asymptomatic infectivity
       update_infection(
         variables$state,
@@ -505,6 +577,22 @@ schedule_infections <- function(
         to_infect_asym
       )
     }
+    if(to_infect_subpatent$size() > 0){
+      # p.v subpatent recovery rate is immunity dependent
+      update_infection(
+        variables$state,
+        'U',
+        variables$infectivity,
+        parameters$cu,
+        variables$recovery_rates,
+        1/anti_parasite_immunity(
+          parameters$dpcr_min, parameters$dpcr_max, parameters$apcr50, parameters$kpcr,
+          variables$iaa$get_values(to_infect_subpatent),
+          variables$iam$get_values(to_infect_subpatent)
+        ),
+        to_infect_subpatent
+      )
+    }
   }
 }
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index d399a527..17d0e958 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -157,6 +157,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     schedule_mock,
     1,
     variables,
+    NULL,
     clinical,
     treated,
     infected,
@@ -659,6 +660,7 @@ test_that('schedule_infections correctly schedules new infections', {
   
   schedule_infections(
     variables,
+    NULL,
     clinical_infections,
     treated,
     infections,
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 47427d99..8dc62d26 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -174,3 +174,68 @@ test_that('that vivax patent prevalence rendering works', {
   )
   
 })
+
+test_that('Test default vivax incidence rendering works', {
+  
+  timestep <- 0
+  year <- 365
+  birth <- individual::IntegerVariable$new(
+    -c(2, 5, 10, 11) * year
+  )
+  vivax_parameters <- get_parameters(
+    parasite = "vivax")
+  
+  renderer <- mock_render(1)
+  incidence_renderer(
+    birth,
+    renderer,
+    individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    'inc_patent_',
+    c(0, 2) * year,
+    c(5, 10) * year,
+    timestep
+  )
+  
+  incidence_probability_renderer(
+    birth,
+    renderer,
+    individual::Bitset$new(4)$insert(seq(4)),
+    c(.1, .2, .3, .4),
+    'inc_patent_',
+    c(0, 2) * year,
+    c(5, 10) * year,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_inc_patent_0_1825',
+    2,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    2,
+    'n_inc_patent_730_3650',
+    2,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    3,
+    'p_inc_patent_0_1825',
+    0.3,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    4,
+    'p_inc_patent_730_3650',
+    .6,
+    timestep
+  )
+})

From 8af7dafad1b52ee9b8d741dd665916f82a8c2976 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Thu, 11 Jul 2024 08:58:43 +0100
Subject: [PATCH 03/14] I've renamed patent infections as light microscopy
 detectable infections or lm_det, to avoid the confusion of what exactly
 patent means.

I've also removed the patent/lm_det infections <- NULL line from the falciparum code, but in doing so, have replicated some infection functions for falciparum and vivax so that schedule_infections may (for vivax) or may not (for falciparum) use the lm_det subset. I've also reordered the arguments for that function so that the lm_det_infections can be at the end, but make sense in the context of other arguments.
---
 R/human_infection.R                         | 117 ++++++++++++--------
 R/parameters.R                              |  14 +--
 tests/testthat/test-infection-integration.R |  18 ++-
 3 files changed, 84 insertions(+), 65 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index 6540e7ff..56fe0340 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -197,7 +197,24 @@ infection_outcome_process <- function(
         renderer
       )
       
-      patent_infections <- NULL
+      treated <- calculate_treated(
+        variables,
+        clinical_infections,
+        parameters,
+        timestep,
+        renderer
+      )
+      
+      renderer$render('n_infections', infected_humans$size(), timestep)
+      
+      schedule_infections(
+        parameters,
+        variables,
+        timestep,
+        infected_humans,
+        treated,
+        clinical_infections
+      )
       
     } else if (parameters$parasite == "vivax"){
       boost_immunity(
@@ -208,8 +225,8 @@ infection_outcome_process <- function(
         parameters$ua
       )
       
-      ## Only S and U infections need to be split using the patent infection function
-      patent_infections <- calculate_patent_infections(
+      ## Only S and U infections need to be split using the lm-detectable infection function
+      lm_det_infections <- calculate_lm_det_infections(
         variables,
         variables$state$get_index_of(c("S","U"))$and(infected_humans),
         parameters,
@@ -217,48 +234,48 @@ infection_outcome_process <- function(
         timestep
       )
       
-      # Patent level infected S and U, and all A infections to get clinical infections
+      # Lm-detectable level infected S and U, and all A infections may receive clinical infections
       clinical_infections <- calculate_clinical_infections(
         variables,
-        variables$state$get_index_of("A")$and(infected_humans)$or(patent_infections),
+        variables$state$get_index_of("A")$and(infected_humans)$or(lm_det_infections),
         parameters,
         renderer,
         timestep
       )
+      
+      treated <- calculate_treated(
+        variables,
+        clinical_infections,
+        parameters,
+        timestep,
+        renderer
+      )
+      
+      renderer$render('n_infections', infected_humans$size(), timestep)
+      
+      schedule_infections(
+        parameters,
+        variables,
+        timestep,
+        infected_humans,
+        treated,
+        clinical_infections,
+        lm_det_infections
+      )
     }
   }
-  
-  treated <- calculate_treated(
-    variables,
-    clinical_infections,
-    parameters,
-    timestep,
-    renderer
-  )
-  
-  renderer$render('n_infections', infected_humans$size(), timestep)
-  
-  schedule_infections(
-    variables,
-    patent_infections,
-    clinical_infections,
-    treated,
-    infected_humans,
-    parameters,
-    timestep
-  )
 }
 
-#' @title Calculate patent infections (p.v only)
+#' @title Calculate light microscopy detectable infections (p.v only)
 #' @description
-#' Sample patent infections from all infections
+#' Sample light microscopy detectable infections from all infections
 #' @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_patent_infections <- function(
+calculate_lm_det_infections <- function(
     variables,
     infections,
     parameters,
@@ -272,15 +289,15 @@ calculate_patent_infections <- function(
   philm <- anti_parasite_immunity(
     min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
     k = parameters$klm, iaa = iaa, iam = iam)
-  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  lm_det_infections <- bitset_at(infections, bernoulli_multi_p(philm))
   
   incidence_renderer(
     variables$birth,
     renderer,
-    patent_infections,
-    'inc_patent_',
-    parameters$patent_incidence_rendering_min_ages,
-    parameters$patent_incidence_rendering_max_ages,
+    lm_det_infections,
+    'inc_lm_det_',
+    parameters$lm_det_incidence_rendering_min_ages,
+    parameters$lm_det_incidence_rendering_max_ages,
     timestep
   )
   incidence_probability_renderer(
@@ -288,12 +305,12 @@ calculate_patent_infections <- function(
     renderer,
     infections,
     philm,
-    'inc_patent_',
-    parameters$patent_incidence_rendering_min_ages,
-    parameters$patent_incidence_rendering_max_ages,
+    'inc_lm_det_',
+    parameters$lm_det_incidence_rendering_min_ages,
+    parameters$lm_det_incidence_rendering_max_ages,
     timestep
   )
-  patent_infections
+  lm_det_infections
 }
 
 #' @title Calculate clinical infections
@@ -518,20 +535,24 @@ calculate_treated <- function(
 #' @title Schedule infections
 #' @description
 #' Schedule infections in humans after the incubation period
-#' @param events a list of all of the model events
-#' @param clinical_infections bitset of clinically infected humans
-#' @param treated bitset of treated humans
-#' @param infections bitset of infected humans
 #' @param parameters model parameters
+#' @param variables a list of all of the model variables
+#' @param timestep current timestep
+#' @param infections bitset of infected humans
+#' @param treated bitset of treated humans
+#' @param clinical_infections bitset of clinically infected humans
+#' @param lm_det_infections bitset of lm-detectable infected humans (p.v only:
+#'  lm_det infections are modelled as a human state, rather than a subset of  
+#'  asymptomatic infections as in the p.f model)
 #' @noRd
 schedule_infections <- function(
+    parameters,
     variables,
-    patent_infections,
-    clinical_infections,
-    treated,
+    timestep,
     infections,
-    parameters,
-    timestep
+    treated,
+    clinical_infections,
+    lm_det_infections = NULL
 ) {
   
   included <- treated$not(TRUE)
@@ -562,8 +583,8 @@ schedule_infections <- function(
       )
     }
   } else if (parameters$parasite == "vivax"){
-    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(patent_infections$not(FALSE))
-    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(patent_infections)$and(clinical_infections$not(FALSE))
+    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(lm_det_infections$not(FALSE))
+    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(lm_det_infections)$and(clinical_infections$not(FALSE))
     
     if(to_infect_asym$size() > 0){
       # p.v has constant asymptomatic infectivity
diff --git a/R/parameters.R b/R/parameters.R
index 95e7f776..5af67575 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -114,7 +114,7 @@
 #' * id0 - scale parameter; default = 1.577533
 #' * kd - shape parameter; default = 0.476614
 #'
-#' probability of patent infection due to anti-parasite immunity (p.v only):
+#' probability of light miscroscopy detectable infection due to anti-parasite immunity (p.v only):
 #'
 #' * phi0lm - maximum probability due to no immunity; default = 0.8918
 #' * phi1lm - maximum reduction due to immunity; default = 0.00482170890334156
@@ -143,7 +143,7 @@
 #' * cd - infectivity of clinically diseased humans towards mosquitoes; default = 0.068
 #' * ca - infectivity of asymptomatic humans towards mosquitoes (p.v only); default = 0.1
 #' * gamma1 - parameter for infectivity of asymptomatic humans; default = 1.82425
-#' * cu - infectivity of sub-patent infection; default = 0.0062
+#' * cu - infectivity of subpatent infection; default = 0.0062
 #' * ct - infectivity of treated infection; default = 0.021896
 #'
 #' mosquito fixed state transitions (including mortality):
@@ -280,8 +280,8 @@
 #' * incidence_rendering_min_ages - the minimum ages for incidence
 #' outputs (includes asymptomatic microscopy +); default = turned off
 #' * incidence_rendering_max_ages - the corresponding max ages; default = turned off
-#' * patent_incidence_rendering_min_ages - the minimum ages for patent incidence outputs (LM detectable), (p.v only); default = turned off
-#' * patent_incidence_rendering_max_ages - the corresponding max ages (p.v only); default = turned off
+#' * lm_det_incidence_rendering_min_ages - the minimum ages for light miscroscopy detectable incidence outputs, (p.v only); default = turned off
+#' * lm_det_incidence_rendering_max_ages - the corresponding max ages (p.v only); 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
 #' * severe_incidence_rendering_min_ages - the minimum ages for severe incidence
@@ -341,7 +341,7 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") {
     # maternal immunity parameters
     # probability of pre-erythrocytic infection/blood immunity
     # probability of asymptomatic detection (p.f only)
-    # probability of patent infection (due to anti-parasite immunity, p.v only)
+    # probability of l infection (due to anti-parasite immunity, p.v only)
     # probability of clinical infection
     # probability of severe infection (p.f only)
     # infectivity towards mosquitos
@@ -467,8 +467,8 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") {
       age_group_rendering_max_ages = numeric(0),
       incidence_rendering_min_ages = numeric(0),
       incidence_rendering_max_ages = numeric(0),
-      patent_incidence_rendering_min_ages = numeric(0),
-      patent_incidence_rendering_max_ages = numeric(0),
+      lm_det_incidence_rendering_min_ages = numeric(0),
+      lm_det_incidence_rendering_max_ages = numeric(0),
       clinical_incidence_rendering_min_ages = numeric(0),
       clinical_incidence_rendering_max_ages = numeric(0),
       severe_incidence_rendering_min_ages = numeric(0),
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 17d0e958..2d99c488 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -156,13 +156,12 @@ test_that('simulate_infection integrates different types of infection and schedu
   mockery::expect_args(
     schedule_mock,
     1,
+    parameters,
     variables,
-    NULL,
-    clinical,
-    treated,
+    timestep,
     infected,
-    parameters,
-    timestep
+    treated,
+    clinical
   )
 })
 
@@ -659,13 +658,12 @@ test_that('schedule_infections correctly schedules new infections', {
   mockery::stub(schedule_infections, 'update_to_asymptomatic_infection', asymp_mock)
   
   schedule_infections(
+    parameters,
     variables,
-    NULL,
-    clinical_infections,
-    treated,
+    42,
     infections,
-    parameters,
-    42 
+    treated,
+    clinical_infections
   )
 
   actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector()

From 6e893d2cd7dc38dd357b3e416d11429fe3f0f9bc Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Fri, 28 Jun 2024 16:40:21 +0100
Subject: [PATCH 04/14] p.v allows subpatent infections to occur, so we now
 calculate patent/subpatent infections, render these and schedule them.

---
 R/human_infection.R                         | 144 ++++++++++++++++----
 tests/testthat/test-infection-integration.R |   2 +
 tests/testthat/test-vivax.R                 |  65 +++++++++
 3 files changed, 183 insertions(+), 28 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index d4450cf2..6540e7ff 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -180,6 +180,25 @@ infection_outcome_process <- function(
         timestep,
         parameters$ud
       )
+      
+      clinical_infections <- calculate_clinical_infections(
+        variables,
+        infected_humans,
+        parameters,
+        renderer,
+        timestep
+      )
+      
+      update_severe_disease(
+        timestep,
+        infected_humans,
+        variables,
+        parameters,
+        renderer
+      )
+      
+      patent_infections <- NULL
+      
     } else if (parameters$parasite == "vivax"){
       boost_immunity(
         variables$iaa,
@@ -188,27 +207,27 @@ infection_outcome_process <- function(
         timestep,
         parameters$ua
       )
+      
+      ## Only S and U infections need to be split using the patent infection function
+      patent_infections <- calculate_patent_infections(
+        variables,
+        variables$state$get_index_of(c("S","U"))$and(infected_humans),
+        parameters,
+        renderer,
+        timestep
+      )
+      
+      # Patent level infected S and U, and all A infections to get clinical infections
+      clinical_infections <- calculate_clinical_infections(
+        variables,
+        variables$state$get_index_of("A")$and(infected_humans)$or(patent_infections),
+        parameters,
+        renderer,
+        timestep
+      )
     }
   }
   
-  clinical_infections <- calculate_clinical_infections(
-    variables,
-    infected_humans,
-    parameters,
-    renderer,
-    timestep
-  )
-  
-  if(parameters$parasite == "falciparum"){
-    update_severe_disease(
-      timestep,
-      infected_humans,
-      variables,
-      parameters,
-      renderer
-    )
-  }
-  
   treated <- calculate_treated(
     variables,
     clinical_infections,
@@ -221,6 +240,7 @@ infection_outcome_process <- function(
   
   schedule_infections(
     variables,
+    patent_infections,
     clinical_infections,
     treated,
     infected_humans,
@@ -229,6 +249,53 @@ infection_outcome_process <- function(
   )
 }
 
+#' @title Calculate patent infections (p.v only)
+#' @description
+#' Sample patent infections from all infections
+#' @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_patent_infections <- function(
+    variables,
+    infections,
+    parameters,
+    renderer,
+    timestep
+) {
+  
+  iaa <- variables$iaa$get_values(infections)
+  iam <- variables$iam$get_values(infections)
+  
+  philm <- anti_parasite_immunity(
+    min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
+    k = parameters$klm, iaa = iaa, iam = iam)
+  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    patent_infections,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  incidence_probability_renderer(
+    variables$birth,
+    renderer,
+    infections,
+    philm,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  patent_infections
+}
+
 #' @title Calculate clinical infections
 #' @description
 #' Sample clinical infections from all infections
@@ -459,20 +526,18 @@ calculate_treated <- function(
 #' @noRd
 schedule_infections <- function(
     variables,
+    patent_infections,
     clinical_infections,
     treated,
     infections,
     parameters,
     timestep
 ) {
+  
   included <- treated$not(TRUE)
+  to_infect_clinical <- clinical_infections$and(included)
 
-  to_infect <- clinical_infections$and(included)
-  to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(
-    included
-  )
-
-  if(to_infect$size() > 0) {
+  if(to_infect_clinical$size() > 0) {
     update_infection(
       variables$state,
       'D',
@@ -480,12 +545,14 @@ schedule_infections <- function(
       parameters$cd,
       variables$recovery_rates,
       1/parameters$dd,
-      to_infect
+      to_infect_clinical
     )
   }
 
-  if(to_infect_asym$size() > 0) {
-    if(parameters$parasite == "falciparum"){
+  if(parameters$parasite == "falciparum"){
+    to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(included)
+    
+    if(to_infect_asym$size() > 0){
       # p.f has immunity-determined asymptomatic infectivity
       update_to_asymptomatic_infection(
         variables,
@@ -493,7 +560,12 @@ schedule_infections <- function(
         timestep,
         to_infect_asym
       )
-    } else if (parameters$parasite == "vivax"){
+    }
+  } else if (parameters$parasite == "vivax"){
+    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(patent_infections$not(FALSE))
+    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(patent_infections)$and(clinical_infections$not(FALSE))
+    
+    if(to_infect_asym$size() > 0){
       # p.v has constant asymptomatic infectivity
       update_infection(
         variables$state,
@@ -505,6 +577,22 @@ schedule_infections <- function(
         to_infect_asym
       )
     }
+    if(to_infect_subpatent$size() > 0){
+      # p.v subpatent recovery rate is immunity dependent
+      update_infection(
+        variables$state,
+        'U',
+        variables$infectivity,
+        parameters$cu,
+        variables$recovery_rates,
+        1/anti_parasite_immunity(
+          parameters$dpcr_min, parameters$dpcr_max, parameters$apcr50, parameters$kpcr,
+          variables$iaa$get_values(to_infect_subpatent),
+          variables$iam$get_values(to_infect_subpatent)
+        ),
+        to_infect_subpatent
+      )
+    }
   }
 }
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index d399a527..17d0e958 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -157,6 +157,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     schedule_mock,
     1,
     variables,
+    NULL,
     clinical,
     treated,
     infected,
@@ -659,6 +660,7 @@ test_that('schedule_infections correctly schedules new infections', {
   
   schedule_infections(
     variables,
+    NULL,
     clinical_infections,
     treated,
     infections,
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 47427d99..8dc62d26 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -174,3 +174,68 @@ test_that('that vivax patent prevalence rendering works', {
   )
   
 })
+
+test_that('Test default vivax incidence rendering works', {
+  
+  timestep <- 0
+  year <- 365
+  birth <- individual::IntegerVariable$new(
+    -c(2, 5, 10, 11) * year
+  )
+  vivax_parameters <- get_parameters(
+    parasite = "vivax")
+  
+  renderer <- mock_render(1)
+  incidence_renderer(
+    birth,
+    renderer,
+    individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    'inc_patent_',
+    c(0, 2) * year,
+    c(5, 10) * year,
+    timestep
+  )
+  
+  incidence_probability_renderer(
+    birth,
+    renderer,
+    individual::Bitset$new(4)$insert(seq(4)),
+    c(.1, .2, .3, .4),
+    'inc_patent_',
+    c(0, 2) * year,
+    c(5, 10) * year,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_inc_patent_0_1825',
+    2,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    2,
+    'n_inc_patent_730_3650',
+    2,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    3,
+    'p_inc_patent_0_1825',
+    0.3,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    4,
+    'p_inc_patent_730_3650',
+    .6,
+    timestep
+  )
+})

From 8755a92eecc15c9e50274af19810b9fe98225b3c Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Thu, 11 Jul 2024 08:58:43 +0100
Subject: [PATCH 05/14] I've renamed patent infections as light microscopy
 detectable infections or lm_det, to avoid the confusion of what exactly
 patent means.

I've also removed the patent/lm_det infections <- NULL line from the falciparum code, but in doing so, have replicated some infection functions for falciparum and vivax so that schedule_infections may (for vivax) or may not (for falciparum) use the lm_det subset. I've also reordered the arguments for that function so that the lm_det_infections can be at the end, but make sense in the context of other arguments.
---
 R/human_infection.R                         | 117 ++++++++++++--------
 R/parameters.R                              |  14 +--
 tests/testthat/test-infection-integration.R |  18 ++-
 3 files changed, 84 insertions(+), 65 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index 6540e7ff..56fe0340 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -197,7 +197,24 @@ infection_outcome_process <- function(
         renderer
       )
       
-      patent_infections <- NULL
+      treated <- calculate_treated(
+        variables,
+        clinical_infections,
+        parameters,
+        timestep,
+        renderer
+      )
+      
+      renderer$render('n_infections', infected_humans$size(), timestep)
+      
+      schedule_infections(
+        parameters,
+        variables,
+        timestep,
+        infected_humans,
+        treated,
+        clinical_infections
+      )
       
     } else if (parameters$parasite == "vivax"){
       boost_immunity(
@@ -208,8 +225,8 @@ infection_outcome_process <- function(
         parameters$ua
       )
       
-      ## Only S and U infections need to be split using the patent infection function
-      patent_infections <- calculate_patent_infections(
+      ## Only S and U infections need to be split using the lm-detectable infection function
+      lm_det_infections <- calculate_lm_det_infections(
         variables,
         variables$state$get_index_of(c("S","U"))$and(infected_humans),
         parameters,
@@ -217,48 +234,48 @@ infection_outcome_process <- function(
         timestep
       )
       
-      # Patent level infected S and U, and all A infections to get clinical infections
+      # Lm-detectable level infected S and U, and all A infections may receive clinical infections
       clinical_infections <- calculate_clinical_infections(
         variables,
-        variables$state$get_index_of("A")$and(infected_humans)$or(patent_infections),
+        variables$state$get_index_of("A")$and(infected_humans)$or(lm_det_infections),
         parameters,
         renderer,
         timestep
       )
+      
+      treated <- calculate_treated(
+        variables,
+        clinical_infections,
+        parameters,
+        timestep,
+        renderer
+      )
+      
+      renderer$render('n_infections', infected_humans$size(), timestep)
+      
+      schedule_infections(
+        parameters,
+        variables,
+        timestep,
+        infected_humans,
+        treated,
+        clinical_infections,
+        lm_det_infections
+      )
     }
   }
-  
-  treated <- calculate_treated(
-    variables,
-    clinical_infections,
-    parameters,
-    timestep,
-    renderer
-  )
-  
-  renderer$render('n_infections', infected_humans$size(), timestep)
-  
-  schedule_infections(
-    variables,
-    patent_infections,
-    clinical_infections,
-    treated,
-    infected_humans,
-    parameters,
-    timestep
-  )
 }
 
-#' @title Calculate patent infections (p.v only)
+#' @title Calculate light microscopy detectable infections (p.v only)
 #' @description
-#' Sample patent infections from all infections
+#' Sample light microscopy detectable infections from all infections
 #' @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_patent_infections <- function(
+calculate_lm_det_infections <- function(
     variables,
     infections,
     parameters,
@@ -272,15 +289,15 @@ calculate_patent_infections <- function(
   philm <- anti_parasite_immunity(
     min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
     k = parameters$klm, iaa = iaa, iam = iam)
-  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  lm_det_infections <- bitset_at(infections, bernoulli_multi_p(philm))
   
   incidence_renderer(
     variables$birth,
     renderer,
-    patent_infections,
-    'inc_patent_',
-    parameters$patent_incidence_rendering_min_ages,
-    parameters$patent_incidence_rendering_max_ages,
+    lm_det_infections,
+    'inc_lm_det_',
+    parameters$lm_det_incidence_rendering_min_ages,
+    parameters$lm_det_incidence_rendering_max_ages,
     timestep
   )
   incidence_probability_renderer(
@@ -288,12 +305,12 @@ calculate_patent_infections <- function(
     renderer,
     infections,
     philm,
-    'inc_patent_',
-    parameters$patent_incidence_rendering_min_ages,
-    parameters$patent_incidence_rendering_max_ages,
+    'inc_lm_det_',
+    parameters$lm_det_incidence_rendering_min_ages,
+    parameters$lm_det_incidence_rendering_max_ages,
     timestep
   )
-  patent_infections
+  lm_det_infections
 }
 
 #' @title Calculate clinical infections
@@ -518,20 +535,24 @@ calculate_treated <- function(
 #' @title Schedule infections
 #' @description
 #' Schedule infections in humans after the incubation period
-#' @param events a list of all of the model events
-#' @param clinical_infections bitset of clinically infected humans
-#' @param treated bitset of treated humans
-#' @param infections bitset of infected humans
 #' @param parameters model parameters
+#' @param variables a list of all of the model variables
+#' @param timestep current timestep
+#' @param infections bitset of infected humans
+#' @param treated bitset of treated humans
+#' @param clinical_infections bitset of clinically infected humans
+#' @param lm_det_infections bitset of lm-detectable infected humans (p.v only:
+#'  lm_det infections are modelled as a human state, rather than a subset of  
+#'  asymptomatic infections as in the p.f model)
 #' @noRd
 schedule_infections <- function(
+    parameters,
     variables,
-    patent_infections,
-    clinical_infections,
-    treated,
+    timestep,
     infections,
-    parameters,
-    timestep
+    treated,
+    clinical_infections,
+    lm_det_infections = NULL
 ) {
   
   included <- treated$not(TRUE)
@@ -562,8 +583,8 @@ schedule_infections <- function(
       )
     }
   } else if (parameters$parasite == "vivax"){
-    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(patent_infections$not(FALSE))
-    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(patent_infections)$and(clinical_infections$not(FALSE))
+    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(lm_det_infections$not(FALSE))
+    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(lm_det_infections)$and(clinical_infections$not(FALSE))
     
     if(to_infect_asym$size() > 0){
       # p.v has constant asymptomatic infectivity
diff --git a/R/parameters.R b/R/parameters.R
index 95e7f776..5af67575 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -114,7 +114,7 @@
 #' * id0 - scale parameter; default = 1.577533
 #' * kd - shape parameter; default = 0.476614
 #'
-#' probability of patent infection due to anti-parasite immunity (p.v only):
+#' probability of light miscroscopy detectable infection due to anti-parasite immunity (p.v only):
 #'
 #' * phi0lm - maximum probability due to no immunity; default = 0.8918
 #' * phi1lm - maximum reduction due to immunity; default = 0.00482170890334156
@@ -143,7 +143,7 @@
 #' * cd - infectivity of clinically diseased humans towards mosquitoes; default = 0.068
 #' * ca - infectivity of asymptomatic humans towards mosquitoes (p.v only); default = 0.1
 #' * gamma1 - parameter for infectivity of asymptomatic humans; default = 1.82425
-#' * cu - infectivity of sub-patent infection; default = 0.0062
+#' * cu - infectivity of subpatent infection; default = 0.0062
 #' * ct - infectivity of treated infection; default = 0.021896
 #'
 #' mosquito fixed state transitions (including mortality):
@@ -280,8 +280,8 @@
 #' * incidence_rendering_min_ages - the minimum ages for incidence
 #' outputs (includes asymptomatic microscopy +); default = turned off
 #' * incidence_rendering_max_ages - the corresponding max ages; default = turned off
-#' * patent_incidence_rendering_min_ages - the minimum ages for patent incidence outputs (LM detectable), (p.v only); default = turned off
-#' * patent_incidence_rendering_max_ages - the corresponding max ages (p.v only); default = turned off
+#' * lm_det_incidence_rendering_min_ages - the minimum ages for light miscroscopy detectable incidence outputs, (p.v only); default = turned off
+#' * lm_det_incidence_rendering_max_ages - the corresponding max ages (p.v only); 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
 #' * severe_incidence_rendering_min_ages - the minimum ages for severe incidence
@@ -341,7 +341,7 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") {
     # maternal immunity parameters
     # probability of pre-erythrocytic infection/blood immunity
     # probability of asymptomatic detection (p.f only)
-    # probability of patent infection (due to anti-parasite immunity, p.v only)
+    # probability of l infection (due to anti-parasite immunity, p.v only)
     # probability of clinical infection
     # probability of severe infection (p.f only)
     # infectivity towards mosquitos
@@ -467,8 +467,8 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") {
       age_group_rendering_max_ages = numeric(0),
       incidence_rendering_min_ages = numeric(0),
       incidence_rendering_max_ages = numeric(0),
-      patent_incidence_rendering_min_ages = numeric(0),
-      patent_incidence_rendering_max_ages = numeric(0),
+      lm_det_incidence_rendering_min_ages = numeric(0),
+      lm_det_incidence_rendering_max_ages = numeric(0),
       clinical_incidence_rendering_min_ages = numeric(0),
       clinical_incidence_rendering_max_ages = numeric(0),
       severe_incidence_rendering_min_ages = numeric(0),
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 17d0e958..2d99c488 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -156,13 +156,12 @@ test_that('simulate_infection integrates different types of infection and schedu
   mockery::expect_args(
     schedule_mock,
     1,
+    parameters,
     variables,
-    NULL,
-    clinical,
-    treated,
+    timestep,
     infected,
-    parameters,
-    timestep
+    treated,
+    clinical
   )
 })
 
@@ -659,13 +658,12 @@ test_that('schedule_infections correctly schedules new infections', {
   mockery::stub(schedule_infections, 'update_to_asymptomatic_infection', asymp_mock)
   
   schedule_infections(
+    parameters,
     variables,
-    NULL,
-    clinical_infections,
-    treated,
+    42,
     infections,
-    parameters,
-    42 
+    treated,
+    clinical_infections
   )
 
   actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector()

From 246df82e05bda5b5834757bca7038acecaaef089 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 15 Jul 2024 12:24:06 +0100
Subject: [PATCH 06/14] Assigning new infections to new states has been
 rewritten such that bitsets for individuals to be moved are created in the
 infection_outcome_process function (to_D, to_A, and to_U), and less
 parasite-specific detail is required in schedule infections. Tests adjusted
 to reflect this change.

---
 R/human_infection.R                         | 123 ++++++++++----------
 tests/testthat/test-infection-integration.R |  27 +++--
 tests/testthat/test-vivax.R                 |  53 +++++++++
 3 files changed, 129 insertions(+), 74 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index 56fe0340..b2807e71 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -165,6 +165,9 @@ infection_outcome_process <- function(
   )
   
   if (infected_humans$size() > 0) {
+    
+    renderer$render('n_infections', infected_humans$size(), timestep)
+    
     boost_immunity(
       variables$ica,
       infected_humans,
@@ -172,6 +175,7 @@ infection_outcome_process <- function(
       timestep,
       parameters$uc
     )
+    
     if(parameters$parasite == "falciparum"){
       boost_immunity(
         variables$id,
@@ -181,7 +185,7 @@ infection_outcome_process <- function(
         parameters$ud
       )
       
-      clinical_infections <- calculate_clinical_infections(
+      clinical <- calculate_clinical_infections(
         variables,
         infected_humans,
         parameters,
@@ -189,34 +193,28 @@ infection_outcome_process <- function(
         timestep
       )
       
-      update_severe_disease(
-        timestep,
-        infected_humans,
-        variables,
-        parameters,
-        renderer
-      )
-      
       treated <- calculate_treated(
         variables,
-        clinical_infections,
+        clinical,
         parameters,
         timestep,
         renderer
       )
       
-      renderer$render('n_infections', infected_humans$size(), timestep)
-      
-      schedule_infections(
-        parameters,
-        variables,
+      update_severe_disease(
         timestep,
         infected_humans,
-        treated,
-        clinical_infections
+        variables,
+        parameters,
+        renderer
       )
       
+      to_D <- treated$not(FALSE)$and(clinical)
+      to_A <- infected_humans$and(clinical$not(FALSE))
+      to_U <- NULL
+      
     } else if (parameters$parasite == "vivax"){
+      
       boost_immunity(
         variables$iaa,
         infected_humans,
@@ -225,8 +223,8 @@ infection_outcome_process <- function(
         parameters$ua
       )
       
-      ## Only S and U infections need to be split using the lm-detectable infection function
-      lm_det_infections <- calculate_lm_det_infections(
+      ## Only S and U infections are considered in generating lm-det infections
+      lm_detectable <- calculate_lm_det_infections(
         variables,
         variables$state$get_index_of(c("S","U"))$and(infected_humans),
         parameters,
@@ -235,9 +233,13 @@ infection_outcome_process <- function(
       )
       
       # Lm-detectable level infected S and U, and all A infections may receive clinical infections
-      clinical_infections <- calculate_clinical_infections(
+      # There is a different calculation to generate clinical infections, based on current infection level
+      # LM infections must only pass through the clinical calculation, therefore all "A" infections are included
+      # "S" and "U" infections must pass through the lm-detectable calculation prior to and in addition to the clinical
+      # calculation. We therefore consider all "A" infections and only the "S" and "U" infections that are now lm-detectable.
+      clinical <- calculate_clinical_infections(
         variables,
-        variables$state$get_index_of("A")$and(infected_humans)$or(lm_det_infections),
+        variables$state$get_index_of("A")$and(infected_humans)$or(lm_detectable),
         parameters,
         renderer,
         timestep
@@ -245,24 +247,25 @@ infection_outcome_process <- function(
       
       treated <- calculate_treated(
         variables,
-        clinical_infections,
+        clinical,
         parameters,
         timestep,
         renderer
       )
       
-      renderer$render('n_infections', infected_humans$size(), timestep)
-      
-      schedule_infections(
-        parameters,
-        variables,
-        timestep,
-        infected_humans,
-        treated,
-        clinical_infections,
-        lm_det_infections
-      )
+      to_U <- infected_humans$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
+      to_A <- lm_detectable$and(clinical$not(F))
+      to_D <- clinical$and(treated$not(F))
     }
+    
+    schedule_infections(
+      parameters,
+      variables,
+      timestep,
+      to_D,
+      to_A,
+      to_U
+    )
   }
 }
 
@@ -532,33 +535,28 @@ calculate_treated <- function(
   
 }
 
+
 #' @title Schedule infections
 #' @description
 #' Schedule infections in humans after the incubation period
 #' @param parameters model parameters
 #' @param variables a list of all of the model variables
 #' @param timestep current timestep
-#' @param infections bitset of infected humans
-#' @param treated bitset of treated humans
-#' @param clinical_infections bitset of clinically infected humans
-#' @param lm_det_infections bitset of lm-detectable infected humans (p.v only:
-#'  lm_det infections are modelled as a human state, rather than a subset of  
-#'  asymptomatic infections as in the p.f model)
+#' @param to_D bitset of humans to move to state D
+#' @param to_A bitset of humans to move to state A
+#' @param to_U bitset of humans to move to state U
+#' @param to_T bitset of humans to move to state T
 #' @noRd
 schedule_infections <- function(
     parameters,
     variables,
     timestep,
-    infections,
-    treated,
-    clinical_infections,
-    lm_det_infections = NULL
+    to_D,
+    to_A,
+    to_U
 ) {
   
-  included <- treated$not(TRUE)
-  to_infect_clinical <- clinical_infections$and(included)
-
-  if(to_infect_clinical$size() > 0) {
+  if(to_D$size() > 0) {
     update_infection(
       variables$state,
       'D',
@@ -566,27 +564,20 @@ schedule_infections <- function(
       parameters$cd,
       variables$recovery_rates,
       1/parameters$dd,
-      to_infect_clinical
+      to_D
     )
   }
 
-  if(parameters$parasite == "falciparum"){
-    to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(included)
-    
-    if(to_infect_asym$size() > 0){
+  if(to_A$size() > 0) {
+    if(parameters$parasite == "falciparum"){
       # p.f has immunity-determined asymptomatic infectivity
       update_to_asymptomatic_infection(
         variables,
         parameters,
         timestep,
-        to_infect_asym
+        to_A
       )
-    }
-  } else if (parameters$parasite == "vivax"){
-    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(lm_det_infections$not(FALSE))
-    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(lm_det_infections)$and(clinical_infections$not(FALSE))
-    
-    if(to_infect_asym$size() > 0){
+    } else if (parameters$parasite == "vivax"){
       # p.v has constant asymptomatic infectivity
       update_infection(
         variables$state,
@@ -595,10 +586,14 @@ schedule_infections <- function(
         parameters$ca,
         variables$recovery_rates,
         1/parameters$da,
-        to_infect_asym
+        to_A
       )
     }
-    if(to_infect_subpatent$size() > 0){
+  }
+  
+  if(parameters$parasite == "vivax"){
+    # new p.v infections may be subpatent
+    if(to_U$size() > 0){
       # p.v subpatent recovery rate is immunity dependent
       update_infection(
         variables$state,
@@ -608,10 +603,10 @@ schedule_infections <- function(
         variables$recovery_rates,
         1/anti_parasite_immunity(
           parameters$dpcr_min, parameters$dpcr_max, parameters$apcr50, parameters$kpcr,
-          variables$iaa$get_values(to_infect_subpatent),
-          variables$iam$get_values(to_infect_subpatent)
+          variables$iaa$get_values(to_U),
+          variables$iam$get_values(to_U)
         ),
-        to_infect_subpatent
+        to_U
       )
     }
   }
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 2d99c488..bf72fedd 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -95,8 +95,11 @@ test_that('simulate_infection integrates different types of infection and schedu
   treated <- individual::Bitset$new(population)$insert(3)
   treated_mock <- mockery::mock(treated)
   schedule_mock <- mockery::mock()
-  
-  
+
+  to_D <- treated$not(FALSE)$and(clinical)
+  to_A <- infected$and(clinical$not(FALSE))
+  to_U <- NULL
+
   mockery::stub(infection_outcome_process, 'incidence_renderer', mockery::mock())
   mockery::stub(infection_outcome_process, 'boost_immunity', boost_immunity_mock)
   mockery::stub(infection_outcome_process, 'calculate_clinical_infections', clinical_infection_mock)
@@ -111,8 +114,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     renderer,
     parameters,
     prob)
-  
-  
+
   mockery::expect_args(
     boost_immunity_mock,
     1,
@@ -143,6 +145,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     renderer
   )
   
+  mockery::mock_args(treated_mock)[[1]][[2]]$to_vector()
   mockery::expect_args(
     treated_mock,
     1,
@@ -159,9 +162,9 @@ test_that('simulate_infection integrates different types of infection and schedu
     parameters,
     variables,
     timestep,
-    infected,
-    treated,
-    clinical
+    to_D,
+    to_A,
+    to_U
   )
 })
 
@@ -654,6 +657,10 @@ test_that('schedule_infections correctly schedules new infections', {
   infection_mock <- mockery::mock()
   asymp_mock <- mockery::mock()
   
+  to_D <- treated$not(FALSE)$and(clinical_infections)
+  to_A <- clinical_infections$not(FALSE)$and(infections)
+  to_U <- NULL
+  
   mockery::stub(schedule_infections, 'update_infection', infection_mock)
   mockery::stub(schedule_infections, 'update_to_asymptomatic_infection', asymp_mock)
   
@@ -661,9 +668,9 @@ test_that('schedule_infections correctly schedules new infections', {
     parameters,
     variables,
     42,
-    infections,
-    treated,
-    clinical_infections
+    to_D,
+    to_A,
+    to_U
   )
 
   actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector()
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 8dc62d26..37fa0d74 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -239,3 +239,56 @@ test_that('Test default vivax incidence rendering works', {
     timestep
   )
 })
+
+test_that('vivax schedule_infections correctly schedules new infections', {
+  parameters <- get_parameters(list(human_population = 20), parasite = "vivax")
+  variables <- create_variables(parameters)
+  
+  variables$state <- individual::CategoricalVariable$new(
+    c('U', 'A', 'D', 'S', 'Tr'),
+    rep(c('S','U','A','D','Tr'), times = 4)
+  )
+  
+  infections <- individual::Bitset$new(20)$insert(1:20)
+  lm_detectable <- individual::Bitset$new(20)$insert(6:20)$and(variables$state$get_index_of(c("S", "U")))
+  clinical <- individual::Bitset$new(20)$insert(11:20)$and(variables$state$get_index_of(c("A"))$or(lm_detectable))
+  treated <- individual::Bitset$new(20)$insert(16:20)$and(clinical)
+  # Only S can be a subpatent infection (1)
+  # Only S and U can be a patent infection (6, 7)
+  # S, U and A can be clinical infections (11, 12, 13), but the model re-infects everyone
+  # Treated only looks at new clinical infections (from SAU, not from D)
+  to_U <- infections$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
+  to_A <- lm_detectable$and(clinical$not(F))
+  to_D <- clinical$and(treated$not(F))
+  
+  infection_mock <- mockery::mock()
+  mockery::stub(schedule_infections, 'update_infection', infection_mock)
+  
+  schedule_infections(
+    parameters,
+    variables,
+    timestep = 42,
+    to_D,
+    to_A,
+    to_U
+  )
+  
+  actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector()
+  actual_asymp_infected <- mockery::mock_args(infection_mock)[[2]][[7]]$to_vector()
+  actual_subpatent_infected <- mockery::mock_args(infection_mock)[[3]][[7]]$to_vector()
+  
+  expect_equal(
+    actual_infected,
+    c(11:13)
+  )
+  
+  expect_equal(
+    actual_asymp_infected,
+    c(6, 7)
+  )
+  
+  expect_equal(
+    actual_subpatent_infected,
+    c(1)
+  )
+})

From 64c9ad76748a079a1474c022a7d24f3c81269450 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Tue, 16 Jul 2024 14:24:55 +0100
Subject: [PATCH 07/14] Added a warning to human_infections to notify
 developers that bitsets are rewritten and subsetted in new state assignments.

---
 R/human_infection.R | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/R/human_infection.R b/R/human_infection.R
index b2807e71..57d301f4 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -209,6 +209,7 @@ infection_outcome_process <- function(
         renderer
       )
       
+      ## The treated and infected_humans bitsets are re-written so be cautious!
       to_D <- treated$not(FALSE)$and(clinical)
       to_A <- infected_humans$and(clinical$not(FALSE))
       to_U <- NULL
@@ -253,6 +254,7 @@ infection_outcome_process <- function(
         renderer
       )
       
+      ## The infected_humans,lm_detectable and clinical bitsets are re-written so be cautious!
       to_U <- infected_humans$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
       to_A <- lm_detectable$and(clinical$not(F))
       to_D <- clinical$and(treated$not(F))

From 220f00f17b4ecbd49cc86b328f3ab785978de463 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Fri, 28 Jun 2024 16:40:21 +0100
Subject: [PATCH 08/14] p.v allows subpatent infections to occur, so we now
 calculate patent/subpatent infections, render these and schedule them.

---
 R/human_infection.R                         | 144 ++++++++++++++++----
 tests/testthat/test-infection-integration.R |   2 +
 tests/testthat/test-vivax.R                 |  65 +++++++++
 3 files changed, 183 insertions(+), 28 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index d4450cf2..6540e7ff 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -180,6 +180,25 @@ infection_outcome_process <- function(
         timestep,
         parameters$ud
       )
+      
+      clinical_infections <- calculate_clinical_infections(
+        variables,
+        infected_humans,
+        parameters,
+        renderer,
+        timestep
+      )
+      
+      update_severe_disease(
+        timestep,
+        infected_humans,
+        variables,
+        parameters,
+        renderer
+      )
+      
+      patent_infections <- NULL
+      
     } else if (parameters$parasite == "vivax"){
       boost_immunity(
         variables$iaa,
@@ -188,27 +207,27 @@ infection_outcome_process <- function(
         timestep,
         parameters$ua
       )
+      
+      ## Only S and U infections need to be split using the patent infection function
+      patent_infections <- calculate_patent_infections(
+        variables,
+        variables$state$get_index_of(c("S","U"))$and(infected_humans),
+        parameters,
+        renderer,
+        timestep
+      )
+      
+      # Patent level infected S and U, and all A infections to get clinical infections
+      clinical_infections <- calculate_clinical_infections(
+        variables,
+        variables$state$get_index_of("A")$and(infected_humans)$or(patent_infections),
+        parameters,
+        renderer,
+        timestep
+      )
     }
   }
   
-  clinical_infections <- calculate_clinical_infections(
-    variables,
-    infected_humans,
-    parameters,
-    renderer,
-    timestep
-  )
-  
-  if(parameters$parasite == "falciparum"){
-    update_severe_disease(
-      timestep,
-      infected_humans,
-      variables,
-      parameters,
-      renderer
-    )
-  }
-  
   treated <- calculate_treated(
     variables,
     clinical_infections,
@@ -221,6 +240,7 @@ infection_outcome_process <- function(
   
   schedule_infections(
     variables,
+    patent_infections,
     clinical_infections,
     treated,
     infected_humans,
@@ -229,6 +249,53 @@ infection_outcome_process <- function(
   )
 }
 
+#' @title Calculate patent infections (p.v only)
+#' @description
+#' Sample patent infections from all infections
+#' @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_patent_infections <- function(
+    variables,
+    infections,
+    parameters,
+    renderer,
+    timestep
+) {
+  
+  iaa <- variables$iaa$get_values(infections)
+  iam <- variables$iam$get_values(infections)
+  
+  philm <- anti_parasite_immunity(
+    min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
+    k = parameters$klm, iaa = iaa, iam = iam)
+  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    patent_infections,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  incidence_probability_renderer(
+    variables$birth,
+    renderer,
+    infections,
+    philm,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  patent_infections
+}
+
 #' @title Calculate clinical infections
 #' @description
 #' Sample clinical infections from all infections
@@ -459,20 +526,18 @@ calculate_treated <- function(
 #' @noRd
 schedule_infections <- function(
     variables,
+    patent_infections,
     clinical_infections,
     treated,
     infections,
     parameters,
     timestep
 ) {
+  
   included <- treated$not(TRUE)
+  to_infect_clinical <- clinical_infections$and(included)
 
-  to_infect <- clinical_infections$and(included)
-  to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(
-    included
-  )
-
-  if(to_infect$size() > 0) {
+  if(to_infect_clinical$size() > 0) {
     update_infection(
       variables$state,
       'D',
@@ -480,12 +545,14 @@ schedule_infections <- function(
       parameters$cd,
       variables$recovery_rates,
       1/parameters$dd,
-      to_infect
+      to_infect_clinical
     )
   }
 
-  if(to_infect_asym$size() > 0) {
-    if(parameters$parasite == "falciparum"){
+  if(parameters$parasite == "falciparum"){
+    to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(included)
+    
+    if(to_infect_asym$size() > 0){
       # p.f has immunity-determined asymptomatic infectivity
       update_to_asymptomatic_infection(
         variables,
@@ -493,7 +560,12 @@ schedule_infections <- function(
         timestep,
         to_infect_asym
       )
-    } else if (parameters$parasite == "vivax"){
+    }
+  } else if (parameters$parasite == "vivax"){
+    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(patent_infections$not(FALSE))
+    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(patent_infections)$and(clinical_infections$not(FALSE))
+    
+    if(to_infect_asym$size() > 0){
       # p.v has constant asymptomatic infectivity
       update_infection(
         variables$state,
@@ -505,6 +577,22 @@ schedule_infections <- function(
         to_infect_asym
       )
     }
+    if(to_infect_subpatent$size() > 0){
+      # p.v subpatent recovery rate is immunity dependent
+      update_infection(
+        variables$state,
+        'U',
+        variables$infectivity,
+        parameters$cu,
+        variables$recovery_rates,
+        1/anti_parasite_immunity(
+          parameters$dpcr_min, parameters$dpcr_max, parameters$apcr50, parameters$kpcr,
+          variables$iaa$get_values(to_infect_subpatent),
+          variables$iam$get_values(to_infect_subpatent)
+        ),
+        to_infect_subpatent
+      )
+    }
   }
 }
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index d399a527..17d0e958 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -157,6 +157,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     schedule_mock,
     1,
     variables,
+    NULL,
     clinical,
     treated,
     infected,
@@ -659,6 +660,7 @@ test_that('schedule_infections correctly schedules new infections', {
   
   schedule_infections(
     variables,
+    NULL,
     clinical_infections,
     treated,
     infections,
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index db1edf19..80f6fa31 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -137,3 +137,68 @@ test_that('that vivax patent prevalence rendering works', {
   )
   
 })
+
+test_that('Test default vivax incidence rendering works', {
+  
+  timestep <- 0
+  year <- 365
+  birth <- individual::IntegerVariable$new(
+    -c(2, 5, 10, 11) * year
+  )
+  vivax_parameters <- get_parameters(
+    parasite = "vivax")
+  
+  renderer <- mock_render(1)
+  incidence_renderer(
+    birth,
+    renderer,
+    individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    'inc_patent_',
+    c(0, 2) * year,
+    c(5, 10) * year,
+    timestep
+  )
+  
+  incidence_probability_renderer(
+    birth,
+    renderer,
+    individual::Bitset$new(4)$insert(seq(4)),
+    c(.1, .2, .3, .4),
+    'inc_patent_',
+    c(0, 2) * year,
+    c(5, 10) * year,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_inc_patent_0_1825',
+    2,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    2,
+    'n_inc_patent_730_3650',
+    2,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    3,
+    'p_inc_patent_0_1825',
+    0.3,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    4,
+    'p_inc_patent_730_3650',
+    .6,
+    timestep
+  )
+})

From a75c8bc9fd7c12df5be34c250e9de57692609fb2 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Thu, 11 Jul 2024 08:58:43 +0100
Subject: [PATCH 09/14] I've renamed patent infections as light microscopy
 detectable infections or lm_det, to avoid the confusion of what exactly
 patent means.

I've also removed the patent/lm_det infections <- NULL line from the falciparum code, but in doing so, have replicated some infection functions for falciparum and vivax so that schedule_infections may (for vivax) or may not (for falciparum) use the lm_det subset. I've also reordered the arguments for that function so that the lm_det_infections can be at the end, but make sense in the context of other arguments.
---
 R/human_infection.R                         | 117 ++++++++++++--------
 R/parameters.R                              |  14 +--
 tests/testthat/test-infection-integration.R |  18 ++-
 3 files changed, 84 insertions(+), 65 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index 6540e7ff..56fe0340 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -197,7 +197,24 @@ infection_outcome_process <- function(
         renderer
       )
       
-      patent_infections <- NULL
+      treated <- calculate_treated(
+        variables,
+        clinical_infections,
+        parameters,
+        timestep,
+        renderer
+      )
+      
+      renderer$render('n_infections', infected_humans$size(), timestep)
+      
+      schedule_infections(
+        parameters,
+        variables,
+        timestep,
+        infected_humans,
+        treated,
+        clinical_infections
+      )
       
     } else if (parameters$parasite == "vivax"){
       boost_immunity(
@@ -208,8 +225,8 @@ infection_outcome_process <- function(
         parameters$ua
       )
       
-      ## Only S and U infections need to be split using the patent infection function
-      patent_infections <- calculate_patent_infections(
+      ## Only S and U infections need to be split using the lm-detectable infection function
+      lm_det_infections <- calculate_lm_det_infections(
         variables,
         variables$state$get_index_of(c("S","U"))$and(infected_humans),
         parameters,
@@ -217,48 +234,48 @@ infection_outcome_process <- function(
         timestep
       )
       
-      # Patent level infected S and U, and all A infections to get clinical infections
+      # Lm-detectable level infected S and U, and all A infections may receive clinical infections
       clinical_infections <- calculate_clinical_infections(
         variables,
-        variables$state$get_index_of("A")$and(infected_humans)$or(patent_infections),
+        variables$state$get_index_of("A")$and(infected_humans)$or(lm_det_infections),
         parameters,
         renderer,
         timestep
       )
+      
+      treated <- calculate_treated(
+        variables,
+        clinical_infections,
+        parameters,
+        timestep,
+        renderer
+      )
+      
+      renderer$render('n_infections', infected_humans$size(), timestep)
+      
+      schedule_infections(
+        parameters,
+        variables,
+        timestep,
+        infected_humans,
+        treated,
+        clinical_infections,
+        lm_det_infections
+      )
     }
   }
-  
-  treated <- calculate_treated(
-    variables,
-    clinical_infections,
-    parameters,
-    timestep,
-    renderer
-  )
-  
-  renderer$render('n_infections', infected_humans$size(), timestep)
-  
-  schedule_infections(
-    variables,
-    patent_infections,
-    clinical_infections,
-    treated,
-    infected_humans,
-    parameters,
-    timestep
-  )
 }
 
-#' @title Calculate patent infections (p.v only)
+#' @title Calculate light microscopy detectable infections (p.v only)
 #' @description
-#' Sample patent infections from all infections
+#' Sample light microscopy detectable infections from all infections
 #' @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_patent_infections <- function(
+calculate_lm_det_infections <- function(
     variables,
     infections,
     parameters,
@@ -272,15 +289,15 @@ calculate_patent_infections <- function(
   philm <- anti_parasite_immunity(
     min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
     k = parameters$klm, iaa = iaa, iam = iam)
-  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  lm_det_infections <- bitset_at(infections, bernoulli_multi_p(philm))
   
   incidence_renderer(
     variables$birth,
     renderer,
-    patent_infections,
-    'inc_patent_',
-    parameters$patent_incidence_rendering_min_ages,
-    parameters$patent_incidence_rendering_max_ages,
+    lm_det_infections,
+    'inc_lm_det_',
+    parameters$lm_det_incidence_rendering_min_ages,
+    parameters$lm_det_incidence_rendering_max_ages,
     timestep
   )
   incidence_probability_renderer(
@@ -288,12 +305,12 @@ calculate_patent_infections <- function(
     renderer,
     infections,
     philm,
-    'inc_patent_',
-    parameters$patent_incidence_rendering_min_ages,
-    parameters$patent_incidence_rendering_max_ages,
+    'inc_lm_det_',
+    parameters$lm_det_incidence_rendering_min_ages,
+    parameters$lm_det_incidence_rendering_max_ages,
     timestep
   )
-  patent_infections
+  lm_det_infections
 }
 
 #' @title Calculate clinical infections
@@ -518,20 +535,24 @@ calculate_treated <- function(
 #' @title Schedule infections
 #' @description
 #' Schedule infections in humans after the incubation period
-#' @param events a list of all of the model events
-#' @param clinical_infections bitset of clinically infected humans
-#' @param treated bitset of treated humans
-#' @param infections bitset of infected humans
 #' @param parameters model parameters
+#' @param variables a list of all of the model variables
+#' @param timestep current timestep
+#' @param infections bitset of infected humans
+#' @param treated bitset of treated humans
+#' @param clinical_infections bitset of clinically infected humans
+#' @param lm_det_infections bitset of lm-detectable infected humans (p.v only:
+#'  lm_det infections are modelled as a human state, rather than a subset of  
+#'  asymptomatic infections as in the p.f model)
 #' @noRd
 schedule_infections <- function(
+    parameters,
     variables,
-    patent_infections,
-    clinical_infections,
-    treated,
+    timestep,
     infections,
-    parameters,
-    timestep
+    treated,
+    clinical_infections,
+    lm_det_infections = NULL
 ) {
   
   included <- treated$not(TRUE)
@@ -562,8 +583,8 @@ schedule_infections <- function(
       )
     }
   } else if (parameters$parasite == "vivax"){
-    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(patent_infections$not(FALSE))
-    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(patent_infections)$and(clinical_infections$not(FALSE))
+    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(lm_det_infections$not(FALSE))
+    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(lm_det_infections)$and(clinical_infections$not(FALSE))
     
     if(to_infect_asym$size() > 0){
       # p.v has constant asymptomatic infectivity
diff --git a/R/parameters.R b/R/parameters.R
index 95e7f776..5af67575 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -114,7 +114,7 @@
 #' * id0 - scale parameter; default = 1.577533
 #' * kd - shape parameter; default = 0.476614
 #'
-#' probability of patent infection due to anti-parasite immunity (p.v only):
+#' probability of light miscroscopy detectable infection due to anti-parasite immunity (p.v only):
 #'
 #' * phi0lm - maximum probability due to no immunity; default = 0.8918
 #' * phi1lm - maximum reduction due to immunity; default = 0.00482170890334156
@@ -143,7 +143,7 @@
 #' * cd - infectivity of clinically diseased humans towards mosquitoes; default = 0.068
 #' * ca - infectivity of asymptomatic humans towards mosquitoes (p.v only); default = 0.1
 #' * gamma1 - parameter for infectivity of asymptomatic humans; default = 1.82425
-#' * cu - infectivity of sub-patent infection; default = 0.0062
+#' * cu - infectivity of subpatent infection; default = 0.0062
 #' * ct - infectivity of treated infection; default = 0.021896
 #'
 #' mosquito fixed state transitions (including mortality):
@@ -280,8 +280,8 @@
 #' * incidence_rendering_min_ages - the minimum ages for incidence
 #' outputs (includes asymptomatic microscopy +); default = turned off
 #' * incidence_rendering_max_ages - the corresponding max ages; default = turned off
-#' * patent_incidence_rendering_min_ages - the minimum ages for patent incidence outputs (LM detectable), (p.v only); default = turned off
-#' * patent_incidence_rendering_max_ages - the corresponding max ages (p.v only); default = turned off
+#' * lm_det_incidence_rendering_min_ages - the minimum ages for light miscroscopy detectable incidence outputs, (p.v only); default = turned off
+#' * lm_det_incidence_rendering_max_ages - the corresponding max ages (p.v only); 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
 #' * severe_incidence_rendering_min_ages - the minimum ages for severe incidence
@@ -341,7 +341,7 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") {
     # maternal immunity parameters
     # probability of pre-erythrocytic infection/blood immunity
     # probability of asymptomatic detection (p.f only)
-    # probability of patent infection (due to anti-parasite immunity, p.v only)
+    # probability of l infection (due to anti-parasite immunity, p.v only)
     # probability of clinical infection
     # probability of severe infection (p.f only)
     # infectivity towards mosquitos
@@ -467,8 +467,8 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") {
       age_group_rendering_max_ages = numeric(0),
       incidence_rendering_min_ages = numeric(0),
       incidence_rendering_max_ages = numeric(0),
-      patent_incidence_rendering_min_ages = numeric(0),
-      patent_incidence_rendering_max_ages = numeric(0),
+      lm_det_incidence_rendering_min_ages = numeric(0),
+      lm_det_incidence_rendering_max_ages = numeric(0),
       clinical_incidence_rendering_min_ages = numeric(0),
       clinical_incidence_rendering_max_ages = numeric(0),
       severe_incidence_rendering_min_ages = numeric(0),
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 17d0e958..2d99c488 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -156,13 +156,12 @@ test_that('simulate_infection integrates different types of infection and schedu
   mockery::expect_args(
     schedule_mock,
     1,
+    parameters,
     variables,
-    NULL,
-    clinical,
-    treated,
+    timestep,
     infected,
-    parameters,
-    timestep
+    treated,
+    clinical
   )
 })
 
@@ -659,13 +658,12 @@ test_that('schedule_infections correctly schedules new infections', {
   mockery::stub(schedule_infections, 'update_to_asymptomatic_infection', asymp_mock)
   
   schedule_infections(
+    parameters,
     variables,
-    NULL,
-    clinical_infections,
-    treated,
+    42,
     infections,
-    parameters,
-    42 
+    treated,
+    clinical_infections
   )
 
   actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector()

From ccc60ec6ed3dd08bacab2e4288572ffea2a2c316 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 29 Jul 2024 15:32:07 +0100
Subject: [PATCH 10/14] Adding calculate patent infection function.

---
 R/human_infection.R | 47 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)

diff --git a/R/human_infection.R b/R/human_infection.R
index 56fe0340..e150a365 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -313,6 +313,53 @@ calculate_lm_det_infections <- function(
   lm_det_infections
 }
 
+#' @title Calculate patent infections (p.v only)
+#' @description
+#' Sample patent infections from all infections
+#' @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_patent_infections <- function(
+    variables,
+    infections,
+    parameters,
+    renderer,
+    timestep
+) {
+  
+  iaa <- variables$iaa$get_values(infections)
+  iam <- variables$iam$get_values(infections)
+  
+  philm <- anti_parasite_immunity(
+    min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
+    k = parameters$klm, iaa = iaa, iam = iam)
+  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    patent_infections,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  incidence_probability_renderer(
+    variables$birth,
+    renderer,
+    infections,
+    philm,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  patent_infections
+}
+
 #' @title Calculate clinical infections
 #' @description
 #' Sample clinical infections from all infections

From 9592c4b9b15882855f4a7d4d7c92fda0c99d637e Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Fri, 28 Jun 2024 16:40:21 +0100
Subject: [PATCH 11/14] p.v allows subpatent infections to occur, so we now
 calculate patent/subpatent infections, render these and schedule them.

---
 R/human_infection.R | 49 ++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 48 insertions(+), 1 deletion(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index e150a365..409cfe07 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -215,7 +215,7 @@ infection_outcome_process <- function(
         treated,
         clinical_infections
       )
-      
+
     } else if (parameters$parasite == "vivax"){
       boost_immunity(
         variables$iaa,
@@ -360,6 +360,53 @@ calculate_patent_infections <- function(
   patent_infections
 }
 
+#' @title Calculate patent infections (p.v only)
+#' @description
+#' Sample patent infections from all infections
+#' @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_patent_infections <- function(
+    variables,
+    infections,
+    parameters,
+    renderer,
+    timestep
+) {
+  
+  iaa <- variables$iaa$get_values(infections)
+  iam <- variables$iam$get_values(infections)
+  
+  philm <- anti_parasite_immunity(
+    min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
+    k = parameters$klm, iaa = iaa, iam = iam)
+  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    patent_infections,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  incidence_probability_renderer(
+    variables$birth,
+    renderer,
+    infections,
+    philm,
+    'inc_patent_',
+    parameters$patent_incidence_rendering_min_ages,
+    parameters$patent_incidence_rendering_max_ages,
+    timestep
+  )
+  patent_infections
+}
+
 #' @title Calculate clinical infections
 #' @description
 #' Sample clinical infections from all infections

From 7fe4709296086a30b3785c3ee31445d8613cd824 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Thu, 11 Jul 2024 08:58:43 +0100
Subject: [PATCH 12/14] I've renamed patent infections as light microscopy
 detectable infections or lm_det, to avoid the confusion of what exactly
 patent means.

I've also removed the patent/lm_det infections <- NULL line from the falciparum code, but in doing so, have replicated some infection functions for falciparum and vivax so that schedule_infections may (for vivax) or may not (for falciparum) use the lm_det subset. I've also reordered the arguments for that function so that the lm_det_infections can be at the end, but make sense in the context of other arguments.
---
 R/human_infection.R | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index 409cfe07..6195e7b7 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -360,16 +360,16 @@ calculate_patent_infections <- function(
   patent_infections
 }
 
-#' @title Calculate patent infections (p.v only)
+#' @title Calculate light microscopy detectable infections (p.v only)
 #' @description
-#' Sample patent infections from all infections
+#' Sample light microscopy detectable infections from all infections
 #' @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_patent_infections <- function(
+calculate_lm_det_infections <- function(
     variables,
     infections,
     parameters,
@@ -383,15 +383,15 @@ calculate_patent_infections <- function(
   philm <- anti_parasite_immunity(
     min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
     k = parameters$klm, iaa = iaa, iam = iam)
-  patent_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  lm_det_infections <- bitset_at(infections, bernoulli_multi_p(philm))
   
   incidence_renderer(
     variables$birth,
     renderer,
-    patent_infections,
-    'inc_patent_',
-    parameters$patent_incidence_rendering_min_ages,
-    parameters$patent_incidence_rendering_max_ages,
+    lm_det_infections,
+    'inc_lm_det_',
+    parameters$lm_det_incidence_rendering_min_ages,
+    parameters$lm_det_incidence_rendering_max_ages,
     timestep
   )
   incidence_probability_renderer(
@@ -399,12 +399,12 @@ calculate_patent_infections <- function(
     renderer,
     infections,
     philm,
-    'inc_patent_',
-    parameters$patent_incidence_rendering_min_ages,
-    parameters$patent_incidence_rendering_max_ages,
+    'inc_lm_det_',
+    parameters$lm_det_incidence_rendering_min_ages,
+    parameters$lm_det_incidence_rendering_max_ages,
     timestep
   )
-  patent_infections
+  lm_det_infections
 }
 
 #' @title Calculate clinical infections

From e8dffbd03d27d1fe251925d1791c86b9b56c53a8 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 15 Jul 2024 12:24:06 +0100
Subject: [PATCH 13/14] Assigning new infections to new states has been
 rewritten such that bitsets for individuals to be moved are created in the
 infection_outcome_process function (to_D, to_A, and to_U), and less
 parasite-specific detail is required in schedule infections. Tests adjusted
 to reflect this change.

---
 R/human_infection.R                         | 123 ++++++++++----------
 tests/testthat/test-infection-integration.R |  27 +++--
 tests/testthat/test-vivax.R                 |  53 +++++++++
 3 files changed, 129 insertions(+), 74 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index 6195e7b7..985844ac 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -165,6 +165,9 @@ infection_outcome_process <- function(
   )
   
   if (infected_humans$size() > 0) {
+    
+    renderer$render('n_infections', infected_humans$size(), timestep)
+    
     boost_immunity(
       variables$ica,
       infected_humans,
@@ -172,6 +175,7 @@ infection_outcome_process <- function(
       timestep,
       parameters$uc
     )
+    
     if(parameters$parasite == "falciparum"){
       boost_immunity(
         variables$id,
@@ -181,7 +185,7 @@ infection_outcome_process <- function(
         parameters$ud
       )
       
-      clinical_infections <- calculate_clinical_infections(
+      clinical <- calculate_clinical_infections(
         variables,
         infected_humans,
         parameters,
@@ -189,34 +193,28 @@ infection_outcome_process <- function(
         timestep
       )
       
-      update_severe_disease(
-        timestep,
-        infected_humans,
+      treated <- calculate_treated(
         variables,
+        clinical,
         parameters,
+        timestep,
         renderer
       )
       
-      treated <- calculate_treated(
+      update_severe_disease(
+        timestep,
+        infected_humans,
         variables,
-        clinical_infections,
         parameters,
-        timestep,
         renderer
       )
       
-      renderer$render('n_infections', infected_humans$size(), timestep)
+      to_D <- treated$not(FALSE)$and(clinical)
+      to_A <- infected_humans$and(clinical$not(FALSE))
+      to_U <- NULL
       
-      schedule_infections(
-        parameters,
-        variables,
-        timestep,
-        infected_humans,
-        treated,
-        clinical_infections
-      )
-
     } else if (parameters$parasite == "vivax"){
+      
       boost_immunity(
         variables$iaa,
         infected_humans,
@@ -225,8 +223,8 @@ infection_outcome_process <- function(
         parameters$ua
       )
       
-      ## Only S and U infections need to be split using the lm-detectable infection function
-      lm_det_infections <- calculate_lm_det_infections(
+      ## Only S and U infections are considered in generating lm-det infections
+      lm_detectable <- calculate_lm_det_infections(
         variables,
         variables$state$get_index_of(c("S","U"))$and(infected_humans),
         parameters,
@@ -235,9 +233,13 @@ infection_outcome_process <- function(
       )
       
       # Lm-detectable level infected S and U, and all A infections may receive clinical infections
-      clinical_infections <- calculate_clinical_infections(
+      # There is a different calculation to generate clinical infections, based on current infection level
+      # LM infections must only pass through the clinical calculation, therefore all "A" infections are included
+      # "S" and "U" infections must pass through the lm-detectable calculation prior to and in addition to the clinical
+      # calculation. We therefore consider all "A" infections and only the "S" and "U" infections that are now lm-detectable.
+      clinical <- calculate_clinical_infections(
         variables,
-        variables$state$get_index_of("A")$and(infected_humans)$or(lm_det_infections),
+        variables$state$get_index_of("A")$and(infected_humans)$or(lm_detectable),
         parameters,
         renderer,
         timestep
@@ -245,24 +247,25 @@ infection_outcome_process <- function(
       
       treated <- calculate_treated(
         variables,
-        clinical_infections,
+        clinical,
         parameters,
         timestep,
         renderer
       )
       
-      renderer$render('n_infections', infected_humans$size(), timestep)
-      
-      schedule_infections(
-        parameters,
-        variables,
-        timestep,
-        infected_humans,
-        treated,
-        clinical_infections,
-        lm_det_infections
-      )
+      to_U <- infected_humans$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
+      to_A <- lm_detectable$and(clinical$not(F))
+      to_D <- clinical$and(treated$not(F))
     }
+    
+    schedule_infections(
+      parameters,
+      variables,
+      timestep,
+      to_D,
+      to_A,
+      to_U
+    )
   }
 }
 
@@ -626,33 +629,28 @@ calculate_treated <- function(
   
 }
 
+
 #' @title Schedule infections
 #' @description
 #' Schedule infections in humans after the incubation period
 #' @param parameters model parameters
 #' @param variables a list of all of the model variables
 #' @param timestep current timestep
-#' @param infections bitset of infected humans
-#' @param treated bitset of treated humans
-#' @param clinical_infections bitset of clinically infected humans
-#' @param lm_det_infections bitset of lm-detectable infected humans (p.v only:
-#'  lm_det infections are modelled as a human state, rather than a subset of  
-#'  asymptomatic infections as in the p.f model)
+#' @param to_D bitset of humans to move to state D
+#' @param to_A bitset of humans to move to state A
+#' @param to_U bitset of humans to move to state U
+#' @param to_T bitset of humans to move to state T
 #' @noRd
 schedule_infections <- function(
     parameters,
     variables,
     timestep,
-    infections,
-    treated,
-    clinical_infections,
-    lm_det_infections = NULL
+    to_D,
+    to_A,
+    to_U
 ) {
   
-  included <- treated$not(TRUE)
-  to_infect_clinical <- clinical_infections$and(included)
-
-  if(to_infect_clinical$size() > 0) {
+  if(to_D$size() > 0) {
     update_infection(
       variables$state,
       'D',
@@ -660,27 +658,20 @@ schedule_infections <- function(
       parameters$cd,
       variables$recovery_rates,
       1/parameters$dd,
-      to_infect_clinical
+      to_D
     )
   }
 
-  if(parameters$parasite == "falciparum"){
-    to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(included)
-    
-    if(to_infect_asym$size() > 0){
+  if(to_A$size() > 0) {
+    if(parameters$parasite == "falciparum"){
       # p.f has immunity-determined asymptomatic infectivity
       update_to_asymptomatic_infection(
         variables,
         parameters,
         timestep,
-        to_infect_asym
+        to_A
       )
-    }
-  } else if (parameters$parasite == "vivax"){
-    to_infect_subpatent <- variables$state$get_index_of(c('S'))$and(included)$and(infections)$and(lm_det_infections$not(FALSE))
-    to_infect_asym <- variables$state$get_index_of(c('S',"U"))$and(included)$and(lm_det_infections)$and(clinical_infections$not(FALSE))
-    
-    if(to_infect_asym$size() > 0){
+    } else if (parameters$parasite == "vivax"){
       # p.v has constant asymptomatic infectivity
       update_infection(
         variables$state,
@@ -689,10 +680,14 @@ schedule_infections <- function(
         parameters$ca,
         variables$recovery_rates,
         1/parameters$da,
-        to_infect_asym
+        to_A
       )
     }
-    if(to_infect_subpatent$size() > 0){
+  }
+  
+  if(parameters$parasite == "vivax"){
+    # new p.v infections may be subpatent
+    if(to_U$size() > 0){
       # p.v subpatent recovery rate is immunity dependent
       update_infection(
         variables$state,
@@ -702,10 +697,10 @@ schedule_infections <- function(
         variables$recovery_rates,
         1/anti_parasite_immunity(
           parameters$dpcr_min, parameters$dpcr_max, parameters$apcr50, parameters$kpcr,
-          variables$iaa$get_values(to_infect_subpatent),
-          variables$iam$get_values(to_infect_subpatent)
+          variables$iaa$get_values(to_U),
+          variables$iam$get_values(to_U)
         ),
-        to_infect_subpatent
+        to_U
       )
     }
   }
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 2d99c488..bf72fedd 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -95,8 +95,11 @@ test_that('simulate_infection integrates different types of infection and schedu
   treated <- individual::Bitset$new(population)$insert(3)
   treated_mock <- mockery::mock(treated)
   schedule_mock <- mockery::mock()
-  
-  
+
+  to_D <- treated$not(FALSE)$and(clinical)
+  to_A <- infected$and(clinical$not(FALSE))
+  to_U <- NULL
+
   mockery::stub(infection_outcome_process, 'incidence_renderer', mockery::mock())
   mockery::stub(infection_outcome_process, 'boost_immunity', boost_immunity_mock)
   mockery::stub(infection_outcome_process, 'calculate_clinical_infections', clinical_infection_mock)
@@ -111,8 +114,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     renderer,
     parameters,
     prob)
-  
-  
+
   mockery::expect_args(
     boost_immunity_mock,
     1,
@@ -143,6 +145,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     renderer
   )
   
+  mockery::mock_args(treated_mock)[[1]][[2]]$to_vector()
   mockery::expect_args(
     treated_mock,
     1,
@@ -159,9 +162,9 @@ test_that('simulate_infection integrates different types of infection and schedu
     parameters,
     variables,
     timestep,
-    infected,
-    treated,
-    clinical
+    to_D,
+    to_A,
+    to_U
   )
 })
 
@@ -654,6 +657,10 @@ test_that('schedule_infections correctly schedules new infections', {
   infection_mock <- mockery::mock()
   asymp_mock <- mockery::mock()
   
+  to_D <- treated$not(FALSE)$and(clinical_infections)
+  to_A <- clinical_infections$not(FALSE)$and(infections)
+  to_U <- NULL
+  
   mockery::stub(schedule_infections, 'update_infection', infection_mock)
   mockery::stub(schedule_infections, 'update_to_asymptomatic_infection', asymp_mock)
   
@@ -661,9 +668,9 @@ test_that('schedule_infections correctly schedules new infections', {
     parameters,
     variables,
     42,
-    infections,
-    treated,
-    clinical_infections
+    to_D,
+    to_A,
+    to_U
   )
 
   actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector()
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 80f6fa31..283e76e5 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -202,3 +202,56 @@ test_that('Test default vivax incidence rendering works', {
     timestep
   )
 })
+
+test_that('vivax schedule_infections correctly schedules new infections', {
+  parameters <- get_parameters(list(human_population = 20), parasite = "vivax")
+  variables <- create_variables(parameters)
+  
+  variables$state <- individual::CategoricalVariable$new(
+    c('U', 'A', 'D', 'S', 'Tr'),
+    rep(c('S','U','A','D','Tr'), times = 4)
+  )
+  
+  infections <- individual::Bitset$new(20)$insert(1:20)
+  lm_detectable <- individual::Bitset$new(20)$insert(6:20)$and(variables$state$get_index_of(c("S", "U")))
+  clinical <- individual::Bitset$new(20)$insert(11:20)$and(variables$state$get_index_of(c("A"))$or(lm_detectable))
+  treated <- individual::Bitset$new(20)$insert(16:20)$and(clinical)
+  # Only S can be a subpatent infection (1)
+  # Only S and U can be a patent infection (6, 7)
+  # S, U and A can be clinical infections (11, 12, 13), but the model re-infects everyone
+  # Treated only looks at new clinical infections (from SAU, not from D)
+  to_U <- infections$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
+  to_A <- lm_detectable$and(clinical$not(F))
+  to_D <- clinical$and(treated$not(F))
+  
+  infection_mock <- mockery::mock()
+  mockery::stub(schedule_infections, 'update_infection', infection_mock)
+  
+  schedule_infections(
+    parameters,
+    variables,
+    timestep = 42,
+    to_D,
+    to_A,
+    to_U
+  )
+  
+  actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector()
+  actual_asymp_infected <- mockery::mock_args(infection_mock)[[2]][[7]]$to_vector()
+  actual_subpatent_infected <- mockery::mock_args(infection_mock)[[3]][[7]]$to_vector()
+  
+  expect_equal(
+    actual_infected,
+    c(11:13)
+  )
+  
+  expect_equal(
+    actual_asymp_infected,
+    c(6, 7)
+  )
+  
+  expect_equal(
+    actual_subpatent_infected,
+    c(1)
+  )
+})

From 78202b9232cb2bc2faec8ef22cb0f1392b140e76 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Tue, 16 Jul 2024 14:24:55 +0100
Subject: [PATCH 14/14] Added a warning to human_infections to notify
 developers that bitsets are rewritten and subsetted in new state assignments.

---
 R/human_infection.R | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/R/human_infection.R b/R/human_infection.R
index 985844ac..d8fb3cc5 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -209,6 +209,7 @@ infection_outcome_process <- function(
         renderer
       )
       
+      ## The treated and infected_humans bitsets are re-written so be cautious!
       to_D <- treated$not(FALSE)$and(clinical)
       to_A <- infected_humans$and(clinical$not(FALSE))
       to_U <- NULL
@@ -253,6 +254,7 @@ infection_outcome_process <- function(
         renderer
       )
       
+      ## The infected_humans,lm_detectable and clinical bitsets are re-written so be cautious!
       to_U <- infected_humans$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
       to_A <- lm_detectable$and(clinical$not(F))
       to_D <- clinical$and(treated$not(F))