From 74b51aa2973253f3e6f98f51908c0bf64fcb6ddd Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Tue, 17 Sep 2024 17:15:53 +0100
Subject: [PATCH 1/3] Hynozoite batches are now modelled, including formation
 via bite infections, decay of hypnozoite batches, new infections via relapse
 and the reseting of batch number to 0 for new births.

A hypnozoite batch variable is created.

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, as currently in the P. falciparum model) 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 relapse 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 the 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 and expanded the p.v vignette.
---
 R/biting_process.R                          |  30 ++-
 R/competing_hazards.R                       |  11 +-
 R/human_infection.R                         | 265 ++++++++++++++------
 R/mortality_processes.R                     |   1 +
 R/parameters.R                              |   4 +-
 R/processes.R                               |  45 +++-
 R/render.R                                  |  40 +++
 R/variables.R                               |   9 +-
 man/get_parameters.Rd                       |   4 +-
 tests/testthat/test-biting-integration.R    |  10 +-
 tests/testthat/test-infection-integration.R |  11 +-
 tests/testthat/test-pev.R                   |   1 +
 tests/testthat/test-vivax.R                 |  60 +++++
 vignettes/Plasmodium_vivax.Rmd              |   6 +
 14 files changed, 391 insertions(+), 106 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index df494604..17b75366 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -34,7 +34,7 @@ create_biting_process <- function(
   function(timestep) {
     # Calculate combined EIR
     age <- get_age(variables$birth$get_values(), timestep)
-    bitten_humans <- simulate_bites(
+    bitten <- simulate_bites(
       renderer,
       solvers,
       models,
@@ -52,7 +52,8 @@ create_biting_process <- function(
     simulate_infection(
       variables,
       events,
-      bitten_humans,
+      bitten$bitten_humans,
+      bitten$n_bites_per_person,
       age,
       parameters,
       timestep,
@@ -78,6 +79,7 @@ simulate_bites <- function(
   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) {
@@ -146,13 +148,24 @@ 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)
+        renderer$render('n_bitten', bitten_humans$size(), timestep)
+        if(parameters$parasite == "vivax"){
+          # p.v must pass through the number of bites per person
+          n_bites_per_person <- tabulate(bitten, nbins = length(lambda))
+        }
       }
     }
 
@@ -202,9 +215,8 @@ simulate_bites <- function(
       )
     }
   }
-  
-  renderer$render('n_bitten', bitten_humans$size(), timestep)
-  bitten_humans
+
+  list(bitten_humans = bitten_humans, n_bites_per_person = n_bites_per_person)
 }
 
 
diff --git a/R/competing_hazards.R b/R/competing_hazards.R
index 13c522ad..5a44e3e6 100644
--- a/R/competing_hazards.R
+++ b/R/competing_hazards.R
@@ -16,6 +16,7 @@ CompetingOutcome <- R6::R6Class(
 
       self$target <- individual::Bitset$new(size)
       self$rates <- NULL
+      self$relative_rates <- NULL
     },
     set_rates = function(target, rates){
       stopifnot(target$size() == length(rates))
@@ -23,15 +24,23 @@ CompetingOutcome <- R6::R6Class(
       self$target$copy_from(target)
       self$rates <- rates
     },
+    set_relative_rates = function(target, relative_rates){
+      stopifnot(target$size() == length(relative_rates))
+      
+      self$target$copy_from(target)
+      self$relative_rates <- relative_rates
+    },
     execute = function(t, target){
       private$targeted_process(t, target)
     },
     reset = function() {
       self$target$clear()
       self$rates <- NULL
+      self$relative_rates <- NULL
     },
     target = NULL,
-    rates = NULL
+    rates = NULL,
+    relative_rates = NULL
   )
 )
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 36c07897..1f5b8705 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,84 +59,116 @@ simulate_infection <- function(
 calculate_infections <- function(
     variables,
     bitten_humans,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
     infection_outcome
 ) {
-  source_humans <- variables$state$get_index_of(
-    c('S', 'A', 'U'))$and(bitten_humans)
 
-  if(parameters$parasite == "falciparum"){
-    ## p.f models blood immunity
-    b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+    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
+      hypnozoites_humans <- variables$hypnozoites$get_index_of(set = 0)$not(T)
+      source_humans <- bitten_humans$copy()$or(hypnozoites_humans)
+    }
+  
+  if(source_humans$size() > 0){
+   
+    source_vector <- source_humans$to_vector()
     
-  } else if (parameters$parasite == "vivax"){
-    ## p.v does not model blood immunity
-    b <- parameters$b
-  }
-
-  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]
+    # 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]
+      )
+    }
+    
+    # infection_rates <- rep(0, length = parameters$human_population)
+    if(parameters$parasite == "falciparum"){
+      
+      ## p.f models blood immunity
+      prob <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+      
+    } 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()]
+      
+      ## get infection rates for all bitten or with hypnozoites
+      infection_rates <- rep(0, length = source_humans$size())
+      infection_rates[bitset_index(source_humans, bitten_humans)] <- prob_to_rate(b)
+      
+      # Add relapse rates for individuals with hypnozoites
+      if(hypnozoites_humans$size()>0){
+        relapse_rates <- variables$hypnozoites$get_values(hypnozoites_humans) * parameters$f
+        hyp_source_index <- bitset_index(source_humans, hypnozoites_humans)
+        infection_rates[hyp_source_index] <- infection_rates[hyp_source_index] + relapse_rates
+        # Get and store relative rates for bite/relapse competing hazards resolution
+        relative_rates <- relapse_rates/infection_rates[hyp_source_index]
+        infection_outcome$set_relative_rates(
+          hypnozoites_humans, 
+          relative_rates)
+      }
+      
+      prob <- rate_to_prob(infection_rates)
+    }
+    
+    prob <- prob * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    infection_rates <- prob_to_rate(prob)
+    
+    ## 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(
+      source_humans, 
+      infection_rates)
   }
-
-  prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
-
-  ## 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(
-    source_humans,
-    prob_to_rate(prob))
 }
 
 #' @title Assigns infections to appropriate human states
@@ -145,6 +181,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,7 +189,8 @@ infection_outcome_process <- function(
     variables,
     renderer,
     parameters,
-    prob){
+    prob,
+    relative_rates = NULL){
   
   if (infected_humans$size() > 0) {
     
@@ -215,6 +253,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,
@@ -227,9 +274,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
@@ -270,31 +315,91 @@ 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
+){
+  
+  if(variables$hypnozoites$get_index_of(0)$not(T)$and(infected_humans)$size()>0){
+    
+    hypnozoite_humans <- variables$hypnozoites$get_index_of(0)$not(T)
+    potential_relapse_index <- bitset_index(hypnozoite_humans, infected_humans)
+    relapse_infections <- bitset_at(
+      hypnozoite_humans$and(infected_humans),
+      bernoulli_multi_p(relative_rates[potential_relapse_index])
+    )
+    
+    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))
+    
+  } else {
+    bite_infections <- infected_humans
+  }
+  
+  ## all bitten humans with an infectious bite (incorporating prophylaxis) get a new batch of hypnozoites
+  if(bite_infections$size()>0){
+    
+    # make sure batches are capped
+    current_batches <- variables$hypnozoites$get_values(bite_infections)
+    new_batch_number <- pmin(current_batches + 1, parameters$kmax)
+    
+    variables$hypnozoites$queue_update(
+      new_batch_number,
+      bite_infections
+    )
+  }
+}
+
 #' @title Calculate light microscopy detectable infections (p.v only)
 #' @description
 #' Sample light microscopy detectable infections from all infections
 #' @param variables a list of all of the model variables
 #' @param infections bitset of infected humans
 #' @param parameters model parameters
-#' @param renderer model render
-#' @param timestep current timestep
 #' @noRd
 calculate_lm_det_infections <- function(
     variables,
     infections,
-    parameters,
-    renderer,
-    timestep
+    parameters
 ) {
   
   iaa <- variables$iaa$get_values(infections)
   iam <- variables$iam$get_values(infections)
-
+  
   philm <- anti_parasite_immunity(
     min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
     k = parameters$klm, iaa = iaa, iam = iam)
 
-  lm_det_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  bitset_at(infections, bernoulli_multi_p(philm))
 }
 
 #' @title Calculate clinical infections
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index b41aca32..9b530c56 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/parameters.R b/R/parameters.R
index 79ef7629..f113d908 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -189,9 +189,9 @@
 #' * delay_gam - lag from parasites to infectious gametocytes; default = 12.5
 #' * dem - extrinsic incubation period in mosquito population model; default = 10
 #' 
-#' hypnozoite parameters (p.v only):
+#' hypnozoite batch parameters (p.v only):
 #'
-#' * f - hypnozoite relapse rate; default = 0.02439024
+#' * f - hypnozoite batch relapse rate; default = 0.02439024
 #' * gammal - hypnozoite batch clearance rate; default = 0.002610966
 #' * init_hyp - initial hypnozoite batch number; default = 0
 #' * kmax - maximum number of hypnozoite batches for use in the equilibrium solution; default = 10
diff --git a/R/processes.R b/R/processes.R
index b67dc559..f8b3b1be 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -47,7 +47,7 @@ create_processes <- function(
     immunity_process = create_exponential_decay_process(variables$ica,
                                                         parameters$rc)
   )
-  
+
   if(parameters$parasite == "falciparum"){
     processes <- c(
       processes,
@@ -71,7 +71,9 @@ create_processes <- function(
       immunity_process = create_exponential_decay_process(variables$iam, 
                                                           parameters$rm),
       immunity_process = create_exponential_decay_process(variables$iaa, 
-                                                          parameters$ra)
+                                                          parameters$ra),
+      hypnozoite_process = create_hypnozoite_batch_decay_process(variables$hypnozoites,
+                                                                 parameters$gammal)
     )
   }
 
@@ -96,7 +98,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
   )
@@ -202,7 +205,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,
+      hypnozoite_prevalence_renderer = create_n_with_hypnozoites_renderer_process(
+        renderer,
+        variables$hypnozoites,
+        parameters
+      ),
+      hypnozoite_prevalence_age_group_renderer = create_n_with_hypnozoites_age_renderer_process(
+        variables$hypnozoites,
+        variables$birth,
+        parameters,
+        renderer
+      )
+    )
   }
 
   processes <- c(
@@ -379,3 +398,21 @@ 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))
+    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/R/render.R b/R/render.R
index 34f4bd99..1357decf 100644
--- a/R/render.R
+++ b/R/render.R
@@ -266,6 +266,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,
@@ -303,3 +308,38 @@ populate_metapopulation_incidence_rendering_columns <- function(renderer, parame
     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
+      )
+    }
+  }
+}
diff --git a/R/variables.R b/R/variables.R
index 024b695d..b6bd1f80 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,7 +99,10 @@ 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(
     initial_immunity(
@@ -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/man/get_parameters.Rd b/man/get_parameters.Rd
index 76e2db6a..76ce99a8 100644
--- a/man/get_parameters.Rd
+++ b/man/get_parameters.Rd
@@ -209,9 +209,9 @@ parasite incubation periods:
 \item dem - extrinsic incubation period in mosquito population model; default = 10
 }
 
-hypnozoite parameters (p.v only):
+hypnozoite batch parameters (p.v only):
 \itemize{
-\item f - hypnozoite relapse rate; default = 0.02439024
+\item f - hypnozoite batch relapse rate; default = 0.02439024
 \item gammal - hypnozoite batch clearance rate; default = 0.002610966
 \item init_hyp - initial hypnozoite batch number; default = 0
 \item kmax - maximum number of hypnozoite batches for use in the equilibrium solution; default = 10
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index fbfd1fd4..9ab38a76 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -33,8 +33,9 @@ test_that('biting_process integrates mosquito effects and human infection', {
     infection_outcome=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)
@@ -63,7 +64,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,
@@ -125,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 ad83c35f..c29420c9 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,
@@ -64,7 +67,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
@@ -156,6 +158,8 @@ test_that('simulate_infection integrates different types of infection and schedu
     renderer
   )
   
+  mockery::mock_args(schedule_mock)
+  
   mockery::expect_args(
     schedule_mock,
     1,
@@ -220,6 +224,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 +236,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,
@@ -699,10 +705,10 @@ 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)
   
-  
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
       infection_outcome_process(timestep, target, variables, renderer, parameters)
@@ -713,6 +719,7 @@ test_that('prophylaxis is considered for medicated humans', {
   infection_rates <- calculate_infections(
     variables,
     bitten,
+    n_bites_per_person,
     parameters,
     mock_render(timestep),
     timestep,
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index f3c3d793..ce41ac20 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -174,6 +174,7 @@ test_that('Infection considers pev efficacy', {
   infection_rates <- calculate_infections(
     variables = variables,
     bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
+    n_bites_per_person = numeric(0),
     parameters = parameters,
     renderer = mock_render(timestep),
     timestep = timestep,
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 8fc77956..4f762728 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)),
+    progression_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(1, 3), 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, 3))
+  )
+
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_infections',
+    4,
+    timestep
+  )
+
+  mockery::expect_args(
+    renderer$render_mock(),
+    2,
+    'n_relapses',
+    2,
+    timestep
+  )
+  
+})
diff --git a/vignettes/Plasmodium_vivax.Rmd b/vignettes/Plasmodium_vivax.Rmd
index 87c56671..67c0227e 100644
--- a/vignettes/Plasmodium_vivax.Rmd
+++ b/vignettes/Plasmodium_vivax.Rmd
@@ -55,6 +55,8 @@ The *P. vivax* model follows a similar structure to the *P. falciparum* model, a
 
 ### New Infections
 
+The rate of infection for an individual who has been bitten increases with the number of bites they have received in the *P. vivax* model. In contrast, the rate of infection in the *P. falciparum* model is independent of the number of bites an individual has received.
+
 Newly infected individuals in the  *P. falciparum* model can move into either the clinically diseased or asymptomatic infection compartment. In addition to these compartments, the *P. vivax* model allows infection to the PCR-detectable compartment (U), where the assignment of light-miscroscopy detectable infections and PCR-detectable infections is mediated by anti-parasite immunity.
 
 ### Immunity
@@ -69,6 +71,10 @@ Anti-parasite immunity has effects in two places. The first is in the separation
 
 While the *P. falciparum* model calculates the onward infectivity of asymptomatic infections (`ca`) using the age and detectability immunity of each individual, the *P. vivax* model uses a constant onwards infectivity for LM-detectable infections (`ca = 0.1`).
 
+### Hypnozoites
+
+The *P. vivax* model tracks the number of hypnozoite batches of each individual which contribute to the rate of new infections. Acquisition of new batches comes from bite infections and are capped (where `kmax = 10` by default). All individuals can acquire new hypnozoite batches via bites, even when these do not result in new blood-stage infection (as in the clinically diseased and treated). Hypnozoite batches are lost probabilistically where an individual's number of batches determines the loss rate of a single batch.
+
 ### Key Model References
 
 White, M. T., Walker, P., Karl, S., Hetzel, M. W., Freeman, T., Waltmann, A., Laman, M., Robinson, L. J., Ghani, A., & Mueller, I. (2018). Mathematical modelling of the impact of expanding levels of malaria control interventions on Plasmodium vivax. Nature Communications, 9(1), 1–10. <https://doi.org/10.1038/s41467-018-05860-8>

From 9a78c5b2760a964b2fff7351af15981f0ad6d783 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Mon, 14 Oct 2024 15:42:48 +0100
Subject: [PATCH 2/3] Rewritten competing hazards resoltion so that it does not
 call itself within the function. This also removed the "set_relative_rates"
 part of the competing hazards, and now we have args, such that the competing
 hazards object can still store the relative rates it needs, but retains
 flexibility.

Calculate_infections has been rewritten as four functions: calculate_falciparum_infections, calculate_vivax_infections, treatment_prophylaxis_efficacy and pev_efficacy.

Infection_outcome_process has been split into falciparum and vivax _infection_outcome_process.

Tests have been adjusted to reflect these changes.
---
 R/competing_hazards.R                       |  18 +-
 R/human_infection.R                         | 518 +++++++++++++-------
 R/processes.R                               |  32 +-
 tests/testthat/test-competing-hazards.R     |  36 +-
 tests/testthat/test-infection-integration.R |  70 ++-
 tests/testthat/test-pev.R                   |  24 +-
 tests/testthat/test-vivax.R                 |   7 +-
 7 files changed, 428 insertions(+), 277 deletions(-)

diff --git a/R/competing_hazards.R b/R/competing_hazards.R
index 5a44e3e6..0ef30e0d 100644
--- a/R/competing_hazards.R
+++ b/R/competing_hazards.R
@@ -16,31 +16,25 @@ CompetingOutcome <- R6::R6Class(
 
       self$target <- individual::Bitset$new(size)
       self$rates <- NULL
-      self$relative_rates <- NULL
+      self$args <- NULL
     },
-    set_rates = function(target, rates){
+    set_rates = function(target, rates, ...){
       stopifnot(target$size() == length(rates))
-
       self$target$copy_from(target)
       self$rates <- rates
-    },
-    set_relative_rates = function(target, relative_rates){
-      stopifnot(target$size() == length(relative_rates))
-      
-      self$target$copy_from(target)
-      self$relative_rates <- relative_rates
+      self$args <- list(...)
     },
     execute = function(t, target){
-      private$targeted_process(t, target)
+      do.call(private$targeted_process, list(t, target, self$args))
     },
     reset = function() {
       self$target$clear()
       self$rates <- NULL
-      self$relative_rates <- NULL
+      self$args <- NULL
     },
     target = NULL,
     rates = NULL,
-    relative_rates = NULL
+    args = NULL
   )
 )
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 1f5b8705..c779efc7 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -23,8 +23,10 @@ simulate_infection <- function(
     renderer,
     infection_outcome
 ) {
+  
   if (bitten_humans$size() > 0) {
     if(parameters$parasite == "falciparum"){
+      
       boost_immunity(
         variables$ib,
         bitten_humans,
@@ -32,22 +34,97 @@ simulate_infection <- function(
         timestep,
         parameters$ub
       )
+      
+      # Calculate Infected
+      calculate_falciparum_infections(
+        variables,
+        bitten_humans,
+        parameters,
+        renderer,
+        timestep,
+        infection_outcome
+      )
+      
+    } else if(parameters$parasite == "vivax"){
+      
+      # Calculate Infected
+      calculate_vivax_infections(
+        variables,
+        bitten_humans,
+        n_bites_per_person,
+        parameters,
+        renderer,
+        timestep,
+        infection_outcome
+      )
     }
   }
+  
+}
 
-  # Calculate Infected
-  calculate_infections(
+#' @title Calculate overall falciparum infections for bitten humans
+#' @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 parameters model parameters
+#' @param renderer model render object
+#' @param timestep current timestep
+#' @noRd
+calculate_falciparum_infections <- function(
     variables,
     bitten_humans,
-    n_bites_per_person,
     parameters,
     renderer,
     timestep,
     infection_outcome
-  )
+) {
+  
+  source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
+  
+  if(source_humans$size() > 0){
+    
+    source_vector <- source_humans$to_vector()
+    
+    # calculate prophylaxis
+    prophylaxis <- treatment_prophylaxis_efficacy(
+      source_vector,
+      variables,
+      parameters,
+      timestep
+    )
+    
+    # calculate vaccine efficacy
+    vaccine_efficacy <- pev_efficacy(
+      source_vector,
+      variables,
+      parameters,
+      timestep)
+    
+    ## p.f models blood immunity
+    prob <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+    prob <- prob * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    infection_rates <- prob_to_rate(prob)
+    
+    ## 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(
+      source_humans, 
+      infection_rates)
+  }
 }
 
-#' @title Calculate overall infections for bitten humans
+#' @title Calculate overall vivax infections for bitten humans
 #' @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
@@ -56,7 +133,7 @@ simulate_infection <- function(
 #' @param renderer model render object
 #' @param timestep current timestep
 #' @noRd
-calculate_infections <- function(
+calculate_vivax_infections <- function(
     variables,
     bitten_humans,
     n_bites_per_person,
@@ -65,91 +142,51 @@ calculate_infections <- function(
     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
-      hypnozoites_humans <- variables$hypnozoites$get_index_of(set = 0)$not(T)
-      source_humans <- bitten_humans$copy()$or(hypnozoites_humans)
-    }
+  
+  ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
+  hypnozoites_humans <- variables$hypnozoites$get_index_of(set = 0)$not(T)
+  source_humans <- bitten_humans$copy()$or(hypnozoites_humans)
   
   if(source_humans$size() > 0){
-   
+    
     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]
-      )
-    }
+    prophylaxis <- treatment_prophylaxis_efficacy(
+      source_vector,
+      variables,
+      parameters,
+      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]
-      )
-    }
+    vaccine_efficacy <- pev_efficacy(
+      source_vector,
+      variables,
+      parameters,
+      timestep)
     
-    # infection_rates <- rep(0, length = parameters$human_population)
-    if(parameters$parasite == "falciparum"){
-      
-      ## p.f models blood immunity
-      prob <- blood_immunity(variables$ib$get_values(source_humans), parameters)
-      
-    } 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()]
-      
-      ## get infection rates for all bitten or with hypnozoites
-      infection_rates <- rep(0, length = source_humans$size())
-      infection_rates[bitset_index(source_humans, bitten_humans)] <- prob_to_rate(b)
-      
-      # Add relapse rates for individuals with hypnozoites
-      if(hypnozoites_humans$size()>0){
-        relapse_rates <- variables$hypnozoites$get_values(hypnozoites_humans) * parameters$f
-        hyp_source_index <- bitset_index(source_humans, hypnozoites_humans)
-        infection_rates[hyp_source_index] <- infection_rates[hyp_source_index] + relapse_rates
-        # Get and store relative rates for bite/relapse competing hazards resolution
-        relative_rates <- relapse_rates/infection_rates[hyp_source_index]
-        infection_outcome$set_relative_rates(
-          hypnozoites_humans, 
-          relative_rates)
-      }
-      
-      prob <- rate_to_prob(infection_rates)
+    ## 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()]
+    
+    ## get infection rates for all bitten or with hypnozoites
+    infection_rates <- rep(0, length = source_humans$size())
+    infection_rates[bitset_index(source_humans, bitten_humans)] <- prob_to_rate(b)
+    
+    # Add relapse rates for individuals with hypnozoites
+    relative_rates <- NULL
+    if(hypnozoites_humans$size()>0){
+      relapse_rates <- variables$hypnozoites$get_values(hypnozoites_humans) * parameters$f
+      hyp_source_index <- bitset_index(source_humans, hypnozoites_humans)
+      infection_rates[hyp_source_index] <- infection_rates[hyp_source_index] + relapse_rates
+      # Get and store relative rates for bite/relapse competing hazards resolution
+      relative_rates <- relapse_rates/infection_rates[hyp_source_index]
+      # infection_outcome$set_relative_rates(
+      # hypnozoites_humans, 
+      # relative_rates)
     }
     
-    prob <- prob * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    prob <- rate_to_prob(infection_rates) * (1 - prophylaxis) * (1 - vaccine_efficacy)
     infection_rates <- prob_to_rate(prob)
     
     ## probability of incidence must be rendered at each timestep
@@ -166,12 +203,85 @@ calculate_infections <- function(
     
     ## capture infection rates to resolve in competing hazards
     infection_outcome$set_rates(
-      source_humans, 
-      infection_rates)
+      source_humans,
+      infection_rates,
+      relative_rates = relative_rates
+    )
+  }
+}
+
+#' @title Calculate protection from infection due to drug prophylaxis
+#' @description This function calculates the probability that drug prophylaxis will 
+#' protect each individual from infection
+#' @param source_vector a vector of individuals with prospective new infections
+#' @param variables a list of all of the model variables
+#' @param parameters model parameters
+#' @param timestep current timestep
+#' @noRd
+treatment_prophylaxis_efficacy <- function(
+    source_vector,
+    variables,
+    parameters,
+    timestep
+){
+  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]
+    )
+  }
+  prophylaxis
+}
+
+#' @title Calculate protection from infection due to vaccination
+#' @description This function calculates the probability that vaccination will 
+#' protect each individual from infection
+#' @param source_vector a vector of individuals with prospective new infections
+#' @param variables a list of all of the model variables
+#' @param parameters model parameters
+#' @param timestep current timestep
+#' @noRd
+pev_efficacy <- function(
+    source_vector,
+    variables,
+    parameters,
+    timestep
+){
+  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]
+    )
   }
+  vaccine_efficacy
 }
 
-#' @title Assigns infections to appropriate human states
+#' @title Assigns falciparum infections to appropriate human states
 #' @description
 #' Updates human states and variables to represent asymptomatic/clinical/severe
 #' and treated malaria; and resulting boosts in immunity
@@ -181,16 +291,101 @@ calculate_infections <- function(
 #' @param renderer model render object
 #' @param parameters model parameters
 #' @param prob vector of population probabilities of infection
+#' @noRd
+falciparum_infection_outcome_process <- function(
+    timestep,
+    infected_humans,
+    variables,
+    renderer,
+    parameters){
+  
+  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,
+      infected_humans,
+      variables$last_boosted_ica,
+      timestep,
+      parameters$uc
+    )
+    
+    boost_immunity(
+      variables$id,
+      infected_humans,
+      variables$last_boosted_id,
+      timestep,
+      parameters$ud
+    )
+    
+    clinical <- calculate_clinical_infections(
+      variables,
+      infected_humans,
+      parameters,
+      renderer,
+      timestep
+    )
+    
+    treated <- calculate_treated(
+      variables,
+      clinical,
+      parameters,
+      timestep,
+      renderer
+    )
+    
+    update_severe_disease(
+      timestep,
+      infected_humans,
+      variables,
+      parameters,
+      renderer
+    )
+    
+    ## The treated and infected_humans bitsets are re-written so be cautious!
+    to_D <- treated$not(FALSE)$and(clinical)
+    to_A <- infected_humans$and(clinical$not(FALSE))
+    to_U <- NULL
+    
+    schedule_infections(
+      parameters,
+      variables,
+      timestep,
+      to_D,
+      to_A,
+      to_U
+    )
+  }
+}
+
+#' @title Assigns vivax infections to appropriate human states
+#' @description
+#' Updates human states and variables to represent asymptomatic/clinical/severe
+#' and treated malaria; and resulting boosts in immunity
+#' @param timestep current timestep
+#' @param infected_humans bitset of infected humans
+#' @param variables a list of all of the model variables
+#' @param renderer model render object
+#' @param parameters model parameters
 #' @param relative_rates relative rates of hypnozoite relapse relative to total infection rate
 #' @noRd
-infection_outcome_process <- function(
+vivax_infection_outcome_process <- function(
     timestep,
     infected_humans,
     variables,
     renderer,
     parameters,
-    prob,
-    relative_rates = NULL){
+    relative_rates){
   
   if (infected_humans$size() > 0) {
     
@@ -205,6 +400,14 @@ infection_outcome_process <- function(
       timestep
     )
     
+    boost_immunity(
+      variables$iaa,
+      infected_humans,
+      variables$last_boosted_iaa,
+      timestep,
+      parameters$ua
+    )
+    
     boost_immunity(
       variables$ica,
       infected_humans,
@@ -212,97 +415,48 @@ infection_outcome_process <- function(
       timestep,
       parameters$uc
     )
-
-    if(parameters$parasite == "falciparum"){
-      boost_immunity(
-        variables$id,
-        infected_humans,
-        variables$last_boosted_id,
-        timestep,
-        parameters$ud
-      )
-      
-      clinical <- calculate_clinical_infections(
-        variables,
-        infected_humans,
-        parameters,
-        renderer,
-        timestep
-      )
-      
-      treated <- calculate_treated(
-        variables,
-        clinical,
-        parameters,
-        timestep,
-        renderer
-      )
-      
-      update_severe_disease(
-        timestep,
-        infected_humans,
-        variables,
-        parameters,
-        renderer
-      )
-      
-      ## The treated and infected_humans bitsets are re-written so be cautious!
-      to_D <- treated$not(FALSE)$and(clinical)
-      to_A <- infected_humans$and(clinical$not(FALSE))
-      to_U <- NULL
-      
-    } else if (parameters$parasite == "vivax"){
-      
-      relapse_bite_infection_hazard_resolution(
-        infected_humans,
-        relative_rates,
-        variables,
-        parameters,
-        renderer,
-        timestep
-      )
-
-      boost_immunity(
-        variables$iaa,
-        infected_humans,
-        variables$last_boosted_iaa,
-        timestep,
-        parameters$ua
-      )
-      
-      ## Only S and U infections are considered in generating lm-det infections
-      lm_detectable <- calculate_lm_det_infections(
-        variables,
-        variables$state$get_index_of(c("S","U"))$and(infected_humans),
-        parameters
-      )
-      
-      # Lm-detectable level infected S and U, and all A infections may receive clinical infections
-      # There is a different calculation to generate clinical infections, based on current infection level
-      # LM infections must only pass through the clinical calculation, therefore all "A" infections are included
-      # "S" and "U" infections must pass through the lm-detectable calculation prior to and in addition to the clinical
-      # calculation. We therefore consider all "A" infections and only the "S" and "U" infections that are now lm-detectable.
-      clinical <- calculate_clinical_infections(
-        variables,
-        variables$state$get_index_of("A")$and(infected_humans)$or(lm_detectable),
-        parameters,
-        renderer,
-        timestep
-      )
-      
-      treated <- calculate_treated(
-        variables,
-        clinical,
-        parameters,
-        timestep,
-        renderer
-      )
-      
-      ## The infected_humans,lm_detectable and clinical bitsets are re-written so be cautious!
-      to_U <- infected_humans$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
-      to_A <- lm_detectable$and(clinical$not(F))
-      to_D <- clinical$and(treated$not(F))
-    }
+    
+    relapse_bite_infection_hazard_resolution(
+      infected_humans,
+      relative_rates,
+      variables,
+      parameters,
+      renderer,
+      timestep
+    )
+    
+    ## Only S and U infections are considered in generating lm-det infections
+    lm_detectable <- calculate_lm_det_infections(
+      variables,
+      variables$state$get_index_of(c("S","U"))$and(infected_humans),
+      parameters
+    )
+    
+    # Lm-detectable level infected S and U, and all A infections may receive clinical infections
+    # There is a different calculation to generate clinical infections, based on current infection level
+    # LM infections must only pass through the clinical calculation, therefore all "A" infections are included
+    # "S" and "U" infections must pass through the lm-detectable calculation prior to and in addition to the clinical
+    # calculation. We therefore consider all "A" infections and only the "S" and "U" infections that are now lm-detectable.
+    clinical <- calculate_clinical_infections(
+      variables,
+      variables$state$get_index_of("A")$and(infected_humans)$or(lm_detectable),
+      parameters,
+      renderer,
+      timestep
+    )
+    
+    treated <- calculate_treated(
+      variables,
+      clinical,
+      parameters,
+      timestep,
+      renderer
+    )
+    
+    ## The infected_humans,lm_detectable and clinical bitsets are re-written so be cautious!
+    to_U <- infected_humans$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
+    to_A <- lm_detectable$and(clinical$not(F))
+    to_D <- clinical$and(treated$not(F))
     
     schedule_infections(
       parameters,
@@ -398,7 +552,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)
-
+  
   bitset_at(infections, bernoulli_multi_p(philm))
 }
 
@@ -679,7 +833,7 @@ schedule_infections <- function(
     to_A,
     to_U
 ) {
-
+  
   if(to_D$size() > 0) {
     update_infection(
       variables$state,
@@ -691,7 +845,7 @@ schedule_infections <- function(
       to_D
     )
   }
-
+  
   if(to_A$size() > 0) {
     if(parameters$parasite == "falciparum"){
       # p.f has immunity-determined asymptomatic infectivity
@@ -790,7 +944,7 @@ severe_immunity <- function(age, acquired_immunity, maternal_immunity, parameter
   parameters$theta0 * (parameters$theta1 + (1 - parameters$theta1) / (
     1 + fv * (
       (acquired_immunity + maternal_immunity) / parameters$iv0) ** parameters$kv
-    )
+  )
   )
 }
 
diff --git a/R/processes.R b/R/processes.R
index f8b3b1be..6ac4f290 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -94,18 +94,30 @@ create_processes <- function(
   # Competing Hazard Outcomes (infections and disease progression)
   # =====================================================
   
-  infection_outcome <- CompetingOutcome$new(
-    targeted_process = function(timestep, target){
-      infection_outcome_process(timestep, target, 
-                                variables, renderer, parameters, 
-                                prob = rate_to_prob(infection_outcome$rates),
-                                relative_rates = infection_outcome$relative_rates)
-    },
-    size = parameters$human_population
-  )
+  if(parameters$parasite == "falciparum"){
+    infection_outcome <- CompetingOutcome$new(
+      targeted_process = function(timestep, target, args){
+        falciparum_infection_outcome_process(timestep, target, 
+                                  variables, renderer, parameters
+        )
+      },
+      size = parameters$human_population
+    )
+    
+  } else if (parameters$parasite == "vivax"){
+    infection_outcome <- CompetingOutcome$new(
+      targeted_process = function(timestep, target, args){
+        vivax_infection_outcome_process(timestep, target, 
+                                  variables, renderer, parameters, 
+                                  args$relative_rates
+        )
+      },
+      size = parameters$human_population
+    )
+  }
   
   progression_outcome <- CompetingOutcome$new(
-    targeted_process = function(timestep, target){
+    targeted_process = function(timestep, target, ...){
       progression_outcome_process(timestep, target, variables, parameters, renderer)
     },
     size = parameters$human_population
diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R
index f7ada3fc..fc817499 100644
--- a/tests/testthat/test-competing-hazards.R
+++ b/tests/testthat/test-competing-hazards.R
@@ -28,11 +28,12 @@ test_that("hazard resolves two disjoint outcomes", {
   outcome_2$set_rates(population, c(0, 10, 0, 10))
   
   hazard$resolve(0)
-   
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(1, 3)))
+                       individual::Bitset$new(size)$insert(c(1, 3)),
+                       list())
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(2, 4)))
+                       individual::Bitset$new(size)$insert(c(2, 4)),
+                       list())
 })
 
 test_that("hazard resolves two competing outcomes", {
@@ -62,9 +63,11 @@ test_that("hazard resolves two competing outcomes", {
   hazard$resolve(0)
    
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(2, 3)))
+                       individual::Bitset$new(size)$insert(c(2, 3)),
+                       list())
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(1, 4)))
+                       individual::Bitset$new(size)$insert(c(1, 4)),
+                       list())
 })
 
 test_that("hazard may resolve to neither outcome", {
@@ -94,9 +97,11 @@ test_that("hazard may resolve to neither outcome", {
   hazard$resolve(0)
    
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(3)))
+                       individual::Bitset$new(size)$insert(c(3)),
+                       list())
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(2, 4)))
+                       individual::Bitset$new(size)$insert(c(2, 4)),
+                       list())
 })
 
 test_that("outcomes can define a partial set of rates", {
@@ -128,10 +133,11 @@ test_that("outcomes can define a partial set of rates", {
   hazard$resolve(0)
    
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(1, 3)))
+                       individual::Bitset$new(size)$insert(c(1, 3)),
+                       list())
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(4))
-  )
+                       individual::Bitset$new(size)$insert(c(4)),
+                       list())
 })
 
 test_that("hazard resolves three competing outcomes", {
@@ -167,10 +173,12 @@ test_that("hazard resolves three competing outcomes", {
   hazard$resolve(0)
    
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(1)))
+                       individual::Bitset$new(size)$insert(c(1)),
+                       list())
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(2)))
+                       individual::Bitset$new(size)$insert(c(2)),
+                       list())
   mockery::expect_args(outcome_3_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(3, 4)))
+                       individual::Bitset$new(size)$insert(c(3, 4)),
+                       list())
 })
-
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index c29420c9..a4bcd1ab 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -17,20 +17,20 @@ 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)
+  n_bites_per_person <- NULL
   
   infection_outcome <- CompetingOutcome$new(
-    targeted_process = function(timestep, target){
-      infection_outcome_process(timestep, target, variables, renderer, parameters)
+    targeted_process = function(timestep, target, args){
+      falciparum_infection_outcome_process(timestep, target, variables, renderer, parameters)
     },
     size = parameters$human_population
   )
   
   mockery::stub(simulate_infection, 'boost_immunity', boost_immunity_mock)
-  mockery::stub(simulate_infection, 'calculate_infections', infection_mock)
+  mockery::stub(simulate_infection, 'calculate_falciparum_infections', infection_mock)
   
   simulate_infection(
     variables,
@@ -59,7 +59,6 @@ test_that('simulate_infection integrates different types of infection and schedu
     1,
     variables,
     bitten,
-    n_bites_per_person,
     parameters,
     renderer,
     timestep,
@@ -85,7 +84,6 @@ test_that('simulate_infection integrates different types of infection and schedu
     # state = list(get_index_of = mockery::mock(asymptomatics, cycle = T))
     state = individual::CategoricalVariable$new(categories = c("S","A","U","D","Tr"), initial_values = rep("S", population))#list(get_index_of = mockery::mock(asymptomatics, cycle = T))
   )
-  prob <- rep(0.5,population)
   
   source_humans <- individual::Bitset$new(population)$insert(c(1, 2, 3, 5))
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
@@ -102,20 +100,19 @@ test_that('simulate_infection integrates different types of infection and schedu
   to_A <- infected$and(clinical$not(FALSE))
   to_U <- NULL
 
-  mockery::stub(infection_outcome_process, 'incidence_renderer', mockery::mock())
-  mockery::stub(infection_outcome_process, 'boost_immunity', boost_immunity_mock)
-  mockery::stub(infection_outcome_process, 'calculate_clinical_infections', clinical_infection_mock)
-  mockery::stub(infection_outcome_process, 'update_severe_disease', severe_infection_mock)
-  mockery::stub(infection_outcome_process, 'calculate_treated', treated_mock)
-  mockery::stub(infection_outcome_process, 'schedule_infections', schedule_mock)
-  
-  infection_outcome_process(
+  mockery::stub(falciparum_infection_outcome_process, 'incidence_renderer', mockery::mock())
+  mockery::stub(falciparum_infection_outcome_process, 'boost_immunity', boost_immunity_mock)
+  mockery::stub(falciparum_infection_outcome_process, 'calculate_clinical_infections', clinical_infection_mock)
+  mockery::stub(falciparum_infection_outcome_process, 'update_severe_disease', severe_infection_mock)
+  mockery::stub(falciparum_infection_outcome_process, 'calculate_treated', treated_mock)
+  mockery::stub(falciparum_infection_outcome_process, 'schedule_infections', schedule_mock)
+  
+  falciparum_infection_outcome_process(
     timestep,
     infected,
     variables,
     renderer,
-    parameters,
-    prob)
+    parameters)
 
   mockery::expect_args(
     boost_immunity_mock,
@@ -207,43 +204,38 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   vaccine_antibodies_mock <- mockery::mock(c(2, 3))
   vaccine_efficacy_mock <- mockery::mock(c(.2, .3))
   bernoulli_mock <- mockery::mock(2)
-  mockery::stub(calculate_infections, 'blood_immunity', immunity_mock)
-  mockery::stub(calculate_infections, 'weibull_survival', weibull_mock)
-  mockery::stub(calculate_infections, 'calculate_pev_antibodies', vaccine_antibodies_mock)
-  mockery::stub(calculate_infections, 'calculate_pev_efficacy', vaccine_efficacy_mock)
-  mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock)
+  mockery::stub(calculate_falciparum_infections, 'blood_immunity', immunity_mock)
+  mockery::stub(calculate_falciparum_infections, 'bernoulli_multi_p', bernoulli_mock)
+  
+  local_mocked_bindings(weibull_survival = weibull_mock)
+  local_mocked_bindings(calculate_pev_antibodies = vaccine_antibodies_mock)
+  local_mocked_bindings(calculate_pev_efficacy = vaccine_efficacy_mock)
   
   # remove randomness from vaccine parameters
-  mockery::stub(
-    calculate_infections,
-    'sample_pev_param',
-    function(index, profiles, name) {
-      vnapply(index, function(i) profiles[[i]][[name]][[1]]) # return mu
-    },
-    depth = 4
-  )
+  local_mocked_bindings(sample_pev_param = function(index, profiles, name) {
+    vnapply(index, function(i) profiles[[i]][[name]][[1]]) # return mu
+  })
   
   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){
-      infection_outcome_process(timestep, target, variables, renderer, parameters)
+      falciparum_infection_outcome_process(timestep, target, variables, renderer, parameters)
     },
     size = 4
   )
   
-  infections <- calculate_infections(
+  infections <- calculate_falciparum_infections(
     variables,
     bitten_humans, 
-    n_bites_per_person,
     parameters,
     mock_render(timestep),
     timestep,
     infection_outcome
   )
-  
-  expect_equal(sum(infections!=0), 3)
+
+  expect_equal(sum(infection_outcome$rates!=0), 3)
   
   mockery::expect_args(immunity_mock, 1, c(.3, .5, .9), parameters)
   mockery::expect_args(
@@ -705,21 +697,19 @@ 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)
+  mockery::stub(calculate_falciparum_infections, 'bernoulli_multi_p', m)
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
-      infection_outcome_process(timestep, target, variables, renderer, parameters)
+      falciparum_infection_outcome_process(timestep, target, variables, renderer, parameters)
     },
     size = 4
   )
   
-  infection_rates <- calculate_infections(
+  infection_rates <- calculate_falciparum_infections(
     variables,
     bitten,
-    n_bites_per_person,
     parameters,
     mock_render(timestep),
     timestep,
@@ -727,7 +717,7 @@ test_that('prophylaxis is considered for medicated humans', {
   )
   
   expect_equal(
-    rate_to_prob(infection_rates[infection_rates!=0]),
+    rate_to_prob(infection_outcome$rates[infection_outcome$rates!=0]),
     c(2.491951e-07, 2.384032e-01, 5.899334e-01),
     tolerance = 1e-3
   )
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index ce41ac20..6abc67b1 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -152,37 +152,31 @@ test_that('Infection considers pev efficacy', {
 
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
-      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters)
+      falciparum_infection_outcome_process(timestep, target, variables, renderer, parameters)
     },
     size = 4
   )
   
   # remove randomness from infection sampling
   bernoulli_mock <- mockery::mock(c(1, 2))
-  mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock)
+  mockery::stub(calculate_falciparum_infections, 'bernoulli_multi_p', bernoulli_mock)
 
   # remove randomness from pev parameters
-  mockery::stub(
-    calculate_infections,
-    'sample_pev_param',
-    function(index, profiles, name) {
-      vnapply(index, function(i) profiles[[i]][[name]][[1]]) # return mu
-    },
-    depth = 4
-  )
-
-  infection_rates <- calculate_infections(
+  local_mocked_bindings(sample_pev_param = function(index, profiles, name) {
+    vnapply(index, function(i) profiles[[i]][[name]][[1]]) # return mu
+  })
+  
+  infection_rates <- calculate_falciparum_infections(
     variables = variables,
     bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
-    n_bites_per_person = numeric(0),
     parameters = parameters,
     renderer = mock_render(timestep),
     timestep = timestep,
     infection_outcome = infection_outcome
   )
-
+  
   expect_equal(
-    rate_to_prob(infection_rates[infection_rates!=0]),
+    rate_to_prob(infection_outcome$rates),
     c(0.590, 0.590, 0.215, 0.244),
     tolerance=1e-3
   )
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 4f762728..58ee404b 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -282,19 +282,18 @@ test_that('relapses are recognised with division between bite infections and rel
   
   bernoulli_mock <- mockery::mock(c(1, 3), 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)
+  mockery::stub(vivax_infection_outcome_process, 'bernoulli_multi_p', bernoulli_mock, depth = 2)
+  mockery::stub(vivax_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(
+  vivax_infection_outcome_process(
     timestep = timestep,
     infected_humans,
     variables,
     renderer,
     parameters,
-    prob = c(rep(0.5, 4)),
     relative_rate = c(rep(0.5, 3))
   )
 

From da381b948207c9c00a371533592b4a164449f596 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Tue, 15 Oct 2024 12:15:32 +0100
Subject: [PATCH 3/3] Removed code redundancies. Adjusted extra arguments
 inputted into competing hazards object and reverting to previous competing
 hazard tests.

---
 R/competing_hazards.R                   |  2 +-
 R/human_infection.R                     |  5 +---
 R/processes.R                           |  8 +++---
 tests/testthat/test-competing-hazards.R | 38 ++++++++++---------------
 4 files changed, 21 insertions(+), 32 deletions(-)

diff --git a/R/competing_hazards.R b/R/competing_hazards.R
index 0ef30e0d..ea7ffd97 100644
--- a/R/competing_hazards.R
+++ b/R/competing_hazards.R
@@ -25,7 +25,7 @@ CompetingOutcome <- R6::R6Class(
       self$args <- list(...)
     },
     execute = function(t, target){
-      do.call(private$targeted_process, list(t, target, self$args))
+      do.call(private$targeted_process, c(t, list(target), self$args))
     },
     reset = function() {
       self$target$clear()
diff --git a/R/human_infection.R b/R/human_infection.R
index c779efc7..36aa7b2a 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -181,9 +181,6 @@ calculate_vivax_infections <- function(
       infection_rates[hyp_source_index] <- infection_rates[hyp_source_index] + relapse_rates
       # Get and store relative rates for bite/relapse competing hazards resolution
       relative_rates <- relapse_rates/infection_rates[hyp_source_index]
-      # infection_outcome$set_relative_rates(
-      # hypnozoites_humans, 
-      # relative_rates)
     }
     
     prob <- rate_to_prob(infection_rates) * (1 - prophylaxis) * (1 - vaccine_efficacy)
@@ -944,7 +941,7 @@ severe_immunity <- function(age, acquired_immunity, maternal_immunity, parameter
   parameters$theta0 * (parameters$theta1 + (1 - parameters$theta1) / (
     1 + fv * (
       (acquired_immunity + maternal_immunity) / parameters$iv0) ** parameters$kv
-  )
+    )
   )
 }
 
diff --git a/R/processes.R b/R/processes.R
index 6ac4f290..3048ba6e 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -96,7 +96,7 @@ create_processes <- function(
   
   if(parameters$parasite == "falciparum"){
     infection_outcome <- CompetingOutcome$new(
-      targeted_process = function(timestep, target, args){
+      targeted_process = function(timestep, target){
         falciparum_infection_outcome_process(timestep, target, 
                                   variables, renderer, parameters
         )
@@ -106,10 +106,10 @@ create_processes <- function(
     
   } else if (parameters$parasite == "vivax"){
     infection_outcome <- CompetingOutcome$new(
-      targeted_process = function(timestep, target, args){
+      targeted_process = function(timestep, target, relative_rates){
         vivax_infection_outcome_process(timestep, target, 
                                   variables, renderer, parameters, 
-                                  args$relative_rates
+                                  relative_rates
         )
       },
       size = parameters$human_population
@@ -117,7 +117,7 @@ create_processes <- function(
   }
   
   progression_outcome <- CompetingOutcome$new(
-    targeted_process = function(timestep, target, ...){
+    targeted_process = function(timestep, target){
       progression_outcome_process(timestep, target, variables, parameters, renderer)
     },
     size = parameters$human_population
diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R
index fc817499..21ddefcb 100644
--- a/tests/testthat/test-competing-hazards.R
+++ b/tests/testthat/test-competing-hazards.R
@@ -26,14 +26,15 @@ test_that("hazard resolves two disjoint outcomes", {
   
   outcome_1$set_rates(population, c(10, 0, 10, 0))
   outcome_2$set_rates(population, c(0, 10, 0, 10))
-  
+
   hazard$resolve(0)
+
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(1, 3)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(1, 3)))
+
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(2, 4)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(2, 4)))
+  
 })
 
 test_that("hazard resolves two competing outcomes", {
@@ -63,11 +64,9 @@ test_that("hazard resolves two competing outcomes", {
   hazard$resolve(0)
    
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(2, 3)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(2, 3)))
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(1, 4)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(1, 4)))
 })
 
 test_that("hazard may resolve to neither outcome", {
@@ -97,11 +96,9 @@ test_that("hazard may resolve to neither outcome", {
   hazard$resolve(0)
    
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(3)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(3)))
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(2, 4)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(2, 4)))
 })
 
 test_that("outcomes can define a partial set of rates", {
@@ -133,11 +130,9 @@ test_that("outcomes can define a partial set of rates", {
   hazard$resolve(0)
    
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(1, 3)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(1, 3)))
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(4)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(4)))
 })
 
 test_that("hazard resolves three competing outcomes", {
@@ -173,12 +168,9 @@ test_that("hazard resolves three competing outcomes", {
   hazard$resolve(0)
    
   mockery::expect_args(outcome_1_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(1)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(1)))
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(2)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(2)))
   mockery::expect_args(outcome_3_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(3, 4)),
-                       list())
+                       individual::Bitset$new(size)$insert(c(3, 4)))
 })