From ee08b7aa4658bf1e6e0ba47a459a71cc3e3efa30 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 1 Jul 2024 10:23:53 +0100
Subject: [PATCH 01/11] We must model hynozoite batch formation via bite
 infections, decay of hypnozoite batches through, new infections via relapse
 and the reseting of batch number to 0 for new births.

We create a hypnozoite variable.

The vivax model considers the number of new bites an individual receives (in contrast to whether an individual has been bitten regardless of number of bites) to calculate the infection rates. This is now captured in bitten_humans, which is now a list that also includes n_bites_each.

The vivax model does not factor the eir by the population level biting heterogeneity (although it does use the heterogeneity to distribute biting, but retains the overall total biting according to the EIR).
biting_rate.R

In vivax, we track bites in all individuals, not just SAU as in falciparum. This is because new hypnozoite batches can occur in any human state, regardless of whether this results in a new infection.

The number of bites is considered in calculating probability of a new infection.

To calculate the infection rates for vivax, we take the biting infection rates and add them to the relase rate. We also save the relative rates in the competing hazards object to resolve these competing hazards later.

A new function to resolve between bites and relapse infections is created. For each bite we increase hte number of batches - with capacity for prophylaxis to be added at a future step. We also cap the total number of batches.

We add a hypnozoite batch decay process. Note that this process differs from the immunity decay processes. We work in discrete outputs (you can't have a fractional batch). We also work in the probability that an individual experiences a single loss, rather than loss at the hypnozoite level.

When an individual dies, we reset the hypnozoite batches back to 0.

We now render number of relapses, average hypnozoite batch number and the number of individuals with hypnozoite batches. We may also output these by specified age groups.

Added a test to check relapses occur as expected.
---
 R/biting_process.R                          |  28 +++-
 R/competing_hazards.R                       |   8 +-
 R/human_infection.R                         | 146 ++++++++++++++++----
 R/mortality_processes.R                     |   1 +
 R/processes.R                               |  40 +++++-
 R/render.R                                  |  35 +++++
 R/variables.R                               |   7 +-
 tests/testthat/test-biting-integration.R    |   2 +-
 tests/testthat/test-infection-integration.R |  10 +-
 tests/testthat/test-pev.R                   |   2 +-
 tests/testthat/test-vivax.R                 |  60 ++++++++
 11 files changed, 293 insertions(+), 46 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index 594b9429..a2c16e68 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -78,6 +78,7 @@ simulate_bites <- function(
 ) {
   
   bitten_humans <- individual::Bitset$new(parameters$human_population)
+  n_bites_each <- NULL
   
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
@@ -144,13 +145,23 @@ simulate_bites <- function(
     
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
-    expected_bites <- species_eir * mean(psi)
+    if(parameters$parasite == "falciparum"){
+      # p.f model factors eir by psi
+      expected_bites <- species_eir * mean(psi)
+    } else if (parameters$parasite == "vivax"){
+      # p.v model standardises biting rate het to eir
+      expected_bites <- species_eir
+    }
+    
     if (expected_bites > 0) {
       n_bites <- rpois(1, expected_bites)
       if (n_bites > 0) {
-        bitten_humans$insert(
-          fast_weighted_sample(n_bites, lambda)
-        )
+        bitten <- fast_weighted_sample(n_bites, lambda)
+        bitten_humans$insert(bitten)
+        if(parameters$parasite == "vivax"){
+          # p.v must pass through the number of bites each individual receives
+          n_bites_each <- table(bitten)
+        }
       }
     }
     
@@ -200,9 +211,14 @@ simulate_bites <- function(
       )
     }
   }
-  
+
   renderer$render('n_bitten', bitten_humans$size(), timestep)
-  bitten_humans
+  
+  if(parameters$parasite == "falciparum"){
+    list(bitten_humans = bitten_humans)
+  } else if (parameters$parasite == "vivax"){
+    list(bitten_humans = bitten_humans, n_bites_each = n_bites_each)
+  }
 }
 
 
diff --git a/R/competing_hazards.R b/R/competing_hazards.R
index cdc6aea0..4f429824 100644
--- a/R/competing_hazards.R
+++ b/R/competing_hazards.R
@@ -14,15 +14,21 @@ CompetingOutcome <- R6::R6Class(
       }
       private$targeted_process <- targeted_process
       self$rates <- rep(0, size)
+      self$relative_rates <- rep(0, size)
     },
     set_rates = function(rates){
       self$rates <- rates
     },
+    set_relative_rates = function(relative_rates){
+      self$relative_rates <- relative_rates
+    },
     execute = function(t, target){
       private$targeted_process(t, target)
       self$rates <- rep(0, length(self$rates))
+      self$relative_rates <- rep(0, length(self$relative_rates))
     },
-    rates = NULL
+    rates = NULL,
+    relative_rates = NULL
   )
 )
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 6540e7ff..1d510320 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -21,11 +21,11 @@ simulate_infection <- function(
     renderer,
     infection_outcome
 ) {
-  if (bitten_humans$size() > 0) {
+  if (bitten_humans$bitten_humans$size() > 0) {
     if(parameters$parasite == "falciparum"){
       boost_immunity(
         variables$ib,
-        bitten_humans,
+        bitten_humans$bitten_humans,
         variables$last_boosted_ib,
         timestep,
         parameters$ub
@@ -60,18 +60,22 @@ calculate_infections <- function(
     timestep,
     infection_outcome
 ) {
-  source_humans <- variables$state$get_index_of(
-    c('S', 'A', 'U'))$and(bitten_humans)
-
+  
+  if(bitten_humans$bitten_humans$size() == 0){return(bitten_humans$bitten_humans)}
+  
   if(parameters$parasite == "falciparum"){
+    source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans$bitten_humans)
+    
     ## p.f models blood immunity
     b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
     
   } else if (parameters$parasite == "vivax"){
-    ## p.v does not model blood immunity
-    b <- parameters$b
+    ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
+    source_humans <- bitten_humans$bitten_humans$copy()$or(variables$hypnozoites$get_index_of(0)$not(TRUE))
+    bitten_vector <- bitten_humans$bitten_humans$to_vector()
+    ## p.v does not model blood immunity but must take into account multiple bites per person
+    b <- 1-(1-parameters$b)^bitten_humans$n_bites_each
   }
-
   source_vector <- source_humans$to_vector()
 
   # calculate prophylaxis
@@ -114,8 +118,26 @@ calculate_infections <- function(
       alpha[pev_profile]
     )
   }
-
-  prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+  
+  infection_rates <- rep(0, length = parameters$human_population)
+  if(parameters$parasite == "falciparum"){
+    prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    infection_rates[source_vector] <- prob_to_rate(prob)
+    
+  } else if (parameters$parasite == "vivax"){
+    ## calculated rate of infection for all bitten or with hypnozoites
+    relapse_rates <- variables$hypnozoites$get_values() * parameters$f
+    infection_rates[bitten_vector] <- infection_rates[bitten_vector] + prob_to_rate(b)
+    relative_rates <- relapse_rates/infection_rates
+    relative_rates[is.nan(relative_rates)] <- 0
+    
+    infection_outcome$set_relative_rates(relative_rates)
+    
+    ## get relative rates to get probability bitten over relapse
+    prob <- rate_to_prob(infection_rates)
+    prob[source_vector] <- prob[source_vector] * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    infection_rates <- prob_to_rate(prob)
+  }
 
   ## probability of incidence must be rendered at each timestep
   incidence_probability_renderer(
@@ -130,8 +152,6 @@ calculate_infections <- function(
   )
   
   ## capture infection rates to resolve in competing hazards
-  infection_rates <- rep(0, length = parameters$human_population)
-  infection_rates[source_vector] <- prob_to_rate(prob)
   infection_outcome$set_rates(infection_rates)
 }
 
@@ -145,6 +165,7 @@ calculate_infections <- function(
 #' @param renderer model render object
 #' @param parameters model parameters
 #' @param prob vector of population probabilities of infection
+#' @param relative_rates relative rates of hypnozoite relapse relative to total infection rate
 #' @noRd
 infection_outcome_process <- function(
     timestep,
@@ -152,8 +173,10 @@ infection_outcome_process <- function(
     variables,
     renderer,
     parameters,
-    prob){
+    prob,
+    relative_rates = NULL){
   
+  renderer$render('n_infections', infected_humans$size(), timestep)
   incidence_renderer(
     variables$birth,
     renderer,
@@ -165,6 +188,8 @@ infection_outcome_process <- function(
   )
   
   if (infected_humans$size() > 0) {
+    # p.f SAU infections get boosted (already subset)
+    # p.v SAUdTr infections all get boosted (therefore not subset)
     boost_immunity(
       variables$ica,
       infected_humans,
@@ -200,6 +225,15 @@ infection_outcome_process <- function(
       patent_infections <- NULL
       
     } else if (parameters$parasite == "vivax"){
+      relapse_bite_infection_hazard_resolution(
+        infected_humans,
+        relative_rates,
+        variables,
+        parameters,
+        renderer,
+        timestep
+      )
+      
       boost_immunity(
         variables$iaa,
         infected_humans,
@@ -226,27 +260,85 @@ infection_outcome_process <- function(
         timestep
       )
     }
+    
+    treated <- calculate_treated(
+      variables,
+      clinical_infections,
+      parameters,
+      timestep,
+      renderer
+    )
+    
+    schedule_infections(
+      variables,
+      patent_infections,
+      clinical_infections,
+      treated,
+      infected_humans,
+      parameters,
+      timestep
+    )
   }
-  
-  treated <- calculate_treated(
+}
+
+#' @title Relapse/bite infection competing hazard resolution (p.v only)
+#' @description
+#' Resolves competing hazards of bite and hypnozoite relapse infections. 
+#' For bite infections we increase the batch number and factor in drug prophylaxis.
+#' 
+#' @param variables a list of all of the model variables
+#' @param infected_humans bitset of infected humans
+#' @param relative_rates relative rate of relapse infection
+#' @param variables model variables
+#' @param parameters model parameters
+#' @param renderer model renderer
+#' @param timestep current timestep
+#' @noRd
+relapse_bite_infection_hazard_resolution <- function(
+    infected_humans,
+    relative_rates,
     variables,
-    clinical_infections,
     parameters,
-    timestep,
-    renderer
-  )
-  
-  renderer$render('n_infections', infected_humans$size(), timestep)
+    renderer,
+    timestep
+){
   
-  schedule_infections(
-    variables,
-    patent_infections,
-    clinical_infections,
-    treated,
+  # draw relapses from total infections
+  relapse_infections <- bitset_at(
     infected_humans,
-    parameters,
+    bernoulli_multi_p(relative_rates[infected_humans$to_vector()])
+  )
+  
+  renderer$render('n_relapses', relapse_infections$size(), timestep)
+  # render relapse infections by age
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    relapse_infections,
+    'inc_relapse_',
+    parameters$incidence_relapse_rendering_min_ages,
+    parameters$incidence_relapse_rendering_max_ages,
     timestep
   )
+  
+  # get bite infections
+  bite_infections <- infected_humans$copy()$and(relapse_infections$not(inplace = F))
+  
+  ## all bitten humans with an infectious bite (incorporating prophylaxis) get a new batch of hypnozoites
+  if(bite_infections$size()>0){
+    new_hypnozoite_batch_formed <- bite_infections # bitset_at(bite_infections, bernoulli_multi_p(1-ls_prophylaxis))
+    
+    # make sure batches are capped
+    current_batches <- variables$hypnozoites$get_values(new_hypnozoite_batch_formed)
+    new_batch_number <- ifelse(current_batches == parameters$kmax,
+                               current_batches,
+                               current_batches + 1)
+    
+    variables$hypnozoites$queue_update(
+      new_batch_number,
+      new_hypnozoite_batch_formed$and(bite_infections)
+    )
+  }
 }
 
 #' @title Calculate patent infections (p.v only)
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index e4fa323f..a66c3eb2 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -112,6 +112,7 @@ reset_target <- function(variables, events, target, state, parameters, timestep)
     } else if (parameters$parasite == "vivax"){
       variables$last_boosted_iaa$queue_update(-1, target)
       variables$iaa$queue_update(0, target)
+      variables$hypnozoites$queue_update(0, target)
     }
 
     # treatment
diff --git a/R/processes.R b/R/processes.R
index b1b3a9f4..b1959c8b 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -59,7 +59,8 @@ create_processes <- function(
       processes,
       # Anti-parasite immunity
       create_exponential_decay_process(variables$iam, parameters$rm),
-      create_exponential_decay_process(variables$iaa, parameters$ra)
+      create_exponential_decay_process(variables$iaa, parameters$ra),
+      create_hypnozoite_batch_decay_process(variables$hypnozoites, parameters$gammal)
     )
   }
 
@@ -84,7 +85,8 @@ create_processes <- function(
     targeted_process = function(timestep, target){
       infection_outcome_process(timestep, target, 
                                 variables, renderer, parameters, 
-                                prob = rate_to_prob(infection_outcome$rates))
+                                prob = rate_to_prob(infection_outcome$rates),
+                                relative_rates = infection_outcome$relative_rates)
     },
     size = parameters$human_population
   )
@@ -189,7 +191,23 @@ create_processes <- function(
   if(parameters$parasite == "falciparum"){
     imm_var_names <- c(imm_var_names,'ib','iva','ivm','id')
   } else if (parameters$parasite == "vivax"){
-    imm_var_names <- c(imm_var_names,'iaa','iam')
+    imm_var_names <- c(imm_var_names,'iaa','iam','hypnozoites')
+    
+    ## hypnozoite infection prevalence rendering
+    processes <- c(
+      processes,
+      create_n_with_hypnozoites_renderer_process(
+        renderer,
+        variables$hypnozoites,
+        parameters
+      ),
+      create_n_with_hypnozoites_age_renderer_process(
+        variables$hypnozoites,
+        variables$birth,
+        parameters,
+        renderer
+      )
+    )
   }
 
   processes <- c(
@@ -344,3 +362,19 @@ create_lagged_eir <- function(variables, solvers, parameters) {
     }
   )
 }
+#' @title Hypnozoite decay function
+#' @description
+#' calulates the number of individuals in whom a batch decay occurs
+#'
+#' @param variable the hypnozoite variable to update
+#' @param rate the hypnozoite decay rate
+#' @noRd
+create_hypnozoite_batch_decay_process <- function(hypnozoites, gammal){
+  function(timestep){
+    to_decay <- bernoulli_multi_p(p = rate_to_prob(hypnozoites$get_values() * gammal))
+    hypnozoites$queue_update(
+      hypnozoites$get_values(to_decay) - 1,
+      to_decay
+    )
+  }
+}
\ No newline at end of file
diff --git a/R/render.R b/R/render.R
index 68d551ad..5943345e 100644
--- a/R/render.R
+++ b/R/render.R
@@ -278,4 +278,39 @@ populate_metapopulation_incidence_rendering_columns <- function(renderer, parame
   for(i in length(parameters)){
     populate_incidence_rendering_columns(renderer[[i]], parameters[[i]])
   }
+}
+
+create_n_with_hypnozoites_renderer_process <- function(
+    renderer,
+    hypnozoites,
+    parameters
+) {
+  function(timestep) {
+    renderer$render(
+      paste0("n_with_hypnozoites"),
+      hypnozoites$size() - hypnozoites$get_size_of(0),
+      timestep
+    )
+  }
+}
+
+create_n_with_hypnozoites_age_renderer_process <- function(
+    hypnozoites,
+    birth,
+    parameters,
+    renderer
+) {
+  function(timestep) {
+    for (i in seq_along(parameters$n_with_hypnozoites_rendering_min_ages)) {
+      lower <- parameters$n_with_hypnozoites_rendering_min_ages[[i]]
+      upper <- parameters$n_with_hypnozoites_rendering_max_ages[[i]]
+      in_age <- in_age_range(birth, timestep, lower, upper)
+      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+      renderer$render(
+        paste0('n_with_hypnozoites_', lower, '_', upper),
+        sum(hypnozoites$get_values(index = in_age)!=0),
+        timestep
+      )
+    }
+  }
 }
\ No newline at end of file
diff --git a/R/variables.R b/R/variables.R
index 7e954a40..757807f3 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -18,6 +18,7 @@
 #' * ICA  - Acquired immunity to clinical disease
 #' * IVA  - Acquired immunity to severe disease (p.f only)
 #' * ID - Acquired immunity to detectability (p.f only)
+#' * hypnozoites - Hypnozoite batch number (p.v only)
 #' * zeta - Heterogeneity of human individuals
 #' * zeta_group - Discretised heterogeneity of human individuals
 #' * last_pev_timestep - The timestep of the last pev vaccination (-1 if there
@@ -98,6 +99,9 @@ create_variables <- function(parameters) {
   initial_states <- initial_state(parameters, initial_age, groups, eq, states)
   state <- individual::CategoricalVariable$new(states, initial_states)
   birth <- individual::IntegerVariable$new(-initial_age)
+  if(parameters$parasite == "vivax"){
+    hypnozoites <- individual::IntegerVariable$new(rep(parameters$init_hyp, parameters$human_population))
+  }
   
   # Maternal immunity to clinical disease
   icm <- individual::DoubleVariable$new(
@@ -299,7 +303,8 @@ create_variables <- function(parameters) {
     variables <- c(variables,
                    last_boosted_iaa = last_boosted_iaa,
                    iaa = iaa,
-                   iam = iam
+                   iam = iam,
+                   hypnozoites = hypnozoites
     )
   }
   
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 4be80329..713fa83a 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -127,7 +127,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     lagged_eir
   )
 
-  expect_equal(bitten$to_vector(), c(2, 3))
+  expect_equal(bitten$bitten_humans$to_vector(), c(2, 3))
 
   f <- parameters$blood_meal_rates[[1]]
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 17d0e958..54289385 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -34,7 +34,7 @@ test_that('simulate_infection integrates different types of infection and schedu
   simulate_infection(
     variables,
     events,
-    bitten,
+    list(bitten_humans = bitten),
     age,
     parameters,
     timestep,
@@ -56,7 +56,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     infection_mock,
     1,
     variables,
-    bitten,
+    list(bitten_humans = bitten),
     parameters,
     renderer,
     timestep,
@@ -64,7 +64,6 @@ test_that('simulate_infection integrates different types of infection and schedu
   )
 })
 
-
 test_that('simulate_infection integrates different types of infection and scheduling', {
   population <- 8
   timestep <- 5
@@ -217,7 +216,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     depth = 4
   )
   
-  bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  bitten_humans <- list(bitten_humans = individual::Bitset$new(4)$insert(c(1, 2, 3, 4)))
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
@@ -700,11 +699,10 @@ test_that('prophylaxis is considered for medicated humans', {
     ib = individual::DoubleVariable$new(c(.2, .3, .5, .9))
   )
   
-  bitten <- individual::Bitset$new(4)$insert(seq(4))
+  bitten <- list(bitten_humans = individual::Bitset$new(4)$insert(seq(4)))
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
   
-  
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
       infection_outcome_process(timestep, target, variables, renderer, parameters)
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index 9a5f1175..2f173edf 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -173,7 +173,7 @@ test_that('Infection considers pev efficacy', {
 
   infection_rates <- calculate_infections(
     variables = variables,
-    bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
+    bitten_humans = list(bitten_humans = individual::Bitset$new(4)$insert(seq(4))),
     parameters = parameters,
     renderer = mock_render(timestep),
     timestep = timestep,
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 289848ce..b45af003 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -216,3 +216,63 @@ test_that('Test default vivax incidence rendering works', {
     timestep
   )
 })
+
+test_that('relapses are recognised with division between bite infections and relapses', {
+  timestep <- 50
+  parameters <- get_parameters(parasite = "vivax", overrides = list(human_population = 4))
+  
+  variables <- list(
+    state = individual::CategoricalVariable$new(
+      c('D', 'S', 'A', 'U', 'Tr'),
+      c('D', 'S', 'A', 'U')
+    ),
+    infectivity = individual::DoubleVariable$new(rep(0, 4)),
+    recovery_rates = individual::DoubleVariable$new(rep(0, 4)),
+    drug = individual::DoubleVariable$new(rep(0, 4)),
+    drug_time = individual::DoubleVariable$new(rep(-1, 4)),
+    iaa = individual::DoubleVariable$new(rep(0, 4)),
+    iam = individual::DoubleVariable$new(rep(0, 4)),
+    ica = individual::DoubleVariable$new(rep(0, 4)),
+    icm = individual::DoubleVariable$new(rep(0, 4)),
+    last_boosted_iaa = individual::DoubleVariable$new(rep(-1, 4)),
+    last_boosted_ica = individual::DoubleVariable$new(rep(-1, 4)),
+    last_eff_pev_timestep = individual::DoubleVariable$new(rep(-1, 4)),
+    pev_profile = individual::IntegerVariable$new(rep(-1, 4)),
+    hypnozoites = individual::IntegerVariable$new(c(0, 1, 2, 3))
+  )
+  
+  bernoulli_mock <- mockery::mock(c(2, 4), 2, cycle = TRUE)
+  calc_mock <- mockery::mock(individual::Bitset$new(4)$insert(2))
+  mockery::stub(infection_outcome_process, 'bernoulli_multi_p', bernoulli_mock, depth = 2)
+  mockery::stub(infection_outcome_process, 'calculate_clinical_infections', calc_mock)
+  
+  renderer <- mock_render(1)
+  infected_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  
+  infection_outcome_process(
+    timestep = timestep,
+    infected_humans,
+    variables,
+    renderer,
+    parameters,
+    prob = c(rep(0.5, 4)),
+    relative_rate = c(rep(0.5, 4))
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_infections',
+    4,
+    timestep
+  )
+
+  mockery::expect_args(
+    renderer$render_mock(),
+    2,
+    'n_relapses',
+    2,
+    timestep
+  )
+  
+})

From 94a1ce2210afbc14321fa78d6ba57a5735487017 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 1 Jul 2024 10:23:53 +0100
Subject: [PATCH 02/11] We must model hynozoite batch formation via bite
 infections, decay of hypnozoite batches through, new infections via relapse
 and the reseting of batch number to 0 for new births.

We create a hypnozoite variable.

The vivax model considers the number of new bites an individual receives (in contrast to whether an individual has been bitten regardless of number of bites) to calculate the infection rates. This is now captured in bitten_humans, which is now a list that also includes n_bites_each.

The vivax model does not factor the eir by the population level biting heterogeneity (although it does use the heterogeneity to distribute biting, but retains the overall total biting according to the EIR).
biting_rate.R

In vivax, we track bites in all individuals, not just SAU as in falciparum. This is because new hypnozoite batches can occur in any human state, regardless of whether this results in a new infection.

The number of bites is considered in calculating probability of a new infection.

To calculate the infection rates for vivax, we take the biting infection rates and add them to the relase rate. We also save the relative rates in the competing hazards object to resolve these competing hazards later.

A new function to resolve between bites and relapse infections is created. For each bite we increase hte number of batches - with capacity for prophylaxis to be added at a future step. We also cap the total number of batches.

We add a hypnozoite batch decay process. Note that this process differs from the immunity decay processes. We work in discrete outputs (you can't have a fractional batch). We also work in the probability that an individual experiences a single loss, rather than loss at the hypnozoite level.

When an individual dies, we reset the hypnozoite batches back to 0.

We now render number of relapses, average hypnozoite batch number and the number of individuals with hypnozoite batches. We may also output these by specified age groups.

Added a test to check relapses occur as expected.
---
 R/biting_process.R                          |  28 +++-
 R/competing_hazards.R                       |   8 +-
 R/human_infection.R                         | 139 ++++++++++++++++++--
 R/mortality_processes.R                     |   1 +
 R/processes.R                               |  40 +++++-
 R/render.R                                  |  35 +++++
 R/variables.R                               |   7 +-
 tests/testthat/test-biting-integration.R    |   2 +-
 tests/testthat/test-infection-integration.R |  10 +-
 tests/testthat/test-pev.R                   |   2 +-
 tests/testthat/test-vivax.R                 |  60 +++++++++
 11 files changed, 300 insertions(+), 32 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index 594b9429..a2c16e68 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -78,6 +78,7 @@ simulate_bites <- function(
 ) {
   
   bitten_humans <- individual::Bitset$new(parameters$human_population)
+  n_bites_each <- NULL
   
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
@@ -144,13 +145,23 @@ simulate_bites <- function(
     
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
-    expected_bites <- species_eir * mean(psi)
+    if(parameters$parasite == "falciparum"){
+      # p.f model factors eir by psi
+      expected_bites <- species_eir * mean(psi)
+    } else if (parameters$parasite == "vivax"){
+      # p.v model standardises biting rate het to eir
+      expected_bites <- species_eir
+    }
+    
     if (expected_bites > 0) {
       n_bites <- rpois(1, expected_bites)
       if (n_bites > 0) {
-        bitten_humans$insert(
-          fast_weighted_sample(n_bites, lambda)
-        )
+        bitten <- fast_weighted_sample(n_bites, lambda)
+        bitten_humans$insert(bitten)
+        if(parameters$parasite == "vivax"){
+          # p.v must pass through the number of bites each individual receives
+          n_bites_each <- table(bitten)
+        }
       }
     }
     
@@ -200,9 +211,14 @@ simulate_bites <- function(
       )
     }
   }
-  
+
   renderer$render('n_bitten', bitten_humans$size(), timestep)
-  bitten_humans
+  
+  if(parameters$parasite == "falciparum"){
+    list(bitten_humans = bitten_humans)
+  } else if (parameters$parasite == "vivax"){
+    list(bitten_humans = bitten_humans, n_bites_each = n_bites_each)
+  }
 }
 
 
diff --git a/R/competing_hazards.R b/R/competing_hazards.R
index cdc6aea0..4f429824 100644
--- a/R/competing_hazards.R
+++ b/R/competing_hazards.R
@@ -14,15 +14,21 @@ CompetingOutcome <- R6::R6Class(
       }
       private$targeted_process <- targeted_process
       self$rates <- rep(0, size)
+      self$relative_rates <- rep(0, size)
     },
     set_rates = function(rates){
       self$rates <- rates
     },
+    set_relative_rates = function(relative_rates){
+      self$relative_rates <- relative_rates
+    },
     execute = function(t, target){
       private$targeted_process(t, target)
       self$rates <- rep(0, length(self$rates))
+      self$relative_rates <- rep(0, length(self$relative_rates))
     },
-    rates = NULL
+    rates = NULL,
+    relative_rates = NULL
   )
 )
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 56fe0340..8de44753 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -21,11 +21,11 @@ simulate_infection <- function(
     renderer,
     infection_outcome
 ) {
-  if (bitten_humans$size() > 0) {
+  if (bitten_humans$bitten_humans$size() > 0) {
     if(parameters$parasite == "falciparum"){
       boost_immunity(
         variables$ib,
-        bitten_humans,
+        bitten_humans$bitten_humans,
         variables$last_boosted_ib,
         timestep,
         parameters$ub
@@ -60,18 +60,22 @@ calculate_infections <- function(
     timestep,
     infection_outcome
 ) {
-  source_humans <- variables$state$get_index_of(
-    c('S', 'A', 'U'))$and(bitten_humans)
-
+  
+  if(bitten_humans$bitten_humans$size() == 0){return(bitten_humans$bitten_humans)}
+  
   if(parameters$parasite == "falciparum"){
+    source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans$bitten_humans)
+    
     ## p.f models blood immunity
     b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
     
   } else if (parameters$parasite == "vivax"){
-    ## p.v does not model blood immunity
-    b <- parameters$b
+    ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
+    source_humans <- bitten_humans$bitten_humans$copy()$or(variables$hypnozoites$get_index_of(0)$not(TRUE))
+    bitten_vector <- bitten_humans$bitten_humans$to_vector()
+    ## p.v does not model blood immunity but must take into account multiple bites per person
+    b <- 1-(1-parameters$b)^bitten_humans$n_bites_each
   }
-
   source_vector <- source_humans$to_vector()
 
   # calculate prophylaxis
@@ -114,8 +118,26 @@ calculate_infections <- function(
       alpha[pev_profile]
     )
   }
-
-  prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+  
+  infection_rates <- rep(0, length = parameters$human_population)
+  if(parameters$parasite == "falciparum"){
+    prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    infection_rates[source_vector] <- prob_to_rate(prob)
+    
+  } else if (parameters$parasite == "vivax"){
+    ## calculated rate of infection for all bitten or with hypnozoites
+    relapse_rates <- variables$hypnozoites$get_values() * parameters$f
+    infection_rates[bitten_vector] <- infection_rates[bitten_vector] + prob_to_rate(b)
+    relative_rates <- relapse_rates/infection_rates
+    relative_rates[is.nan(relative_rates)] <- 0
+    
+    infection_outcome$set_relative_rates(relative_rates)
+    
+    ## get relative rates to get probability bitten over relapse
+    prob <- rate_to_prob(infection_rates)
+    prob[source_vector] <- prob[source_vector] * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    infection_rates <- prob_to_rate(prob)
+  }
 
   ## probability of incidence must be rendered at each timestep
   incidence_probability_renderer(
@@ -130,8 +152,6 @@ calculate_infections <- function(
   )
   
   ## capture infection rates to resolve in competing hazards
-  infection_rates <- rep(0, length = parameters$human_population)
-  infection_rates[source_vector] <- prob_to_rate(prob)
   infection_outcome$set_rates(infection_rates)
 }
 
@@ -145,6 +165,7 @@ calculate_infections <- function(
 #' @param renderer model render object
 #' @param parameters model parameters
 #' @param prob vector of population probabilities of infection
+#' @param relative_rates relative rates of hypnozoite relapse relative to total infection rate
 #' @noRd
 infection_outcome_process <- function(
     timestep,
@@ -152,8 +173,10 @@ infection_outcome_process <- function(
     variables,
     renderer,
     parameters,
-    prob){
+    prob,
+    relative_rates = NULL){
   
+  renderer$render('n_infections', infected_humans$size(), timestep)
   incidence_renderer(
     variables$birth,
     renderer,
@@ -165,6 +188,8 @@ infection_outcome_process <- function(
   )
   
   if (infected_humans$size() > 0) {
+    # p.f SAU infections get boosted (already subset)
+    # p.v SAUdTr infections all get boosted (therefore not subset)
     boost_immunity(
       variables$ica,
       infected_humans,
@@ -217,6 +242,15 @@ infection_outcome_process <- function(
       )
       
     } else if (parameters$parasite == "vivax"){
+      relapse_bite_infection_hazard_resolution(
+        infected_humans,
+        relative_rates,
+        variables,
+        parameters,
+        renderer,
+        timestep
+      )
+      
       boost_immunity(
         variables$iaa,
         infected_humans,
@@ -262,10 +296,89 @@ infection_outcome_process <- function(
         clinical_infections,
         lm_det_infections
       )
+      
+      treated <- calculate_treated(
+        variables,
+        clinical_infections,
+        parameters,
+        timestep,
+        renderer
+      )
+      
+      schedule_infections(
+        parameters,
+        variables,
+        timestep,
+        infected_humans,
+        treated,
+        clinical_infections,
+        lm_det_infections
+      )
+      
     }
   }
 }
 
+#' @title Relapse/bite infection competing hazard resolution (p.v only)
+#' @description
+#' Resolves competing hazards of bite and hypnozoite relapse infections. 
+#' For bite infections we increase the batch number and factor in drug prophylaxis.
+#' 
+#' @param variables a list of all of the model variables
+#' @param infected_humans bitset of infected humans
+#' @param relative_rates relative rate of relapse infection
+#' @param variables model variables
+#' @param parameters model parameters
+#' @param renderer model renderer
+#' @param timestep current timestep
+#' @noRd
+relapse_bite_infection_hazard_resolution <- function(
+    infected_humans,
+    relative_rates,
+    variables,
+    parameters,
+    renderer,
+    timestep
+){
+  
+  # draw relapses from total infections
+  relapse_infections <- bitset_at(
+    infected_humans,
+    bernoulli_multi_p(relative_rates[infected_humans$to_vector()])
+  )
+  
+  renderer$render('n_relapses', relapse_infections$size(), timestep)
+  # render relapse infections by age
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    relapse_infections,
+    'inc_relapse_',
+    parameters$incidence_relapse_rendering_min_ages,
+    parameters$incidence_relapse_rendering_max_ages,
+    timestep
+  )
+  
+  # get bite infections
+  bite_infections <- infected_humans$copy()$and(relapse_infections$not(inplace = F))
+  
+  ## all bitten humans with an infectious bite (incorporating prophylaxis) get a new batch of hypnozoites
+  if(bite_infections$size()>0){
+    new_hypnozoite_batch_formed <- bite_infections # bitset_at(bite_infections, bernoulli_multi_p(1-ls_prophylaxis))
+    
+    # make sure batches are capped
+    current_batches <- variables$hypnozoites$get_values(new_hypnozoite_batch_formed)
+    new_batch_number <- ifelse(current_batches == parameters$kmax,
+                               current_batches,
+                               current_batches + 1)
+    
+    variables$hypnozoites$queue_update(
+      new_batch_number,
+      new_hypnozoite_batch_formed$and(bite_infections)
+    )
+  }
+}
+
 #' @title Calculate light microscopy detectable infections (p.v only)
 #' @description
 #' Sample light microscopy detectable infections from all infections
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index e4fa323f..a66c3eb2 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -112,6 +112,7 @@ reset_target <- function(variables, events, target, state, parameters, timestep)
     } else if (parameters$parasite == "vivax"){
       variables$last_boosted_iaa$queue_update(-1, target)
       variables$iaa$queue_update(0, target)
+      variables$hypnozoites$queue_update(0, target)
     }
 
     # treatment
diff --git a/R/processes.R b/R/processes.R
index afa74f50..2b0adef7 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -59,7 +59,8 @@ create_processes <- function(
       processes,
       # Anti-parasite immunity
       create_exponential_decay_process(variables$iam, parameters$rm),
-      create_exponential_decay_process(variables$iaa, parameters$ra)
+      create_exponential_decay_process(variables$iaa, parameters$ra),
+      create_hypnozoite_batch_decay_process(variables$hypnozoites, parameters$gammal)
     )
   }
 
@@ -84,7 +85,8 @@ create_processes <- function(
     targeted_process = function(timestep, target){
       infection_outcome_process(timestep, target, 
                                 variables, renderer, parameters, 
-                                prob = rate_to_prob(infection_outcome$rates))
+                                prob = rate_to_prob(infection_outcome$rates),
+                                relative_rates = infection_outcome$relative_rates)
     },
     size = parameters$human_population
   )
@@ -189,7 +191,23 @@ create_processes <- function(
   if(parameters$parasite == "falciparum"){
     imm_var_names <- c(imm_var_names, 'ib', 'iva', 'ivm', 'id')
   } else if (parameters$parasite == "vivax"){
-    imm_var_names <- c(imm_var_names, 'iaa', 'iam')
+    imm_var_names <- c(imm_var_names, 'iaa', 'iam','hypnozoites')
+    
+    ## hypnozoite infection prevalence rendering
+    processes <- c(
+      processes,
+      create_n_with_hypnozoites_renderer_process(
+        renderer,
+        variables$hypnozoites,
+        parameters
+      ),
+      create_n_with_hypnozoites_age_renderer_process(
+        variables$hypnozoites,
+        variables$birth,
+        parameters,
+        renderer
+      )
+    )
   }
 
   processes <- c(
@@ -344,3 +362,19 @@ create_lagged_eir <- function(variables, solvers, parameters) {
     }
   )
 }
+#' @title Hypnozoite decay function
+#' @description
+#' calulates the number of individuals in whom a batch decay occurs
+#'
+#' @param variable the hypnozoite variable to update
+#' @param rate the hypnozoite decay rate
+#' @noRd
+create_hypnozoite_batch_decay_process <- function(hypnozoites, gammal){
+  function(timestep){
+    to_decay <- bernoulli_multi_p(p = rate_to_prob(hypnozoites$get_values() * gammal))
+    hypnozoites$queue_update(
+      hypnozoites$get_values(to_decay) - 1,
+      to_decay
+    )
+  }
+}
\ No newline at end of file
diff --git a/R/render.R b/R/render.R
index b219cdb9..0d066b46 100644
--- a/R/render.R
+++ b/R/render.R
@@ -278,4 +278,39 @@ populate_metapopulation_incidence_rendering_columns <- function(renderer, parame
   for(i in length(parameters)){
     populate_incidence_rendering_columns(renderer[[i]], parameters[[i]])
   }
+}
+
+create_n_with_hypnozoites_renderer_process <- function(
+    renderer,
+    hypnozoites,
+    parameters
+) {
+  function(timestep) {
+    renderer$render(
+      paste0("n_with_hypnozoites"),
+      hypnozoites$size() - hypnozoites$get_size_of(0),
+      timestep
+    )
+  }
+}
+
+create_n_with_hypnozoites_age_renderer_process <- function(
+    hypnozoites,
+    birth,
+    parameters,
+    renderer
+) {
+  function(timestep) {
+    for (i in seq_along(parameters$n_with_hypnozoites_rendering_min_ages)) {
+      lower <- parameters$n_with_hypnozoites_rendering_min_ages[[i]]
+      upper <- parameters$n_with_hypnozoites_rendering_max_ages[[i]]
+      in_age <- in_age_range(birth, timestep, lower, upper)
+      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+      renderer$render(
+        paste0('n_with_hypnozoites_', lower, '_', upper),
+        sum(hypnozoites$get_values(index = in_age)!=0),
+        timestep
+      )
+    }
+  }
 }
\ No newline at end of file
diff --git a/R/variables.R b/R/variables.R
index 7e954a40..757807f3 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -18,6 +18,7 @@
 #' * ICA  - Acquired immunity to clinical disease
 #' * IVA  - Acquired immunity to severe disease (p.f only)
 #' * ID - Acquired immunity to detectability (p.f only)
+#' * hypnozoites - Hypnozoite batch number (p.v only)
 #' * zeta - Heterogeneity of human individuals
 #' * zeta_group - Discretised heterogeneity of human individuals
 #' * last_pev_timestep - The timestep of the last pev vaccination (-1 if there
@@ -98,6 +99,9 @@ create_variables <- function(parameters) {
   initial_states <- initial_state(parameters, initial_age, groups, eq, states)
   state <- individual::CategoricalVariable$new(states, initial_states)
   birth <- individual::IntegerVariable$new(-initial_age)
+  if(parameters$parasite == "vivax"){
+    hypnozoites <- individual::IntegerVariable$new(rep(parameters$init_hyp, parameters$human_population))
+  }
   
   # Maternal immunity to clinical disease
   icm <- individual::DoubleVariable$new(
@@ -299,7 +303,8 @@ create_variables <- function(parameters) {
     variables <- c(variables,
                    last_boosted_iaa = last_boosted_iaa,
                    iaa = iaa,
-                   iam = iam
+                   iam = iam,
+                   hypnozoites = hypnozoites
     )
   }
   
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 4be80329..713fa83a 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -127,7 +127,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     lagged_eir
   )
 
-  expect_equal(bitten$to_vector(), c(2, 3))
+  expect_equal(bitten$bitten_humans$to_vector(), c(2, 3))
 
   f <- parameters$blood_meal_rates[[1]]
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 2d99c488..ef4cf03c 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -34,7 +34,7 @@ test_that('simulate_infection integrates different types of infection and schedu
   simulate_infection(
     variables,
     events,
-    bitten,
+    list(bitten_humans = bitten),
     age,
     parameters,
     timestep,
@@ -56,7 +56,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     infection_mock,
     1,
     variables,
-    bitten,
+    list(bitten_humans = bitten),
     parameters,
     renderer,
     timestep,
@@ -64,7 +64,6 @@ test_that('simulate_infection integrates different types of infection and schedu
   )
 })
 
-
 test_that('simulate_infection integrates different types of infection and scheduling', {
   population <- 8
   timestep <- 5
@@ -216,7 +215,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     depth = 4
   )
   
-  bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  bitten_humans <- list(bitten_humans = individual::Bitset$new(4)$insert(c(1, 2, 3, 4)))
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
@@ -698,11 +697,10 @@ test_that('prophylaxis is considered for medicated humans', {
     ib = individual::DoubleVariable$new(c(.2, .3, .5, .9))
   )
   
-  bitten <- individual::Bitset$new(4)$insert(seq(4))
+  bitten <- list(bitten_humans = individual::Bitset$new(4)$insert(seq(4)))
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
   
-  
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
       infection_outcome_process(timestep, target, variables, renderer, parameters)
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index 9a5f1175..2f173edf 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -173,7 +173,7 @@ test_that('Infection considers pev efficacy', {
 
   infection_rates <- calculate_infections(
     variables = variables,
-    bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
+    bitten_humans = list(bitten_humans = individual::Bitset$new(4)$insert(seq(4))),
     parameters = parameters,
     renderer = mock_render(timestep),
     timestep = timestep,
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 8dc62d26..39513006 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -239,3 +239,63 @@ test_that('Test default vivax incidence rendering works', {
     timestep
   )
 })
+
+test_that('relapses are recognised with division between bite infections and relapses', {
+  timestep <- 50
+  parameters <- get_parameters(parasite = "vivax", overrides = list(human_population = 4))
+  
+  variables <- list(
+    state = individual::CategoricalVariable$new(
+      c('D', 'S', 'A', 'U', 'Tr'),
+      c('D', 'S', 'A', 'U')
+    ),
+    infectivity = individual::DoubleVariable$new(rep(0, 4)),
+    recovery_rates = individual::DoubleVariable$new(rep(0, 4)),
+    drug = individual::DoubleVariable$new(rep(0, 4)),
+    drug_time = individual::DoubleVariable$new(rep(-1, 4)),
+    iaa = individual::DoubleVariable$new(rep(0, 4)),
+    iam = individual::DoubleVariable$new(rep(0, 4)),
+    ica = individual::DoubleVariable$new(rep(0, 4)),
+    icm = individual::DoubleVariable$new(rep(0, 4)),
+    last_boosted_iaa = individual::DoubleVariable$new(rep(-1, 4)),
+    last_boosted_ica = individual::DoubleVariable$new(rep(-1, 4)),
+    last_eff_pev_timestep = individual::DoubleVariable$new(rep(-1, 4)),
+    pev_profile = individual::IntegerVariable$new(rep(-1, 4)),
+    hypnozoites = individual::IntegerVariable$new(c(0, 1, 2, 3))
+  )
+  
+  bernoulli_mock <- mockery::mock(c(2, 4), 2, cycle = TRUE)
+  calc_mock <- mockery::mock(individual::Bitset$new(4)$insert(2))
+  mockery::stub(infection_outcome_process, 'bernoulli_multi_p', bernoulli_mock, depth = 2)
+  mockery::stub(infection_outcome_process, 'calculate_clinical_infections', calc_mock)
+  
+  renderer <- mock_render(1)
+  infected_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  
+  infection_outcome_process(
+    timestep = timestep,
+    infected_humans,
+    variables,
+    renderer,
+    parameters,
+    prob = c(rep(0.5, 4)),
+    relative_rate = c(rep(0.5, 4))
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_infections',
+    4,
+    timestep
+  )
+
+  mockery::expect_args(
+    renderer$render_mock(),
+    2,
+    'n_relapses',
+    2,
+    timestep
+  )
+  
+})

From 8600743dd2e1d01a8cd3fa30e925e3e3084b03dc Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Thu, 11 Jul 2024 10:21:53 +0100
Subject: [PATCH 03/11] Removed duplicated code from human_infections.R

---
 R/human_infection.R | 37 -------------------------------------
 1 file changed, 37 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index ccbff2d0..d6266d6a 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -296,44 +296,7 @@ infection_outcome_process <- function(
         clinical_infections,
         lm_det_infections
       )
-      
-      treated <- calculate_treated(
-        variables,
-        clinical_infections,
-        parameters,
-        timestep,
-        renderer
-      )
-      
-      schedule_infections(
-        parameters,
-        variables,
-        timestep,
-        infected_humans,
-        treated,
-        clinical_infections,
-        lm_det_infections
-      )
-      
     }
-    
-    treated <- calculate_treated(
-      variables,
-      clinical_infections,
-      parameters,
-      timestep,
-      renderer
-    )
-    
-    schedule_infections(
-      variables,
-      patent_infections,
-      clinical_infections,
-      treated,
-      infected_humans,
-      parameters,
-      timestep
-    )
   }
 }
 

From ef475cdb56138ec1358864e24f1b34dd9aaaaed8 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Wed, 17 Jul 2024 12:08:42 +0100
Subject: [PATCH 04/11] Removed n_bites_each, opting instead to create a
 variable  that contains the vector of number of bites per individual that can
 be passed in place of the individuals bitten bitset. I think this simplifies
 a lot of the code, but there may still be ways of improving this.

The relapse bite hazard resolution function was replicated - this has been removed. lm-det incidence rendering has been removed. An if statement has been added to the hypnozoite decay process, and tests adjusted to reflect bitten_humans changes.
---
 R/biting_process.R                          |  24 ++--
 R/human_infection.R                         | 123 ++++----------------
 R/parameters.R                              |   4 -
 R/processes.R                               |  10 +-
 tests/testthat/test-biting-integration.R    |   2 +-
 tests/testthat/test-infection-integration.R |   8 +-
 tests/testthat/test-pev.R                   |   2 +-
 7 files changed, 43 insertions(+), 130 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index a2c16e68..0149f3f6 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -77,9 +77,6 @@ simulate_bites <- function(
     mixing_index = 1
 ) {
   
-  bitten_humans <- individual::Bitset$new(parameters$human_population)
-  n_bites_each <- NULL
-  
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
     human_infectivity <- account_for_tbv(
@@ -146,9 +143,11 @@ simulate_bites <- function(
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     if(parameters$parasite == "falciparum"){
+      bitten_humans <- individual::Bitset$new(parameters$human_population)
       # p.f model factors eir by psi
       expected_bites <- species_eir * mean(psi)
     } else if (parameters$parasite == "vivax"){
+      bitten_humans <- individual::IntegerVariable$new(initial_values = rep(0, parameters$human_population))
       # p.v model standardises biting rate het to eir
       expected_bites <- species_eir
     }
@@ -157,10 +156,15 @@ simulate_bites <- function(
       n_bites <- rpois(1, expected_bites)
       if (n_bites > 0) {
         bitten <- fast_weighted_sample(n_bites, lambda)
-        bitten_humans$insert(bitten)
+        if(parameters$parasite == "falciparum"){
+          bitten_humans$insert(bitten)
+          renderer$render('n_bitten', bitten_humans$size(), timestep)
+        }
         if(parameters$parasite == "vivax"){
-          # p.v must pass through the number of bites each individual receives
-          n_bites_each <- table(bitten)
+          # p.v must pass through the number of bites per person
+          bites_per_person <- tabulate(bitten, nbins = length(lambda))
+          bitten_humans <- individual::IntegerVariable$new(initial_values = bites_per_person)
+          renderer$render('n_bitten', bitten_humans$get_size_of(!0), timestep)
         }
       }
     }
@@ -212,13 +216,7 @@ simulate_bites <- function(
     }
   }
 
-  renderer$render('n_bitten', bitten_humans$size(), timestep)
-  
-  if(parameters$parasite == "falciparum"){
-    list(bitten_humans = bitten_humans)
-  } else if (parameters$parasite == "vivax"){
-    list(bitten_humans = bitten_humans, n_bites_each = n_bites_each)
-  }
+  bitten_humans
 }
 
 
diff --git a/R/human_infection.R b/R/human_infection.R
index d6266d6a..9eea52cc 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -21,11 +21,11 @@ simulate_infection <- function(
     renderer,
     infection_outcome
 ) {
-  if (bitten_humans$bitten_humans$size() > 0) {
-    if(parameters$parasite == "falciparum"){
+  if(parameters$parasite == "falciparum"){
+    if (bitten_humans$size() > 0) {
       boost_immunity(
         variables$ib,
-        bitten_humans$bitten_humans,
+        bitten_humans,
         variables$last_boosted_ib,
         timestep,
         parameters$ub
@@ -61,20 +61,25 @@ calculate_infections <- function(
     infection_outcome
 ) {
   
-  if(bitten_humans$bitten_humans$size() == 0){return(bitten_humans$bitten_humans)}
-  
   if(parameters$parasite == "falciparum"){
-    source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans$bitten_humans)
+    
+    if(bitten_humans$size() == 0){return(bitten_humans)}
+    
+    source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
     
     ## p.f models blood immunity
     b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
     
   } else if (parameters$parasite == "vivax"){
     ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
-    source_humans <- bitten_humans$bitten_humans$copy()$or(variables$hypnozoites$get_index_of(0)$not(TRUE))
-    bitten_vector <- bitten_humans$bitten_humans$to_vector()
+    bitten_humans_index <- bitten_humans$get_index_of(!0)
+    source_humans <- variables$hypnozoites$get_index_of(0)$not(T)$or(bitten_humans_index)
+
+    if(source_humans$size() == 0){return(source_humans)}
+    
+    # bitten_vector <- bitten_humans$to_vector()
     ## p.v does not model blood immunity but must take into account multiple bites per person
-    b <- 1-(1-parameters$b)^bitten_humans$n_bites_each
+    b <- 1 - (1 - parameters$b)^bitten_humans$get_values(bitten_humans_index)
   }
   source_vector <- source_humans$to_vector()
 
@@ -126,8 +131,9 @@ calculate_infections <- function(
     
   } else if (parameters$parasite == "vivax"){
     ## calculated rate of infection for all bitten or with hypnozoites
+    infection_rates[bitten_humans_index$to_vector()] <- prob_to_rate(b)
     relapse_rates <- variables$hypnozoites$get_values() * parameters$f
-    infection_rates[bitten_vector] <- infection_rates[bitten_vector] + prob_to_rate(b)
+    infection_rates <- infection_rates + relapse_rates
     relative_rates <- relapse_rates/infection_rates
     relative_rates[is.nan(relative_rates)] <- 0
     
@@ -263,9 +269,7 @@ infection_outcome_process <- function(
       lm_det_infections <- calculate_lm_det_infections(
         variables,
         variables$state$get_index_of(c("S","U"))$and(infected_humans),
-        parameters,
-        renderer,
-        timestep
+        parameters
       )
       
       # Lm-detectable level infected S and U, and all A infections may receive clinical infections
@@ -349,69 +353,7 @@ relapse_bite_infection_hazard_resolution <- function(
     
     # make sure batches are capped
     current_batches <- variables$hypnozoites$get_values(new_hypnozoite_batch_formed)
-    new_batch_number <- ifelse(current_batches == parameters$kmax,
-                               current_batches,
-                               current_batches + 1)
-    
-    variables$hypnozoites$queue_update(
-      new_batch_number,
-      new_hypnozoite_batch_formed$and(bite_infections)
-    )
-  }
-}
-
-#' @title Relapse/bite infection competing hazard resolution (p.v only)
-#' @description
-#' Resolves competing hazards of bite and hypnozoite relapse infections. 
-#' For bite infections we increase the batch number and factor in drug prophylaxis.
-#' 
-#' @param variables a list of all of the model variables
-#' @param infected_humans bitset of infected humans
-#' @param relative_rates relative rate of relapse infection
-#' @param variables model variables
-#' @param parameters model parameters
-#' @param renderer model renderer
-#' @param timestep current timestep
-#' @noRd
-relapse_bite_infection_hazard_resolution <- function(
-    infected_humans,
-    relative_rates,
-    variables,
-    parameters,
-    renderer,
-    timestep
-){
-  
-  # draw relapses from total infections
-  relapse_infections <- bitset_at(
-    infected_humans,
-    bernoulli_multi_p(relative_rates[infected_humans$to_vector()])
-  )
-  
-  renderer$render('n_relapses', relapse_infections$size(), timestep)
-  # render relapse infections by age
-  incidence_renderer(
-    variables$birth,
-    renderer,
-    relapse_infections,
-    'inc_relapse_',
-    parameters$incidence_relapse_rendering_min_ages,
-    parameters$incidence_relapse_rendering_max_ages,
-    timestep
-  )
-  
-  # get bite infections
-  bite_infections <- infected_humans$copy()$and(relapse_infections$not(inplace = F))
-  
-  ## all bitten humans with an infectious bite (incorporating prophylaxis) get a new batch of hypnozoites
-  if(bite_infections$size()>0){
-    new_hypnozoite_batch_formed <- bite_infections # bitset_at(bite_infections, bernoulli_multi_p(1-ls_prophylaxis))
-    
-    # make sure batches are capped
-    current_batches <- variables$hypnozoites$get_values(new_hypnozoite_batch_formed)
-    new_batch_number <- ifelse(current_batches == parameters$kmax,
-                               current_batches,
-                               current_batches + 1)
+    new_batch_number <- pmin(current_batches + 1, parameters$kmax)
     
     variables$hypnozoites$queue_update(
       new_batch_number,
@@ -426,15 +368,11 @@ relapse_bite_infection_hazard_resolution <- function(
 #' @param variables a list of all of the model variables
 #' @param infections bitset of infected humans
 #' @param parameters model parameters
-#' @param renderer model render
-#' @param timestep current timestep
 #' @noRd
 calculate_lm_det_infections <- function(
     variables,
     infections,
-    parameters,
-    renderer,
-    timestep
+    parameters
 ) {
   
   iaa <- variables$iaa$get_values(infections)
@@ -443,28 +381,7 @@ calculate_lm_det_infections <- function(
   philm <- anti_parasite_immunity(
     min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
     k = parameters$klm, iaa = iaa, iam = iam)
-  lm_det_infections <- bitset_at(infections, bernoulli_multi_p(philm))
-  
-  incidence_renderer(
-    variables$birth,
-    renderer,
-    lm_det_infections,
-    'inc_lm_det_',
-    parameters$lm_det_incidence_rendering_min_ages,
-    parameters$lm_det_incidence_rendering_max_ages,
-    timestep
-  )
-  incidence_probability_renderer(
-    variables$birth,
-    renderer,
-    infections,
-    philm,
-    'inc_lm_det_',
-    parameters$lm_det_incidence_rendering_min_ages,
-    parameters$lm_det_incidence_rendering_max_ages,
-    timestep
-  )
-  lm_det_infections
+  bitset_at(infections, bernoulli_multi_p(philm))
 }
 
 #' @title Calculate clinical infections
diff --git a/R/parameters.R b/R/parameters.R
index 5af67575..8bcd8ef3 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -280,8 +280,6 @@
 #' * 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
-#' * 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
@@ -467,8 +465,6 @@ 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),
-      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/R/processes.R b/R/processes.R
index 4bec896f..f5e5a7b7 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -372,9 +372,11 @@ create_lagged_eir <- function(variables, solvers, parameters) {
 create_hypnozoite_batch_decay_process <- function(hypnozoites, gammal){
   function(timestep){
     to_decay <- bernoulli_multi_p(p = rate_to_prob(hypnozoites$get_values() * gammal))
-    hypnozoites$queue_update(
-      hypnozoites$get_values(to_decay) - 1,
-      to_decay
-    )
+    if(length(to_decay) > 0){
+      hypnozoites$queue_update(
+        hypnozoites$get_values(to_decay) - 1,
+        to_decay
+      )
+    }
   }
 }
\ No newline at end of file
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 713fa83a..4be80329 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -127,7 +127,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     lagged_eir
   )
 
-  expect_equal(bitten$bitten_humans$to_vector(), c(2, 3))
+  expect_equal(bitten$to_vector(), c(2, 3))
 
   f <- parameters$blood_meal_rates[[1]]
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index ef4cf03c..50c54bf2 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -34,7 +34,7 @@ test_that('simulate_infection integrates different types of infection and schedu
   simulate_infection(
     variables,
     events,
-    list(bitten_humans = bitten),
+    bitten,
     age,
     parameters,
     timestep,
@@ -56,7 +56,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     infection_mock,
     1,
     variables,
-    list(bitten_humans = bitten),
+    bitten,
     parameters,
     renderer,
     timestep,
@@ -215,7 +215,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     depth = 4
   )
   
-  bitten_humans <- list(bitten_humans = individual::Bitset$new(4)$insert(c(1, 2, 3, 4)))
+  bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
@@ -697,7 +697,7 @@ test_that('prophylaxis is considered for medicated humans', {
     ib = individual::DoubleVariable$new(c(.2, .3, .5, .9))
   )
   
-  bitten <- list(bitten_humans = individual::Bitset$new(4)$insert(seq(4)))
+  bitten <- individual::Bitset$new(4)$insert(seq(4))
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
   
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index 2f173edf..9a5f1175 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -173,7 +173,7 @@ test_that('Infection considers pev efficacy', {
 
   infection_rates <- calculate_infections(
     variables = variables,
-    bitten_humans = list(bitten_humans = individual::Bitset$new(4)$insert(seq(4))),
+    bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
     parameters = parameters,
     renderer = mock_render(timestep),
     timestep = timestep,

From 75a48a67ce2a44d4463b64f60c565d54dcfa83da Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 22 Jul 2024 10:12:27 +0100
Subject: [PATCH 05/11] Went back to two biting process outputs: bitset of
 bitten humans and a vector that contains the number of bites each person
 receives.

Adjusted tests to reflect this.
---
 R/biting_process.R                          |  22 ++-
 R/human_infection.R                         | 185 ++++++++++----------
 tests/testthat/test-biting-integration.R    |  10 +-
 tests/testthat/test-infection-integration.R |   7 +
 4 files changed, 117 insertions(+), 107 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index 0149f3f6..30804ba2 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -33,7 +33,7 @@ create_biting_process <- function(
 ) {
   function(timestep){
     age <- get_age(variables$birth$get_values(), timestep)
-    bitten_humans <- simulate_bites(
+    bitten <- simulate_bites(
       renderer,
       solvers,
       models,
@@ -51,7 +51,8 @@ create_biting_process <- function(
     simulate_infection(
       variables,
       events,
-      bitten_humans,
+      bitten$bitten_humans,
+      bitten$n_bites_per_person,
       age,
       parameters,
       timestep,
@@ -76,6 +77,9 @@ simulate_bites <- function(
     mixing = 1,
     mixing_index = 1
 ) {
+
+  bitten_humans <- individual::Bitset$new(parameters$human_population)
+  n_bites_per_person <- numeric(0)
   
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
@@ -143,11 +147,9 @@ simulate_bites <- function(
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     if(parameters$parasite == "falciparum"){
-      bitten_humans <- individual::Bitset$new(parameters$human_population)
       # p.f model factors eir by psi
       expected_bites <- species_eir * mean(psi)
     } else if (parameters$parasite == "vivax"){
-      bitten_humans <- individual::IntegerVariable$new(initial_values = rep(0, parameters$human_population))
       # p.v model standardises biting rate het to eir
       expected_bites <- species_eir
     }
@@ -156,15 +158,11 @@ simulate_bites <- function(
       n_bites <- rpois(1, expected_bites)
       if (n_bites > 0) {
         bitten <- fast_weighted_sample(n_bites, lambda)
-        if(parameters$parasite == "falciparum"){
-          bitten_humans$insert(bitten)
-          renderer$render('n_bitten', bitten_humans$size(), timestep)
-        }
+        bitten_humans$insert(bitten)
+        renderer$render('n_bitten', bitten_humans$size(), timestep)
         if(parameters$parasite == "vivax"){
           # p.v must pass through the number of bites per person
-          bites_per_person <- tabulate(bitten, nbins = length(lambda))
-          bitten_humans <- individual::IntegerVariable$new(initial_values = bites_per_person)
-          renderer$render('n_bitten', bitten_humans$get_size_of(!0), timestep)
+          n_bites_per_person <- tabulate(bitten, nbins = length(lambda))
         }
       }
     }
@@ -216,7 +214,7 @@ simulate_bites <- function(
     }
   }
 
-  bitten_humans
+  list(bitten_humans = bitten_humans, n_bites_per_person = n_bites_per_person)
 }
 
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 9eea52cc..24835bf1 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -5,6 +5,7 @@
 #' @param variables a list of all of the model variables
 #' @param events a list of all of the model events
 #' @param bitten_humans a bitset of bitten humans
+#' @param n_bites_per_person vector of number of bites each person receives (p.v only)
 #' @param age of each human (timesteps)
 #' @param parameters of the model
 #' @param timestep current timestep
@@ -15,6 +16,7 @@ simulate_infection <- function(
     variables,
     events,
     bitten_humans,
+    n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -37,6 +39,7 @@ simulate_infection <- function(
   calculate_infections(
     variables,
     bitten_humans,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
@@ -48,6 +51,7 @@ simulate_infection <- function(
 #' @description Infection rates are stored in the infection outcome competing hazards object
 #' @param variables a list of all of the model variables
 #' @param bitten_humans bitset of bitten humans
+#' @param n_bites_per_person vector of number of bites each person receives (p.v only)
 #' @param parameters model parameters
 #' @param renderer model render object
 #' @param timestep current timestep
@@ -55,110 +59,109 @@ simulate_infection <- function(
 calculate_infections <- function(
     variables,
     bitten_humans,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
     infection_outcome
 ) {
+
+    if(parameters$parasite == "falciparum"){
+      source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
+    } else if (parameters$parasite == "vivax"){
+      ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
+      source_humans <- variables$hypnozoites$get_index_of(0)$not(T)$or(bitten_humans)
+    }
   
-  if(parameters$parasite == "falciparum"){
-    
-    if(bitten_humans$size() == 0){return(bitten_humans)}
+  if(source_humans$size() > 0){
+   
+    source_vector <- source_humans$to_vector()
     
-    source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
+    # calculate prophylaxis
+    prophylaxis <- rep(0, length(source_vector))
+    drug <- variables$drug$get_values(source_vector)
+    medicated <- (drug > 0)
+    if (any(medicated)) {
+      drug <- drug[medicated]
+      drug_time <- variables$drug_time$get_values(source_vector[medicated])
+      prophylaxis[medicated] <- weibull_survival(
+        timestep - drug_time,
+        parameters$drug_prophylaxis_shape[drug],
+        parameters$drug_prophylaxis_scale[drug]
+      )
+    }
     
-    ## p.f models blood immunity
-    b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+    # calculate vaccine efficacy
+    vaccine_efficacy <- rep(0, length(source_vector))
+    vaccine_times <- variables$last_eff_pev_timestep$get_values(source_vector)
+    pev_profile <- variables$pev_profile$get_values(source_vector)
+    # get vector of individuals who have received their 3rd dose
+    vaccinated <- vaccine_times > -1
+    pev_profile <- pev_profile[vaccinated]
+    if (length(vaccinated) > 0) {
+      antibodies <- calculate_pev_antibodies(
+        timestep - vaccine_times[vaccinated],
+        exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'cs')),
+        invlogit(sample_pev_param(pev_profile, parameters$pev_profiles, 'rho')),
+        exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'ds')),
+        exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'dl')),
+        parameters
+      )
+      vmax <- vnapply(parameters$pev_profiles, function(p) p$vmax)
+      beta <- vnapply(parameters$pev_profiles, function(p) p$beta)
+      alpha <- vnapply(parameters$pev_profiles, function(p) p$alpha)
+      vaccine_efficacy[vaccinated] <- calculate_pev_efficacy(
+        antibodies,
+        vmax[pev_profile],
+        beta[pev_profile],
+        alpha[pev_profile]
+      )
+    }
     
-  } else if (parameters$parasite == "vivax"){
-    ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
-    bitten_humans_index <- bitten_humans$get_index_of(!0)
-    source_humans <- variables$hypnozoites$get_index_of(0)$not(T)$or(bitten_humans_index)
-
-    if(source_humans$size() == 0){return(source_humans)}
+    infection_rates <- rep(0, length = parameters$human_population)
+    if(parameters$parasite == "falciparum"){
+      
+      ## p.f models blood immunity
+      b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+      
+      prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+      infection_rates[source_vector] <- prob_to_rate(prob)
+      
+    } else if (parameters$parasite == "vivax"){
+      
+      ## p.v does not model blood immunity but must take into account multiple bites per person
+      b <- 1 - (1 - parameters$b)^n_bites_per_person[bitten_humans$to_vector()]
+      
+      ## calculated rate of infection for all bitten or with hypnozoites
+      infection_rates[bitten_humans$to_vector()] <- prob_to_rate(b)
+      relapse_rates <- variables$hypnozoites$get_values() * parameters$f
+      infection_rates <- infection_rates + relapse_rates
+      relative_rates <- relapse_rates/infection_rates
+      relative_rates[is.nan(relative_rates)] <- 0
+      
+      infection_outcome$set_relative_rates(relative_rates)
+      
+      ## get relative rates to get probability bitten over relapse
+      prob <- rate_to_prob(infection_rates)
+      prob[source_vector] <- prob[source_vector] * (1 - prophylaxis) * (1 - vaccine_efficacy)
+      infection_rates <- prob_to_rate(prob)
+    }
     
-    # bitten_vector <- bitten_humans$to_vector()
-    ## p.v does not model blood immunity but must take into account multiple bites per person
-    b <- 1 - (1 - parameters$b)^bitten_humans$get_values(bitten_humans_index)
-  }
-  source_vector <- source_humans$to_vector()
-
-  # calculate prophylaxis
-  prophylaxis <- rep(0, length(source_vector))
-  drug <- variables$drug$get_values(source_vector)
-  medicated <- (drug > 0)
-  if (any(medicated)) {
-    drug <- drug[medicated]
-    drug_time <- variables$drug_time$get_values(source_vector[medicated])
-    prophylaxis[medicated] <- weibull_survival(
-      timestep - drug_time,
-      parameters$drug_prophylaxis_shape[drug],
-      parameters$drug_prophylaxis_scale[drug]
+    ## probability of incidence must be rendered at each timestep
+    incidence_probability_renderer(
+      variables$birth,
+      renderer,
+      source_humans,
+      prob,
+      "inc_",
+      parameters$incidence_min_ages,
+      parameters$incidence_max_ages,
+      timestep
     )
-  }
-
-  # calculate vaccine efficacy
-  vaccine_efficacy <- rep(0, length(source_vector))
-  vaccine_times <- variables$last_eff_pev_timestep$get_values(source_vector)
-  pev_profile <- variables$pev_profile$get_values(source_vector)
-  # get vector of individuals who have received their 3rd dose
-  vaccinated <- vaccine_times > -1
-  pev_profile <- pev_profile[vaccinated]
-  if (length(vaccinated) > 0) {
-    antibodies <- calculate_pev_antibodies(
-      timestep - vaccine_times[vaccinated],
-      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'cs')),
-      invlogit(sample_pev_param(pev_profile, parameters$pev_profiles, 'rho')),
-      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'ds')),
-      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'dl')),
-      parameters
-    )
-    vmax <- vnapply(parameters$pev_profiles, function(p) p$vmax)
-    beta <- vnapply(parameters$pev_profiles, function(p) p$beta)
-    alpha <- vnapply(parameters$pev_profiles, function(p) p$alpha)
-    vaccine_efficacy[vaccinated] <- calculate_pev_efficacy(
-      antibodies,
-      vmax[pev_profile],
-      beta[pev_profile],
-      alpha[pev_profile]
-    )
-  }
-  
-  infection_rates <- rep(0, length = parameters$human_population)
-  if(parameters$parasite == "falciparum"){
-    prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
-    infection_rates[source_vector] <- prob_to_rate(prob)
     
-  } else if (parameters$parasite == "vivax"){
-    ## calculated rate of infection for all bitten or with hypnozoites
-    infection_rates[bitten_humans_index$to_vector()] <- prob_to_rate(b)
-    relapse_rates <- variables$hypnozoites$get_values() * parameters$f
-    infection_rates <- infection_rates + relapse_rates
-    relative_rates <- relapse_rates/infection_rates
-    relative_rates[is.nan(relative_rates)] <- 0
-    
-    infection_outcome$set_relative_rates(relative_rates)
-    
-    ## get relative rates to get probability bitten over relapse
-    prob <- rate_to_prob(infection_rates)
-    prob[source_vector] <- prob[source_vector] * (1 - prophylaxis) * (1 - vaccine_efficacy)
-    infection_rates <- prob_to_rate(prob)
+    ## capture infection rates to resolve in competing hazards
+    infection_outcome$set_rates(infection_rates)
   }
-
-  ## probability of incidence must be rendered at each timestep
-  incidence_probability_renderer(
-    variables$birth,
-    renderer,
-    source_humans,
-    prob,
-    "inc_",
-    parameters$incidence_min_ages,
-    parameters$incidence_max_ages,
-    timestep
-  )
-  
-  ## capture infection rates to resolve in competing hazards
-  infection_outcome$set_rates(infection_rates)
 }
 
 #' @title Assigns infections to appropriate human states
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 4be80329..fe73d1db 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -35,8 +35,9 @@ test_that('biting_process integrates mosquito effects and human infection', {
     infection_outcome
   )
   
-  bitten <- individual::Bitset$new(parameters$human_population)
-  bites_mock <- mockery::mock(bitten)
+  bitten <- list(bitten_humans = individual::Bitset$new(parameters$human_population),
+                 n_bites_per_person = numeric(0))
+  bites_mock <- mockery::mock(bitten, cycle = T)
   infection_mock <- mockery::mock()
 
   mockery::stub(biting_process, 'simulate_bites', bites_mock)
@@ -65,7 +66,8 @@ test_that('biting_process integrates mosquito effects and human infection', {
     1,
     variables,
     events,
-    bitten,
+    bitten$bitten,
+    bitten$n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -127,7 +129,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     lagged_eir
   )
 
-  expect_equal(bitten$to_vector(), c(2, 3))
+  expect_equal(bitten$bitten_humans$to_vector(), c(2, 3))
 
   f <- parameters$blood_meal_rates[[1]]
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 50c54bf2..8554f9b6 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -17,6 +17,7 @@ test_that('simulate_infection integrates different types of infection and schedu
   )
   
   bitten <- individual::Bitset$new(population)$insert(c(1, 3, 5, 7))
+  n_bites_per_person <- numeric(0)
   boost_immunity_mock <- mockery::mock()
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
   infection_mock <- mockery::mock(infected)
@@ -35,6 +36,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     variables,
     events,
     bitten,
+    n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -57,6 +59,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     1,
     variables,
     bitten,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
@@ -216,6 +219,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   )
   
   bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  n_bites_per_person <- numeric(0)
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
@@ -227,6 +231,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   infections <- calculate_infections(
     variables,
     bitten_humans, 
+    n_bites_per_person,
     parameters,
     mock_render(timestep),
     timestep,
@@ -698,6 +703,7 @@ test_that('prophylaxis is considered for medicated humans', {
   )
   
   bitten <- individual::Bitset$new(4)$insert(seq(4))
+  n_bites_per_person <- numeric(0)
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
   
@@ -711,6 +717,7 @@ test_that('prophylaxis is considered for medicated humans', {
   infection_rates <- calculate_infections(
     variables,
     bitten,
+    n_bites_per_person,
     parameters,
     mock_render(timestep),
     timestep,

From a450963c0f76c4937301c896f155d13f2ec27442 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Thu, 25 Jul 2024 10:18:56 +0100
Subject: [PATCH 06/11] Relapses can sometimes be 0, so must be included in the
 initial rendering to avoid NAs in the render object.

---
 R/render.R | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/R/render.R b/R/render.R
index 0d066b46..e0b03087 100644
--- a/R/render.R
+++ b/R/render.R
@@ -242,6 +242,11 @@ populate_incidence_rendering_columns <- function(renderer, parameters){
     renderer$set_default('n_early_treatment_failure', 0)
     renderer$set_default('n_slow_parasite_clearance', 0)
   }
+
+  # relapses only render for the vivax model
+  if(parameters$parasite == "vivax"){
+    renderer$set_default('n_relapses', 0)
+  }
   
   if(length(parameters$incidence_rendering_min_ages)>0){
     render_initial_incidence(renderer,

From 11a2c2c4ef348ad5eec87b8a0ae8c2403800bc80 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 1 Jul 2024 10:23:53 +0100
Subject: [PATCH 07/11] We must model hynozoite batch formation via bite
 infections, decay of hypnozoite batches through, new infections via relapse
 and the reseting of batch number to 0 for new births.

We create a hypnozoite variable.

The vivax model considers the number of new bites an individual receives (in contrast to whether an individual has been bitten regardless of number of bites) to calculate the infection rates. This is now captured in bitten_humans, which is now a list that also includes n_bites_each.

The vivax model does not factor the eir by the population level biting heterogeneity (although it does use the heterogeneity to distribute biting, but retains the overall total biting according to the EIR).
biting_rate.R

In vivax, we track bites in all individuals, not just SAU as in falciparum. This is because new hypnozoite batches can occur in any human state, regardless of whether this results in a new infection.

The number of bites is considered in calculating probability of a new infection.

To calculate the infection rates for vivax, we take the biting infection rates and add them to the relase rate. We also save the relative rates in the competing hazards object to resolve these competing hazards later.

A new function to resolve between bites and relapse infections is created. For each bite we increase hte number of batches - with capacity for prophylaxis to be added at a future step. We also cap the total number of batches.

We add a hypnozoite batch decay process. Note that this process differs from the immunity decay processes. We work in discrete outputs (you can't have a fractional batch). We also work in the probability that an individual experiences a single loss, rather than loss at the hypnozoite level.

When an individual dies, we reset the hypnozoite batches back to 0.

We now render number of relapses, average hypnozoite batch number and the number of individuals with hypnozoite batches. We may also output these by specified age groups.

Added a test to check relapses occur as expected.
---
 R/biting_process.R                          |  28 +++-
 R/competing_hazards.R                       |   8 +-
 R/human_infection.R                         | 138 ++++++++++++++++----
 R/mortality_processes.R                     |   1 +
 R/processes.R                               |  40 +++++-
 R/render.R                                  |  35 +++++
 R/variables.R                               |   7 +-
 tests/testthat/test-biting-integration.R    |   2 +-
 tests/testthat/test-infection-integration.R |  13 +-
 tests/testthat/test-pev.R                   |   2 +-
 tests/testthat/test-vivax.R                 |  60 +++++++++
 11 files changed, 290 insertions(+), 44 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index 594b9429..a2c16e68 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -78,6 +78,7 @@ simulate_bites <- function(
 ) {
   
   bitten_humans <- individual::Bitset$new(parameters$human_population)
+  n_bites_each <- NULL
   
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
@@ -144,13 +145,23 @@ simulate_bites <- function(
     
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
-    expected_bites <- species_eir * mean(psi)
+    if(parameters$parasite == "falciparum"){
+      # p.f model factors eir by psi
+      expected_bites <- species_eir * mean(psi)
+    } else if (parameters$parasite == "vivax"){
+      # p.v model standardises biting rate het to eir
+      expected_bites <- species_eir
+    }
+    
     if (expected_bites > 0) {
       n_bites <- rpois(1, expected_bites)
       if (n_bites > 0) {
-        bitten_humans$insert(
-          fast_weighted_sample(n_bites, lambda)
-        )
+        bitten <- fast_weighted_sample(n_bites, lambda)
+        bitten_humans$insert(bitten)
+        if(parameters$parasite == "vivax"){
+          # p.v must pass through the number of bites each individual receives
+          n_bites_each <- table(bitten)
+        }
       }
     }
     
@@ -200,9 +211,14 @@ simulate_bites <- function(
       )
     }
   }
-  
+
   renderer$render('n_bitten', bitten_humans$size(), timestep)
-  bitten_humans
+  
+  if(parameters$parasite == "falciparum"){
+    list(bitten_humans = bitten_humans)
+  } else if (parameters$parasite == "vivax"){
+    list(bitten_humans = bitten_humans, n_bites_each = n_bites_each)
+  }
 }
 
 
diff --git a/R/competing_hazards.R b/R/competing_hazards.R
index cdc6aea0..4f429824 100644
--- a/R/competing_hazards.R
+++ b/R/competing_hazards.R
@@ -14,15 +14,21 @@ CompetingOutcome <- R6::R6Class(
       }
       private$targeted_process <- targeted_process
       self$rates <- rep(0, size)
+      self$relative_rates <- rep(0, size)
     },
     set_rates = function(rates){
       self$rates <- rates
     },
+    set_relative_rates = function(relative_rates){
+      self$relative_rates <- relative_rates
+    },
     execute = function(t, target){
       private$targeted_process(t, target)
       self$rates <- rep(0, length(self$rates))
+      self$relative_rates <- rep(0, length(self$relative_rates))
     },
-    rates = NULL
+    rates = NULL,
+    relative_rates = NULL
   )
 )
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 57d301f4..9976a198 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -21,11 +21,11 @@ simulate_infection <- function(
     renderer,
     infection_outcome
 ) {
-  if (bitten_humans$size() > 0) {
+  if (bitten_humans$bitten_humans$size() > 0) {
     if(parameters$parasite == "falciparum"){
       boost_immunity(
         variables$ib,
-        bitten_humans,
+        bitten_humans$bitten_humans,
         variables$last_boosted_ib,
         timestep,
         parameters$ub
@@ -60,18 +60,22 @@ calculate_infections <- function(
     timestep,
     infection_outcome
 ) {
-  source_humans <- variables$state$get_index_of(
-    c('S', 'A', 'U'))$and(bitten_humans)
-
+  
+  if(bitten_humans$bitten_humans$size() == 0){return(bitten_humans$bitten_humans)}
+  
   if(parameters$parasite == "falciparum"){
+    source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans$bitten_humans)
+    
     ## p.f models blood immunity
     b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
     
   } else if (parameters$parasite == "vivax"){
-    ## p.v does not model blood immunity
-    b <- parameters$b
+    ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
+    source_humans <- bitten_humans$bitten_humans$copy()$or(variables$hypnozoites$get_index_of(0)$not(TRUE))
+    bitten_vector <- bitten_humans$bitten_humans$to_vector()
+    ## p.v does not model blood immunity but must take into account multiple bites per person
+    b <- 1-(1-parameters$b)^bitten_humans$n_bites_each
   }
-
   source_vector <- source_humans$to_vector()
 
   # calculate prophylaxis
@@ -114,8 +118,26 @@ calculate_infections <- function(
       alpha[pev_profile]
     )
   }
-
-  prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+  
+  infection_rates <- rep(0, length = parameters$human_population)
+  if(parameters$parasite == "falciparum"){
+    prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    infection_rates[source_vector] <- prob_to_rate(prob)
+    
+  } else if (parameters$parasite == "vivax"){
+    ## calculated rate of infection for all bitten or with hypnozoites
+    relapse_rates <- variables$hypnozoites$get_values() * parameters$f
+    infection_rates[bitten_vector] <- infection_rates[bitten_vector] + prob_to_rate(b)
+    relative_rates <- relapse_rates/infection_rates
+    relative_rates[is.nan(relative_rates)] <- 0
+    
+    infection_outcome$set_relative_rates(relative_rates)
+    
+    ## get relative rates to get probability bitten over relapse
+    prob <- rate_to_prob(infection_rates)
+    prob[source_vector] <- prob[source_vector] * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    infection_rates <- prob_to_rate(prob)
+  }
 
   ## probability of incidence must be rendered at each timestep
   incidence_probability_renderer(
@@ -130,8 +152,6 @@ calculate_infections <- function(
   )
   
   ## capture infection rates to resolve in competing hazards
-  infection_rates <- rep(0, length = parameters$human_population)
-  infection_rates[source_vector] <- prob_to_rate(prob)
   infection_outcome$set_rates(infection_rates)
 }
 
@@ -145,6 +165,7 @@ calculate_infections <- function(
 #' @param renderer model render object
 #' @param parameters model parameters
 #' @param prob vector of population probabilities of infection
+#' @param relative_rates relative rates of hypnozoite relapse relative to total infection rate
 #' @noRd
 infection_outcome_process <- function(
     timestep,
@@ -152,21 +173,21 @@ infection_outcome_process <- function(
     variables,
     renderer,
     parameters,
-    prob){
-  
-  incidence_renderer(
-    variables$birth,
-    renderer,
-    infected_humans,
-    'inc_',
-    parameters$incidence_rendering_min_ages,
-    parameters$incidence_rendering_max_ages,
-    timestep
-  )
+    prob,
+    relative_rates = NULL){
   
   if (infected_humans$size() > 0) {
-    
+
     renderer$render('n_infections', infected_humans$size(), timestep)
+    incidence_renderer(
+      variables$birth,
+      renderer,
+      infected_humans,
+      'inc_',
+      parameters$incidence_rendering_min_ages,
+      parameters$incidence_rendering_max_ages,
+      timestep
+    )
     
     boost_immunity(
       variables$ica,
@@ -215,6 +236,14 @@ infection_outcome_process <- function(
       to_U <- NULL
       
     } else if (parameters$parasite == "vivax"){
+      relapse_bite_infection_hazard_resolution(
+        infected_humans,
+        relative_rates,
+        variables,
+        parameters,
+        renderer,
+        timestep
+      )
       
       boost_immunity(
         variables$iaa,
@@ -258,6 +287,7 @@ infection_outcome_process <- function(
       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(
@@ -271,6 +301,66 @@ infection_outcome_process <- function(
   }
 }
 
+#' @title Relapse/bite infection competing hazard resolution (p.v only)
+#' @description
+#' Resolves competing hazards of bite and hypnozoite relapse infections. 
+#' For bite infections we increase the batch number and factor in drug prophylaxis.
+#' 
+#' @param variables a list of all of the model variables
+#' @param infected_humans bitset of infected humans
+#' @param relative_rates relative rate of relapse infection
+#' @param variables model variables
+#' @param parameters model parameters
+#' @param renderer model renderer
+#' @param timestep current timestep
+#' @noRd
+relapse_bite_infection_hazard_resolution <- function(
+    infected_humans,
+    relative_rates,
+    variables,
+    parameters,
+    renderer,
+    timestep
+){
+  
+  # draw relapses from total infections
+  relapse_infections <- bitset_at(
+    infected_humans,
+    bernoulli_multi_p(relative_rates[infected_humans$to_vector()])
+  )
+  
+  renderer$render('n_relapses', relapse_infections$size(), timestep)
+  # render relapse infections by age
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    relapse_infections,
+    'inc_relapse_',
+    parameters$incidence_relapse_rendering_min_ages,
+    parameters$incidence_relapse_rendering_max_ages,
+    timestep
+  )
+  
+  # get bite infections
+  bite_infections <- infected_humans$copy()$and(relapse_infections$not(inplace = F))
+  
+  ## all bitten humans with an infectious bite (incorporating prophylaxis) get a new batch of hypnozoites
+  if(bite_infections$size()>0){
+    new_hypnozoite_batch_formed <- bite_infections # bitset_at(bite_infections, bernoulli_multi_p(1-ls_prophylaxis))
+    
+    # make sure batches are capped
+    current_batches <- variables$hypnozoites$get_values(new_hypnozoite_batch_formed)
+    new_batch_number <- ifelse(current_batches == parameters$kmax,
+                               current_batches,
+                               current_batches + 1)
+    
+    variables$hypnozoites$queue_update(
+      new_batch_number,
+      new_hypnozoite_batch_formed$and(bite_infections)
+    )
+  }
+}
+
 #' @title Calculate light microscopy detectable infections (p.v only)
 #' @description
 #' Sample light microscopy detectable infections from all infections
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index e4fa323f..a66c3eb2 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -112,6 +112,7 @@ reset_target <- function(variables, events, target, state, parameters, timestep)
     } else if (parameters$parasite == "vivax"){
       variables$last_boosted_iaa$queue_update(-1, target)
       variables$iaa$queue_update(0, target)
+      variables$hypnozoites$queue_update(0, target)
     }
 
     # treatment
diff --git a/R/processes.R b/R/processes.R
index afa74f50..2b0adef7 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -59,7 +59,8 @@ create_processes <- function(
       processes,
       # Anti-parasite immunity
       create_exponential_decay_process(variables$iam, parameters$rm),
-      create_exponential_decay_process(variables$iaa, parameters$ra)
+      create_exponential_decay_process(variables$iaa, parameters$ra),
+      create_hypnozoite_batch_decay_process(variables$hypnozoites, parameters$gammal)
     )
   }
 
@@ -84,7 +85,8 @@ create_processes <- function(
     targeted_process = function(timestep, target){
       infection_outcome_process(timestep, target, 
                                 variables, renderer, parameters, 
-                                prob = rate_to_prob(infection_outcome$rates))
+                                prob = rate_to_prob(infection_outcome$rates),
+                                relative_rates = infection_outcome$relative_rates)
     },
     size = parameters$human_population
   )
@@ -189,7 +191,23 @@ create_processes <- function(
   if(parameters$parasite == "falciparum"){
     imm_var_names <- c(imm_var_names, 'ib', 'iva', 'ivm', 'id')
   } else if (parameters$parasite == "vivax"){
-    imm_var_names <- c(imm_var_names, 'iaa', 'iam')
+    imm_var_names <- c(imm_var_names, 'iaa', 'iam','hypnozoites')
+    
+    ## hypnozoite infection prevalence rendering
+    processes <- c(
+      processes,
+      create_n_with_hypnozoites_renderer_process(
+        renderer,
+        variables$hypnozoites,
+        parameters
+      ),
+      create_n_with_hypnozoites_age_renderer_process(
+        variables$hypnozoites,
+        variables$birth,
+        parameters,
+        renderer
+      )
+    )
   }
 
   processes <- c(
@@ -344,3 +362,19 @@ create_lagged_eir <- function(variables, solvers, parameters) {
     }
   )
 }
+#' @title Hypnozoite decay function
+#' @description
+#' calulates the number of individuals in whom a batch decay occurs
+#'
+#' @param variable the hypnozoite variable to update
+#' @param rate the hypnozoite decay rate
+#' @noRd
+create_hypnozoite_batch_decay_process <- function(hypnozoites, gammal){
+  function(timestep){
+    to_decay <- bernoulli_multi_p(p = rate_to_prob(hypnozoites$get_values() * gammal))
+    hypnozoites$queue_update(
+      hypnozoites$get_values(to_decay) - 1,
+      to_decay
+    )
+  }
+}
\ No newline at end of file
diff --git a/R/render.R b/R/render.R
index b219cdb9..0d066b46 100644
--- a/R/render.R
+++ b/R/render.R
@@ -278,4 +278,39 @@ populate_metapopulation_incidence_rendering_columns <- function(renderer, parame
   for(i in length(parameters)){
     populate_incidence_rendering_columns(renderer[[i]], parameters[[i]])
   }
+}
+
+create_n_with_hypnozoites_renderer_process <- function(
+    renderer,
+    hypnozoites,
+    parameters
+) {
+  function(timestep) {
+    renderer$render(
+      paste0("n_with_hypnozoites"),
+      hypnozoites$size() - hypnozoites$get_size_of(0),
+      timestep
+    )
+  }
+}
+
+create_n_with_hypnozoites_age_renderer_process <- function(
+    hypnozoites,
+    birth,
+    parameters,
+    renderer
+) {
+  function(timestep) {
+    for (i in seq_along(parameters$n_with_hypnozoites_rendering_min_ages)) {
+      lower <- parameters$n_with_hypnozoites_rendering_min_ages[[i]]
+      upper <- parameters$n_with_hypnozoites_rendering_max_ages[[i]]
+      in_age <- in_age_range(birth, timestep, lower, upper)
+      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+      renderer$render(
+        paste0('n_with_hypnozoites_', lower, '_', upper),
+        sum(hypnozoites$get_values(index = in_age)!=0),
+        timestep
+      )
+    }
+  }
 }
\ No newline at end of file
diff --git a/R/variables.R b/R/variables.R
index 7e954a40..757807f3 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -18,6 +18,7 @@
 #' * ICA  - Acquired immunity to clinical disease
 #' * IVA  - Acquired immunity to severe disease (p.f only)
 #' * ID - Acquired immunity to detectability (p.f only)
+#' * hypnozoites - Hypnozoite batch number (p.v only)
 #' * zeta - Heterogeneity of human individuals
 #' * zeta_group - Discretised heterogeneity of human individuals
 #' * last_pev_timestep - The timestep of the last pev vaccination (-1 if there
@@ -98,6 +99,9 @@ create_variables <- function(parameters) {
   initial_states <- initial_state(parameters, initial_age, groups, eq, states)
   state <- individual::CategoricalVariable$new(states, initial_states)
   birth <- individual::IntegerVariable$new(-initial_age)
+  if(parameters$parasite == "vivax"){
+    hypnozoites <- individual::IntegerVariable$new(rep(parameters$init_hyp, parameters$human_population))
+  }
   
   # Maternal immunity to clinical disease
   icm <- individual::DoubleVariable$new(
@@ -299,7 +303,8 @@ create_variables <- function(parameters) {
     variables <- c(variables,
                    last_boosted_iaa = last_boosted_iaa,
                    iaa = iaa,
-                   iam = iam
+                   iam = iam,
+                   hypnozoites = hypnozoites
     )
   }
   
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 4be80329..713fa83a 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -127,7 +127,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     lagged_eir
   )
 
-  expect_equal(bitten$to_vector(), c(2, 3))
+  expect_equal(bitten$bitten_humans$to_vector(), c(2, 3))
 
   f <- parameters$blood_meal_rates[[1]]
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index bf72fedd..1628798f 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -34,7 +34,7 @@ test_that('simulate_infection integrates different types of infection and schedu
   simulate_infection(
     variables,
     events,
-    bitten,
+    list(bitten_humans = bitten),
     age,
     parameters,
     timestep,
@@ -56,7 +56,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     infection_mock,
     1,
     variables,
-    bitten,
+    list(bitten_humans = bitten),
     parameters,
     renderer,
     timestep,
@@ -64,7 +64,6 @@ test_that('simulate_infection integrates different types of infection and schedu
   )
 })
 
-
 test_that('simulate_infection integrates different types of infection and scheduling', {
   population <- 8
   timestep <- 5
@@ -145,7 +144,6 @@ 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,
@@ -156,6 +154,8 @@ test_that('simulate_infection integrates different types of infection and schedu
     renderer
   )
   
+  mockery::mock_args(schedule_mock)
+  
   mockery::expect_args(
     schedule_mock,
     1,
@@ -219,7 +219,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     depth = 4
   )
   
-  bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  bitten_humans <- list(bitten_humans = individual::Bitset$new(4)$insert(c(1, 2, 3, 4)))
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
@@ -705,11 +705,10 @@ test_that('prophylaxis is considered for medicated humans', {
     ib = individual::DoubleVariable$new(c(.2, .3, .5, .9))
   )
   
-  bitten <- individual::Bitset$new(4)$insert(seq(4))
+  bitten <- list(bitten_humans = individual::Bitset$new(4)$insert(seq(4)))
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
   
-  
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
       infection_outcome_process(timestep, target, variables, renderer, parameters)
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index 9a5f1175..2f173edf 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -173,7 +173,7 @@ test_that('Infection considers pev efficacy', {
 
   infection_rates <- calculate_infections(
     variables = variables,
-    bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
+    bitten_humans = list(bitten_humans = individual::Bitset$new(4)$insert(seq(4))),
     parameters = parameters,
     renderer = mock_render(timestep),
     timestep = timestep,
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 283e76e5..f73e5bb3 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -255,3 +255,63 @@ test_that('vivax schedule_infections correctly schedules new infections', {
     c(1)
   )
 })
+
+test_that('relapses are recognised with division between bite infections and relapses', {
+  timestep <- 50
+  parameters <- get_parameters(parasite = "vivax", overrides = list(human_population = 4))
+  
+  variables <- list(
+    state = individual::CategoricalVariable$new(
+      c('D', 'S', 'A', 'U', 'Tr'),
+      c('D', 'S', 'A', 'U')
+    ),
+    infectivity = individual::DoubleVariable$new(rep(0, 4)),
+    recovery_rates = individual::DoubleVariable$new(rep(0, 4)),
+    drug = individual::DoubleVariable$new(rep(0, 4)),
+    drug_time = individual::DoubleVariable$new(rep(-1, 4)),
+    iaa = individual::DoubleVariable$new(rep(0, 4)),
+    iam = individual::DoubleVariable$new(rep(0, 4)),
+    ica = individual::DoubleVariable$new(rep(0, 4)),
+    icm = individual::DoubleVariable$new(rep(0, 4)),
+    last_boosted_iaa = individual::DoubleVariable$new(rep(-1, 4)),
+    last_boosted_ica = individual::DoubleVariable$new(rep(-1, 4)),
+    last_eff_pev_timestep = individual::DoubleVariable$new(rep(-1, 4)),
+    pev_profile = individual::IntegerVariable$new(rep(-1, 4)),
+    hypnozoites = individual::IntegerVariable$new(c(0, 1, 2, 3))
+  )
+  
+  bernoulli_mock <- mockery::mock(c(2, 4), 2, cycle = TRUE)
+  calc_mock <- mockery::mock(individual::Bitset$new(4)$insert(2))
+  mockery::stub(infection_outcome_process, 'bernoulli_multi_p', bernoulli_mock, depth = 2)
+  mockery::stub(infection_outcome_process, 'calculate_clinical_infections', calc_mock)
+  
+  renderer <- mock_render(1)
+  infected_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  
+  infection_outcome_process(
+    timestep = timestep,
+    infected_humans,
+    variables,
+    renderer,
+    parameters,
+    prob = c(rep(0.5, 4)),
+    relative_rate = c(rep(0.5, 4))
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_infections',
+    4,
+    timestep
+  )
+
+  mockery::expect_args(
+    renderer$render_mock(),
+    2,
+    'n_relapses',
+    2,
+    timestep
+  )
+  
+})

From 820d28087115e3ad6bd749eec5af200db9cc2a8f Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 1 Jul 2024 10:23:53 +0100
Subject: [PATCH 08/11] We must model hynozoite batch formation via bite
 infections, decay of hypnozoite batches through, new infections via relapse
 and the reseting of batch number to 0 for new births.

We create a hypnozoite variable.

The vivax model considers the number of new bites an individual receives (in contrast to whether an individual has been bitten regardless of number of bites) to calculate the infection rates. This is now captured in bitten_humans, which is now a list that also includes n_bites_each.

The vivax model does not factor the eir by the population level biting heterogeneity (although it does use the heterogeneity to distribute biting, but retains the overall total biting according to the EIR).
biting_rate.R

In vivax, we track bites in all individuals, not just SAU as in falciparum. This is because new hypnozoite batches can occur in any human state, regardless of whether this results in a new infection.

The number of bites is considered in calculating probability of a new infection.

To calculate the infection rates for vivax, we take the biting infection rates and add them to the relase rate. We also save the relative rates in the competing hazards object to resolve these competing hazards later.

A new function to resolve between bites and relapse infections is created. For each bite we increase hte number of batches - with capacity for prophylaxis to be added at a future step. We also cap the total number of batches.

We add a hypnozoite batch decay process. Note that this process differs from the immunity decay processes. We work in discrete outputs (you can't have a fractional batch). We also work in the probability that an individual experiences a single loss, rather than loss at the hypnozoite level.

When an individual dies, we reset the hypnozoite batches back to 0.

We now render number of relapses, average hypnozoite batch number and the number of individuals with hypnozoite batches. We may also output these by specified age groups.

Added a test to check relapses occur as expected.
---
 R/human_infection.R | 2 +-
 R/processes.R       | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/human_infection.R b/R/human_infection.R
index 9976a198..92abe126 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -175,7 +175,7 @@ infection_outcome_process <- function(
     parameters,
     prob,
     relative_rates = NULL){
-  
+
   if (infected_humans$size() > 0) {
 
     renderer$render('n_infections', infected_humans$size(), timestep)
diff --git a/R/processes.R b/R/processes.R
index 2b0adef7..80fa6912 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -192,7 +192,7 @@ create_processes <- function(
     imm_var_names <- c(imm_var_names, 'ib', 'iva', 'ivm', 'id')
   } else if (parameters$parasite == "vivax"){
     imm_var_names <- c(imm_var_names, 'iaa', 'iam','hypnozoites')
-    
+
     ## hypnozoite infection prevalence rendering
     processes <- c(
       processes,

From eb8308a4972d946ed1b62a4efe59aae1731c4e78 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Wed, 17 Jul 2024 12:08:42 +0100
Subject: [PATCH 09/11] Removed n_bites_each, opting instead to create a
 variable  that contains the vector of number of bites per individual that can
 be passed in place of the individuals bitten bitset. I think this simplifies
 a lot of the code, but there may still be ways of improving this.

The relapse bite hazard resolution function was replicated - this has been removed. lm-det incidence rendering has been removed. An if statement has been added to the hypnozoite decay process, and tests adjusted to reflect bitten_humans changes.
---
 R/biting_process.R                          | 24 ++++----
 R/human_infection.R                         | 63 +++++++--------------
 R/parameters.R                              |  4 --
 R/processes.R                               | 10 ++--
 tests/testthat/test-biting-integration.R    |  2 +-
 tests/testthat/test-infection-integration.R |  8 +--
 tests/testthat/test-pev.R                   |  2 +-
 7 files changed, 43 insertions(+), 70 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index a2c16e68..0149f3f6 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -77,9 +77,6 @@ simulate_bites <- function(
     mixing_index = 1
 ) {
   
-  bitten_humans <- individual::Bitset$new(parameters$human_population)
-  n_bites_each <- NULL
-  
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
     human_infectivity <- account_for_tbv(
@@ -146,9 +143,11 @@ simulate_bites <- function(
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     if(parameters$parasite == "falciparum"){
+      bitten_humans <- individual::Bitset$new(parameters$human_population)
       # p.f model factors eir by psi
       expected_bites <- species_eir * mean(psi)
     } else if (parameters$parasite == "vivax"){
+      bitten_humans <- individual::IntegerVariable$new(initial_values = rep(0, parameters$human_population))
       # p.v model standardises biting rate het to eir
       expected_bites <- species_eir
     }
@@ -157,10 +156,15 @@ simulate_bites <- function(
       n_bites <- rpois(1, expected_bites)
       if (n_bites > 0) {
         bitten <- fast_weighted_sample(n_bites, lambda)
-        bitten_humans$insert(bitten)
+        if(parameters$parasite == "falciparum"){
+          bitten_humans$insert(bitten)
+          renderer$render('n_bitten', bitten_humans$size(), timestep)
+        }
         if(parameters$parasite == "vivax"){
-          # p.v must pass through the number of bites each individual receives
-          n_bites_each <- table(bitten)
+          # p.v must pass through the number of bites per person
+          bites_per_person <- tabulate(bitten, nbins = length(lambda))
+          bitten_humans <- individual::IntegerVariable$new(initial_values = bites_per_person)
+          renderer$render('n_bitten', bitten_humans$get_size_of(!0), timestep)
         }
       }
     }
@@ -212,13 +216,7 @@ simulate_bites <- function(
     }
   }
 
-  renderer$render('n_bitten', bitten_humans$size(), timestep)
-  
-  if(parameters$parasite == "falciparum"){
-    list(bitten_humans = bitten_humans)
-  } else if (parameters$parasite == "vivax"){
-    list(bitten_humans = bitten_humans, n_bites_each = n_bites_each)
-  }
+  bitten_humans
 }
 
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 92abe126..0df45355 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -21,11 +21,11 @@ simulate_infection <- function(
     renderer,
     infection_outcome
 ) {
-  if (bitten_humans$bitten_humans$size() > 0) {
-    if(parameters$parasite == "falciparum"){
+  if(parameters$parasite == "falciparum"){
+    if (bitten_humans$size() > 0) {
       boost_immunity(
         variables$ib,
-        bitten_humans$bitten_humans,
+        bitten_humans,
         variables$last_boosted_ib,
         timestep,
         parameters$ub
@@ -61,20 +61,25 @@ calculate_infections <- function(
     infection_outcome
 ) {
   
-  if(bitten_humans$bitten_humans$size() == 0){return(bitten_humans$bitten_humans)}
-  
   if(parameters$parasite == "falciparum"){
-    source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans$bitten_humans)
+    
+    if(bitten_humans$size() == 0){return(bitten_humans)}
+    
+    source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
     
     ## p.f models blood immunity
     b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
     
   } else if (parameters$parasite == "vivax"){
     ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
-    source_humans <- bitten_humans$bitten_humans$copy()$or(variables$hypnozoites$get_index_of(0)$not(TRUE))
-    bitten_vector <- bitten_humans$bitten_humans$to_vector()
+    bitten_humans_index <- bitten_humans$get_index_of(!0)
+    source_humans <- variables$hypnozoites$get_index_of(0)$not(T)$or(bitten_humans_index)
+
+    if(source_humans$size() == 0){return(source_humans)}
+    
+    # bitten_vector <- bitten_humans$to_vector()
     ## p.v does not model blood immunity but must take into account multiple bites per person
-    b <- 1-(1-parameters$b)^bitten_humans$n_bites_each
+    b <- 1 - (1 - parameters$b)^bitten_humans$get_values(bitten_humans_index)
   }
   source_vector <- source_humans$to_vector()
 
@@ -126,8 +131,9 @@ calculate_infections <- function(
     
   } else if (parameters$parasite == "vivax"){
     ## calculated rate of infection for all bitten or with hypnozoites
+    infection_rates[bitten_humans_index$to_vector()] <- prob_to_rate(b)
     relapse_rates <- variables$hypnozoites$get_values() * parameters$f
-    infection_rates[bitten_vector] <- infection_rates[bitten_vector] + prob_to_rate(b)
+    infection_rates <- infection_rates + relapse_rates
     relative_rates <- relapse_rates/infection_rates
     relative_rates[is.nan(relative_rates)] <- 0
     
@@ -257,9 +263,7 @@ infection_outcome_process <- function(
       lm_detectable <- calculate_lm_det_infections(
         variables,
         variables$state$get_index_of(c("S","U"))$and(infected_humans),
-        parameters,
-        renderer,
-        timestep
+        parameters
       )
       
       # Lm-detectable level infected S and U, and all A infections may receive clinical infections
@@ -350,9 +354,7 @@ relapse_bite_infection_hazard_resolution <- function(
     
     # make sure batches are capped
     current_batches <- variables$hypnozoites$get_values(new_hypnozoite_batch_formed)
-    new_batch_number <- ifelse(current_batches == parameters$kmax,
-                               current_batches,
-                               current_batches + 1)
+    new_batch_number <- pmin(current_batches + 1, parameters$kmax)
     
     variables$hypnozoites$queue_update(
       new_batch_number,
@@ -367,15 +369,11 @@ relapse_bite_infection_hazard_resolution <- function(
 #' @param variables a list of all of the model variables
 #' @param infections bitset of infected humans
 #' @param parameters model parameters
-#' @param renderer model render
-#' @param timestep current timestep
 #' @noRd
 calculate_lm_det_infections <- function(
     variables,
     infections,
-    parameters,
-    renderer,
-    timestep
+    parameters
 ) {
   
   iaa <- variables$iaa$get_values(infections)
@@ -384,28 +382,7 @@ calculate_lm_det_infections <- function(
   philm <- anti_parasite_immunity(
     min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
     k = parameters$klm, iaa = iaa, iam = iam)
-  lm_det_infections <- bitset_at(infections, bernoulli_multi_p(philm))
-  
-  incidence_renderer(
-    variables$birth,
-    renderer,
-    lm_det_infections,
-    'inc_lm_det_',
-    parameters$lm_det_incidence_rendering_min_ages,
-    parameters$lm_det_incidence_rendering_max_ages,
-    timestep
-  )
-  incidence_probability_renderer(
-    variables$birth,
-    renderer,
-    infections,
-    philm,
-    'inc_lm_det_',
-    parameters$lm_det_incidence_rendering_min_ages,
-    parameters$lm_det_incidence_rendering_max_ages,
-    timestep
-  )
-  lm_det_infections
+  bitset_at(infections, bernoulli_multi_p(philm))
 }
 
 #' @title Calculate clinical infections
diff --git a/R/parameters.R b/R/parameters.R
index 5af67575..8bcd8ef3 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -280,8 +280,6 @@
 #' * 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
-#' * 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
@@ -467,8 +465,6 @@ 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),
-      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/R/processes.R b/R/processes.R
index 80fa6912..5efadc88 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -372,9 +372,11 @@ create_lagged_eir <- function(variables, solvers, parameters) {
 create_hypnozoite_batch_decay_process <- function(hypnozoites, gammal){
   function(timestep){
     to_decay <- bernoulli_multi_p(p = rate_to_prob(hypnozoites$get_values() * gammal))
-    hypnozoites$queue_update(
-      hypnozoites$get_values(to_decay) - 1,
-      to_decay
-    )
+    if(length(to_decay) > 0){
+      hypnozoites$queue_update(
+        hypnozoites$get_values(to_decay) - 1,
+        to_decay
+      )
+    }
   }
 }
\ No newline at end of file
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 713fa83a..4be80329 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -127,7 +127,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     lagged_eir
   )
 
-  expect_equal(bitten$bitten_humans$to_vector(), c(2, 3))
+  expect_equal(bitten$to_vector(), c(2, 3))
 
   f <- parameters$blood_meal_rates[[1]]
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 1628798f..99f890ec 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -34,7 +34,7 @@ test_that('simulate_infection integrates different types of infection and schedu
   simulate_infection(
     variables,
     events,
-    list(bitten_humans = bitten),
+    bitten,
     age,
     parameters,
     timestep,
@@ -56,7 +56,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     infection_mock,
     1,
     variables,
-    list(bitten_humans = bitten),
+    bitten,
     parameters,
     renderer,
     timestep,
@@ -219,7 +219,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     depth = 4
   )
   
-  bitten_humans <- list(bitten_humans = individual::Bitset$new(4)$insert(c(1, 2, 3, 4)))
+  bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
@@ -705,7 +705,7 @@ test_that('prophylaxis is considered for medicated humans', {
     ib = individual::DoubleVariable$new(c(.2, .3, .5, .9))
   )
   
-  bitten <- list(bitten_humans = individual::Bitset$new(4)$insert(seq(4)))
+  bitten <- individual::Bitset$new(4)$insert(seq(4))
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
   
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index 2f173edf..9a5f1175 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -173,7 +173,7 @@ test_that('Infection considers pev efficacy', {
 
   infection_rates <- calculate_infections(
     variables = variables,
-    bitten_humans = list(bitten_humans = individual::Bitset$new(4)$insert(seq(4))),
+    bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
     parameters = parameters,
     renderer = mock_render(timestep),
     timestep = timestep,

From c3bed712cef6e5a896b6d42638a9f847ad52f145 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 22 Jul 2024 10:12:27 +0100
Subject: [PATCH 10/11] Went back to two biting process outputs: bitset of
 bitten humans and a vector that contains the number of bites each person
 receives.

Adjusted tests to reflect this.
---
 R/biting_process.R                          |  22 ++-
 R/human_infection.R                         | 185 ++++++++++----------
 tests/testthat/test-biting-integration.R    |  10 +-
 tests/testthat/test-infection-integration.R |   7 +
 4 files changed, 117 insertions(+), 107 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index 0149f3f6..30804ba2 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -33,7 +33,7 @@ create_biting_process <- function(
 ) {
   function(timestep){
     age <- get_age(variables$birth$get_values(), timestep)
-    bitten_humans <- simulate_bites(
+    bitten <- simulate_bites(
       renderer,
       solvers,
       models,
@@ -51,7 +51,8 @@ create_biting_process <- function(
     simulate_infection(
       variables,
       events,
-      bitten_humans,
+      bitten$bitten_humans,
+      bitten$n_bites_per_person,
       age,
       parameters,
       timestep,
@@ -76,6 +77,9 @@ simulate_bites <- function(
     mixing = 1,
     mixing_index = 1
 ) {
+
+  bitten_humans <- individual::Bitset$new(parameters$human_population)
+  n_bites_per_person <- numeric(0)
   
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
@@ -143,11 +147,9 @@ simulate_bites <- function(
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     if(parameters$parasite == "falciparum"){
-      bitten_humans <- individual::Bitset$new(parameters$human_population)
       # p.f model factors eir by psi
       expected_bites <- species_eir * mean(psi)
     } else if (parameters$parasite == "vivax"){
-      bitten_humans <- individual::IntegerVariable$new(initial_values = rep(0, parameters$human_population))
       # p.v model standardises biting rate het to eir
       expected_bites <- species_eir
     }
@@ -156,15 +158,11 @@ simulate_bites <- function(
       n_bites <- rpois(1, expected_bites)
       if (n_bites > 0) {
         bitten <- fast_weighted_sample(n_bites, lambda)
-        if(parameters$parasite == "falciparum"){
-          bitten_humans$insert(bitten)
-          renderer$render('n_bitten', bitten_humans$size(), timestep)
-        }
+        bitten_humans$insert(bitten)
+        renderer$render('n_bitten', bitten_humans$size(), timestep)
         if(parameters$parasite == "vivax"){
           # p.v must pass through the number of bites per person
-          bites_per_person <- tabulate(bitten, nbins = length(lambda))
-          bitten_humans <- individual::IntegerVariable$new(initial_values = bites_per_person)
-          renderer$render('n_bitten', bitten_humans$get_size_of(!0), timestep)
+          n_bites_per_person <- tabulate(bitten, nbins = length(lambda))
         }
       }
     }
@@ -216,7 +214,7 @@ simulate_bites <- function(
     }
   }
 
-  bitten_humans
+  list(bitten_humans = bitten_humans, n_bites_per_person = n_bites_per_person)
 }
 
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 0df45355..52ddfeb4 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -5,6 +5,7 @@
 #' @param variables a list of all of the model variables
 #' @param events a list of all of the model events
 #' @param bitten_humans a bitset of bitten humans
+#' @param n_bites_per_person vector of number of bites each person receives (p.v only)
 #' @param age of each human (timesteps)
 #' @param parameters of the model
 #' @param timestep current timestep
@@ -15,6 +16,7 @@ simulate_infection <- function(
     variables,
     events,
     bitten_humans,
+    n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -37,6 +39,7 @@ simulate_infection <- function(
   calculate_infections(
     variables,
     bitten_humans,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
@@ -48,6 +51,7 @@ simulate_infection <- function(
 #' @description Infection rates are stored in the infection outcome competing hazards object
 #' @param variables a list of all of the model variables
 #' @param bitten_humans bitset of bitten humans
+#' @param n_bites_per_person vector of number of bites each person receives (p.v only)
 #' @param parameters model parameters
 #' @param renderer model render object
 #' @param timestep current timestep
@@ -55,110 +59,109 @@ simulate_infection <- function(
 calculate_infections <- function(
     variables,
     bitten_humans,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
     infection_outcome
 ) {
+
+    if(parameters$parasite == "falciparum"){
+      source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
+    } else if (parameters$parasite == "vivax"){
+      ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
+      source_humans <- variables$hypnozoites$get_index_of(0)$not(T)$or(bitten_humans)
+    }
   
-  if(parameters$parasite == "falciparum"){
-    
-    if(bitten_humans$size() == 0){return(bitten_humans)}
+  if(source_humans$size() > 0){
+   
+    source_vector <- source_humans$to_vector()
     
-    source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
+    # calculate prophylaxis
+    prophylaxis <- rep(0, length(source_vector))
+    drug <- variables$drug$get_values(source_vector)
+    medicated <- (drug > 0)
+    if (any(medicated)) {
+      drug <- drug[medicated]
+      drug_time <- variables$drug_time$get_values(source_vector[medicated])
+      prophylaxis[medicated] <- weibull_survival(
+        timestep - drug_time,
+        parameters$drug_prophylaxis_shape[drug],
+        parameters$drug_prophylaxis_scale[drug]
+      )
+    }
     
-    ## p.f models blood immunity
-    b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+    # calculate vaccine efficacy
+    vaccine_efficacy <- rep(0, length(source_vector))
+    vaccine_times <- variables$last_eff_pev_timestep$get_values(source_vector)
+    pev_profile <- variables$pev_profile$get_values(source_vector)
+    # get vector of individuals who have received their 3rd dose
+    vaccinated <- vaccine_times > -1
+    pev_profile <- pev_profile[vaccinated]
+    if (length(vaccinated) > 0) {
+      antibodies <- calculate_pev_antibodies(
+        timestep - vaccine_times[vaccinated],
+        exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'cs')),
+        invlogit(sample_pev_param(pev_profile, parameters$pev_profiles, 'rho')),
+        exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'ds')),
+        exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'dl')),
+        parameters
+      )
+      vmax <- vnapply(parameters$pev_profiles, function(p) p$vmax)
+      beta <- vnapply(parameters$pev_profiles, function(p) p$beta)
+      alpha <- vnapply(parameters$pev_profiles, function(p) p$alpha)
+      vaccine_efficacy[vaccinated] <- calculate_pev_efficacy(
+        antibodies,
+        vmax[pev_profile],
+        beta[pev_profile],
+        alpha[pev_profile]
+      )
+    }
     
-  } else if (parameters$parasite == "vivax"){
-    ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
-    bitten_humans_index <- bitten_humans$get_index_of(!0)
-    source_humans <- variables$hypnozoites$get_index_of(0)$not(T)$or(bitten_humans_index)
-
-    if(source_humans$size() == 0){return(source_humans)}
+    infection_rates <- rep(0, length = parameters$human_population)
+    if(parameters$parasite == "falciparum"){
+      
+      ## p.f models blood immunity
+      b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+      
+      prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+      infection_rates[source_vector] <- prob_to_rate(prob)
+      
+    } else if (parameters$parasite == "vivax"){
+      
+      ## p.v does not model blood immunity but must take into account multiple bites per person
+      b <- 1 - (1 - parameters$b)^n_bites_per_person[bitten_humans$to_vector()]
+      
+      ## calculated rate of infection for all bitten or with hypnozoites
+      infection_rates[bitten_humans$to_vector()] <- prob_to_rate(b)
+      relapse_rates <- variables$hypnozoites$get_values() * parameters$f
+      infection_rates <- infection_rates + relapse_rates
+      relative_rates <- relapse_rates/infection_rates
+      relative_rates[is.nan(relative_rates)] <- 0
+      
+      infection_outcome$set_relative_rates(relative_rates)
+      
+      ## get relative rates to get probability bitten over relapse
+      prob <- rate_to_prob(infection_rates)
+      prob[source_vector] <- prob[source_vector] * (1 - prophylaxis) * (1 - vaccine_efficacy)
+      infection_rates <- prob_to_rate(prob)
+    }
     
-    # bitten_vector <- bitten_humans$to_vector()
-    ## p.v does not model blood immunity but must take into account multiple bites per person
-    b <- 1 - (1 - parameters$b)^bitten_humans$get_values(bitten_humans_index)
-  }
-  source_vector <- source_humans$to_vector()
-
-  # calculate prophylaxis
-  prophylaxis <- rep(0, length(source_vector))
-  drug <- variables$drug$get_values(source_vector)
-  medicated <- (drug > 0)
-  if (any(medicated)) {
-    drug <- drug[medicated]
-    drug_time <- variables$drug_time$get_values(source_vector[medicated])
-    prophylaxis[medicated] <- weibull_survival(
-      timestep - drug_time,
-      parameters$drug_prophylaxis_shape[drug],
-      parameters$drug_prophylaxis_scale[drug]
-    )
-  }
-
-  # calculate vaccine efficacy
-  vaccine_efficacy <- rep(0, length(source_vector))
-  vaccine_times <- variables$last_eff_pev_timestep$get_values(source_vector)
-  pev_profile <- variables$pev_profile$get_values(source_vector)
-  # get vector of individuals who have received their 3rd dose
-  vaccinated <- vaccine_times > -1
-  pev_profile <- pev_profile[vaccinated]
-  if (length(vaccinated) > 0) {
-    antibodies <- calculate_pev_antibodies(
-      timestep - vaccine_times[vaccinated],
-      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'cs')),
-      invlogit(sample_pev_param(pev_profile, parameters$pev_profiles, 'rho')),
-      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'ds')),
-      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'dl')),
-      parameters
-    )
-    vmax <- vnapply(parameters$pev_profiles, function(p) p$vmax)
-    beta <- vnapply(parameters$pev_profiles, function(p) p$beta)
-    alpha <- vnapply(parameters$pev_profiles, function(p) p$alpha)
-    vaccine_efficacy[vaccinated] <- calculate_pev_efficacy(
-      antibodies,
-      vmax[pev_profile],
-      beta[pev_profile],
-      alpha[pev_profile]
+    ## probability of incidence must be rendered at each timestep
+    incidence_probability_renderer(
+      variables$birth,
+      renderer,
+      source_humans,
+      prob,
+      "inc_",
+      parameters$incidence_min_ages,
+      parameters$incidence_max_ages,
+      timestep
     )
-  }
-  
-  infection_rates <- rep(0, length = parameters$human_population)
-  if(parameters$parasite == "falciparum"){
-    prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
-    infection_rates[source_vector] <- prob_to_rate(prob)
-    
-  } else if (parameters$parasite == "vivax"){
-    ## calculated rate of infection for all bitten or with hypnozoites
-    infection_rates[bitten_humans_index$to_vector()] <- prob_to_rate(b)
-    relapse_rates <- variables$hypnozoites$get_values() * parameters$f
-    infection_rates <- infection_rates + relapse_rates
-    relative_rates <- relapse_rates/infection_rates
-    relative_rates[is.nan(relative_rates)] <- 0
-    
-    infection_outcome$set_relative_rates(relative_rates)
     
-    ## get relative rates to get probability bitten over relapse
-    prob <- rate_to_prob(infection_rates)
-    prob[source_vector] <- prob[source_vector] * (1 - prophylaxis) * (1 - vaccine_efficacy)
-    infection_rates <- prob_to_rate(prob)
+    ## capture infection rates to resolve in competing hazards
+    infection_outcome$set_rates(infection_rates)
   }
-
-  ## probability of incidence must be rendered at each timestep
-  incidence_probability_renderer(
-    variables$birth,
-    renderer,
-    source_humans,
-    prob,
-    "inc_",
-    parameters$incidence_min_ages,
-    parameters$incidence_max_ages,
-    timestep
-  )
-  
-  ## capture infection rates to resolve in competing hazards
-  infection_outcome$set_rates(infection_rates)
 }
 
 #' @title Assigns infections to appropriate human states
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 4be80329..fe73d1db 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -35,8 +35,9 @@ test_that('biting_process integrates mosquito effects and human infection', {
     infection_outcome
   )
   
-  bitten <- individual::Bitset$new(parameters$human_population)
-  bites_mock <- mockery::mock(bitten)
+  bitten <- list(bitten_humans = individual::Bitset$new(parameters$human_population),
+                 n_bites_per_person = numeric(0))
+  bites_mock <- mockery::mock(bitten, cycle = T)
   infection_mock <- mockery::mock()
 
   mockery::stub(biting_process, 'simulate_bites', bites_mock)
@@ -65,7 +66,8 @@ test_that('biting_process integrates mosquito effects and human infection', {
     1,
     variables,
     events,
-    bitten,
+    bitten$bitten,
+    bitten$n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -127,7 +129,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     lagged_eir
   )
 
-  expect_equal(bitten$to_vector(), c(2, 3))
+  expect_equal(bitten$bitten_humans$to_vector(), c(2, 3))
 
   f <- parameters$blood_meal_rates[[1]]
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 99f890ec..072d3695 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -17,6 +17,7 @@ test_that('simulate_infection integrates different types of infection and schedu
   )
   
   bitten <- individual::Bitset$new(population)$insert(c(1, 3, 5, 7))
+  n_bites_per_person <- numeric(0)
   boost_immunity_mock <- mockery::mock()
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
   infection_mock <- mockery::mock(infected)
@@ -35,6 +36,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     variables,
     events,
     bitten,
+    n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -57,6 +59,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     1,
     variables,
     bitten,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
@@ -220,6 +223,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   )
   
   bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  n_bites_per_person <- numeric(0)
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
@@ -231,6 +235,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   infections <- calculate_infections(
     variables,
     bitten_humans, 
+    n_bites_per_person,
     parameters,
     mock_render(timestep),
     timestep,
@@ -706,6 +711,7 @@ test_that('prophylaxis is considered for medicated humans', {
   )
   
   bitten <- individual::Bitset$new(4)$insert(seq(4))
+  n_bites_per_person <- numeric(0)
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
   
@@ -719,6 +725,7 @@ test_that('prophylaxis is considered for medicated humans', {
   infection_rates <- calculate_infections(
     variables,
     bitten,
+    n_bites_per_person,
     parameters,
     mock_render(timestep),
     timestep,

From 8b03d441425210d8c7deefd7aec6749ad826b89b Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Thu, 25 Jul 2024 10:18:56 +0100
Subject: [PATCH 11/11] Relapses can sometimes be 0, so must be included in the
 initial rendering to avoid NAs in the render object.

---
 R/render.R | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/R/render.R b/R/render.R
index 0d066b46..e0b03087 100644
--- a/R/render.R
+++ b/R/render.R
@@ -242,6 +242,11 @@ populate_incidence_rendering_columns <- function(renderer, parameters){
     renderer$set_default('n_early_treatment_failure', 0)
     renderer$set_default('n_slow_parasite_clearance', 0)
   }
+
+  # relapses only render for the vivax model
+  if(parameters$parasite == "vivax"){
+    renderer$set_default('n_relapses', 0)
+  }
   
   if(length(parameters$incidence_rendering_min_ages)>0){
     render_initial_incidence(renderer,