From b36be83c34086b14c724450309791a734f527afd Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 9 May 2024 12:16:02 +0100 Subject: [PATCH 01/35] competing_hazards.R: Giovanni designed two objects: CompetingOutcome (containing an outcome following resolution and rates that that outcome occurs for each individual in the population) and CompetingHazard (contains all CompetingOutcomes and process function to resolve hazards). processes.R: create_biting_process now replaced with create_infection_rates_process and biting_process, that uses essentially the same functional pathway, but which stops when the rates of infection are generated, storing them in the infection_outcome object within human_infection.R Similarly, disease_progression_process(es) are replaced with new functions: create_recovery_rates_process and calculate_recovery_rates, that get the recovery rates for each individual based on disease state. We then add these processes to the process list, including the resolution function create_infection_recovery_hazard_process$resolve. This process (create_infection_recovery_hazard_process) contains infection_outcome and recovery_outcome. infection_outcome calls infection_process_resolved_hazard, which follows the rest of the infection pathway functions: assigning new infections to different human states. Note that this requires probability of infection which gets used in the incidence renderer function. recovery_outcome also updates human states and infectivities. I've adjusted some of the tests to reflect the changes in functions required to generate results. Removed white space deletions. --- R/biting_process.R | 71 +++--- R/competing_hazards.R | 78 ++++++ R/disease_progression.R | 125 ++++++---- R/human_infection.R | 166 +++++++------ R/processes.R | 88 ++++--- R/utils.R | 14 ++ tests/testthat/test-biology.R | 12 +- tests/testthat/test-biting-integration.R | 35 ++- tests/testthat/test-competing-hazards.R | 248 ++++++++++++++++++++ tests/testthat/test-infection-integration.R | 119 ++++++++-- tests/testthat/test-pev.R | 14 +- 11 files changed, 728 insertions(+), 242 deletions(-) create mode 100644 R/competing_hazards.R create mode 100644 tests/testthat/test-competing-hazards.R diff --git a/R/biting_process.R b/R/biting_process.R index 801aea20..8ae1147f 100644 --- a/R/biting_process.R +++ b/R/biting_process.R @@ -15,8 +15,10 @@ #' values (default: 1) #' @param mixing_index an index for this population's position in the #' lagged_infectivity list (default: 1) +#' @param infection_outcome competing hazards object for infection rates +#' @param timestep the current timestep #' @noRd -create_biting_process <- function( +biting_process <- function( renderer, solvers, models, @@ -26,37 +28,38 @@ create_biting_process <- function( lagged_infectivity, lagged_eir, mixing = 1, - mixing_index = 1 + mixing_index = 1, + infection_outcome, + timestep ) { - function(timestep) { - # Calculate combined EIR - age <- get_age(variables$birth$get_values(), timestep) - - bitten_humans <- simulate_bites( - renderer, - solvers, - models, - variables, - events, - age, - parameters, - timestep, - lagged_infectivity, - lagged_eir, - mixing, - mixing_index - ) - - simulate_infection( - variables, - events, - bitten_humans, - age, - parameters, - timestep, - renderer - ) - } + + age <- get_age(variables$birth$get_values(), timestep) + + bitten_humans <- simulate_bites( + renderer, + solvers, + models, + variables, + events, + age, + parameters, + timestep, + lagged_infectivity, + lagged_eir, + mixing, + mixing_index + ) + + simulate_infection( + variables, + events, + bitten_humans, + age, + parameters, + timestep, + renderer, + infection_outcome + ) } #' @importFrom stats rpois @@ -74,6 +77,7 @@ simulate_bites <- function( mixing = 1, mixing_index = 1 ) { + bitten_humans <- individual::Bitset$new(parameters$human_population) human_infectivity <- variables$infectivity$get_values() @@ -210,9 +214,7 @@ simulate_bites <- function( calculate_eir <- function(species, solvers, variables, parameters, timestep) { a <- human_blood_meal_rate(species, variables, parameters, timestep) infectious <- calculate_infectious(species, solvers, variables, parameters) - age <- get_age(variables$birth$get_values(), timestep) - psi <- unique_biting_rate(age, parameters) - infectious * a * mean(psi) + infectious * a } effective_biting_rates <- function(a, .pi, p_bitten) { @@ -246,7 +248,6 @@ calculate_infectious_individual <- function( species_index, parameters ) { - infectious_index$copy()$and(species_index)$size() } diff --git a/R/competing_hazards.R b/R/competing_hazards.R new file mode 100644 index 00000000..cdc6aea0 --- /dev/null +++ b/R/competing_hazards.R @@ -0,0 +1,78 @@ +## Define classes to resolve competing hazards +CompetingOutcome <- R6::R6Class( + "CompetingOutcome", + private = list( + targeted_process = NULL + ), + public = list( + initialize = function(targeted_process, size){ + if (!is.function(targeted_process)){ + stop("targeted_process must be a function") + } + if (!is.numeric(size) || size <= 0){ + stop("size must be positive integer") + } + private$targeted_process <- targeted_process + self$rates <- rep(0, size) + }, + set_rates = function(rates){ + self$rates <- rates + }, + execute = function(t, target){ + private$targeted_process(t, target) + self$rates <- rep(0, length(self$rates)) + }, + rates = NULL + ) +) + +CompetingHazard <- R6::R6Class( + "CompetingHazard", + private = list( + outcomes = list(), + size = NULL, + # RNG is passed in because mockery is not able to stub runif + # TODO: change when fixed + rng = NULL + ), + public = list( + initialize = function(outcomes, rng = runif){ + if (length(outcomes) == 0){ + stop("At least one outcome must be provided") + } + if (!all(sapply(outcomes, function(x) inherits(x, "CompetingOutcome")))){ + stop("All outcomes must be of class CompetingOutcome") + } + private$outcomes <- outcomes + private$size <- length(outcomes[[1]]$rates) + private$rng <- rng + }, + resolve = function(t){ + event_rates <- do.call( + 'cbind', + lapply(private$outcomes, function(x) x$rates) + ) + occur_rates <- rowSums(event_rates) + occur_rng <- private$rng(private$size) + occurs <- occur_rng < rate_to_prob(occur_rates) + norm_probs <- event_rates / occur_rates + norm_probs[is.na(norm_probs)] <- 0 + cum_probs <- apply(norm_probs, 1, cumsum) + event_rng <- private$rng(private$size) + for(o in seq_along(private$outcomes)){ + if (o == 1) { + selected <- (event_rng <= cum_probs[o,]) + } else { + selected <- (event_rng > cum_probs[o - 1,]) & + (event_rng <= cum_probs[o,]) + } + target <- individual::Bitset$new(private$size)$insert( + which(selected & occurs) + ) + if (target$size() > 0){ + private$outcomes[[o]]$execute(t, target) + } + } + } + ) +) diff --git a/R/disease_progression.R b/R/disease_progression.R index d9ac11e4..d72d2c75 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -1,3 +1,79 @@ +#' @title Calculate recovery rates +#' @description Calculates recovery rates for each individual in the population +#' for storage in competing hazards object and future resolution +#' +#' @param variables the available human variables +#' @param parameters model parameters +#' @param recovery_outcome competing hazards object for recovery rates +#' @noRd +calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt_input){ + + # Get correct input for dt depending on spc + if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ + dt_v <- dt_input + } else if (isTRUE(parameters$antimalarial_resistance)) { + dt_v <- dt_input$and(variables$state$get_index_of("Tr")$to_vector())$to_vector() + } + + recovery_rates <- numeric(length = parameters$human_population) + recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd + recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da + recovery_rates[variables$state$get_index_of("U")$to_vector()] <- 1/parameters$du + recovery_rates[variables$state$get_index_of("Tr")$to_vector()] <- 1/dt_v + recovery_outcome$set_rates(recovery_rates) +} + +#' @title Disease progression outcomes (recovery) +#' @description Following resolution of competing hazards, update state and +#' infectivity of sampled individuals +#' +#' @param timestep the current timestep +#' @param target the sampled recovering individuals +#' @param variables the available human variables +#' @param parameters model parameters +#' @param renderer competing hazards object for recovery rates +#' @noRd +recovery_process_resolved_hazard <- function( + timestep, + target, + variables, + parameters, + renderer +){ + + update_to_asymptomatic_infection( + variables, + parameters, + timestep, + variables$state$get_index_of("D")$and(target) + ) + + update_infection( + variables$state, + "U", + variables$infectivity, + parameters$cu, + variables$state$get_index_of("A")$and(target) + ) + + update_infection( + variables$state, + "S", + variables$infectivity, + 0, + variables$state$get_index_of("U")$and(target) + ) + + update_infection( + variables$state, + "S", + variables$infectivity, + 0, + variables$state$get_index_of("Tr")$and(target) + ) + +} + #' @title Update the state of an individual as infection events occur #' @description Randomly moves individuals towards the later stages of disease #' and updates their infectivity @@ -18,55 +94,6 @@ update_infection <- function( infectivity$queue_update(new_infectivity, to_move) } -create_progression_process <- function( - state, - from_state, - to_state, - rate, - infectivity, - new_infectivity -) { - function(timestep) { - - # Retrieve the indices of all individuals in the to_move state: - index <- state$get_index_of(from_state) - - # If the length of rate is greater than 1 (when it's a variable): - if (inherits(rate, "DoubleVariable")) { - rate <- rate$get_values(index) - } - - # Sample the individuals to be moved into a new Bitset using the transition rate(s): - to_move <- index$sample(1/rate) - - # Update the infection status of those individuals who are moving: - update_infection( - state, - to_state, - infectivity, - new_infectivity, - to_move - ) - } -} - -create_asymptomatic_progression_process <- function( - state, - rate, - variables, - parameters - ) { - function(timestep) { - to_move <- state$get_index_of('D')$sample(1/rate) - update_to_asymptomatic_infection( - variables, - parameters, - timestep, - to_move - ) - } -} - #' @title Modelling the progression to asymptomatic disease #' @description Randomly moves individuals to asymptomatic disease and #' calculates the infectivity for their age and immunity diff --git a/R/human_infection.R b/R/human_infection.R index d62b16f0..df2e0b1b 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -1,13 +1,15 @@ #' @title Simulate malaria infection in humans #' @description -#' Updates human states and variables to represent asymptomatic/clinical/severe -#' and treated malaria; and resulting boosts in immunity +#' This function ends with the assignment of rates of infection to the competing +#' hazard resolution object. Boosts immunity given infectious bites. #' @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 age of each human (timesteps) #' @param parameters of the model #' @param timestep current timestep +#' @param renderer the model renderer object +#' @param infection_outcome competing hazards object for infection rates #' @noRd simulate_infection <- function( variables, @@ -16,7 +18,8 @@ simulate_infection <- function( age, parameters, timestep, - renderer + renderer, + infection_outcome ) { if (bitten_humans$size() > 0) { boost_immunity( @@ -29,70 +32,18 @@ simulate_infection <- function( } # Calculate Infected - infected_humans <- calculate_infections( + calculate_infections( variables, bitten_humans, parameters, renderer, - timestep - ) - - if (infected_humans$size() > 0) { - 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_infections <- calculate_clinical_infections( - variables, - infected_humans, - parameters, - renderer, - timestep - ) - - update_severe_disease( - timestep, - infected_humans, - variables, - parameters, - renderer - ) - - treated <- calculate_treated( - variables, - clinical_infections, - parameters, timestep, - renderer - ) - - renderer$render('n_infections', infected_humans$size(), timestep) - - schedule_infections( - variables, - clinical_infections, - treated, - infected_humans, - parameters, - timestep + infection_outcome ) } #' @title Calculate overall infections for bitten humans -#' @description -#' Sample infected humans given prophylaxis and vaccination +#' @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 @@ -100,12 +51,13 @@ simulate_infection <- function( #' @param timestep current timestep #' @noRd calculate_infections <- function( - variables, - bitten_humans, - parameters, - renderer, - timestep - ) { + variables, + bitten_humans, + parameters, + renderer, + timestep, + infection_outcome +) { source_humans <- variables$state$get_index_of( c('S', 'A', 'U'))$and(bitten_humans) @@ -155,21 +107,97 @@ calculate_infections <- function( } prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy) - infected <- bitset_at(source_humans, bernoulli_multi_p(prob)) + ## Capture infection rates to resolve in competing hazards + infection_rates <- numeric(length = parameters$human_population) + infection_rates[source_vector] <- prob_to_rate(prob) + infection_outcome$set_rates(infection_rates) +} + +#' @title Assigns 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 prob vector of population probabilities of infection +#' @noRd +infection_process_resolved_hazard <- function( + timestep, + infected_humans, + variables, + renderer, + parameters, + prob){ + + source_humans <- variables$state$get_index_of(values = c("S","A","U")) + incidence_renderer( variables$birth, renderer, - infected, + infected_humans, source_humans, - prob, + prob[source_humans$to_vector()], 'inc_', parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages, timestep ) - - infected + + if (infected_humans$size() > 0) { + 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_infections <- calculate_clinical_infections( + variables, + infected_humans, + parameters, + renderer, + timestep + ) + + update_severe_disease( + timestep, + infected_humans, + variables, + parameters, + renderer + ) + + treated <- calculate_treated( + variables, + clinical_infections, + parameters, + timestep, + renderer + ) + + renderer$render('n_infections', infected_humans$size(), timestep) + + schedule_infections( + variables, + clinical_infections, + treated, + infected_humans, + parameters, + timestep + ) } #' @title Calculate clinical infections diff --git a/R/processes.R b/R/processes.R index 12f8c3db..322cd560 100644 --- a/R/processes.R +++ b/R/processes.R @@ -61,16 +61,37 @@ create_processes <- function( ) ) } + + # ===================================================== + # Competing Hazard Outcomes (Infections and Recoveries) + # ===================================================== + + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + }, + size = parameters$human_population + ) + + recovery_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + recovery_process_resolved_hazard(timestep, target, variables, parameters, renderer) + }, + size = parameters$human_population + ) + + create_infection_recovery_hazard_process <- CompetingHazard$new( + outcomes = list(infection_outcome, recovery_outcome) + ) # ============================== # Biting and mortality processes # ============================== - # schedule infections for humans and set last_boosted_* + # simulate infections for humans and set last_boosted_* # move mosquitoes into incubating state # kill mosquitoes caught in vector control - processes <- c( - processes, - create_biting_process( + create_infection_rates_process <- function(timestep){ + biting_process( renderer, solvers, models, @@ -80,31 +101,11 @@ create_processes <- function( lagged_infectivity, lagged_eir, mixing, - mixing_index - ), - create_asymptomatic_progression_process( - variables$state, - parameters$dd, - variables, - parameters - ), - create_progression_process( - variables$state, - 'A', - 'U', - parameters$da, - variables$infectivity, - parameters$cu - ), - create_progression_process( - variables$state, - 'U', - 'S', - parameters$du, - variables$infectivity, - 0 + mixing_index, + infection_outcome, + timestep ) - ) + } # ======================= # Antimalarial Resistance @@ -121,18 +122,23 @@ create_processes <- function( dt_input <- variables$dt } - # Create the progression process for Tr --> S specifying dt_input as the rate: - processes <- c( - processes, - create_progression_process( - variables$state, - 'Tr', - 'S', - dt_input, - variables$infectivity, - 0 + # =================== + # Disease Progression + # =================== + + create_recovery_rates_process <- function(timestep){ + calculate_recovery_rates( + variables, + parameters, + recovery_outcome, + dt_input ) - ) + } + + processes <- c(processes, + create_infection_rates_process, + create_recovery_rates_process, + create_infection_recovery_hazard_process$resolve) # =============== # ODE integration @@ -264,7 +270,11 @@ create_processes <- function( ) } + # ====================== # Mortality step + # ====================== + # Mortality is not resolved as a competing hazard + processes <- c( processes, create_mortality_process(variables, events, renderer, parameters)) diff --git a/R/utils.R b/R/utils.R index 2262c939..34de912e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -88,3 +88,17 @@ RandomState <- R6::R6Class( } ) ) + +#'@title Convert probability to a rate +#'@param prob probability +#'@noRd +prob_to_rate <- function(prob){ + -log(1 - prob) +} + +#'@title Convert rate to a probability +#'@param rate rate +#'@noRd +rate_to_prob <- function(rate){ + 1- exp(-rate) +} diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R index c739f7e7..77ee42f6 100644 --- a/tests/testthat/test-biology.R +++ b/tests/testthat/test-biology.R @@ -16,10 +16,7 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR', { vector_models <- parameterise_mosquito_models(parameters, 1) solvers <- parameterise_solvers(vector_models, parameters) estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0) - age <- get_age(variables$birth$get_values(), 0) - psi <- unique_biting_rate(age, parameters) - omega <- mean(psi) - mean(estimated_eir) / omega / population * 365 + mean(estimated_eir) / population * 365 }) expect_equal( @@ -44,10 +41,7 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR (with h vector_models <- parameterise_mosquito_models(parameters, 1) solvers <- parameterise_solvers(vector_models, parameters) estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0) - age <- get_age(variables$birth$get_values(), 0) - psi <- unique_biting_rate(age, parameters) - omega <- mean(psi) - mean(estimated_eir) / omega / population * 365 + mean(estimated_eir) / population * 365 }) expect_equal( @@ -92,7 +86,7 @@ test_that('FOIM is consistent with equilibrium', { psi <- unique_biting_rate(age, parameters) zeta <- variables$zeta$get_values() .pi <- human_pi(psi, zeta) - calculate_foim(a, sum(.pi * variables$infectivity$get_values()), 1.) + calculate_foim(a, sum(.pi * variables$infectivity$get_values()), mixing = 1) } ) expect_equal( diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R index 7da79b66..9712a243 100644 --- a/tests/testthat/test-biting-integration.R +++ b/tests/testthat/test-biting-integration.R @@ -14,7 +14,22 @@ test_that('biting_process integrates mosquito effects and human infection', { models <- parameterise_mosquito_models(parameters, timestep) solvers <- parameterise_solvers(models, parameters) - biting_process <- create_biting_process( + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + + bitten <- individual::Bitset$new(parameters$human_population) + bites_mock <- mockery::mock(bitten) + infection_mock <- mockery::mock() + + mockery::stub(biting_process, 'simulate_bites', bites_mock) + mockery::stub(biting_process, 'simulate_infection', infection_mock) + + + biting_process( renderer, solvers, models, @@ -22,17 +37,12 @@ test_that('biting_process integrates mosquito effects and human infection', { events, parameters, lagged_foim, - lagged_eir + lagged_eir, + mixing = 1, + mixing_index = 1, + infection_outcome, + timestep ) - - bitten <- individual::Bitset$new(parameters$human_population) - bites_mock <- mockery::mock(bitten) - infection_mock <- mockery::mock() - - mockery::stub(biting_process, 'simulate_bites', bites_mock) - mockery::stub(biting_process, 'simulate_infection', infection_mock) - biting_process(timestep) - mockery::expect_args( bites_mock, 1, @@ -59,7 +69,8 @@ test_that('biting_process integrates mosquito effects and human infection', { age, parameters, timestep, - renderer + renderer, + infection_outcome ) }) diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R new file mode 100644 index 00000000..ef67e2d0 --- /dev/null +++ b/tests/testthat/test-competing-hazards.R @@ -0,0 +1,248 @@ +# ============================== +# Test the CompetingHazard class +# ============================== + +test_that('hazard resolves two normalised outcomes when all events occur', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock( + c(0, 0, 0, 0), # all events occur + c(.05, .3, .2, .5) # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4)) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6)) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 3)) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(2, 4)) + ) +}) + +test_that('hazard resolves two unnormalised outcomes when all events occur', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock( + c(0, 0, 0, 0), # all events occur + c(.05, .3, .2, .5) # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 3)) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(2, 4)) + ) +}) + +test_that('hazard resolves two outcomes when some events occur', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock( + c(0, 0, 1, 1), # some events occur + c(.05, .3, .2, .5), # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(1) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(2) + ) +}) + + +test_that('hazard resolves when no rates are set', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2) + ) + + hazard$resolve(0) + + mockery::expect_called( + outcome_1_process, + 0 + ) + mockery::expect_called( + outcome_2_process, + 0 + ) +}) + +test_that('hazard resolves when rates are set to one', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2) + ) + + outcome_1$set_rates(rep(100, size)) + hazard$resolve(0) + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 2, 3, 4)) + ) + mockery::expect_called( + outcome_2_process, + 0 + ) +}) + +test_that('hazard resolves three outcomes', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + outcome_3_process <- mockery::mock() + outcome_3 <- CompetingOutcome$new( + targeted_process = outcome_3_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2, outcome_3), + rng = mockery::mock( + c(0, 0, 0, 0), # all events occur + c(.04, .3, .14, .6), # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + outcome_3$set_rates(c(0.5, 0.5, 0.5, 0.5)) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 3)) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(2) + ) + mockery::expect_args( + outcome_3_process, + 1, + 0, + individual::Bitset$new(size)$insert(4) + ) +}) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index 83d7def9..f6eacb80 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -20,22 +20,17 @@ test_that('simulate_infection integrates different types of infection and schedu boost_immunity_mock <- mockery::mock() infected <- individual::Bitset$new(population)$insert(c(1, 3, 5)) infection_mock <- mockery::mock(infected) - clinical <- individual::Bitset$new(population)$insert(c(1, 3)) - clinical_infection_mock <- mockery::mock(clinical) - severe <- individual::Bitset$new(population)$insert(1) - severe_infection_mock <- mockery::mock(severe) - treated <- individual::Bitset$new(population)$insert(3) - treated_mock <- mockery::mock(treated) - schedule_mock <- mockery::mock() + + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(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_clinical_infections', clinical_infection_mock) - mockery::stub(simulate_infection, 'update_severe_disease', severe_infection_mock) - mockery::stub(simulate_infection, 'calculate_treated', treated_mock) - mockery::stub(simulate_infection, 'schedule_infections', schedule_mock) - mockery::stub(simulate_infection, 'incidence_renderer', mockery::mock()) - mockery::stub(simulate_infection, 'clinical_incidence_renderer', mockery::mock()) + simulate_infection( variables, events, @@ -43,7 +38,8 @@ test_that('simulate_infection integrates different types of infection and schedu age, parameters, timestep, - renderer + renderer, + infection_outcome ) mockery::expect_args( @@ -63,7 +59,67 @@ test_that('simulate_infection integrates different types of infection and schedu bitten, parameters, renderer, - timestep + timestep, + infection_outcome + ) +}) + + +test_that('simulate_infection integrates different types of infection and scheduling', { + population <- 8 + timestep <- 5 + parameters <- get_parameters(list( + human_population = population + )) + events <- create_events(parameters) + renderer <- mock_render(timestep) + + age <- c(20, 24, 5, 39, 20, 24, 5, 39) * 365 + immunity <- c(.2, .3, .5, .9, .2, .3, .5, .9) + # asymptomatics <- mockery::mock() + variables <- list( + ib = individual::DoubleVariable$new(immunity), + id = individual::DoubleVariable$new(immunity), + # 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) + + infected <- individual::Bitset$new(population)$insert(c(1, 3, 5)) + clinical <- individual::Bitset$new(population)$insert(c(1, 3)) + clinical_infection_mock <- mockery::mock(clinical) + boost_immunity_mock <- mockery::mock() + severe <- individual::Bitset$new(population)$insert(1) + severe_infection_mock <- mockery::mock(severe) + treated <- individual::Bitset$new(population)$insert(3) + treated_mock <- mockery::mock(treated) + schedule_mock <- mockery::mock() + + + mockery::stub(infection_process_resolved_hazard, 'incidence_renderer', mockery::mock()) + mockery::stub(infection_process_resolved_hazard, 'boost_immunity', boost_immunity_mock) + mockery::stub(infection_process_resolved_hazard, 'calculate_clinical_infections', clinical_infection_mock) + mockery::stub(infection_process_resolved_hazard, 'update_severe_disease', severe_infection_mock) + mockery::stub(infection_process_resolved_hazard, 'calculate_treated', treated_mock) + mockery::stub(infection_process_resolved_hazard, 'schedule_infections', schedule_mock) + + infection_process_resolved_hazard( + timestep, + infected, + variables, + renderer, + parameters, + prob) + + + mockery::expect_args( + boost_immunity_mock, + 1, + variables$ica, + infected, + variables$last_boosted_ica, + 5, + parameters$uc ) mockery::expect_args( @@ -161,15 +217,23 @@ test_that('calculate_infections works various combinations of drug and vaccinati bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4)) + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + infections <- calculate_infections( variables, bitten_humans, parameters, mock_render(timestep), - timestep + timestep, + infection_outcome ) - expect_equal(infections$to_vector(), 3) + expect_equal(sum(infections!=0), 3) mockery::expect_args(immunity_mock, 1, c(.3, .5, .9), parameters) mockery::expect_args( @@ -198,12 +262,6 @@ test_that('calculate_infections works various combinations of drug and vaccinati c(rtss_profile$beta, rtss_booster_profile$beta), c(rtss_profile$alpha, rtss_booster_profile$alpha) ) - mockery::expect_args( - bernoulli_mock, - 1, - c(.2 * .8 * .8, .3 * .7, .4) - ) - }) test_that('calculate_clinical_infections correctly samples clinically infected', { @@ -642,16 +700,25 @@ test_that('prophylaxis is considered for medicated humans', { m <- mockery::mock(seq(3)) mockery::stub(calculate_infections, 'bernoulli_multi_p', m) - calculate_infections( + + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + + infection_rates <- calculate_infections( variables, bitten, parameters, mock_render(timestep), - timestep + timestep, + infection_outcome ) expect_equal( - mockery::mock_args(m)[[1]][[1]], + rate_to_prob(infection_rates[infection_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 09c4cf7d..9a5f1175 100644 --- a/tests/testthat/test-pev.R +++ b/tests/testthat/test-pev.R @@ -150,6 +150,13 @@ test_that('Infection considers pev efficacy', { rep(.2, 4) ) + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + # remove randomness from infection sampling bernoulli_mock <- mockery::mock(c(1, 2)) mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock) @@ -164,16 +171,17 @@ test_that('Infection considers pev efficacy', { depth = 4 ) - calculate_infections( + infection_rates <- calculate_infections( variables = variables, bitten_humans = individual::Bitset$new(4)$insert(seq(4)), parameters = parameters, renderer = mock_render(timestep), - timestep = timestep + timestep = timestep, + infection_outcome = infection_outcome ) expect_equal( - mockery::mock_args(bernoulli_mock)[[1]][[1]], + rate_to_prob(infection_rates[infection_rates!=0]), c(0.590, 0.590, 0.215, 0.244), tolerance=1e-3 ) From 0e05f75d7470b50b400d6b08c848cdf0d46998bf Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 9 May 2024 16:11:41 +0100 Subject: [PATCH 02/35] Added a stash source_human function to the competing hazards function for use in infection outcome. --- R/competing_hazards.R | 6 +++++- R/human_infection.R | 5 +++-- R/processes.R | 3 ++- tests/testthat/test-infection-integration.R | 2 ++ 4 files changed, 12 insertions(+), 4 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index cdc6aea0..4ad9b8f6 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -18,11 +18,15 @@ CompetingOutcome <- R6::R6Class( set_rates = function(rates){ self$rates <- rates }, + stash_source_humans = function(source_humans){ + self$source_humans <- source_humans + }, execute = function(t, target){ private$targeted_process(t, target) self$rates <- rep(0, length(self$rates)) }, - rates = NULL + rates = NULL, + source_humans = NULL ) ) diff --git a/R/human_infection.R b/R/human_infection.R index df2e0b1b..75534555 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -60,6 +60,7 @@ calculate_infections <- function( ) { source_humans <- variables$state$get_index_of( c('S', 'A', 'U'))$and(bitten_humans) + infection_outcome$stash_source_humans(source_humans) b <- blood_immunity(variables$ib$get_values(source_humans), parameters) @@ -120,6 +121,7 @@ calculate_infections <- function( #' and treated malaria; and resulting boosts in immunity #' @param timestep current timestep #' @param infected_humans bitset of infected humans +#' @param source_humans bitset of humans with infection potential #' @param variables a list of all of the model variables #' @param renderer model render object #' @param parameters model parameters @@ -128,13 +130,12 @@ calculate_infections <- function( infection_process_resolved_hazard <- function( timestep, infected_humans, + source_humans, variables, renderer, parameters, prob){ - source_humans <- variables$state$get_index_of(values = c("S","A","U")) - incidence_renderer( variables$birth, renderer, diff --git a/R/processes.R b/R/processes.R index 322cd560..ac664a8c 100644 --- a/R/processes.R +++ b/R/processes.R @@ -68,7 +68,8 @@ create_processes <- function( infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + infection_process_resolved_hazard(timestep, target, infection_outcome$source_humans, + variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) }, size = parameters$human_population ) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index f6eacb80..eda048a1 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -85,6 +85,7 @@ test_that('simulate_infection integrates different types of infection and schedu ) 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)) clinical <- individual::Bitset$new(population)$insert(c(1, 3)) clinical_infection_mock <- mockery::mock(clinical) @@ -106,6 +107,7 @@ test_that('simulate_infection integrates different types of infection and schedu infection_process_resolved_hazard( timestep, infected, + source_humans, variables, renderer, parameters, From f47bb3e0d261f14ccea58181981e3b806691fb8f Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 9 May 2024 12:16:02 +0100 Subject: [PATCH 03/35] competing_hazards.R: Giovanni designed two objects: CompetingOutcome (containing an outcome following resolution and rates that that outcome occurs for each individual in the population) and CompetingHazard (contains all CompetingOutcomes and process function to resolve hazards). processes.R: create_biting_process now replaced with create_infection_rates_process and biting_process, that uses essentially the same functional pathway, but which stops when the rates of infection are generated, storing them in the infection_outcome object within human_infection.R Similarly, disease_progression_process(es) are replaced with new functions: create_recovery_rates_process and calculate_recovery_rates, that get the recovery rates for each individual based on disease state. We then add these processes to the process list, including the resolution function create_infection_recovery_hazard_process$resolve. This process (create_infection_recovery_hazard_process) contains infection_outcome and recovery_outcome. infection_outcome calls infection_process_resolved_hazard, which follows the rest of the infection pathway functions: assigning new infections to different human states. Note that this requires probability of infection which gets used in the incidence renderer function. recovery_outcome also updates human states and infectivities. I've adjusted some of the tests to reflect the changes in functions required to generate results. Removed white space deletions. --- R/biting_process.R | 67 +++--- R/competing_hazards.R | 78 ++++++ R/disease_progression.R | 125 ++++++---- R/human_infection.R | 166 +++++++------ R/processes.R | 88 ++++--- R/utils.R | 14 ++ tests/testthat/test-biology.R | 2 +- tests/testthat/test-biting-integration.R | 35 ++- tests/testthat/test-competing-hazards.R | 248 ++++++++++++++++++++ tests/testthat/test-infection-integration.R | 119 ++++++++-- tests/testthat/test-pev.R | 14 +- 11 files changed, 725 insertions(+), 231 deletions(-) create mode 100644 R/competing_hazards.R create mode 100644 tests/testthat/test-competing-hazards.R diff --git a/R/biting_process.R b/R/biting_process.R index 279802ea..8ae1147f 100644 --- a/R/biting_process.R +++ b/R/biting_process.R @@ -15,8 +15,10 @@ #' values (default: 1) #' @param mixing_index an index for this population's position in the #' lagged_infectivity list (default: 1) +#' @param infection_outcome competing hazards object for infection rates +#' @param timestep the current timestep #' @noRd -create_biting_process <- function( +biting_process <- function( renderer, solvers, models, @@ -26,37 +28,38 @@ create_biting_process <- function( lagged_infectivity, lagged_eir, mixing = 1, - mixing_index = 1 + mixing_index = 1, + infection_outcome, + timestep ) { - function(timestep) { - # Calculate combined EIR - age <- get_age(variables$birth$get_values(), timestep) - - bitten_humans <- simulate_bites( - renderer, - solvers, - models, - variables, - events, - age, - parameters, - timestep, - lagged_infectivity, - lagged_eir, - mixing, - mixing_index - ) - - simulate_infection( - variables, - events, - bitten_humans, - age, - parameters, - timestep, - renderer - ) - } + + age <- get_age(variables$birth$get_values(), timestep) + + bitten_humans <- simulate_bites( + renderer, + solvers, + models, + variables, + events, + age, + parameters, + timestep, + lagged_infectivity, + lagged_eir, + mixing, + mixing_index + ) + + simulate_infection( + variables, + events, + bitten_humans, + age, + parameters, + timestep, + renderer, + infection_outcome + ) } #' @importFrom stats rpois @@ -74,6 +77,7 @@ simulate_bites <- function( mixing = 1, mixing_index = 1 ) { + bitten_humans <- individual::Bitset$new(parameters$human_population) human_infectivity <- variables$infectivity$get_values() @@ -244,7 +248,6 @@ calculate_infectious_individual <- function( species_index, parameters ) { - infectious_index$copy()$and(species_index)$size() } diff --git a/R/competing_hazards.R b/R/competing_hazards.R new file mode 100644 index 00000000..cdc6aea0 --- /dev/null +++ b/R/competing_hazards.R @@ -0,0 +1,78 @@ +## Define classes to resolve competing hazards +CompetingOutcome <- R6::R6Class( + "CompetingOutcome", + private = list( + targeted_process = NULL + ), + public = list( + initialize = function(targeted_process, size){ + if (!is.function(targeted_process)){ + stop("targeted_process must be a function") + } + if (!is.numeric(size) || size <= 0){ + stop("size must be positive integer") + } + private$targeted_process <- targeted_process + self$rates <- rep(0, size) + }, + set_rates = function(rates){ + self$rates <- rates + }, + execute = function(t, target){ + private$targeted_process(t, target) + self$rates <- rep(0, length(self$rates)) + }, + rates = NULL + ) +) + +CompetingHazard <- R6::R6Class( + "CompetingHazard", + private = list( + outcomes = list(), + size = NULL, + # RNG is passed in because mockery is not able to stub runif + # TODO: change when fixed + rng = NULL + ), + public = list( + initialize = function(outcomes, rng = runif){ + if (length(outcomes) == 0){ + stop("At least one outcome must be provided") + } + if (!all(sapply(outcomes, function(x) inherits(x, "CompetingOutcome")))){ + stop("All outcomes must be of class CompetingOutcome") + } + private$outcomes <- outcomes + private$size <- length(outcomes[[1]]$rates) + private$rng <- rng + }, + resolve = function(t){ + event_rates <- do.call( + 'cbind', + lapply(private$outcomes, function(x) x$rates) + ) + occur_rates <- rowSums(event_rates) + occur_rng <- private$rng(private$size) + occurs <- occur_rng < rate_to_prob(occur_rates) + norm_probs <- event_rates / occur_rates + norm_probs[is.na(norm_probs)] <- 0 + cum_probs <- apply(norm_probs, 1, cumsum) + event_rng <- private$rng(private$size) + for(o in seq_along(private$outcomes)){ + if (o == 1) { + selected <- (event_rng <= cum_probs[o,]) + } else { + selected <- (event_rng > cum_probs[o - 1,]) & + (event_rng <= cum_probs[o,]) + } + target <- individual::Bitset$new(private$size)$insert( + which(selected & occurs) + ) + if (target$size() > 0){ + private$outcomes[[o]]$execute(t, target) + } + } + } + ) +) diff --git a/R/disease_progression.R b/R/disease_progression.R index d9ac11e4..d72d2c75 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -1,3 +1,79 @@ +#' @title Calculate recovery rates +#' @description Calculates recovery rates for each individual in the population +#' for storage in competing hazards object and future resolution +#' +#' @param variables the available human variables +#' @param parameters model parameters +#' @param recovery_outcome competing hazards object for recovery rates +#' @noRd +calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt_input){ + + # Get correct input for dt depending on spc + if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ + dt_v <- dt_input + } else if (isTRUE(parameters$antimalarial_resistance)) { + dt_v <- dt_input$and(variables$state$get_index_of("Tr")$to_vector())$to_vector() + } + + recovery_rates <- numeric(length = parameters$human_population) + recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd + recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da + recovery_rates[variables$state$get_index_of("U")$to_vector()] <- 1/parameters$du + recovery_rates[variables$state$get_index_of("Tr")$to_vector()] <- 1/dt_v + recovery_outcome$set_rates(recovery_rates) +} + +#' @title Disease progression outcomes (recovery) +#' @description Following resolution of competing hazards, update state and +#' infectivity of sampled individuals +#' +#' @param timestep the current timestep +#' @param target the sampled recovering individuals +#' @param variables the available human variables +#' @param parameters model parameters +#' @param renderer competing hazards object for recovery rates +#' @noRd +recovery_process_resolved_hazard <- function( + timestep, + target, + variables, + parameters, + renderer +){ + + update_to_asymptomatic_infection( + variables, + parameters, + timestep, + variables$state$get_index_of("D")$and(target) + ) + + update_infection( + variables$state, + "U", + variables$infectivity, + parameters$cu, + variables$state$get_index_of("A")$and(target) + ) + + update_infection( + variables$state, + "S", + variables$infectivity, + 0, + variables$state$get_index_of("U")$and(target) + ) + + update_infection( + variables$state, + "S", + variables$infectivity, + 0, + variables$state$get_index_of("Tr")$and(target) + ) + +} + #' @title Update the state of an individual as infection events occur #' @description Randomly moves individuals towards the later stages of disease #' and updates their infectivity @@ -18,55 +94,6 @@ update_infection <- function( infectivity$queue_update(new_infectivity, to_move) } -create_progression_process <- function( - state, - from_state, - to_state, - rate, - infectivity, - new_infectivity -) { - function(timestep) { - - # Retrieve the indices of all individuals in the to_move state: - index <- state$get_index_of(from_state) - - # If the length of rate is greater than 1 (when it's a variable): - if (inherits(rate, "DoubleVariable")) { - rate <- rate$get_values(index) - } - - # Sample the individuals to be moved into a new Bitset using the transition rate(s): - to_move <- index$sample(1/rate) - - # Update the infection status of those individuals who are moving: - update_infection( - state, - to_state, - infectivity, - new_infectivity, - to_move - ) - } -} - -create_asymptomatic_progression_process <- function( - state, - rate, - variables, - parameters - ) { - function(timestep) { - to_move <- state$get_index_of('D')$sample(1/rate) - update_to_asymptomatic_infection( - variables, - parameters, - timestep, - to_move - ) - } -} - #' @title Modelling the progression to asymptomatic disease #' @description Randomly moves individuals to asymptomatic disease and #' calculates the infectivity for their age and immunity diff --git a/R/human_infection.R b/R/human_infection.R index d62b16f0..df2e0b1b 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -1,13 +1,15 @@ #' @title Simulate malaria infection in humans #' @description -#' Updates human states and variables to represent asymptomatic/clinical/severe -#' and treated malaria; and resulting boosts in immunity +#' This function ends with the assignment of rates of infection to the competing +#' hazard resolution object. Boosts immunity given infectious bites. #' @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 age of each human (timesteps) #' @param parameters of the model #' @param timestep current timestep +#' @param renderer the model renderer object +#' @param infection_outcome competing hazards object for infection rates #' @noRd simulate_infection <- function( variables, @@ -16,7 +18,8 @@ simulate_infection <- function( age, parameters, timestep, - renderer + renderer, + infection_outcome ) { if (bitten_humans$size() > 0) { boost_immunity( @@ -29,70 +32,18 @@ simulate_infection <- function( } # Calculate Infected - infected_humans <- calculate_infections( + calculate_infections( variables, bitten_humans, parameters, renderer, - timestep - ) - - if (infected_humans$size() > 0) { - 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_infections <- calculate_clinical_infections( - variables, - infected_humans, - parameters, - renderer, - timestep - ) - - update_severe_disease( - timestep, - infected_humans, - variables, - parameters, - renderer - ) - - treated <- calculate_treated( - variables, - clinical_infections, - parameters, timestep, - renderer - ) - - renderer$render('n_infections', infected_humans$size(), timestep) - - schedule_infections( - variables, - clinical_infections, - treated, - infected_humans, - parameters, - timestep + infection_outcome ) } #' @title Calculate overall infections for bitten humans -#' @description -#' Sample infected humans given prophylaxis and vaccination +#' @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 @@ -100,12 +51,13 @@ simulate_infection <- function( #' @param timestep current timestep #' @noRd calculate_infections <- function( - variables, - bitten_humans, - parameters, - renderer, - timestep - ) { + variables, + bitten_humans, + parameters, + renderer, + timestep, + infection_outcome +) { source_humans <- variables$state$get_index_of( c('S', 'A', 'U'))$and(bitten_humans) @@ -155,21 +107,97 @@ calculate_infections <- function( } prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy) - infected <- bitset_at(source_humans, bernoulli_multi_p(prob)) + ## Capture infection rates to resolve in competing hazards + infection_rates <- numeric(length = parameters$human_population) + infection_rates[source_vector] <- prob_to_rate(prob) + infection_outcome$set_rates(infection_rates) +} + +#' @title Assigns 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 prob vector of population probabilities of infection +#' @noRd +infection_process_resolved_hazard <- function( + timestep, + infected_humans, + variables, + renderer, + parameters, + prob){ + + source_humans <- variables$state$get_index_of(values = c("S","A","U")) + incidence_renderer( variables$birth, renderer, - infected, + infected_humans, source_humans, - prob, + prob[source_humans$to_vector()], 'inc_', parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages, timestep ) - - infected + + if (infected_humans$size() > 0) { + 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_infections <- calculate_clinical_infections( + variables, + infected_humans, + parameters, + renderer, + timestep + ) + + update_severe_disease( + timestep, + infected_humans, + variables, + parameters, + renderer + ) + + treated <- calculate_treated( + variables, + clinical_infections, + parameters, + timestep, + renderer + ) + + renderer$render('n_infections', infected_humans$size(), timestep) + + schedule_infections( + variables, + clinical_infections, + treated, + infected_humans, + parameters, + timestep + ) } #' @title Calculate clinical infections diff --git a/R/processes.R b/R/processes.R index 12f8c3db..322cd560 100644 --- a/R/processes.R +++ b/R/processes.R @@ -61,16 +61,37 @@ create_processes <- function( ) ) } + + # ===================================================== + # Competing Hazard Outcomes (Infections and Recoveries) + # ===================================================== + + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + }, + size = parameters$human_population + ) + + recovery_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + recovery_process_resolved_hazard(timestep, target, variables, parameters, renderer) + }, + size = parameters$human_population + ) + + create_infection_recovery_hazard_process <- CompetingHazard$new( + outcomes = list(infection_outcome, recovery_outcome) + ) # ============================== # Biting and mortality processes # ============================== - # schedule infections for humans and set last_boosted_* + # simulate infections for humans and set last_boosted_* # move mosquitoes into incubating state # kill mosquitoes caught in vector control - processes <- c( - processes, - create_biting_process( + create_infection_rates_process <- function(timestep){ + biting_process( renderer, solvers, models, @@ -80,31 +101,11 @@ create_processes <- function( lagged_infectivity, lagged_eir, mixing, - mixing_index - ), - create_asymptomatic_progression_process( - variables$state, - parameters$dd, - variables, - parameters - ), - create_progression_process( - variables$state, - 'A', - 'U', - parameters$da, - variables$infectivity, - parameters$cu - ), - create_progression_process( - variables$state, - 'U', - 'S', - parameters$du, - variables$infectivity, - 0 + mixing_index, + infection_outcome, + timestep ) - ) + } # ======================= # Antimalarial Resistance @@ -121,18 +122,23 @@ create_processes <- function( dt_input <- variables$dt } - # Create the progression process for Tr --> S specifying dt_input as the rate: - processes <- c( - processes, - create_progression_process( - variables$state, - 'Tr', - 'S', - dt_input, - variables$infectivity, - 0 + # =================== + # Disease Progression + # =================== + + create_recovery_rates_process <- function(timestep){ + calculate_recovery_rates( + variables, + parameters, + recovery_outcome, + dt_input ) - ) + } + + processes <- c(processes, + create_infection_rates_process, + create_recovery_rates_process, + create_infection_recovery_hazard_process$resolve) # =============== # ODE integration @@ -264,7 +270,11 @@ create_processes <- function( ) } + # ====================== # Mortality step + # ====================== + # Mortality is not resolved as a competing hazard + processes <- c( processes, create_mortality_process(variables, events, renderer, parameters)) diff --git a/R/utils.R b/R/utils.R index 2262c939..34de912e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -88,3 +88,17 @@ RandomState <- R6::R6Class( } ) ) + +#'@title Convert probability to a rate +#'@param prob probability +#'@noRd +prob_to_rate <- function(prob){ + -log(1 - prob) +} + +#'@title Convert rate to a probability +#'@param rate rate +#'@noRd +rate_to_prob <- function(rate){ + 1- exp(-rate) +} diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R index e429a48f..77ee42f6 100644 --- a/tests/testthat/test-biology.R +++ b/tests/testthat/test-biology.R @@ -86,7 +86,7 @@ test_that('FOIM is consistent with equilibrium', { psi <- unique_biting_rate(age, parameters) zeta <- variables$zeta$get_values() .pi <- human_pi(psi, zeta) - calculate_foim(a, sum(.pi * variables$infectivity$get_values()), 1.) + calculate_foim(a, sum(.pi * variables$infectivity$get_values()), mixing = 1) } ) expect_equal( diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R index 7da79b66..9712a243 100644 --- a/tests/testthat/test-biting-integration.R +++ b/tests/testthat/test-biting-integration.R @@ -14,7 +14,22 @@ test_that('biting_process integrates mosquito effects and human infection', { models <- parameterise_mosquito_models(parameters, timestep) solvers <- parameterise_solvers(models, parameters) - biting_process <- create_biting_process( + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + + bitten <- individual::Bitset$new(parameters$human_population) + bites_mock <- mockery::mock(bitten) + infection_mock <- mockery::mock() + + mockery::stub(biting_process, 'simulate_bites', bites_mock) + mockery::stub(biting_process, 'simulate_infection', infection_mock) + + + biting_process( renderer, solvers, models, @@ -22,17 +37,12 @@ test_that('biting_process integrates mosquito effects and human infection', { events, parameters, lagged_foim, - lagged_eir + lagged_eir, + mixing = 1, + mixing_index = 1, + infection_outcome, + timestep ) - - bitten <- individual::Bitset$new(parameters$human_population) - bites_mock <- mockery::mock(bitten) - infection_mock <- mockery::mock() - - mockery::stub(biting_process, 'simulate_bites', bites_mock) - mockery::stub(biting_process, 'simulate_infection', infection_mock) - biting_process(timestep) - mockery::expect_args( bites_mock, 1, @@ -59,7 +69,8 @@ test_that('biting_process integrates mosquito effects and human infection', { age, parameters, timestep, - renderer + renderer, + infection_outcome ) }) diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R new file mode 100644 index 00000000..ef67e2d0 --- /dev/null +++ b/tests/testthat/test-competing-hazards.R @@ -0,0 +1,248 @@ +# ============================== +# Test the CompetingHazard class +# ============================== + +test_that('hazard resolves two normalised outcomes when all events occur', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock( + c(0, 0, 0, 0), # all events occur + c(.05, .3, .2, .5) # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4)) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6)) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 3)) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(2, 4)) + ) +}) + +test_that('hazard resolves two unnormalised outcomes when all events occur', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock( + c(0, 0, 0, 0), # all events occur + c(.05, .3, .2, .5) # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 3)) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(2, 4)) + ) +}) + +test_that('hazard resolves two outcomes when some events occur', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock( + c(0, 0, 1, 1), # some events occur + c(.05, .3, .2, .5), # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(1) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(2) + ) +}) + + +test_that('hazard resolves when no rates are set', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2) + ) + + hazard$resolve(0) + + mockery::expect_called( + outcome_1_process, + 0 + ) + mockery::expect_called( + outcome_2_process, + 0 + ) +}) + +test_that('hazard resolves when rates are set to one', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2) + ) + + outcome_1$set_rates(rep(100, size)) + hazard$resolve(0) + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 2, 3, 4)) + ) + mockery::expect_called( + outcome_2_process, + 0 + ) +}) + +test_that('hazard resolves three outcomes', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + outcome_3_process <- mockery::mock() + outcome_3 <- CompetingOutcome$new( + targeted_process = outcome_3_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2, outcome_3), + rng = mockery::mock( + c(0, 0, 0, 0), # all events occur + c(.04, .3, .14, .6), # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + outcome_3$set_rates(c(0.5, 0.5, 0.5, 0.5)) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 3)) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(2) + ) + mockery::expect_args( + outcome_3_process, + 1, + 0, + individual::Bitset$new(size)$insert(4) + ) +}) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index 83d7def9..f6eacb80 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -20,22 +20,17 @@ test_that('simulate_infection integrates different types of infection and schedu boost_immunity_mock <- mockery::mock() infected <- individual::Bitset$new(population)$insert(c(1, 3, 5)) infection_mock <- mockery::mock(infected) - clinical <- individual::Bitset$new(population)$insert(c(1, 3)) - clinical_infection_mock <- mockery::mock(clinical) - severe <- individual::Bitset$new(population)$insert(1) - severe_infection_mock <- mockery::mock(severe) - treated <- individual::Bitset$new(population)$insert(3) - treated_mock <- mockery::mock(treated) - schedule_mock <- mockery::mock() + + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(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_clinical_infections', clinical_infection_mock) - mockery::stub(simulate_infection, 'update_severe_disease', severe_infection_mock) - mockery::stub(simulate_infection, 'calculate_treated', treated_mock) - mockery::stub(simulate_infection, 'schedule_infections', schedule_mock) - mockery::stub(simulate_infection, 'incidence_renderer', mockery::mock()) - mockery::stub(simulate_infection, 'clinical_incidence_renderer', mockery::mock()) + simulate_infection( variables, events, @@ -43,7 +38,8 @@ test_that('simulate_infection integrates different types of infection and schedu age, parameters, timestep, - renderer + renderer, + infection_outcome ) mockery::expect_args( @@ -63,7 +59,67 @@ test_that('simulate_infection integrates different types of infection and schedu bitten, parameters, renderer, - timestep + timestep, + infection_outcome + ) +}) + + +test_that('simulate_infection integrates different types of infection and scheduling', { + population <- 8 + timestep <- 5 + parameters <- get_parameters(list( + human_population = population + )) + events <- create_events(parameters) + renderer <- mock_render(timestep) + + age <- c(20, 24, 5, 39, 20, 24, 5, 39) * 365 + immunity <- c(.2, .3, .5, .9, .2, .3, .5, .9) + # asymptomatics <- mockery::mock() + variables <- list( + ib = individual::DoubleVariable$new(immunity), + id = individual::DoubleVariable$new(immunity), + # 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) + + infected <- individual::Bitset$new(population)$insert(c(1, 3, 5)) + clinical <- individual::Bitset$new(population)$insert(c(1, 3)) + clinical_infection_mock <- mockery::mock(clinical) + boost_immunity_mock <- mockery::mock() + severe <- individual::Bitset$new(population)$insert(1) + severe_infection_mock <- mockery::mock(severe) + treated <- individual::Bitset$new(population)$insert(3) + treated_mock <- mockery::mock(treated) + schedule_mock <- mockery::mock() + + + mockery::stub(infection_process_resolved_hazard, 'incidence_renderer', mockery::mock()) + mockery::stub(infection_process_resolved_hazard, 'boost_immunity', boost_immunity_mock) + mockery::stub(infection_process_resolved_hazard, 'calculate_clinical_infections', clinical_infection_mock) + mockery::stub(infection_process_resolved_hazard, 'update_severe_disease', severe_infection_mock) + mockery::stub(infection_process_resolved_hazard, 'calculate_treated', treated_mock) + mockery::stub(infection_process_resolved_hazard, 'schedule_infections', schedule_mock) + + infection_process_resolved_hazard( + timestep, + infected, + variables, + renderer, + parameters, + prob) + + + mockery::expect_args( + boost_immunity_mock, + 1, + variables$ica, + infected, + variables$last_boosted_ica, + 5, + parameters$uc ) mockery::expect_args( @@ -161,15 +217,23 @@ test_that('calculate_infections works various combinations of drug and vaccinati bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4)) + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + infections <- calculate_infections( variables, bitten_humans, parameters, mock_render(timestep), - timestep + timestep, + infection_outcome ) - expect_equal(infections$to_vector(), 3) + expect_equal(sum(infections!=0), 3) mockery::expect_args(immunity_mock, 1, c(.3, .5, .9), parameters) mockery::expect_args( @@ -198,12 +262,6 @@ test_that('calculate_infections works various combinations of drug and vaccinati c(rtss_profile$beta, rtss_booster_profile$beta), c(rtss_profile$alpha, rtss_booster_profile$alpha) ) - mockery::expect_args( - bernoulli_mock, - 1, - c(.2 * .8 * .8, .3 * .7, .4) - ) - }) test_that('calculate_clinical_infections correctly samples clinically infected', { @@ -642,16 +700,25 @@ test_that('prophylaxis is considered for medicated humans', { m <- mockery::mock(seq(3)) mockery::stub(calculate_infections, 'bernoulli_multi_p', m) - calculate_infections( + + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + + infection_rates <- calculate_infections( variables, bitten, parameters, mock_render(timestep), - timestep + timestep, + infection_outcome ) expect_equal( - mockery::mock_args(m)[[1]][[1]], + rate_to_prob(infection_rates[infection_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 09c4cf7d..9a5f1175 100644 --- a/tests/testthat/test-pev.R +++ b/tests/testthat/test-pev.R @@ -150,6 +150,13 @@ test_that('Infection considers pev efficacy', { rep(.2, 4) ) + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + # remove randomness from infection sampling bernoulli_mock <- mockery::mock(c(1, 2)) mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock) @@ -164,16 +171,17 @@ test_that('Infection considers pev efficacy', { depth = 4 ) - calculate_infections( + infection_rates <- calculate_infections( variables = variables, bitten_humans = individual::Bitset$new(4)$insert(seq(4)), parameters = parameters, renderer = mock_render(timestep), - timestep = timestep + timestep = timestep, + infection_outcome = infection_outcome ) expect_equal( - mockery::mock_args(bernoulli_mock)[[1]][[1]], + rate_to_prob(infection_rates[infection_rates!=0]), c(0.590, 0.590, 0.215, 0.244), tolerance=1e-3 ) From c8c32452be47f29ebc29eade02619b9b3176224f Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 9 May 2024 16:11:41 +0100 Subject: [PATCH 04/35] Added a stash source_human function to the competing hazards function for use in infection outcome. --- R/competing_hazards.R | 6 +++++- R/human_infection.R | 5 +++-- R/processes.R | 3 ++- tests/testthat/test-infection-integration.R | 2 ++ 4 files changed, 12 insertions(+), 4 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index cdc6aea0..4ad9b8f6 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -18,11 +18,15 @@ CompetingOutcome <- R6::R6Class( set_rates = function(rates){ self$rates <- rates }, + stash_source_humans = function(source_humans){ + self$source_humans <- source_humans + }, execute = function(t, target){ private$targeted_process(t, target) self$rates <- rep(0, length(self$rates)) }, - rates = NULL + rates = NULL, + source_humans = NULL ) ) diff --git a/R/human_infection.R b/R/human_infection.R index df2e0b1b..75534555 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -60,6 +60,7 @@ calculate_infections <- function( ) { source_humans <- variables$state$get_index_of( c('S', 'A', 'U'))$and(bitten_humans) + infection_outcome$stash_source_humans(source_humans) b <- blood_immunity(variables$ib$get_values(source_humans), parameters) @@ -120,6 +121,7 @@ calculate_infections <- function( #' and treated malaria; and resulting boosts in immunity #' @param timestep current timestep #' @param infected_humans bitset of infected humans +#' @param source_humans bitset of humans with infection potential #' @param variables a list of all of the model variables #' @param renderer model render object #' @param parameters model parameters @@ -128,13 +130,12 @@ calculate_infections <- function( infection_process_resolved_hazard <- function( timestep, infected_humans, + source_humans, variables, renderer, parameters, prob){ - source_humans <- variables$state$get_index_of(values = c("S","A","U")) - incidence_renderer( variables$birth, renderer, diff --git a/R/processes.R b/R/processes.R index 322cd560..ac664a8c 100644 --- a/R/processes.R +++ b/R/processes.R @@ -68,7 +68,8 @@ create_processes <- function( infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + infection_process_resolved_hazard(timestep, target, infection_outcome$source_humans, + variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) }, size = parameters$human_population ) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index f6eacb80..eda048a1 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -85,6 +85,7 @@ test_that('simulate_infection integrates different types of infection and schedu ) 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)) clinical <- individual::Bitset$new(population)$insert(c(1, 3)) clinical_infection_mock <- mockery::mock(clinical) @@ -106,6 +107,7 @@ test_that('simulate_infection integrates different types of infection and schedu infection_process_resolved_hazard( timestep, infected, + source_humans, variables, renderer, parameters, From c5947ef4f37edab9c4eb09c4677b42f89068dd98 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 23 May 2024 10:12:38 +0100 Subject: [PATCH 05/35] To correctly render ages and incidence, given the competing hazards solution (any rendering in the infection outcome does not occur), I've split the incidence renderer into three functions: age rendering, n incidence rendering and p incidence rendering. This allows much more flexibility in when these are called in the code sequence. All age inputs are now consolidated and rendered in a single function. In addition, we intitate and populate these incidence columns with 0 when they are used. --- R/human_infection.R | 102 ++++++++----- R/model.R | 1 + R/parameters.R | 2 +- R/render.R | 158 +++++++++++++++----- tests/testthat/test-infection-integration.R | 14 +- tests/testthat/test-render.R | 92 +++++++++--- vignettes/Vaccines.Rmd | 2 +- 7 files changed, 275 insertions(+), 96 deletions(-) diff --git a/R/human_infection.R b/R/human_infection.R index 75534555..50568a26 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -12,15 +12,15 @@ #' @param infection_outcome competing hazards object for infection rates #' @noRd simulate_infection <- function( - variables, - events, - bitten_humans, - age, - parameters, - timestep, - renderer, - infection_outcome - ) { + variables, + events, + bitten_humans, + age, + parameters, + timestep, + renderer, + infection_outcome +) { if (bitten_humans$size() > 0) { boost_immunity( variables$ib, @@ -109,7 +109,19 @@ calculate_infections <- function( prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy) - ## Capture infection rates to resolve in competing hazards + ## probability of incidence must be rendered at each timestep + incidence_probability_renderer( + variables$birth, + renderer, + source_pop, + prob, + "inc_", + parameters$incidence_min_ages, + parameters$incidence_max_ages, + timestep + ) + + ## capture infection rates to resolve in competing hazards infection_rates <- numeric(length = parameters$human_population) infection_rates[source_vector] <- prob_to_rate(prob) infection_outcome$set_rates(infection_rates) @@ -140,8 +152,6 @@ infection_process_resolved_hazard <- function( variables$birth, renderer, infected_humans, - source_humans, - prob[source_humans$to_vector()], 'inc_', parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages, @@ -211,12 +221,12 @@ infection_process_resolved_hazard <- function( #' @param timestep current timestep #' @noRd calculate_clinical_infections <- function( - variables, - infections, - parameters, - renderer, - timestep - ) { + variables, + infections, + parameters, + renderer, + timestep +) { ica <- variables$ica$get_values(infections) icm <- variables$icm$get_values(infections) phi <- clinical_immunity(ica, icm, parameters) @@ -225,6 +235,14 @@ calculate_clinical_infections <- function( variables$birth, renderer, clinical_infections, + 'inc_clinical_', + parameters$clinical_incidence_rendering_min_ages, + parameters$clinical_incidence_rendering_max_ages, + timestep + ) + incidence_probability_renderer( + variables$birth, + renderer, infections, phi, 'inc_clinical_', @@ -245,12 +263,12 @@ calculate_clinical_infections <- function( #' @param renderer model outputs #' @noRd update_severe_disease <- function( - timestep, - infections, - variables, - parameters, - renderer - ) { + timestep, + infections, + variables, + parameters, + renderer +) { age <- get_age(variables$birth$get_values(infections), timestep) iva <- variables$iva$get_values(infections) ivm <- variables$ivm$get_values(infections) @@ -266,6 +284,14 @@ update_severe_disease <- function( variables$birth, renderer, severe_infections, + 'inc_severe_', + parameters$severe_incidence_rendering_min_ages, + parameters$severe_incidence_rendering_max_ages, + timestep + ) + incidence_probability_renderer( + variables$birth, + renderer, infections, theta, 'inc_severe_', @@ -409,13 +435,13 @@ calculate_treated <- function( #' @param parameters model parameters #' @noRd schedule_infections <- function( - variables, - clinical_infections, - treated, - infections, - parameters, - timestep - ) { + variables, + clinical_infections, + treated, + infections, + parameters, + timestep +) { included <- treated$not(TRUE) to_infect <- clinical_infections$and(included) @@ -447,12 +473,12 @@ schedule_infections <- function( # Utility functions # ================= boost_immunity <- function( - immunity_variable, - exposed_index, - last_boosted_variable, - timestep, - delay - ) { + immunity_variable, + exposed_index, + last_boosted_variable, + timestep, + delay +) { # record who can be boosted exposed_index_vector <- exposed_index$to_vector() last_boosted <- last_boosted_variable$get_values(exposed_index) @@ -497,7 +523,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/model.R b/R/model.R index 67f79758..fb9b2fab 100644 --- a/R/model.R +++ b/R/model.R @@ -122,6 +122,7 @@ run_resumable_simulation <- function( events <- create_events(parameters) initialise_events(events, variables, parameters) renderer <- individual::Render$new(timesteps) + populate_incidence_rendering_columns(renderer, parameters) attach_event_listeners( events, variables, diff --git a/R/parameters.R b/R/parameters.R index 647b78b2..0bb6e6f6 100644 --- a/R/parameters.R +++ b/R/parameters.R @@ -415,7 +415,7 @@ get_parameters <- function(overrides = list()) { incidence_rendering_min_ages = numeric(0), incidence_rendering_max_ages = numeric(0), clinical_incidence_rendering_min_ages = numeric(0), - clinical_incidence_rendering_max_ages = 5 * 365, + clinical_incidence_rendering_max_ages = numeric(0), severe_prevalence_rendering_min_ages = numeric(0), severe_prevalence_rendering_max_ages = numeric(0), severe_incidence_rendering_min_ages = numeric(0), diff --git a/R/render.R b/R/render.R index 60cb5b73..96896741 100644 --- a/R/render.R +++ b/R/render.R @@ -15,12 +15,12 @@ in_age_range <- function(birth, timestep, lower, upper) { #' #' @noRd create_prevelance_renderer <- function( - state, - birth, - immunity, - parameters, - renderer - ) { + state, + birth, + immunity, + parameters, + renderer +) { function(timestep) { asymptomatic <- state$get_index_of('A') prob <- probability_of_detection( @@ -37,11 +37,6 @@ create_prevelance_renderer <- function( lower <- parameters$prevalence_rendering_min_ages[[i]] upper <- parameters$prevalence_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_detect_', lower, '_', upper), in_age$copy()$and(detected)$size(), @@ -60,13 +55,11 @@ create_prevelance_renderer <- function( #' @title Render incidence statistics #' -#' @description renders incidence (new for this timestep) for indivduals +#' @description renders incidence (new for this timestep) for individuals #' #' @param birth variable for birth of the individual #' @param renderer object for model outputs #' @param target incidence population -#' @param source_pop the population which is sampled for infection -#' @param prob probability of infection #' @param prefix for model outputs #' @param lowers age bounds #' @param uppers age bounds @@ -74,27 +67,54 @@ create_prevelance_renderer <- function( #' #' @noRd incidence_renderer <- function( - birth, - renderer, - target, - source_pop, - prob, - prefix, - lowers, - uppers, - timestep - ) { + birth, + renderer, + target, + prefix, + lowers, + uppers, + timestep +) { for (i in seq_along(lowers)) { lower <- lowers[[i]] upper <- uppers[[i]] in_age <- in_age_range(birth, timestep, lower, upper) - renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep) renderer$render( paste0('n_', prefix, lower, '_', upper), in_age$copy()$and(target)$size(), timestep ) + } +} +#' @title Render probability of incidence statistics +#' +#' @description renders probability of incidence (new for this timestep) for individuals +#' +#' @param birth variable for birth of the individual +#' @param renderer object for model outputs +#' @param source_pop the population which is sampled for infection (bitten SAU for incidence, infecte for clinical/severe) +#' @param prob probability of infection +#' @param prefix for model outputs +#' @param lowers age bounds +#' @param uppers age bounds +#' @param timestep current target +#' +#' @noRd +incidence_probability_renderer <- function( + birth, + renderer, + source_pop, + prob, + prefix, + lowers, + uppers, + timestep +) { + for (i in seq_along(lowers)) { + lower <- lowers[[i]] + upper <- uppers[[i]] + in_age <- in_age_range(birth, timestep, lower, upper) renderer$render( paste0('p_', prefix, lower, '_', upper), sum(prob[bitset_index(source_pop, in_age)]), @@ -104,9 +124,9 @@ incidence_renderer <- function( } create_variable_mean_renderer_process <- function( - renderer, - names, - variables + renderer, + names, + variables ) { function(timestep) { for (i in seq_along(variables)) { @@ -120,12 +140,12 @@ create_variable_mean_renderer_process <- function( } create_vector_count_renderer_individual <- function( - mosquito_state, - species, - state, - renderer, - parameters - ) { + mosquito_state, + species, + state, + renderer, + parameters +) { function(timestep) { adult <- mosquito_state$get_index_of('NonExistent')$not(TRUE) for (i in seq_along(parameters$species)) { @@ -160,15 +180,77 @@ create_age_group_renderer <- function( renderer ) { function(timestep) { - for (i in seq_along(parameters$age_group_rendering_min_ages)) { - lower <- parameters$age_group_rendering_min_ages[[i]] - upper <- parameters$age_group_rendering_max_ages[[i]] + + unique_age_combinations <- as.data.frame(unique(rbind(cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages), + cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages), + cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages), + cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages), + cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages)))) + + ordered_unique_age_combinations <- unique_age_combinations[ + with(unique_age_combinations, order(V1, V2)), + ] + + for (i in seq_along(ordered_unique_age_combinations$V1)) { + lower <- ordered_unique_age_combinations$V1[[i]] + upper <- ordered_unique_age_combinations$V2[[i]] in_age <- in_age_range(birth, timestep, lower, upper) renderer$render( - paste0('n_age_', lower, '_', upper), + paste0('n_', lower, '_', upper), in_age$size(), timestep ) } } } + + +populate_incidence_rendering_columns <- function(renderer, parameters){ + + # infections must render in all simulations + renderer$set_default('n_infections', 0) + + # treatment associated only renders when drugs are used + if(sum(unlist(parameters$clinical_treatment_coverages))>0){ + renderer$set_default('ft', 0) + renderer$set_default('n_treated', 0) + renderer$set_default('n_drug_efficacy_failures', 0) + renderer$set_default('n_successfully_treated', 0) + } + + # ETC, SPC only render when antimalarial resistance is on + if(parameters$antimalarial_resistance){ + renderer$set_default('n_early_treatment_failure', 0) + renderer$set_default('n_slow_parasite_clearance', 0) + } + + if(length(parameters$incidence_rendering_min_ages)>0){ + render_initial_incidence(renderer, + parameters$incidence_rendering_min_ages, + parameters$incidence_rendering_max_ages, + "inc") + } + + if(length(parameters$clinical_incidence_rendering_min_ages)>0){ + render_initial_incidence(renderer, + parameters$clinical_incidence_rendering_min_ages, + parameters$clinical_incidence_rendering_max_ages, + "inc_clinical") + } + + if(length(parameters$severe_incidence_rendering_min_ages)>0){ + render_initial_incidence(renderer, + parameters$severe_incidence_rendering_min_ages, + parameters$severe_incidence_rendering_max_ages, + "inc_severe") + } + +} + +render_initial_incidence <- function(renderer, lower_vals, upper_vals, inc_name){ + for (i in seq_along(lower_vals)){ + renderer$set_default(paste0('n_', inc_name, "_", lower_vals[i], '_', upper_vals[i]), 0) + renderer$set_default(paste0('p_', inc_name, "_", lower_vals[i], '_', upper_vals[i]), 0) + } +} + \ No newline at end of file diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index eda048a1..25637f0d 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -843,6 +843,7 @@ test_that('update_severe_disease renders with no infections', { render_function <- mockery::mock() mockery::stub(update_severe_disease, 'incidence_renderer', render_function) + mockery::stub(update_severe_disease, 'incidence_probability_renderer', render_function) empty <- individual::Bitset$new(population) update_severe_disease( @@ -852,13 +853,24 @@ test_that('update_severe_disease renders with no infections', { parameters, renderer ) - + mockery::expect_args( render_function, 1, variables$birth, renderer, empty, + 'inc_severe_', + 0, + 5 * 365, + timestep + ) + + mockery::expect_args( + render_function, + 2, + variables$birth, + renderer, empty, double(0), 'inc_severe_', diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R index 920f09bb..8b7a177b 100644 --- a/tests/testthat/test-render.R +++ b/tests/testthat/test-render.R @@ -26,14 +26,6 @@ test_that('that default rendering works', { mockery::expect_args( renderer$render_mock(), 1, - 'n_730_3650', - 3, - timestep - ) - - mockery::expect_args( - renderer$render_mock(), - 2, 'n_detect_730_3650', 2, timestep @@ -41,7 +33,7 @@ test_that('that default rendering works', { mockery::expect_args( renderer$render_mock(), - 3, + 2, 'p_detect_730_3650', 1.5, timestep @@ -75,7 +67,7 @@ test_that('that default rendering works when no one is in the age range', { mockery::expect_args( renderer$render_mock(), 1, - 'n_730_3650', + 'n_detect_730_3650', 0, timestep ) @@ -93,6 +85,15 @@ test_that('that clinical incidence rendering works', { birth, renderer, individual::Bitset$new(4)$insert(c(1, 2, 4)), + 'inc_clinical_', + c(0, 2) * year, + c(5, 10) * year, + timestep + ) + + incidence_probability_renderer( + birth, + renderer, individual::Bitset$new(4)$insert(seq(4)), c(.1, .2, .3, .4), 'inc_clinical_', @@ -104,7 +105,7 @@ test_that('that clinical incidence rendering works', { mockery::expect_args( renderer$render_mock(), 1, - 'n_0_1825', + 'n_inc_clinical_0_1825', 2, timestep ) @@ -112,7 +113,7 @@ test_that('that clinical incidence rendering works', { mockery::expect_args( renderer$render_mock(), 2, - 'n_inc_clinical_0_1825', + 'n_inc_clinical_730_3650', 2, timestep ) @@ -124,11 +125,68 @@ test_that('that clinical incidence rendering works', { .3, timestep ) + + mockery::expect_args( + renderer$render_mock(), + 4, + 'p_inc_clinical_730_3650', + .6, + timestep + ) +}) + +test_that('that age rendering works', { + timestep <- 0 + year <- 365 + birth <- individual::IntegerVariable$new( + -c(0:101) * year + ) + + parameters <- get_parameters(overrides = list(prevalence_rendering_min_ages = c(0, 1, 5) * year, + prevalence_rendering_max_ages = c(2, 5, 8) * year - 1, + incidence_rendering_min_ages = c(3, 5, 10) * year, + incidence_rendering_max_ages = c(10, 20, 100) * year - 1, + clinical_incidence_rendering_min_ages = c(0, 1, 5) * year, + clinical_incidence_rendering_max_ages = c(2, 5, 8) * year - 1)) + + renderer <- mock_render(1) + process <- create_age_group_renderer( + birth, + parameters, + renderer) + + process(timestep) + + mockery::expect_args( renderer$render_mock(), + 1, + 'n_age_0_729', + 2, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 2, + 'n_age_365_1824', 4, - 'n_730_3650', + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 3, + 'n_age_1095_3649', + 7, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 4, + 'n_age_1825_2919', 3, timestep ) @@ -136,16 +194,16 @@ test_that('that clinical incidence rendering works', { mockery::expect_args( renderer$render_mock(), 5, - 'n_inc_clinical_730_3650', - 2, + 'n_age_1825_7299', + 15, timestep ) mockery::expect_args( renderer$render_mock(), 6, - 'p_inc_clinical_730_3650', - .6, + 'n_age_3650_36499', + 90, timestep ) }) diff --git a/vignettes/Vaccines.Rmd b/vignettes/Vaccines.Rmd index 02eb9358..0d1d91ad 100644 --- a/vignettes/Vaccines.Rmd +++ b/vignettes/Vaccines.Rmd @@ -388,7 +388,7 @@ plot_prev() ```{r, fig.align = 'center', out.width='100%', message=FALSE} # Plot doses output$month <- ceiling(output$timestep / month) -doses <- output[, c(2:6, 38)] +doses <- output[, c(grep(pattern = "pev", names(output), value = T), "month")] doses <- aggregate(cbind(doses[1:5]), by = list(doses$month), FUN = sum) From 4820186dc35c6a4f3580a3a17203bb9fc601a4b3 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 23 May 2024 11:23:28 +0100 Subject: [PATCH 06/35] Removed "age" from age outputs in tests and vignettes. Corrected dt variables recovery rate assignment when antimalarial resistance is modelled. --- R/disease_progression.R | 6 +++--- tests/testthat/test-render.R | 12 ++++++------ vignettes/Demography.Rmd | 2 +- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/R/disease_progression.R b/R/disease_progression.R index d72d2c75..55103ac1 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -7,14 +7,14 @@ #' @param recovery_outcome competing hazards object for recovery rates #' @noRd calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt_input){ - + # Get correct input for dt depending on spc if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ dt_v <- dt_input } else if (isTRUE(parameters$antimalarial_resistance)) { - dt_v <- dt_input$and(variables$state$get_index_of("Tr")$to_vector())$to_vector() + dt_v <- dt_input$get_values(variables$state$get_index_of("Tr")) } - + recovery_rates <- numeric(length = parameters$human_population) recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R index 8b7a177b..d6d4787a 100644 --- a/tests/testthat/test-render.R +++ b/tests/testthat/test-render.R @@ -162,7 +162,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 1, - 'n_age_0_729', + 'n_0_729', 2, timestep ) @@ -170,7 +170,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 2, - 'n_age_365_1824', + 'n_365_1824', 4, timestep ) @@ -178,7 +178,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 3, - 'n_age_1095_3649', + 'n_1095_3649', 7, timestep ) @@ -186,7 +186,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 4, - 'n_age_1825_2919', + 'n_1825_2919', 3, timestep ) @@ -194,7 +194,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 5, - 'n_age_1825_7299', + 'n_1825_7299', 15, timestep ) @@ -202,7 +202,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 6, - 'n_age_3650_36499', + 'n_3650_36499', 90, timestep ) diff --git a/vignettes/Demography.Rmd b/vignettes/Demography.Rmd index d8c96a93..8b53421f 100644 --- a/vignettes/Demography.Rmd +++ b/vignettes/Demography.Rmd @@ -124,7 +124,7 @@ convert_to_long <- function(age_min, age_max, output) { data.frame( age_lower = age_min[[i]], age_upper = age_max[[i]], - n = output[,paste0('n_age_', age_min[[i]], '_',age_max[[i]])], + n = output[,paste0('n_', age_min[[i]], '_',age_max[[i]])], age_plot = age_min[[i]]/365, run = output$run, timestep = output$timestep) From 62791a2c9d4ebd4100972550a5b2c81009038742 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Fri, 24 May 2024 10:32:45 +0100 Subject: [PATCH 07/35] Removed missed merge commit line. --- vignettes/Carrying-capacity.Rmd | 1 - 1 file changed, 1 deletion(-) diff --git a/vignettes/Carrying-capacity.Rmd b/vignettes/Carrying-capacity.Rmd index ed8d7738..e17ac25d 100644 --- a/vignettes/Carrying-capacity.Rmd +++ b/vignettes/Carrying-capacity.Rmd @@ -107,7 +107,6 @@ p_lsm <- p |> set.seed(123) s_lsm <- run_simulation(timesteps = timesteps, parameters = p_lsm) s_lsm$pfpr <- s_lsm$n_detect_lm_730_3650 / s_lsm$n_age_730_3650 -s_lsm$pfpr <- s_lsm$n_detect_lm_730_3650 / s_lsm$n_730_3650 par(mfrow = c(1, 2)) plot(s$EIR_gamb ~ s$timestep, t = "l", ylim = c(0, 150), xlab = "Time", ylab = "EIR") From 1b305628c8aeafb073521a2f7675c814a75a1e95 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Fri, 24 May 2024 11:00:26 +0100 Subject: [PATCH 08/35] Correcting rendering tests (merging error). --- tests/testthat/test-render.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R index 0cd2d0bc..8184c66f 100644 --- a/tests/testthat/test-render.R +++ b/tests/testthat/test-render.R @@ -41,7 +41,7 @@ test_that('that default rendering works', { mockery::expect_args( renderer$render_mock(), - 4, + 3, 'n_detect_pcr_730_3650', 3, timestep @@ -75,7 +75,7 @@ test_that('that default rendering works when no one is in the age range', { mockery::expect_args( renderer$render_mock(), 1, - 'n_detect_730_3650', + 'n_detect_lm_730_3650', 0, timestep ) From 352fbae2b6aaaf506d5fd9013a10926ed7da813e Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Mon, 24 Jun 2024 16:31:00 +0100 Subject: [PATCH 09/35] biting_process renamed as create_biting_process, where functionality now occurs in biting_process.R. Similarly, create_recovery_rates_process functionality is now in disease_progression.R stash_source_humans in competing hazard resolution is now redundant and has been removed from infection processes. recovery rates are now assigned to a variable which is updated as infections, recoveries (via update_infection) and deaths occur. This also takes into account recovery time of treated due to slow parasite clearance and replaced the dt variable. I've also noticed (I think) that Tr->S and U->S recoveries can be done in a single step, so have attempted to combine these at risk of making the code more complicated! Added a metapopulation incidence rendering wrapper function for initialising multiple renderers. Rewrote create_age_group_renderer for clarity following feedback. Some tests have been adjusted to reflect these changes. --- R/biting_process.R | 61 ++++++++-------- R/competing_hazards.R | 6 +- R/disease_progression.R | 77 +++++++++++---------- R/human_infection.R | 28 ++++---- R/model.R | 1 + R/mortality_processes.R | 6 +- R/processes.R | 61 ++++++---------- R/render.R | 30 +++++--- R/utils.R | 2 +- R/variables.R | 20 +++--- tests/testthat/test-biting-integration.R | 26 +++---- tests/testthat/test-infection-integration.R | 50 ++++++------- 12 files changed, 181 insertions(+), 187 deletions(-) diff --git a/R/biting_process.R b/R/biting_process.R index 8ae1147f..594b9429 100644 --- a/R/biting_process.R +++ b/R/biting_process.R @@ -18,7 +18,7 @@ #' @param infection_outcome competing hazards object for infection rates #' @param timestep the current timestep #' @noRd -biting_process <- function( +create_biting_process <- function( renderer, solvers, models, @@ -29,37 +29,36 @@ biting_process <- function( lagged_eir, mixing = 1, mixing_index = 1, - infection_outcome, - timestep -) { - - age <- get_age(variables$birth$get_values(), timestep) - - bitten_humans <- simulate_bites( - renderer, - solvers, - models, - variables, - events, - age, - parameters, - timestep, - lagged_infectivity, - lagged_eir, - mixing, - mixing_index - ) - - simulate_infection( - variables, - events, - bitten_humans, - age, - parameters, - timestep, - renderer, infection_outcome - ) +) { + function(timestep){ + age <- get_age(variables$birth$get_values(), timestep) + bitten_humans <- simulate_bites( + renderer, + solvers, + models, + variables, + events, + age, + parameters, + timestep, + lagged_infectivity, + lagged_eir, + mixing, + mixing_index + ) + + simulate_infection( + variables, + events, + bitten_humans, + age, + parameters, + timestep, + renderer, + infection_outcome + ) + } } #' @importFrom stats rpois diff --git a/R/competing_hazards.R b/R/competing_hazards.R index 4ad9b8f6..cdc6aea0 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -18,15 +18,11 @@ CompetingOutcome <- R6::R6Class( set_rates = function(rates){ self$rates <- rates }, - stash_source_humans = function(source_humans){ - self$source_humans <- source_humans - }, execute = function(t, target){ private$targeted_process(t, target) self$rates <- rep(0, length(self$rates)) }, - rates = NULL, - source_humans = NULL + rates = NULL ) ) diff --git a/R/disease_progression.R b/R/disease_progression.R index 78ac0f7a..f5cef8c5 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -1,28 +1,20 @@ #' @title Calculate recovery rates #' @description Calculates recovery rates for each individual in the population -#' for storage in competing hazards object and future resolution +#' for storage in competing hazards object and subsequent resolution #' #' @param variables the available human variables -#' @param parameters model parameters #' @param recovery_outcome competing hazards object for recovery rates #' @noRd -calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt_input){ - - # Get correct input for dt depending on antimalarial resistance - if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ - dt_v <- dt_input - } else if (isTRUE(parameters$antimalarial_resistance)) { - dt_v <- dt_input$get_values(variables$state$get_index_of("Tr")) +create_recovery_rates_process <- function( + variables, + recovery_outcome +) { + function(timestep){ + recovery_outcome$set_rates(variables$recovery_rates$get_values()) } - - recovery_rates <- numeric(length = parameters$human_population) - recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd - recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da - recovery_rates[variables$state$get_index_of("U")$to_vector()] <- 1/parameters$du - recovery_rates[variables$state$get_index_of("Tr")$to_vector()] <- 1/dt_v - recovery_outcome$set_rates(recovery_rates) } + #' @title Disease progression outcomes (recovery) #' @description Following resolution of competing hazards, update state and #' infectivity of sampled individuals @@ -33,7 +25,7 @@ calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt #' @param parameters model parameters #' @param renderer competing hazards object for recovery rates #' @noRd -recovery_process_resolved_hazard <- function( +recovery_outcome_process <- function( timestep, target, variables, @@ -53,25 +45,31 @@ recovery_process_resolved_hazard <- function( "U", variables$infectivity, parameters$cu, + variables$recovery_rates, + 1/parameters$du, variables$state$get_index_of("A")$and(target) ) + # Is there a reason we aren't doing this in one step...? These all recover to S... update_infection( variables$state, "S", variables$infectivity, 0, - variables$state$get_index_of("U")$and(target) - ) - - update_infection( - variables$state, - "S", - variables$infectivity, + variables$recovery_rates, 0, - variables$state$get_index_of("Tr")$and(target) + variables$state$get_index_of(c("U","Tr"))$and(target) ) + # update_infection( + # variables$state, + # "S", + # variables$infectivity, + # 0, + # variables$recovery_rates, + # 0, + # variables$state$get_index_of("Tr")$and(target) + # ) } #' @title Update the state of an individual as infection events occur @@ -84,14 +82,17 @@ recovery_process_resolved_hazard <- function( #' @param new_infectivity the new infectivity of the progressed individuals #' @noRd update_infection <- function( - state, - to_state, - infectivity, - new_infectivity, - to_move - ) { + state, + to_state, + infectivity, + new_infectivity, + recovery_rates, + new_recovery_rate, + to_move +) { state$queue_update(to_state, to_move) infectivity$queue_update(new_infectivity, to_move) + recovery_rates$queue_update(new_recovery_rate, to_move) } #' @title Modelling the progression to asymptomatic disease @@ -102,11 +103,11 @@ update_infection <- function( #' @param parameters model parameters #' @noRd update_to_asymptomatic_infection <- function( - variables, - parameters, - timestep, - to_move - ) { + variables, + parameters, + timestep, + to_move +) { if (to_move$size() > 0) { variables$state$queue_update('A', to_move) new_infectivity <- asymptomatic_infectivity( @@ -121,5 +122,9 @@ update_to_asymptomatic_infection <- function( new_infectivity, to_move ) + variables$recovery_rates$queue_update( + 1/parameters$da, + to_move + ) } } diff --git a/R/human_infection.R b/R/human_infection.R index 50568a26..d8a39d90 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -1,7 +1,7 @@ #' @title Simulate malaria infection in humans #' @description #' This function ends with the assignment of rates of infection to the competing -#' hazard resolution object. Boosts immunity given infectious bites. +#' hazard resolution object and boosts immunity given infectious bites. #' @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 @@ -60,7 +60,6 @@ calculate_infections <- function( ) { source_humans <- variables$state$get_index_of( c('S', 'A', 'U'))$and(bitten_humans) - infection_outcome$stash_source_humans(source_humans) b <- blood_immunity(variables$ib$get_values(source_humans), parameters) @@ -113,7 +112,7 @@ calculate_infections <- function( incidence_probability_renderer( variables$birth, renderer, - source_pop, + source_humans, prob, "inc_", parameters$incidence_min_ages, @@ -122,7 +121,7 @@ calculate_infections <- function( ) ## capture infection rates to resolve in competing hazards - infection_rates <- numeric(length = parameters$human_population) + infection_rates <- rep(0, length = parameters$human_population) infection_rates[source_vector] <- prob_to_rate(prob) infection_outcome$set_rates(infection_rates) } @@ -133,16 +132,14 @@ calculate_infections <- function( #' and treated malaria; and resulting boosts in immunity #' @param timestep current timestep #' @param infected_humans bitset of infected humans -#' @param source_humans bitset of humans with infection potential #' @param variables a list of all of the model variables #' @param renderer model render object #' @param parameters model parameters #' @param prob vector of population probabilities of infection #' @noRd -infection_process_resolved_hazard <- function( +infection_outcome_process <- function( timestep, infected_humans, - source_humans, variables, renderer, parameters, @@ -410,14 +407,19 @@ calculate_treated <- function( successfully_treated ) if(parameters$antimalarial_resistance) { - variables$dt$queue_update( - parameters$dt, + variables$recovery_rates$queue_update( + 1/parameters$dt, non_slow_parasite_clearance_individuals ) - variables$dt$queue_update( - dt_slow_parasite_clearance, + variables$recovery_rates$queue_update( + 1/dt_slow_parasite_clearance, slow_parasite_clearance_individuals ) + } else { + variables$recovery_rates$queue_update( + 1/parameters$dt, + successfully_treated + ) } } @@ -455,6 +457,8 @@ schedule_infections <- function( 'D', variables$infectivity, parameters$cd, + variables$recovery_rates, + 1/parameters$dd, to_infect ) } @@ -523,7 +527,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/model.R b/R/model.R index de5d001a..21df1173 100644 --- a/R/model.R +++ b/R/model.R @@ -220,6 +220,7 @@ run_metapop_simulation <- function( variables <- lapply(parameters, create_variables) events <- lapply(parameters, create_events) renderer <- lapply(parameters, function(.) individual::Render$new(timesteps)) + populate_metapopulation_incidence_rendering_columns(renderer, parameters) for (i in seq_along(parameters)) { # NOTE: forceAndCall is necessary here to make sure i refers to the current # iteration diff --git a/R/mortality_processes.R b/R/mortality_processes.R index 37792f6c..e67baf0e 100644 --- a/R/mortality_processes.R +++ b/R/mortality_processes.R @@ -113,12 +113,8 @@ reset_target <- function(variables, events, target, state, parameters, timestep) # onwards infectiousness variables$infectivity$queue_update(0, target) + variables$recovery_rates$queue_update(0, target) - # treated compartment residence time: - if(!is.null(variables$dt)) { - variables$dt$queue_update(parameters$dt, target) - } - # zeta and zeta group and vector controls survive rebirth } } diff --git a/R/processes.R b/R/processes.R index ac664a8c..583a9d09 100644 --- a/R/processes.R +++ b/R/processes.R @@ -68,31 +68,29 @@ create_processes <- function( infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, infection_outcome$source_humans, - variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + infection_outcome_process(timestep, target, + variables, renderer, parameters, + prob = rate_to_prob(infection_outcome$rates)) }, size = parameters$human_population ) recovery_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - recovery_process_resolved_hazard(timestep, target, variables, parameters, renderer) + recovery_outcome_process(timestep, target, variables, parameters, renderer) }, size = parameters$human_population ) - create_infection_recovery_hazard_process <- CompetingHazard$new( - outcomes = list(infection_outcome, recovery_outcome) - ) - # ============================== # Biting and mortality processes # ============================== - # simulate infections for humans and set last_boosted_* + # simulate bites, calculates infection rates for bitten humans and set last_boosted_* # move mosquitoes into incubating state # kill mosquitoes caught in vector control - create_infection_rates_process <- function(timestep){ - biting_process( + processes <- c( + processes, + create_biting_process( renderer, solvers, models, @@ -103,44 +101,27 @@ create_processes <- function( lagged_eir, mixing, mixing_index, - infection_outcome, - timestep + infection_outcome ) - } - - # ======================= - # Antimalarial Resistance - # ======================= - # Add an a new process which governs the transition from Tr to S when - # antimalarial resistance is simulated. The rate of transition switches - # from a parameter to a variable when antimalarial resistance == TRUE. - - # Assign the dt input to a separate object with the default single parameter value: - dt_input <- parameters$dt - - # If antimalarial resistance is switched on, assign dt variable values to the - if(parameters$antimalarial_resistance) { - dt_input <- variables$dt - } + ) # =================== # Disease Progression # =================== - create_recovery_rates_process <- function(timestep){ - calculate_recovery_rates( + processes <- c( + processes, + create_recovery_rates_process( variables, - parameters, - recovery_outcome, - dt_input - ) - } + recovery_outcome + ), + + # Resolve competing hazards of infection with disease progression + CompetingHazard$new( + outcomes = list(infection_outcome, recovery_outcome) + )$resolve + ) - processes <- c(processes, - create_infection_rates_process, - create_recovery_rates_process, - create_infection_recovery_hazard_process$resolve) - # =============== # ODE integration # =============== diff --git a/R/render.R b/R/render.R index 0d1c252e..4b6335eb 100644 --- a/R/render.R +++ b/R/render.R @@ -188,18 +188,20 @@ create_age_group_renderer <- function( parameters, renderer ) { + + age_ranges <- rbind( + cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages), + cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages), + cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages), + cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages), + cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages) + ) + + unique_age_combinations <- as.data.frame(unique(age_ranges)) + ordered_unique_age_combinations <- unique_age_combinations[order(unique_age_combinations$V1, unique_age_combinations$V2), ] + function(timestep) { - unique_age_combinations <- as.data.frame(unique(rbind(cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages), - cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages), - cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages), - cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages), - cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages)))) - - ordered_unique_age_combinations <- unique_age_combinations[ - with(unique_age_combinations, order(V1, V2)), - ] - for (i in seq_along(ordered_unique_age_combinations$V1)) { lower <- ordered_unique_age_combinations$V1[[i]] upper <- ordered_unique_age_combinations$V2[[i]] @@ -262,4 +264,10 @@ render_initial_incidence <- function(renderer, lower_vals, upper_vals, inc_name) renderer$set_default(paste0('p_', inc_name, "_", lower_vals[i], '_', upper_vals[i]), 0) } } - \ No newline at end of file + + +populate_metapopulation_incidence_rendering_columns <- function(renderer, parameters){ + for(i in length(parameters)){ + populate_incidence_rendering_columns(renderer[[i]], parameters[[i]]) + } +} \ No newline at end of file diff --git a/R/utils.R b/R/utils.R index 34de912e..8ce50545 100644 --- a/R/utils.R +++ b/R/utils.R @@ -100,5 +100,5 @@ prob_to_rate <- function(prob){ #'@param rate rate #'@noRd rate_to_prob <- function(rate){ - 1- exp(-rate) + 1 - exp(-rate) } diff --git a/R/variables.R b/R/variables.R index 092a6431..b044392b 100644 --- a/R/variables.R +++ b/R/variables.R @@ -178,6 +178,7 @@ create_variables <- function(parameters) { diseased <- state$get_index_of('D')$to_vector() asymptomatic <- state$get_index_of('A')$to_vector() subpatent <- state$get_index_of('U')$to_vector() + treated <- state$get_index_of('Tr')$to_vector() # Set the initial infectivity values for each individual infectivity_values[diseased] <- parameters$cd @@ -191,6 +192,16 @@ create_variables <- function(parameters) { # Initialise the infectivity variable infectivity <- individual::DoubleVariable$new(infectivity_values) + # Set recovery rate for each individual + recovery_values <- rep(0, get_human_population(parameters, 0)) + recovery_values[diseased] <- 1/parameters$dd + recovery_values[asymptomatic] <- 1/parameters$da + recovery_values[subpatent] <- 1/parameters$du + recovery_values[treated] <- 1/parameters$dt + + # Initialise the recovery rate variable + recovery_rates <- individual::DoubleVariable$new(recovery_values) + drug <- individual::IntegerVariable$new(rep(0, size)) drug_time <- individual::IntegerVariable$new(rep(-1, size)) @@ -220,6 +231,7 @@ create_variables <- function(parameters) { zeta = zeta, zeta_group = zeta_group, infectivity = infectivity, + recovery_rates = recovery_rates, drug = drug, drug_time = drug_time, last_pev_timestep = last_pev_timestep, @@ -230,14 +242,6 @@ create_variables <- function(parameters) { spray_time = spray_time ) - # Add variables for antimalarial resistance state residency times (dt) - if(parameters$antimalarial_resistance) { - dt <- individual::DoubleVariable$new(rep(parameters$dt, size)) - variables <- c( - variables, - dt = dt - ) - } # Add variables for individual mosquitoes if (parameters$individual_mosquitoes) { diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R index 9712a243..4be80329 100644 --- a/tests/testthat/test-biting-integration.R +++ b/tests/testthat/test-biting-integration.R @@ -21,15 +21,7 @@ test_that('biting_process integrates mosquito effects and human infection', { size = parameters$human_population ) - bitten <- individual::Bitset$new(parameters$human_population) - bites_mock <- mockery::mock(bitten) - infection_mock <- mockery::mock() - - mockery::stub(biting_process, 'simulate_bites', bites_mock) - mockery::stub(biting_process, 'simulate_infection', infection_mock) - - - biting_process( + biting_process <- create_biting_process( renderer, solvers, models, @@ -38,11 +30,19 @@ test_that('biting_process integrates mosquito effects and human infection', { parameters, lagged_foim, lagged_eir, - mixing = 1, - mixing_index = 1, - infection_outcome, - timestep + 1, + 1, + infection_outcome ) + + bitten <- individual::Bitset$new(parameters$human_population) + bites_mock <- mockery::mock(bitten) + infection_mock <- mockery::mock() + + mockery::stub(biting_process, 'simulate_bites', bites_mock) + mockery::stub(biting_process, 'simulate_infection', infection_mock) + biting_process(timestep) + mockery::expect_args( bites_mock, 1, diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index 25637f0d..d399a527 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -23,7 +23,7 @@ test_that('simulate_infection integrates different types of infection and schedu infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + infection_outcome_process(timestep, target, variables, renderer, parameters) }, size = parameters$human_population ) @@ -97,17 +97,16 @@ test_that('simulate_infection integrates different types of infection and schedu schedule_mock <- mockery::mock() - mockery::stub(infection_process_resolved_hazard, 'incidence_renderer', mockery::mock()) - mockery::stub(infection_process_resolved_hazard, 'boost_immunity', boost_immunity_mock) - mockery::stub(infection_process_resolved_hazard, 'calculate_clinical_infections', clinical_infection_mock) - mockery::stub(infection_process_resolved_hazard, 'update_severe_disease', severe_infection_mock) - mockery::stub(infection_process_resolved_hazard, 'calculate_treated', treated_mock) - mockery::stub(infection_process_resolved_hazard, 'schedule_infections', schedule_mock) + 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_process_resolved_hazard( + infection_outcome_process( timestep, infected, - source_humans, variables, renderer, parameters, @@ -221,7 +220,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + infection_outcome_process(timestep, target, variables, renderer, parameters) }, size = parameters$human_population ) @@ -321,6 +320,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat variables <- list( state = list(queue_update = mockery::mock()), infectivity = list(queue_update = mockery::mock()), + recovery_rates = list(queue_update = mockery::mock()), drug = list(queue_update = mockery::mock()), drug_time = list(queue_update = mockery::mock()) ) @@ -410,9 +410,9 @@ test_that('calculate_treated correctly samples treated and updates the drug stat variables <- list( state = list(queue_update = mockery::mock()), infectivity = list(queue_update = mockery::mock()), + recovery_rates = list(queue_update = mockery::mock()), drug = list(queue_update = mockery::mock()), - drug_time = list(queue_update = mockery::mock()), - dt = list(queue_update = mockery::mock()) + drug_time = list(queue_update = mockery::mock()) ) renderer <- individual::Render$new(timesteps = timestep) @@ -487,8 +487,8 @@ test_that('calculate_treated correctly samples treated and updates the drug stat ) expect_bitset_update(variables$drug$queue_update, c(2, 1, 1, 1, 2, 2, 2), c(1, 2, 3, 4, 5, 6, 7)) expect_bitset_update(variables$drug_time$queue_update, 5, c(1, 2, 3, 4, 5, 6, 7)) - expect_bitset_update(variables$dt$queue_update, 5, c(2, 3, 4, 5, 6, 7), 1) - expect_bitset_update(variables$dt$queue_update, 15, c(1), 2) + expect_bitset_update(variables$recovery_rates$queue_update, 1/5, c(2, 3, 4, 5, 6, 7), 1) + expect_bitset_update(variables$recovery_rates$queue_update, 1/15, c(1), 2) }) test_that('calculate_treated correctly samples treated and updates the drug state when resistance not set for all drugs', { @@ -523,9 +523,9 @@ test_that('calculate_treated correctly samples treated and updates the drug stat variables <- list( state = list(queue_update = mockery::mock()), infectivity = list(queue_update = mockery::mock()), + recovery_rates = list(queue_update = mockery::mock()), drug = list(queue_update = mockery::mock()), - drug_time = list(queue_update = mockery::mock()), - dt = list(queue_update = mockery::mock()) + drug_time = list(queue_update = mockery::mock()) ) # Create a Bitset of individuals seeking treatment individuals: @@ -629,15 +629,15 @@ test_that('calculate_treated correctly samples treated and updates the drug stat # Check that update queued for dt for the non-slow parasite clearance individuals is correct: expect_bitset_update( - variables$dt$queue_update, - parameters$dt, + variables$recovery_rates$queue_update, + 1/parameters$dt, c(2, 3, 4, 5, 6, 7), 1) # Check that update queued for dt for the slow parasite clearance individuals is correct: expect_bitset_update( - variables$dt$queue_update, - unlist(parameters$dt_slow_parasite_clearance), + variables$recovery_rates$queue_update, + 1/unlist(parameters$dt_slow_parasite_clearance), c(1), 2) @@ -665,8 +665,8 @@ test_that('schedule_infections correctly schedules new infections', { parameters, 42 ) - - actual_infected <- mockery::mock_args(infection_mock)[[1]][[5]]$to_vector() + + actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector() actual_asymp_infected <- mockery::mock_args(asymp_mock)[[1]][[4]]$to_vector() expect_equal( @@ -705,7 +705,7 @@ test_that('prophylaxis is considered for medicated humans', { infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + infection_outcome_process(timestep, target, variables, renderer, parameters) }, size = parameters$human_population ) @@ -1053,9 +1053,9 @@ test_that('calculate_treated() successfully adds additional resistance columns t variables <- list( state = list(queue_update = mockery::mock()), infectivity = list(queue_update = mockery::mock()), + recovery_rates = list(queue_update = mockery::mock()), drug = list(queue_update = mockery::mock()), - drug_time = list(queue_update = mockery::mock()), - dt = list(queue_update = mockery::mock()) + drug_time = list(queue_update = mockery::mock()) ) renderer <- individual::Render$new(timesteps = 10) From 8866b5a7bb9e7dc7e44cf0a437f6384a3275f111 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 9 May 2024 12:16:02 +0100 Subject: [PATCH 10/35] competing_hazards.R: Giovanni designed two objects: CompetingOutcome (containing an outcome following resolution and rates that that outcome occurs for each individual in the population) and CompetingHazard (contains all CompetingOutcomes and process function to resolve hazards). processes.R: create_biting_process now replaced with create_infection_rates_process and biting_process, that uses essentially the same functional pathway, but which stops when the rates of infection are generated, storing them in the infection_outcome object within human_infection.R Similarly, disease_progression_process(es) are replaced with new functions: create_recovery_rates_process and calculate_recovery_rates, that get the recovery rates for each individual based on disease state. We then add these processes to the process list, including the resolution function create_infection_recovery_hazard_process$resolve. This process (create_infection_recovery_hazard_process) contains infection_outcome and recovery_outcome. infection_outcome calls infection_process_resolved_hazard, which follows the rest of the infection pathway functions: assigning new infections to different human states. Note that this requires probability of infection which gets used in the incidence renderer function. recovery_outcome also updates human states and infectivities. I've adjusted some of the tests to reflect the changes in functions required to generate results. Removed white space deletions. --- R/biting_process.R | 67 +++--- R/competing_hazards.R | 78 ++++++ R/disease_progression.R | 125 ++++++---- R/human_infection.R | 166 +++++++------ R/processes.R | 88 ++++--- R/utils.R | 14 ++ tests/testthat/test-biology.R | 2 +- tests/testthat/test-biting-integration.R | 35 ++- tests/testthat/test-competing-hazards.R | 248 ++++++++++++++++++++ tests/testthat/test-infection-integration.R | 119 ++++++++-- tests/testthat/test-pev.R | 14 +- 11 files changed, 725 insertions(+), 231 deletions(-) create mode 100644 R/competing_hazards.R create mode 100644 tests/testthat/test-competing-hazards.R diff --git a/R/biting_process.R b/R/biting_process.R index 279802ea..8ae1147f 100644 --- a/R/biting_process.R +++ b/R/biting_process.R @@ -15,8 +15,10 @@ #' values (default: 1) #' @param mixing_index an index for this population's position in the #' lagged_infectivity list (default: 1) +#' @param infection_outcome competing hazards object for infection rates +#' @param timestep the current timestep #' @noRd -create_biting_process <- function( +biting_process <- function( renderer, solvers, models, @@ -26,37 +28,38 @@ create_biting_process <- function( lagged_infectivity, lagged_eir, mixing = 1, - mixing_index = 1 + mixing_index = 1, + infection_outcome, + timestep ) { - function(timestep) { - # Calculate combined EIR - age <- get_age(variables$birth$get_values(), timestep) - - bitten_humans <- simulate_bites( - renderer, - solvers, - models, - variables, - events, - age, - parameters, - timestep, - lagged_infectivity, - lagged_eir, - mixing, - mixing_index - ) - - simulate_infection( - variables, - events, - bitten_humans, - age, - parameters, - timestep, - renderer - ) - } + + age <- get_age(variables$birth$get_values(), timestep) + + bitten_humans <- simulate_bites( + renderer, + solvers, + models, + variables, + events, + age, + parameters, + timestep, + lagged_infectivity, + lagged_eir, + mixing, + mixing_index + ) + + simulate_infection( + variables, + events, + bitten_humans, + age, + parameters, + timestep, + renderer, + infection_outcome + ) } #' @importFrom stats rpois @@ -74,6 +77,7 @@ simulate_bites <- function( mixing = 1, mixing_index = 1 ) { + bitten_humans <- individual::Bitset$new(parameters$human_population) human_infectivity <- variables$infectivity$get_values() @@ -244,7 +248,6 @@ calculate_infectious_individual <- function( species_index, parameters ) { - infectious_index$copy()$and(species_index)$size() } diff --git a/R/competing_hazards.R b/R/competing_hazards.R new file mode 100644 index 00000000..cdc6aea0 --- /dev/null +++ b/R/competing_hazards.R @@ -0,0 +1,78 @@ +## Define classes to resolve competing hazards +CompetingOutcome <- R6::R6Class( + "CompetingOutcome", + private = list( + targeted_process = NULL + ), + public = list( + initialize = function(targeted_process, size){ + if (!is.function(targeted_process)){ + stop("targeted_process must be a function") + } + if (!is.numeric(size) || size <= 0){ + stop("size must be positive integer") + } + private$targeted_process <- targeted_process + self$rates <- rep(0, size) + }, + set_rates = function(rates){ + self$rates <- rates + }, + execute = function(t, target){ + private$targeted_process(t, target) + self$rates <- rep(0, length(self$rates)) + }, + rates = NULL + ) +) + +CompetingHazard <- R6::R6Class( + "CompetingHazard", + private = list( + outcomes = list(), + size = NULL, + # RNG is passed in because mockery is not able to stub runif + # TODO: change when fixed + rng = NULL + ), + public = list( + initialize = function(outcomes, rng = runif){ + if (length(outcomes) == 0){ + stop("At least one outcome must be provided") + } + if (!all(sapply(outcomes, function(x) inherits(x, "CompetingOutcome")))){ + stop("All outcomes must be of class CompetingOutcome") + } + private$outcomes <- outcomes + private$size <- length(outcomes[[1]]$rates) + private$rng <- rng + }, + resolve = function(t){ + event_rates <- do.call( + 'cbind', + lapply(private$outcomes, function(x) x$rates) + ) + occur_rates <- rowSums(event_rates) + occur_rng <- private$rng(private$size) + occurs <- occur_rng < rate_to_prob(occur_rates) + norm_probs <- event_rates / occur_rates + norm_probs[is.na(norm_probs)] <- 0 + cum_probs <- apply(norm_probs, 1, cumsum) + event_rng <- private$rng(private$size) + for(o in seq_along(private$outcomes)){ + if (o == 1) { + selected <- (event_rng <= cum_probs[o,]) + } else { + selected <- (event_rng > cum_probs[o - 1,]) & + (event_rng <= cum_probs[o,]) + } + target <- individual::Bitset$new(private$size)$insert( + which(selected & occurs) + ) + if (target$size() > 0){ + private$outcomes[[o]]$execute(t, target) + } + } + } + ) +) diff --git a/R/disease_progression.R b/R/disease_progression.R index d9ac11e4..d72d2c75 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -1,3 +1,79 @@ +#' @title Calculate recovery rates +#' @description Calculates recovery rates for each individual in the population +#' for storage in competing hazards object and future resolution +#' +#' @param variables the available human variables +#' @param parameters model parameters +#' @param recovery_outcome competing hazards object for recovery rates +#' @noRd +calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt_input){ + + # Get correct input for dt depending on spc + if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ + dt_v <- dt_input + } else if (isTRUE(parameters$antimalarial_resistance)) { + dt_v <- dt_input$and(variables$state$get_index_of("Tr")$to_vector())$to_vector() + } + + recovery_rates <- numeric(length = parameters$human_population) + recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd + recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da + recovery_rates[variables$state$get_index_of("U")$to_vector()] <- 1/parameters$du + recovery_rates[variables$state$get_index_of("Tr")$to_vector()] <- 1/dt_v + recovery_outcome$set_rates(recovery_rates) +} + +#' @title Disease progression outcomes (recovery) +#' @description Following resolution of competing hazards, update state and +#' infectivity of sampled individuals +#' +#' @param timestep the current timestep +#' @param target the sampled recovering individuals +#' @param variables the available human variables +#' @param parameters model parameters +#' @param renderer competing hazards object for recovery rates +#' @noRd +recovery_process_resolved_hazard <- function( + timestep, + target, + variables, + parameters, + renderer +){ + + update_to_asymptomatic_infection( + variables, + parameters, + timestep, + variables$state$get_index_of("D")$and(target) + ) + + update_infection( + variables$state, + "U", + variables$infectivity, + parameters$cu, + variables$state$get_index_of("A")$and(target) + ) + + update_infection( + variables$state, + "S", + variables$infectivity, + 0, + variables$state$get_index_of("U")$and(target) + ) + + update_infection( + variables$state, + "S", + variables$infectivity, + 0, + variables$state$get_index_of("Tr")$and(target) + ) + +} + #' @title Update the state of an individual as infection events occur #' @description Randomly moves individuals towards the later stages of disease #' and updates their infectivity @@ -18,55 +94,6 @@ update_infection <- function( infectivity$queue_update(new_infectivity, to_move) } -create_progression_process <- function( - state, - from_state, - to_state, - rate, - infectivity, - new_infectivity -) { - function(timestep) { - - # Retrieve the indices of all individuals in the to_move state: - index <- state$get_index_of(from_state) - - # If the length of rate is greater than 1 (when it's a variable): - if (inherits(rate, "DoubleVariable")) { - rate <- rate$get_values(index) - } - - # Sample the individuals to be moved into a new Bitset using the transition rate(s): - to_move <- index$sample(1/rate) - - # Update the infection status of those individuals who are moving: - update_infection( - state, - to_state, - infectivity, - new_infectivity, - to_move - ) - } -} - -create_asymptomatic_progression_process <- function( - state, - rate, - variables, - parameters - ) { - function(timestep) { - to_move <- state$get_index_of('D')$sample(1/rate) - update_to_asymptomatic_infection( - variables, - parameters, - timestep, - to_move - ) - } -} - #' @title Modelling the progression to asymptomatic disease #' @description Randomly moves individuals to asymptomatic disease and #' calculates the infectivity for their age and immunity diff --git a/R/human_infection.R b/R/human_infection.R index d62b16f0..df2e0b1b 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -1,13 +1,15 @@ #' @title Simulate malaria infection in humans #' @description -#' Updates human states and variables to represent asymptomatic/clinical/severe -#' and treated malaria; and resulting boosts in immunity +#' This function ends with the assignment of rates of infection to the competing +#' hazard resolution object. Boosts immunity given infectious bites. #' @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 age of each human (timesteps) #' @param parameters of the model #' @param timestep current timestep +#' @param renderer the model renderer object +#' @param infection_outcome competing hazards object for infection rates #' @noRd simulate_infection <- function( variables, @@ -16,7 +18,8 @@ simulate_infection <- function( age, parameters, timestep, - renderer + renderer, + infection_outcome ) { if (bitten_humans$size() > 0) { boost_immunity( @@ -29,70 +32,18 @@ simulate_infection <- function( } # Calculate Infected - infected_humans <- calculate_infections( + calculate_infections( variables, bitten_humans, parameters, renderer, - timestep - ) - - if (infected_humans$size() > 0) { - 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_infections <- calculate_clinical_infections( - variables, - infected_humans, - parameters, - renderer, - timestep - ) - - update_severe_disease( - timestep, - infected_humans, - variables, - parameters, - renderer - ) - - treated <- calculate_treated( - variables, - clinical_infections, - parameters, timestep, - renderer - ) - - renderer$render('n_infections', infected_humans$size(), timestep) - - schedule_infections( - variables, - clinical_infections, - treated, - infected_humans, - parameters, - timestep + infection_outcome ) } #' @title Calculate overall infections for bitten humans -#' @description -#' Sample infected humans given prophylaxis and vaccination +#' @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 @@ -100,12 +51,13 @@ simulate_infection <- function( #' @param timestep current timestep #' @noRd calculate_infections <- function( - variables, - bitten_humans, - parameters, - renderer, - timestep - ) { + variables, + bitten_humans, + parameters, + renderer, + timestep, + infection_outcome +) { source_humans <- variables$state$get_index_of( c('S', 'A', 'U'))$and(bitten_humans) @@ -155,21 +107,97 @@ calculate_infections <- function( } prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy) - infected <- bitset_at(source_humans, bernoulli_multi_p(prob)) + ## Capture infection rates to resolve in competing hazards + infection_rates <- numeric(length = parameters$human_population) + infection_rates[source_vector] <- prob_to_rate(prob) + infection_outcome$set_rates(infection_rates) +} + +#' @title Assigns 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 prob vector of population probabilities of infection +#' @noRd +infection_process_resolved_hazard <- function( + timestep, + infected_humans, + variables, + renderer, + parameters, + prob){ + + source_humans <- variables$state$get_index_of(values = c("S","A","U")) + incidence_renderer( variables$birth, renderer, - infected, + infected_humans, source_humans, - prob, + prob[source_humans$to_vector()], 'inc_', parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages, timestep ) - - infected + + if (infected_humans$size() > 0) { + 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_infections <- calculate_clinical_infections( + variables, + infected_humans, + parameters, + renderer, + timestep + ) + + update_severe_disease( + timestep, + infected_humans, + variables, + parameters, + renderer + ) + + treated <- calculate_treated( + variables, + clinical_infections, + parameters, + timestep, + renderer + ) + + renderer$render('n_infections', infected_humans$size(), timestep) + + schedule_infections( + variables, + clinical_infections, + treated, + infected_humans, + parameters, + timestep + ) } #' @title Calculate clinical infections diff --git a/R/processes.R b/R/processes.R index a8c20efb..fe73b174 100644 --- a/R/processes.R +++ b/R/processes.R @@ -61,16 +61,37 @@ create_processes <- function( ) ) } + + # ===================================================== + # Competing Hazard Outcomes (Infections and Recoveries) + # ===================================================== + + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + }, + size = parameters$human_population + ) + + recovery_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + recovery_process_resolved_hazard(timestep, target, variables, parameters, renderer) + }, + size = parameters$human_population + ) + + create_infection_recovery_hazard_process <- CompetingHazard$new( + outcomes = list(infection_outcome, recovery_outcome) + ) # ============================== # Biting and mortality processes # ============================== - # schedule infections for humans and set last_boosted_* + # simulate infections for humans and set last_boosted_* # move mosquitoes into incubating state # kill mosquitoes caught in vector control - processes <- c( - processes, - create_biting_process( + create_infection_rates_process <- function(timestep){ + biting_process( renderer, solvers, models, @@ -80,31 +101,11 @@ create_processes <- function( lagged_infectivity, lagged_eir, mixing, - mixing_index - ), - create_asymptomatic_progression_process( - variables$state, - parameters$dd, - variables, - parameters - ), - create_progression_process( - variables$state, - 'A', - 'U', - parameters$da, - variables$infectivity, - parameters$cu - ), - create_progression_process( - variables$state, - 'U', - 'S', - parameters$du, - variables$infectivity, - 0 + mixing_index, + infection_outcome, + timestep ) - ) + } # ======================= # Antimalarial Resistance @@ -121,18 +122,23 @@ create_processes <- function( dt_input <- variables$dt } - # Create the progression process for Tr --> S specifying dt_input as the rate: - processes <- c( - processes, - create_progression_process( - variables$state, - 'Tr', - 'S', - dt_input, - variables$infectivity, - 0 + # =================== + # Disease Progression + # =================== + + create_recovery_rates_process <- function(timestep){ + calculate_recovery_rates( + variables, + parameters, + recovery_outcome, + dt_input ) - ) + } + + processes <- c(processes, + create_infection_rates_process, + create_recovery_rates_process, + create_infection_recovery_hazard_process$resolve) # =============== # ODE integration @@ -264,7 +270,11 @@ create_processes <- function( ) } + # ====================== # Mortality step + # ====================== + # Mortality is not resolved as a competing hazard + processes <- c( processes, create_mortality_process(variables, events, renderer, parameters)) diff --git a/R/utils.R b/R/utils.R index 2262c939..34de912e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -88,3 +88,17 @@ RandomState <- R6::R6Class( } ) ) + +#'@title Convert probability to a rate +#'@param prob probability +#'@noRd +prob_to_rate <- function(prob){ + -log(1 - prob) +} + +#'@title Convert rate to a probability +#'@param rate rate +#'@noRd +rate_to_prob <- function(rate){ + 1- exp(-rate) +} diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R index e429a48f..77ee42f6 100644 --- a/tests/testthat/test-biology.R +++ b/tests/testthat/test-biology.R @@ -86,7 +86,7 @@ test_that('FOIM is consistent with equilibrium', { psi <- unique_biting_rate(age, parameters) zeta <- variables$zeta$get_values() .pi <- human_pi(psi, zeta) - calculate_foim(a, sum(.pi * variables$infectivity$get_values()), 1.) + calculate_foim(a, sum(.pi * variables$infectivity$get_values()), mixing = 1) } ) expect_equal( diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R index 7da79b66..9712a243 100644 --- a/tests/testthat/test-biting-integration.R +++ b/tests/testthat/test-biting-integration.R @@ -14,7 +14,22 @@ test_that('biting_process integrates mosquito effects and human infection', { models <- parameterise_mosquito_models(parameters, timestep) solvers <- parameterise_solvers(models, parameters) - biting_process <- create_biting_process( + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + + bitten <- individual::Bitset$new(parameters$human_population) + bites_mock <- mockery::mock(bitten) + infection_mock <- mockery::mock() + + mockery::stub(biting_process, 'simulate_bites', bites_mock) + mockery::stub(biting_process, 'simulate_infection', infection_mock) + + + biting_process( renderer, solvers, models, @@ -22,17 +37,12 @@ test_that('biting_process integrates mosquito effects and human infection', { events, parameters, lagged_foim, - lagged_eir + lagged_eir, + mixing = 1, + mixing_index = 1, + infection_outcome, + timestep ) - - bitten <- individual::Bitset$new(parameters$human_population) - bites_mock <- mockery::mock(bitten) - infection_mock <- mockery::mock() - - mockery::stub(biting_process, 'simulate_bites', bites_mock) - mockery::stub(biting_process, 'simulate_infection', infection_mock) - biting_process(timestep) - mockery::expect_args( bites_mock, 1, @@ -59,7 +69,8 @@ test_that('biting_process integrates mosquito effects and human infection', { age, parameters, timestep, - renderer + renderer, + infection_outcome ) }) diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R new file mode 100644 index 00000000..ef67e2d0 --- /dev/null +++ b/tests/testthat/test-competing-hazards.R @@ -0,0 +1,248 @@ +# ============================== +# Test the CompetingHazard class +# ============================== + +test_that('hazard resolves two normalised outcomes when all events occur', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock( + c(0, 0, 0, 0), # all events occur + c(.05, .3, .2, .5) # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4)) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6)) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 3)) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(2, 4)) + ) +}) + +test_that('hazard resolves two unnormalised outcomes when all events occur', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock( + c(0, 0, 0, 0), # all events occur + c(.05, .3, .2, .5) # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 3)) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(2, 4)) + ) +}) + +test_that('hazard resolves two outcomes when some events occur', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock( + c(0, 0, 1, 1), # some events occur + c(.05, .3, .2, .5), # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(1) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(2) + ) +}) + + +test_that('hazard resolves when no rates are set', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2) + ) + + hazard$resolve(0) + + mockery::expect_called( + outcome_1_process, + 0 + ) + mockery::expect_called( + outcome_2_process, + 0 + ) +}) + +test_that('hazard resolves when rates are set to one', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2) + ) + + outcome_1$set_rates(rep(100, size)) + hazard$resolve(0) + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 2, 3, 4)) + ) + mockery::expect_called( + outcome_2_process, + 0 + ) +}) + +test_that('hazard resolves three outcomes', { + + size <- 4 + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + outcome_3_process <- mockery::mock() + outcome_3 <- CompetingOutcome$new( + targeted_process = outcome_3_process, + size = size + ) + + hazard <- CompetingHazard$new( + outcomes = list(outcome_1, outcome_2, outcome_3), + rng = mockery::mock( + c(0, 0, 0, 0), # all events occur + c(.04, .3, .14, .6), # event_rng + ) + ) + + outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) + outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + outcome_3$set_rates(c(0.5, 0.5, 0.5, 0.5)) + + hazard$resolve(0) + + mockery::expect_args( + outcome_1_process, + 1, + 0, + individual::Bitset$new(size)$insert(c(1, 3)) + ) + mockery::expect_args( + outcome_2_process, + 1, + 0, + individual::Bitset$new(size)$insert(2) + ) + mockery::expect_args( + outcome_3_process, + 1, + 0, + individual::Bitset$new(size)$insert(4) + ) +}) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index 83d7def9..f6eacb80 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -20,22 +20,17 @@ test_that('simulate_infection integrates different types of infection and schedu boost_immunity_mock <- mockery::mock() infected <- individual::Bitset$new(population)$insert(c(1, 3, 5)) infection_mock <- mockery::mock(infected) - clinical <- individual::Bitset$new(population)$insert(c(1, 3)) - clinical_infection_mock <- mockery::mock(clinical) - severe <- individual::Bitset$new(population)$insert(1) - severe_infection_mock <- mockery::mock(severe) - treated <- individual::Bitset$new(population)$insert(3) - treated_mock <- mockery::mock(treated) - schedule_mock <- mockery::mock() + + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(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_clinical_infections', clinical_infection_mock) - mockery::stub(simulate_infection, 'update_severe_disease', severe_infection_mock) - mockery::stub(simulate_infection, 'calculate_treated', treated_mock) - mockery::stub(simulate_infection, 'schedule_infections', schedule_mock) - mockery::stub(simulate_infection, 'incidence_renderer', mockery::mock()) - mockery::stub(simulate_infection, 'clinical_incidence_renderer', mockery::mock()) + simulate_infection( variables, events, @@ -43,7 +38,8 @@ test_that('simulate_infection integrates different types of infection and schedu age, parameters, timestep, - renderer + renderer, + infection_outcome ) mockery::expect_args( @@ -63,7 +59,67 @@ test_that('simulate_infection integrates different types of infection and schedu bitten, parameters, renderer, - timestep + timestep, + infection_outcome + ) +}) + + +test_that('simulate_infection integrates different types of infection and scheduling', { + population <- 8 + timestep <- 5 + parameters <- get_parameters(list( + human_population = population + )) + events <- create_events(parameters) + renderer <- mock_render(timestep) + + age <- c(20, 24, 5, 39, 20, 24, 5, 39) * 365 + immunity <- c(.2, .3, .5, .9, .2, .3, .5, .9) + # asymptomatics <- mockery::mock() + variables <- list( + ib = individual::DoubleVariable$new(immunity), + id = individual::DoubleVariable$new(immunity), + # 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) + + infected <- individual::Bitset$new(population)$insert(c(1, 3, 5)) + clinical <- individual::Bitset$new(population)$insert(c(1, 3)) + clinical_infection_mock <- mockery::mock(clinical) + boost_immunity_mock <- mockery::mock() + severe <- individual::Bitset$new(population)$insert(1) + severe_infection_mock <- mockery::mock(severe) + treated <- individual::Bitset$new(population)$insert(3) + treated_mock <- mockery::mock(treated) + schedule_mock <- mockery::mock() + + + mockery::stub(infection_process_resolved_hazard, 'incidence_renderer', mockery::mock()) + mockery::stub(infection_process_resolved_hazard, 'boost_immunity', boost_immunity_mock) + mockery::stub(infection_process_resolved_hazard, 'calculate_clinical_infections', clinical_infection_mock) + mockery::stub(infection_process_resolved_hazard, 'update_severe_disease', severe_infection_mock) + mockery::stub(infection_process_resolved_hazard, 'calculate_treated', treated_mock) + mockery::stub(infection_process_resolved_hazard, 'schedule_infections', schedule_mock) + + infection_process_resolved_hazard( + timestep, + infected, + variables, + renderer, + parameters, + prob) + + + mockery::expect_args( + boost_immunity_mock, + 1, + variables$ica, + infected, + variables$last_boosted_ica, + 5, + parameters$uc ) mockery::expect_args( @@ -161,15 +217,23 @@ test_that('calculate_infections works various combinations of drug and vaccinati bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4)) + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + infections <- calculate_infections( variables, bitten_humans, parameters, mock_render(timestep), - timestep + timestep, + infection_outcome ) - expect_equal(infections$to_vector(), 3) + expect_equal(sum(infections!=0), 3) mockery::expect_args(immunity_mock, 1, c(.3, .5, .9), parameters) mockery::expect_args( @@ -198,12 +262,6 @@ test_that('calculate_infections works various combinations of drug and vaccinati c(rtss_profile$beta, rtss_booster_profile$beta), c(rtss_profile$alpha, rtss_booster_profile$alpha) ) - mockery::expect_args( - bernoulli_mock, - 1, - c(.2 * .8 * .8, .3 * .7, .4) - ) - }) test_that('calculate_clinical_infections correctly samples clinically infected', { @@ -642,16 +700,25 @@ test_that('prophylaxis is considered for medicated humans', { m <- mockery::mock(seq(3)) mockery::stub(calculate_infections, 'bernoulli_multi_p', m) - calculate_infections( + + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + + infection_rates <- calculate_infections( variables, bitten, parameters, mock_render(timestep), - timestep + timestep, + infection_outcome ) expect_equal( - mockery::mock_args(m)[[1]][[1]], + rate_to_prob(infection_rates[infection_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 09c4cf7d..9a5f1175 100644 --- a/tests/testthat/test-pev.R +++ b/tests/testthat/test-pev.R @@ -150,6 +150,13 @@ test_that('Infection considers pev efficacy', { rep(.2, 4) ) + infection_outcome <- CompetingOutcome$new( + targeted_process = function(timestep, target){ + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + }, + size = parameters$human_population + ) + # remove randomness from infection sampling bernoulli_mock <- mockery::mock(c(1, 2)) mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock) @@ -164,16 +171,17 @@ test_that('Infection considers pev efficacy', { depth = 4 ) - calculate_infections( + infection_rates <- calculate_infections( variables = variables, bitten_humans = individual::Bitset$new(4)$insert(seq(4)), parameters = parameters, renderer = mock_render(timestep), - timestep = timestep + timestep = timestep, + infection_outcome = infection_outcome ) expect_equal( - mockery::mock_args(bernoulli_mock)[[1]][[1]], + rate_to_prob(infection_rates[infection_rates!=0]), c(0.590, 0.590, 0.215, 0.244), tolerance=1e-3 ) From 7702e19dd3916b2491f56cd0dc8bae46b7d8d41a Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 9 May 2024 16:11:41 +0100 Subject: [PATCH 11/35] Added a stash source_human function to the competing hazards function for use in infection outcome. --- R/competing_hazards.R | 6 +++++- R/human_infection.R | 5 +++-- R/processes.R | 3 ++- tests/testthat/test-infection-integration.R | 2 ++ 4 files changed, 12 insertions(+), 4 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index cdc6aea0..4ad9b8f6 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -18,11 +18,15 @@ CompetingOutcome <- R6::R6Class( set_rates = function(rates){ self$rates <- rates }, + stash_source_humans = function(source_humans){ + self$source_humans <- source_humans + }, execute = function(t, target){ private$targeted_process(t, target) self$rates <- rep(0, length(self$rates)) }, - rates = NULL + rates = NULL, + source_humans = NULL ) ) diff --git a/R/human_infection.R b/R/human_infection.R index df2e0b1b..75534555 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -60,6 +60,7 @@ calculate_infections <- function( ) { source_humans <- variables$state$get_index_of( c('S', 'A', 'U'))$and(bitten_humans) + infection_outcome$stash_source_humans(source_humans) b <- blood_immunity(variables$ib$get_values(source_humans), parameters) @@ -120,6 +121,7 @@ calculate_infections <- function( #' and treated malaria; and resulting boosts in immunity #' @param timestep current timestep #' @param infected_humans bitset of infected humans +#' @param source_humans bitset of humans with infection potential #' @param variables a list of all of the model variables #' @param renderer model render object #' @param parameters model parameters @@ -128,13 +130,12 @@ calculate_infections <- function( infection_process_resolved_hazard <- function( timestep, infected_humans, + source_humans, variables, renderer, parameters, prob){ - source_humans <- variables$state$get_index_of(values = c("S","A","U")) - incidence_renderer( variables$birth, renderer, diff --git a/R/processes.R b/R/processes.R index fe73b174..59e0b8c1 100644 --- a/R/processes.R +++ b/R/processes.R @@ -68,7 +68,8 @@ create_processes <- function( infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + infection_process_resolved_hazard(timestep, target, infection_outcome$source_humans, + variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) }, size = parameters$human_population ) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index f6eacb80..eda048a1 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -85,6 +85,7 @@ test_that('simulate_infection integrates different types of infection and schedu ) 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)) clinical <- individual::Bitset$new(population)$insert(c(1, 3)) clinical_infection_mock <- mockery::mock(clinical) @@ -106,6 +107,7 @@ test_that('simulate_infection integrates different types of infection and schedu infection_process_resolved_hazard( timestep, infected, + source_humans, variables, renderer, parameters, From 479c26e85d4f4b158127fac1305011d7cfdc7b69 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 23 May 2024 10:12:38 +0100 Subject: [PATCH 12/35] To correctly render ages and incidence, given the competing hazards solution (any rendering in the infection outcome does not occur), I've split the incidence renderer into three functions: age rendering, n incidence rendering and p incidence rendering. This allows much more flexibility in when these are called in the code sequence. All age inputs are now consolidated and rendered in a single function. In addition, we intitate and populate these incidence columns with 0 when they are used. Rebasing --- R/human_infection.R | 102 ++++++++----- R/model.R | 1 + R/parameters.R | 2 +- R/render.R | 156 +++++++++++++++----- tests/testthat/test-infection-integration.R | 14 +- tests/testthat/test-render.R | 94 +++++++++--- 6 files changed, 274 insertions(+), 95 deletions(-) diff --git a/R/human_infection.R b/R/human_infection.R index 75534555..50568a26 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -12,15 +12,15 @@ #' @param infection_outcome competing hazards object for infection rates #' @noRd simulate_infection <- function( - variables, - events, - bitten_humans, - age, - parameters, - timestep, - renderer, - infection_outcome - ) { + variables, + events, + bitten_humans, + age, + parameters, + timestep, + renderer, + infection_outcome +) { if (bitten_humans$size() > 0) { boost_immunity( variables$ib, @@ -109,7 +109,19 @@ calculate_infections <- function( prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy) - ## Capture infection rates to resolve in competing hazards + ## probability of incidence must be rendered at each timestep + incidence_probability_renderer( + variables$birth, + renderer, + source_pop, + prob, + "inc_", + parameters$incidence_min_ages, + parameters$incidence_max_ages, + timestep + ) + + ## capture infection rates to resolve in competing hazards infection_rates <- numeric(length = parameters$human_population) infection_rates[source_vector] <- prob_to_rate(prob) infection_outcome$set_rates(infection_rates) @@ -140,8 +152,6 @@ infection_process_resolved_hazard <- function( variables$birth, renderer, infected_humans, - source_humans, - prob[source_humans$to_vector()], 'inc_', parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages, @@ -211,12 +221,12 @@ infection_process_resolved_hazard <- function( #' @param timestep current timestep #' @noRd calculate_clinical_infections <- function( - variables, - infections, - parameters, - renderer, - timestep - ) { + variables, + infections, + parameters, + renderer, + timestep +) { ica <- variables$ica$get_values(infections) icm <- variables$icm$get_values(infections) phi <- clinical_immunity(ica, icm, parameters) @@ -225,6 +235,14 @@ calculate_clinical_infections <- function( variables$birth, renderer, clinical_infections, + 'inc_clinical_', + parameters$clinical_incidence_rendering_min_ages, + parameters$clinical_incidence_rendering_max_ages, + timestep + ) + incidence_probability_renderer( + variables$birth, + renderer, infections, phi, 'inc_clinical_', @@ -245,12 +263,12 @@ calculate_clinical_infections <- function( #' @param renderer model outputs #' @noRd update_severe_disease <- function( - timestep, - infections, - variables, - parameters, - renderer - ) { + timestep, + infections, + variables, + parameters, + renderer +) { age <- get_age(variables$birth$get_values(infections), timestep) iva <- variables$iva$get_values(infections) ivm <- variables$ivm$get_values(infections) @@ -266,6 +284,14 @@ update_severe_disease <- function( variables$birth, renderer, severe_infections, + 'inc_severe_', + parameters$severe_incidence_rendering_min_ages, + parameters$severe_incidence_rendering_max_ages, + timestep + ) + incidence_probability_renderer( + variables$birth, + renderer, infections, theta, 'inc_severe_', @@ -409,13 +435,13 @@ calculate_treated <- function( #' @param parameters model parameters #' @noRd schedule_infections <- function( - variables, - clinical_infections, - treated, - infections, - parameters, - timestep - ) { + variables, + clinical_infections, + treated, + infections, + parameters, + timestep +) { included <- treated$not(TRUE) to_infect <- clinical_infections$and(included) @@ -447,12 +473,12 @@ schedule_infections <- function( # Utility functions # ================= boost_immunity <- function( - immunity_variable, - exposed_index, - last_boosted_variable, - timestep, - delay - ) { + immunity_variable, + exposed_index, + last_boosted_variable, + timestep, + delay +) { # record who can be boosted exposed_index_vector <- exposed_index$to_vector() last_boosted <- last_boosted_variable$get_values(exposed_index) @@ -497,7 +523,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/model.R b/R/model.R index d59132f5..de5d001a 100644 --- a/R/model.R +++ b/R/model.R @@ -119,6 +119,7 @@ run_resumable_simulation <- function( events <- create_events(parameters) initialise_events(events, variables, parameters) renderer <- individual::Render$new(timesteps) + populate_incidence_rendering_columns(renderer, parameters) attach_event_listeners( events, variables, diff --git a/R/parameters.R b/R/parameters.R index bb692873..90ad4842 100644 --- a/R/parameters.R +++ b/R/parameters.R @@ -438,7 +438,7 @@ get_parameters <- function(overrides = list()) { incidence_rendering_min_ages = numeric(0), incidence_rendering_max_ages = numeric(0), clinical_incidence_rendering_min_ages = numeric(0), - clinical_incidence_rendering_max_ages = 5 * 365, + clinical_incidence_rendering_max_ages = numeric(0), severe_prevalence_rendering_min_ages = numeric(0), severe_prevalence_rendering_max_ages = numeric(0), severe_incidence_rendering_min_ages = numeric(0), diff --git a/R/render.R b/R/render.R index 4f5816f3..0d1c252e 100644 --- a/R/render.R +++ b/R/render.R @@ -18,12 +18,12 @@ in_age_range <- function(birth, timestep, lower, upper) { #' #' @noRd create_prevelance_renderer <- function( - state, - birth, - immunity, - parameters, - renderer - ) { + state, + birth, + immunity, + parameters, + renderer +) { function(timestep) { asymptomatic <- state$get_index_of('A') prob <- probability_of_detection( @@ -41,11 +41,6 @@ create_prevelance_renderer <- function( lower <- parameters$prevalence_rendering_min_ages[[i]] upper <- parameters$prevalence_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_detect_lm_', lower, '_', upper), in_age$copy()$and(detected)$size(), @@ -69,13 +64,11 @@ create_prevelance_renderer <- function( #' @title Render incidence statistics #' -#' @description renders incidence (new for this timestep) for indivduals +#' @description renders incidence (new for this timestep) for individuals #' #' @param birth variable for birth of the individual #' @param renderer object for model outputs #' @param target incidence population -#' @param source_pop the population which is sampled for infection -#' @param prob probability of infection #' @param prefix for model outputs #' @param lowers age bounds #' @param uppers age bounds @@ -83,27 +76,54 @@ create_prevelance_renderer <- function( #' #' @noRd incidence_renderer <- function( - birth, - renderer, - target, - source_pop, - prob, - prefix, - lowers, - uppers, - timestep - ) { + birth, + renderer, + target, + prefix, + lowers, + uppers, + timestep +) { for (i in seq_along(lowers)) { lower <- lowers[[i]] upper <- uppers[[i]] in_age <- in_age_range(birth, timestep, lower, upper) - renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep) renderer$render( paste0('n_', prefix, lower, '_', upper), in_age$copy()$and(target)$size(), timestep ) + } +} +#' @title Render probability of incidence statistics +#' +#' @description renders probability of incidence (new for this timestep) for individuals +#' +#' @param birth variable for birth of the individual +#' @param renderer object for model outputs +#' @param source_pop the population which is sampled for infection (bitten SAU for incidence, infecte for clinical/severe) +#' @param prob probability of infection +#' @param prefix for model outputs +#' @param lowers age bounds +#' @param uppers age bounds +#' @param timestep current target +#' +#' @noRd +incidence_probability_renderer <- function( + birth, + renderer, + source_pop, + prob, + prefix, + lowers, + uppers, + timestep +) { + for (i in seq_along(lowers)) { + lower <- lowers[[i]] + upper <- uppers[[i]] + in_age <- in_age_range(birth, timestep, lower, upper) renderer$render( paste0('p_', prefix, lower, '_', upper), sum(prob[bitset_index(source_pop, in_age)]), @@ -113,9 +133,9 @@ incidence_renderer <- function( } create_variable_mean_renderer_process <- function( - renderer, - names, - variables + renderer, + names, + variables ) { function(timestep) { for (i in seq_along(variables)) { @@ -129,12 +149,12 @@ create_variable_mean_renderer_process <- function( } create_vector_count_renderer_individual <- function( - mosquito_state, - species, - state, - renderer, - parameters - ) { + mosquito_state, + species, + state, + renderer, + parameters +) { function(timestep) { adult <- mosquito_state$get_index_of('NonExistent')$not(TRUE) for (i in seq_along(parameters$species)) { @@ -169,9 +189,20 @@ create_age_group_renderer <- function( renderer ) { function(timestep) { - for (i in seq_along(parameters$age_group_rendering_min_ages)) { - lower <- parameters$age_group_rendering_min_ages[[i]] - upper <- parameters$age_group_rendering_max_ages[[i]] + + unique_age_combinations <- as.data.frame(unique(rbind(cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages), + cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages), + cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages), + cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages), + cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages)))) + + ordered_unique_age_combinations <- unique_age_combinations[ + with(unique_age_combinations, order(V1, V2)), + ] + + for (i in seq_along(ordered_unique_age_combinations$V1)) { + lower <- ordered_unique_age_combinations$V1[[i]] + upper <- ordered_unique_age_combinations$V2[[i]] in_age <- in_age_range(birth, timestep, lower, upper) renderer$render( paste0('n_age_', lower, '_', upper), @@ -181,3 +212,54 @@ create_age_group_renderer <- function( } } } + + +populate_incidence_rendering_columns <- function(renderer, parameters){ + + # infections must render in all simulations + renderer$set_default('n_infections', 0) + + # treatment associated only renders when drugs are used + if(sum(unlist(parameters$clinical_treatment_coverages))>0){ + renderer$set_default('ft', 0) + renderer$set_default('n_treated', 0) + renderer$set_default('n_drug_efficacy_failures', 0) + renderer$set_default('n_successfully_treated', 0) + } + + # ETC, SPC only render when antimalarial resistance is on + if(parameters$antimalarial_resistance){ + renderer$set_default('n_early_treatment_failure', 0) + renderer$set_default('n_slow_parasite_clearance', 0) + } + + if(length(parameters$incidence_rendering_min_ages)>0){ + render_initial_incidence(renderer, + parameters$incidence_rendering_min_ages, + parameters$incidence_rendering_max_ages, + "inc") + } + + if(length(parameters$clinical_incidence_rendering_min_ages)>0){ + render_initial_incidence(renderer, + parameters$clinical_incidence_rendering_min_ages, + parameters$clinical_incidence_rendering_max_ages, + "inc_clinical") + } + + if(length(parameters$severe_incidence_rendering_min_ages)>0){ + render_initial_incidence(renderer, + parameters$severe_incidence_rendering_min_ages, + parameters$severe_incidence_rendering_max_ages, + "inc_severe") + } + +} + +render_initial_incidence <- function(renderer, lower_vals, upper_vals, inc_name){ + for (i in seq_along(lower_vals)){ + renderer$set_default(paste0('n_', inc_name, "_", lower_vals[i], '_', upper_vals[i]), 0) + renderer$set_default(paste0('p_', inc_name, "_", lower_vals[i], '_', upper_vals[i]), 0) + } +} + \ No newline at end of file diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index eda048a1..25637f0d 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -843,6 +843,7 @@ test_that('update_severe_disease renders with no infections', { render_function <- mockery::mock() mockery::stub(update_severe_disease, 'incidence_renderer', render_function) + mockery::stub(update_severe_disease, 'incidence_probability_renderer', render_function) empty <- individual::Bitset$new(population) update_severe_disease( @@ -852,13 +853,24 @@ test_that('update_severe_disease renders with no infections', { parameters, renderer ) - + mockery::expect_args( render_function, 1, variables$birth, renderer, empty, + 'inc_severe_', + 0, + 5 * 365, + timestep + ) + + mockery::expect_args( + render_function, + 2, + variables$birth, + renderer, empty, double(0), 'inc_severe_', diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R index cc4e9f87..8184c66f 100644 --- a/tests/testthat/test-render.R +++ b/tests/testthat/test-render.R @@ -26,14 +26,6 @@ test_that('that default rendering works', { mockery::expect_args( renderer$render_mock(), 1, - 'n_730_3650', - 3, - timestep - ) - - mockery::expect_args( - renderer$render_mock(), - 2, 'n_detect_lm_730_3650', 2, timestep @@ -41,7 +33,7 @@ test_that('that default rendering works', { mockery::expect_args( renderer$render_mock(), - 3, + 2, 'p_detect_lm_730_3650', 1.5, timestep @@ -49,7 +41,7 @@ test_that('that default rendering works', { mockery::expect_args( renderer$render_mock(), - 4, + 3, 'n_detect_pcr_730_3650', 3, timestep @@ -83,7 +75,7 @@ test_that('that default rendering works when no one is in the age range', { mockery::expect_args( renderer$render_mock(), 1, - 'n_730_3650', + 'n_detect_lm_730_3650', 0, timestep ) @@ -101,6 +93,15 @@ test_that('that clinical incidence rendering works', { birth, renderer, individual::Bitset$new(4)$insert(c(1, 2, 4)), + 'inc_clinical_', + c(0, 2) * year, + c(5, 10) * year, + timestep + ) + + incidence_probability_renderer( + birth, + renderer, individual::Bitset$new(4)$insert(seq(4)), c(.1, .2, .3, .4), 'inc_clinical_', @@ -112,7 +113,7 @@ test_that('that clinical incidence rendering works', { mockery::expect_args( renderer$render_mock(), 1, - 'n_0_1825', + 'n_inc_clinical_0_1825', 2, timestep ) @@ -120,7 +121,7 @@ test_that('that clinical incidence rendering works', { mockery::expect_args( renderer$render_mock(), 2, - 'n_inc_clinical_0_1825', + 'n_inc_clinical_730_3650', 2, timestep ) @@ -132,11 +133,68 @@ test_that('that clinical incidence rendering works', { .3, timestep ) + + mockery::expect_args( + renderer$render_mock(), + 4, + 'p_inc_clinical_730_3650', + .6, + timestep + ) +}) + +test_that('that age rendering works', { + timestep <- 0 + year <- 365 + birth <- individual::IntegerVariable$new( + -c(0:101) * year + ) + + parameters <- get_parameters(overrides = list(prevalence_rendering_min_ages = c(0, 1, 5) * year, + prevalence_rendering_max_ages = c(2, 5, 8) * year - 1, + incidence_rendering_min_ages = c(3, 5, 10) * year, + incidence_rendering_max_ages = c(10, 20, 100) * year - 1, + clinical_incidence_rendering_min_ages = c(0, 1, 5) * year, + clinical_incidence_rendering_max_ages = c(2, 5, 8) * year - 1)) + + renderer <- mock_render(1) + process <- create_age_group_renderer( + birth, + parameters, + renderer) + + process(timestep) + + + mockery::expect_args( + renderer$render_mock(), + 1, + 'n_age_0_729', + 2, + timestep + ) + mockery::expect_args( renderer$render_mock(), + 2, + 'n_age_365_1824', 4, - 'n_730_3650', + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 3, + 'n_age_1095_3649', + 7, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 4, + 'n_age_1825_2919', 3, timestep ) @@ -144,16 +202,16 @@ test_that('that clinical incidence rendering works', { mockery::expect_args( renderer$render_mock(), 5, - 'n_inc_clinical_730_3650', - 2, + 'n_age_1825_7299', + 15, timestep ) mockery::expect_args( renderer$render_mock(), 6, - 'p_inc_clinical_730_3650', - .6, + 'n_age_3650_36499', + 90, timestep ) }) From 6a4311a44e7733ef7f9660f272bcb9e59a78889e Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 23 May 2024 11:23:28 +0100 Subject: [PATCH 13/35] Removed "age" from age outputs in tests and vignettes. Corrected dt variables recovery rate assignment when antimalarial resistance is modelled. --- R/disease_progression.R | 6 +++--- tests/testthat/test-render.R | 12 ++++++------ vignettes/Demography.Rmd | 2 +- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/R/disease_progression.R b/R/disease_progression.R index d72d2c75..55103ac1 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -7,14 +7,14 @@ #' @param recovery_outcome competing hazards object for recovery rates #' @noRd calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt_input){ - + # Get correct input for dt depending on spc if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ dt_v <- dt_input } else if (isTRUE(parameters$antimalarial_resistance)) { - dt_v <- dt_input$and(variables$state$get_index_of("Tr")$to_vector())$to_vector() + dt_v <- dt_input$get_values(variables$state$get_index_of("Tr")) } - + recovery_rates <- numeric(length = parameters$human_population) recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R index 8184c66f..ce3ebcac 100644 --- a/tests/testthat/test-render.R +++ b/tests/testthat/test-render.R @@ -170,7 +170,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 1, - 'n_age_0_729', + 'n_0_729', 2, timestep ) @@ -178,7 +178,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 2, - 'n_age_365_1824', + 'n_365_1824', 4, timestep ) @@ -186,7 +186,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 3, - 'n_age_1095_3649', + 'n_1095_3649', 7, timestep ) @@ -194,7 +194,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 4, - 'n_age_1825_2919', + 'n_1825_2919', 3, timestep ) @@ -202,7 +202,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 5, - 'n_age_1825_7299', + 'n_1825_7299', 15, timestep ) @@ -210,7 +210,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 6, - 'n_age_3650_36499', + 'n_3650_36499', 90, timestep ) diff --git a/vignettes/Demography.Rmd b/vignettes/Demography.Rmd index d8c96a93..8b53421f 100644 --- a/vignettes/Demography.Rmd +++ b/vignettes/Demography.Rmd @@ -124,7 +124,7 @@ convert_to_long <- function(age_min, age_max, output) { data.frame( age_lower = age_min[[i]], age_upper = age_max[[i]], - n = output[,paste0('n_age_', age_min[[i]], '_',age_max[[i]])], + n = output[,paste0('n_', age_min[[i]], '_',age_max[[i]])], age_plot = age_min[[i]]/365, run = output$run, timestep = output$timestep) From 13748fdbe7fa5eaca821a7df60ecebdb89a7b224 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 9 May 2024 12:16:02 +0100 Subject: [PATCH 14/35] competing_hazards.R: Giovanni designed two objects: CompetingOutcome (containing an outcome following resolution and rates that that outcome occurs for each individual in the population) and CompetingHazard (contains all CompetingOutcomes and process function to resolve hazards). processes.R: create_biting_process now replaced with create_infection_rates_process and biting_process, that uses essentially the same functional pathway, but which stops when the rates of infection are generated, storing them in the infection_outcome object within human_infection.R Similarly, disease_progression_process(es) are replaced with new functions: create_recovery_rates_process and calculate_recovery_rates, that get the recovery rates for each individual based on disease state. We then add these processes to the process list, including the resolution function create_infection_recovery_hazard_process$resolve. This process (create_infection_recovery_hazard_process) contains infection_outcome and recovery_outcome. infection_outcome calls infection_process_resolved_hazard, which follows the rest of the infection pathway functions: assigning new infections to different human states. Note that this requires probability of infection which gets used in the incidence renderer function. recovery_outcome also updates human states and infectivities. I've adjusted some of the tests to reflect the changes in functions required to generate results. Removed white space deletions. --- R/competing_hazards.R | 6 +----- R/disease_progression.R | 6 +++--- R/human_infection.R | 2 -- R/processes.R | 3 +-- tests/testthat/test-infection-integration.R | 2 -- 5 files changed, 5 insertions(+), 14 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index 4ad9b8f6..cdc6aea0 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -18,15 +18,11 @@ CompetingOutcome <- R6::R6Class( set_rates = function(rates){ self$rates <- rates }, - stash_source_humans = function(source_humans){ - self$source_humans <- source_humans - }, execute = function(t, target){ private$targeted_process(t, target) self$rates <- rep(0, length(self$rates)) }, - rates = NULL, - source_humans = NULL + rates = NULL ) ) diff --git a/R/disease_progression.R b/R/disease_progression.R index 55103ac1..d72d2c75 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -7,14 +7,14 @@ #' @param recovery_outcome competing hazards object for recovery rates #' @noRd calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt_input){ - + # Get correct input for dt depending on spc if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ dt_v <- dt_input } else if (isTRUE(parameters$antimalarial_resistance)) { - dt_v <- dt_input$get_values(variables$state$get_index_of("Tr")) + dt_v <- dt_input$and(variables$state$get_index_of("Tr")$to_vector())$to_vector() } - + recovery_rates <- numeric(length = parameters$human_population) recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da diff --git a/R/human_infection.R b/R/human_infection.R index 50568a26..9cb9c2c1 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -133,7 +133,6 @@ calculate_infections <- function( #' and treated malaria; and resulting boosts in immunity #' @param timestep current timestep #' @param infected_humans bitset of infected humans -#' @param source_humans bitset of humans with infection potential #' @param variables a list of all of the model variables #' @param renderer model render object #' @param parameters model parameters @@ -142,7 +141,6 @@ calculate_infections <- function( infection_process_resolved_hazard <- function( timestep, infected_humans, - source_humans, variables, renderer, parameters, diff --git a/R/processes.R b/R/processes.R index 59e0b8c1..fe73b174 100644 --- a/R/processes.R +++ b/R/processes.R @@ -68,8 +68,7 @@ create_processes <- function( infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, infection_outcome$source_humans, - variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + infection_process_resolved_hazard(timestep, target, variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) }, size = parameters$human_population ) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index 25637f0d..e2ef644d 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -85,7 +85,6 @@ test_that('simulate_infection integrates different types of infection and schedu ) 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)) clinical <- individual::Bitset$new(population)$insert(c(1, 3)) clinical_infection_mock <- mockery::mock(clinical) @@ -107,7 +106,6 @@ test_that('simulate_infection integrates different types of infection and schedu infection_process_resolved_hazard( timestep, infected, - source_humans, variables, renderer, parameters, From 98c982f1ff17afaadb3dfe2e0092e300763901ec Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 9 May 2024 16:11:41 +0100 Subject: [PATCH 15/35] Added a stash source_human function to the competing hazards function for use in infection outcome. --- R/competing_hazards.R | 6 +++++- R/human_infection.R | 2 ++ R/processes.R | 3 ++- tests/testthat/test-infection-integration.R | 2 ++ 4 files changed, 11 insertions(+), 2 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index cdc6aea0..4ad9b8f6 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -18,11 +18,15 @@ CompetingOutcome <- R6::R6Class( set_rates = function(rates){ self$rates <- rates }, + stash_source_humans = function(source_humans){ + self$source_humans <- source_humans + }, execute = function(t, target){ private$targeted_process(t, target) self$rates <- rep(0, length(self$rates)) }, - rates = NULL + rates = NULL, + source_humans = NULL ) ) diff --git a/R/human_infection.R b/R/human_infection.R index 9cb9c2c1..50568a26 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -133,6 +133,7 @@ calculate_infections <- function( #' and treated malaria; and resulting boosts in immunity #' @param timestep current timestep #' @param infected_humans bitset of infected humans +#' @param source_humans bitset of humans with infection potential #' @param variables a list of all of the model variables #' @param renderer model render object #' @param parameters model parameters @@ -141,6 +142,7 @@ calculate_infections <- function( infection_process_resolved_hazard <- function( timestep, infected_humans, + source_humans, variables, renderer, parameters, diff --git a/R/processes.R b/R/processes.R index fe73b174..59e0b8c1 100644 --- a/R/processes.R +++ b/R/processes.R @@ -68,7 +68,8 @@ create_processes <- function( infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + infection_process_resolved_hazard(timestep, target, infection_outcome$source_humans, + variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) }, size = parameters$human_population ) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index e2ef644d..25637f0d 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -85,6 +85,7 @@ test_that('simulate_infection integrates different types of infection and schedu ) 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)) clinical <- individual::Bitset$new(population)$insert(c(1, 3)) clinical_infection_mock <- mockery::mock(clinical) @@ -106,6 +107,7 @@ test_that('simulate_infection integrates different types of infection and schedu infection_process_resolved_hazard( timestep, infected, + source_humans, variables, renderer, parameters, From 16d83d35cc4a88365db251718881446cb928d40f Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Fri, 24 May 2024 10:32:45 +0100 Subject: [PATCH 16/35] Removed missed merge commit line. --- vignettes/Carrying-capacity.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/Carrying-capacity.Rmd b/vignettes/Carrying-capacity.Rmd index db38e54a..0494e341 100644 --- a/vignettes/Carrying-capacity.Rmd +++ b/vignettes/Carrying-capacity.Rmd @@ -106,7 +106,7 @@ p_lsm <- p |> ) set.seed(123) s_lsm <- run_simulation(timesteps = timesteps, parameters = p_lsm) -s_lsm$pfpr <- s_lsm$n_detect_lm_730_3650 / s_lsm$n_730_3650 +s_lsm$pfpr <- s_lsm$n_detect_lm_730_3650 / s_lsm$n_age_730_3650 par(mfrow = c(1, 2)) plot(s$EIR_gamb ~ s$timestep, t = "l", ylim = c(0, 150), xlab = "Time", ylab = "EIR") From c649292c2ea3b5ac0cfe3feef126578d5c5d4332 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Mon, 24 Jun 2024 16:31:00 +0100 Subject: [PATCH 17/35] biting_process renamed as create_biting_process, where functionality now occurs in biting_process.R. Similarly, create_recovery_rates_process functionality is now in disease_progression.R stash_source_humans in competing hazard resolution is now redundant and has been removed from infection processes. recovery rates are now assigned to a variable which is updated as infections, recoveries (via update_infection) and deaths occur. This also takes into account recovery time of treated due to slow parasite clearance and replaced the dt variable. I've also noticed (I think) that Tr->S and U->S recoveries can be done in a single step, so have attempted to combine these at risk of making the code more complicated! Added a metapopulation incidence rendering wrapper function for initialising multiple renderers. Rewrote create_age_group_renderer for clarity following feedback. Some tests have been adjusted to reflect these changes. --- R/biting_process.R | 61 ++++++++-------- R/competing_hazards.R | 6 +- R/disease_progression.R | 77 +++++++++++---------- R/human_infection.R | 28 ++++---- R/model.R | 1 + R/mortality_processes.R | 6 +- R/processes.R | 61 ++++++---------- R/render.R | 30 +++++--- R/utils.R | 2 +- R/variables.R | 20 +++--- tests/testthat/test-biting-integration.R | 26 +++---- tests/testthat/test-infection-integration.R | 50 ++++++------- 12 files changed, 181 insertions(+), 187 deletions(-) diff --git a/R/biting_process.R b/R/biting_process.R index 8ae1147f..594b9429 100644 --- a/R/biting_process.R +++ b/R/biting_process.R @@ -18,7 +18,7 @@ #' @param infection_outcome competing hazards object for infection rates #' @param timestep the current timestep #' @noRd -biting_process <- function( +create_biting_process <- function( renderer, solvers, models, @@ -29,37 +29,36 @@ biting_process <- function( lagged_eir, mixing = 1, mixing_index = 1, - infection_outcome, - timestep -) { - - age <- get_age(variables$birth$get_values(), timestep) - - bitten_humans <- simulate_bites( - renderer, - solvers, - models, - variables, - events, - age, - parameters, - timestep, - lagged_infectivity, - lagged_eir, - mixing, - mixing_index - ) - - simulate_infection( - variables, - events, - bitten_humans, - age, - parameters, - timestep, - renderer, infection_outcome - ) +) { + function(timestep){ + age <- get_age(variables$birth$get_values(), timestep) + bitten_humans <- simulate_bites( + renderer, + solvers, + models, + variables, + events, + age, + parameters, + timestep, + lagged_infectivity, + lagged_eir, + mixing, + mixing_index + ) + + simulate_infection( + variables, + events, + bitten_humans, + age, + parameters, + timestep, + renderer, + infection_outcome + ) + } } #' @importFrom stats rpois diff --git a/R/competing_hazards.R b/R/competing_hazards.R index 4ad9b8f6..cdc6aea0 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -18,15 +18,11 @@ CompetingOutcome <- R6::R6Class( set_rates = function(rates){ self$rates <- rates }, - stash_source_humans = function(source_humans){ - self$source_humans <- source_humans - }, execute = function(t, target){ private$targeted_process(t, target) self$rates <- rep(0, length(self$rates)) }, - rates = NULL, - source_humans = NULL + rates = NULL ) ) diff --git a/R/disease_progression.R b/R/disease_progression.R index d72d2c75..f5cef8c5 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -1,28 +1,20 @@ #' @title Calculate recovery rates #' @description Calculates recovery rates for each individual in the population -#' for storage in competing hazards object and future resolution +#' for storage in competing hazards object and subsequent resolution #' #' @param variables the available human variables -#' @param parameters model parameters #' @param recovery_outcome competing hazards object for recovery rates #' @noRd -calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt_input){ - - # Get correct input for dt depending on spc - if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ - dt_v <- dt_input - } else if (isTRUE(parameters$antimalarial_resistance)) { - dt_v <- dt_input$and(variables$state$get_index_of("Tr")$to_vector())$to_vector() +create_recovery_rates_process <- function( + variables, + recovery_outcome +) { + function(timestep){ + recovery_outcome$set_rates(variables$recovery_rates$get_values()) } - - recovery_rates <- numeric(length = parameters$human_population) - recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd - recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da - recovery_rates[variables$state$get_index_of("U")$to_vector()] <- 1/parameters$du - recovery_rates[variables$state$get_index_of("Tr")$to_vector()] <- 1/dt_v - recovery_outcome$set_rates(recovery_rates) } + #' @title Disease progression outcomes (recovery) #' @description Following resolution of competing hazards, update state and #' infectivity of sampled individuals @@ -33,7 +25,7 @@ calculate_recovery_rates <- function(variables, parameters, recovery_outcome, dt #' @param parameters model parameters #' @param renderer competing hazards object for recovery rates #' @noRd -recovery_process_resolved_hazard <- function( +recovery_outcome_process <- function( timestep, target, variables, @@ -53,25 +45,31 @@ recovery_process_resolved_hazard <- function( "U", variables$infectivity, parameters$cu, + variables$recovery_rates, + 1/parameters$du, variables$state$get_index_of("A")$and(target) ) + # Is there a reason we aren't doing this in one step...? These all recover to S... update_infection( variables$state, "S", variables$infectivity, 0, - variables$state$get_index_of("U")$and(target) - ) - - update_infection( - variables$state, - "S", - variables$infectivity, + variables$recovery_rates, 0, - variables$state$get_index_of("Tr")$and(target) + variables$state$get_index_of(c("U","Tr"))$and(target) ) + # update_infection( + # variables$state, + # "S", + # variables$infectivity, + # 0, + # variables$recovery_rates, + # 0, + # variables$state$get_index_of("Tr")$and(target) + # ) } #' @title Update the state of an individual as infection events occur @@ -84,14 +82,17 @@ recovery_process_resolved_hazard <- function( #' @param new_infectivity the new infectivity of the progressed individuals #' @noRd update_infection <- function( - state, - to_state, - infectivity, - new_infectivity, - to_move - ) { + state, + to_state, + infectivity, + new_infectivity, + recovery_rates, + new_recovery_rate, + to_move +) { state$queue_update(to_state, to_move) infectivity$queue_update(new_infectivity, to_move) + recovery_rates$queue_update(new_recovery_rate, to_move) } #' @title Modelling the progression to asymptomatic disease @@ -102,11 +103,11 @@ update_infection <- function( #' @param parameters model parameters #' @noRd update_to_asymptomatic_infection <- function( - variables, - parameters, - timestep, - to_move - ) { + variables, + parameters, + timestep, + to_move +) { if (to_move$size() > 0) { variables$state$queue_update('A', to_move) new_infectivity <- asymptomatic_infectivity( @@ -121,5 +122,9 @@ update_to_asymptomatic_infection <- function( new_infectivity, to_move ) + variables$recovery_rates$queue_update( + 1/parameters$da, + to_move + ) } } diff --git a/R/human_infection.R b/R/human_infection.R index 50568a26..d8a39d90 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -1,7 +1,7 @@ #' @title Simulate malaria infection in humans #' @description #' This function ends with the assignment of rates of infection to the competing -#' hazard resolution object. Boosts immunity given infectious bites. +#' hazard resolution object and boosts immunity given infectious bites. #' @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 @@ -60,7 +60,6 @@ calculate_infections <- function( ) { source_humans <- variables$state$get_index_of( c('S', 'A', 'U'))$and(bitten_humans) - infection_outcome$stash_source_humans(source_humans) b <- blood_immunity(variables$ib$get_values(source_humans), parameters) @@ -113,7 +112,7 @@ calculate_infections <- function( incidence_probability_renderer( variables$birth, renderer, - source_pop, + source_humans, prob, "inc_", parameters$incidence_min_ages, @@ -122,7 +121,7 @@ calculate_infections <- function( ) ## capture infection rates to resolve in competing hazards - infection_rates <- numeric(length = parameters$human_population) + infection_rates <- rep(0, length = parameters$human_population) infection_rates[source_vector] <- prob_to_rate(prob) infection_outcome$set_rates(infection_rates) } @@ -133,16 +132,14 @@ calculate_infections <- function( #' and treated malaria; and resulting boosts in immunity #' @param timestep current timestep #' @param infected_humans bitset of infected humans -#' @param source_humans bitset of humans with infection potential #' @param variables a list of all of the model variables #' @param renderer model render object #' @param parameters model parameters #' @param prob vector of population probabilities of infection #' @noRd -infection_process_resolved_hazard <- function( +infection_outcome_process <- function( timestep, infected_humans, - source_humans, variables, renderer, parameters, @@ -410,14 +407,19 @@ calculate_treated <- function( successfully_treated ) if(parameters$antimalarial_resistance) { - variables$dt$queue_update( - parameters$dt, + variables$recovery_rates$queue_update( + 1/parameters$dt, non_slow_parasite_clearance_individuals ) - variables$dt$queue_update( - dt_slow_parasite_clearance, + variables$recovery_rates$queue_update( + 1/dt_slow_parasite_clearance, slow_parasite_clearance_individuals ) + } else { + variables$recovery_rates$queue_update( + 1/parameters$dt, + successfully_treated + ) } } @@ -455,6 +457,8 @@ schedule_infections <- function( 'D', variables$infectivity, parameters$cd, + variables$recovery_rates, + 1/parameters$dd, to_infect ) } @@ -523,7 +527,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/model.R b/R/model.R index de5d001a..21df1173 100644 --- a/R/model.R +++ b/R/model.R @@ -220,6 +220,7 @@ run_metapop_simulation <- function( variables <- lapply(parameters, create_variables) events <- lapply(parameters, create_events) renderer <- lapply(parameters, function(.) individual::Render$new(timesteps)) + populate_metapopulation_incidence_rendering_columns(renderer, parameters) for (i in seq_along(parameters)) { # NOTE: forceAndCall is necessary here to make sure i refers to the current # iteration diff --git a/R/mortality_processes.R b/R/mortality_processes.R index 37792f6c..e67baf0e 100644 --- a/R/mortality_processes.R +++ b/R/mortality_processes.R @@ -113,12 +113,8 @@ reset_target <- function(variables, events, target, state, parameters, timestep) # onwards infectiousness variables$infectivity$queue_update(0, target) + variables$recovery_rates$queue_update(0, target) - # treated compartment residence time: - if(!is.null(variables$dt)) { - variables$dt$queue_update(parameters$dt, target) - } - # zeta and zeta group and vector controls survive rebirth } } diff --git a/R/processes.R b/R/processes.R index 59e0b8c1..c416a334 100644 --- a/R/processes.R +++ b/R/processes.R @@ -68,31 +68,29 @@ create_processes <- function( infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, infection_outcome$source_humans, - variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates)) + infection_outcome_process(timestep, target, + variables, renderer, parameters, + prob = rate_to_prob(infection_outcome$rates)) }, size = parameters$human_population ) recovery_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - recovery_process_resolved_hazard(timestep, target, variables, parameters, renderer) + recovery_outcome_process(timestep, target, variables, parameters, renderer) }, size = parameters$human_population ) - create_infection_recovery_hazard_process <- CompetingHazard$new( - outcomes = list(infection_outcome, recovery_outcome) - ) - # ============================== # Biting and mortality processes # ============================== - # simulate infections for humans and set last_boosted_* + # simulate bites, calculates infection rates for bitten humans and set last_boosted_* # move mosquitoes into incubating state # kill mosquitoes caught in vector control - create_infection_rates_process <- function(timestep){ - biting_process( + processes <- c( + processes, + create_biting_process( renderer, solvers, models, @@ -103,44 +101,27 @@ create_processes <- function( lagged_eir, mixing, mixing_index, - infection_outcome, - timestep + infection_outcome ) - } - - # ======================= - # Antimalarial Resistance - # ======================= - # Add an a new process which governs the transition from Tr to S when - # antimalarial resistance is simulated. The rate of transition switches - # from a parameter to a variable when antimalarial resistance == TRUE. - - # Assign the dt input to a separate object with the default single parameter value: - dt_input <- parameters$dt - - # If antimalarial resistance is switched on, assign dt variable values to the - if(parameters$antimalarial_resistance) { - dt_input <- variables$dt - } + ) # =================== # Disease Progression # =================== - create_recovery_rates_process <- function(timestep){ - calculate_recovery_rates( + processes <- c( + processes, + create_recovery_rates_process( variables, - parameters, - recovery_outcome, - dt_input - ) - } + recovery_outcome + ), + + # Resolve competing hazards of infection with disease progression + CompetingHazard$new( + outcomes = list(infection_outcome, recovery_outcome) + )$resolve + ) - processes <- c(processes, - create_infection_rates_process, - create_recovery_rates_process, - create_infection_recovery_hazard_process$resolve) - # =============== # ODE integration # =============== diff --git a/R/render.R b/R/render.R index 0d1c252e..4b6335eb 100644 --- a/R/render.R +++ b/R/render.R @@ -188,18 +188,20 @@ create_age_group_renderer <- function( parameters, renderer ) { + + age_ranges <- rbind( + cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages), + cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages), + cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages), + cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages), + cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages) + ) + + unique_age_combinations <- as.data.frame(unique(age_ranges)) + ordered_unique_age_combinations <- unique_age_combinations[order(unique_age_combinations$V1, unique_age_combinations$V2), ] + function(timestep) { - unique_age_combinations <- as.data.frame(unique(rbind(cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages), - cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages), - cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages), - cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages), - cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages)))) - - ordered_unique_age_combinations <- unique_age_combinations[ - with(unique_age_combinations, order(V1, V2)), - ] - for (i in seq_along(ordered_unique_age_combinations$V1)) { lower <- ordered_unique_age_combinations$V1[[i]] upper <- ordered_unique_age_combinations$V2[[i]] @@ -262,4 +264,10 @@ render_initial_incidence <- function(renderer, lower_vals, upper_vals, inc_name) renderer$set_default(paste0('p_', inc_name, "_", lower_vals[i], '_', upper_vals[i]), 0) } } - \ No newline at end of file + + +populate_metapopulation_incidence_rendering_columns <- function(renderer, parameters){ + for(i in length(parameters)){ + populate_incidence_rendering_columns(renderer[[i]], parameters[[i]]) + } +} \ No newline at end of file diff --git a/R/utils.R b/R/utils.R index 34de912e..8ce50545 100644 --- a/R/utils.R +++ b/R/utils.R @@ -100,5 +100,5 @@ prob_to_rate <- function(prob){ #'@param rate rate #'@noRd rate_to_prob <- function(rate){ - 1- exp(-rate) + 1 - exp(-rate) } diff --git a/R/variables.R b/R/variables.R index 86a13185..5a6dc749 100644 --- a/R/variables.R +++ b/R/variables.R @@ -178,6 +178,7 @@ create_variables <- function(parameters) { diseased <- state$get_index_of('D')$to_vector() asymptomatic <- state$get_index_of('A')$to_vector() subpatent <- state$get_index_of('U')$to_vector() + treated <- state$get_index_of('Tr')$to_vector() # Set the initial infectivity values for each individual infectivity_values[diseased] <- parameters$cd @@ -191,6 +192,16 @@ create_variables <- function(parameters) { # Initialise the infectivity variable infectivity <- individual::DoubleVariable$new(infectivity_values) + # Set recovery rate for each individual + recovery_values <- rep(0, get_human_population(parameters, 0)) + recovery_values[diseased] <- 1/parameters$dd + recovery_values[asymptomatic] <- 1/parameters$da + recovery_values[subpatent] <- 1/parameters$du + recovery_values[treated] <- 1/parameters$dt + + # Initialise the recovery rate variable + recovery_rates <- individual::DoubleVariable$new(recovery_values) + drug <- individual::IntegerVariable$new(rep(0, size)) drug_time <- individual::IntegerVariable$new(rep(-1, size)) @@ -220,6 +231,7 @@ create_variables <- function(parameters) { zeta = zeta, zeta_group = zeta_group, infectivity = infectivity, + recovery_rates = recovery_rates, drug = drug, drug_time = drug_time, last_pev_timestep = last_pev_timestep, @@ -230,14 +242,6 @@ create_variables <- function(parameters) { spray_time = spray_time ) - # Add variables for antimalarial resistance state residency times (dt) - if(parameters$antimalarial_resistance) { - dt <- individual::DoubleVariable$new(rep(parameters$dt, size)) - variables <- c( - variables, - dt = dt - ) - } # Add variables for individual mosquitoes if (parameters$individual_mosquitoes) { diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R index 9712a243..4be80329 100644 --- a/tests/testthat/test-biting-integration.R +++ b/tests/testthat/test-biting-integration.R @@ -21,15 +21,7 @@ test_that('biting_process integrates mosquito effects and human infection', { size = parameters$human_population ) - bitten <- individual::Bitset$new(parameters$human_population) - bites_mock <- mockery::mock(bitten) - infection_mock <- mockery::mock() - - mockery::stub(biting_process, 'simulate_bites', bites_mock) - mockery::stub(biting_process, 'simulate_infection', infection_mock) - - - biting_process( + biting_process <- create_biting_process( renderer, solvers, models, @@ -38,11 +30,19 @@ test_that('biting_process integrates mosquito effects and human infection', { parameters, lagged_foim, lagged_eir, - mixing = 1, - mixing_index = 1, - infection_outcome, - timestep + 1, + 1, + infection_outcome ) + + bitten <- individual::Bitset$new(parameters$human_population) + bites_mock <- mockery::mock(bitten) + infection_mock <- mockery::mock() + + mockery::stub(biting_process, 'simulate_bites', bites_mock) + mockery::stub(biting_process, 'simulate_infection', infection_mock) + biting_process(timestep) + mockery::expect_args( bites_mock, 1, diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index 25637f0d..d399a527 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -23,7 +23,7 @@ test_that('simulate_infection integrates different types of infection and schedu infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + infection_outcome_process(timestep, target, variables, renderer, parameters) }, size = parameters$human_population ) @@ -97,17 +97,16 @@ test_that('simulate_infection integrates different types of infection and schedu schedule_mock <- mockery::mock() - mockery::stub(infection_process_resolved_hazard, 'incidence_renderer', mockery::mock()) - mockery::stub(infection_process_resolved_hazard, 'boost_immunity', boost_immunity_mock) - mockery::stub(infection_process_resolved_hazard, 'calculate_clinical_infections', clinical_infection_mock) - mockery::stub(infection_process_resolved_hazard, 'update_severe_disease', severe_infection_mock) - mockery::stub(infection_process_resolved_hazard, 'calculate_treated', treated_mock) - mockery::stub(infection_process_resolved_hazard, 'schedule_infections', schedule_mock) + 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_process_resolved_hazard( + infection_outcome_process( timestep, infected, - source_humans, variables, renderer, parameters, @@ -221,7 +220,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + infection_outcome_process(timestep, target, variables, renderer, parameters) }, size = parameters$human_population ) @@ -321,6 +320,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat variables <- list( state = list(queue_update = mockery::mock()), infectivity = list(queue_update = mockery::mock()), + recovery_rates = list(queue_update = mockery::mock()), drug = list(queue_update = mockery::mock()), drug_time = list(queue_update = mockery::mock()) ) @@ -410,9 +410,9 @@ test_that('calculate_treated correctly samples treated and updates the drug stat variables <- list( state = list(queue_update = mockery::mock()), infectivity = list(queue_update = mockery::mock()), + recovery_rates = list(queue_update = mockery::mock()), drug = list(queue_update = mockery::mock()), - drug_time = list(queue_update = mockery::mock()), - dt = list(queue_update = mockery::mock()) + drug_time = list(queue_update = mockery::mock()) ) renderer <- individual::Render$new(timesteps = timestep) @@ -487,8 +487,8 @@ test_that('calculate_treated correctly samples treated and updates the drug stat ) expect_bitset_update(variables$drug$queue_update, c(2, 1, 1, 1, 2, 2, 2), c(1, 2, 3, 4, 5, 6, 7)) expect_bitset_update(variables$drug_time$queue_update, 5, c(1, 2, 3, 4, 5, 6, 7)) - expect_bitset_update(variables$dt$queue_update, 5, c(2, 3, 4, 5, 6, 7), 1) - expect_bitset_update(variables$dt$queue_update, 15, c(1), 2) + expect_bitset_update(variables$recovery_rates$queue_update, 1/5, c(2, 3, 4, 5, 6, 7), 1) + expect_bitset_update(variables$recovery_rates$queue_update, 1/15, c(1), 2) }) test_that('calculate_treated correctly samples treated and updates the drug state when resistance not set for all drugs', { @@ -523,9 +523,9 @@ test_that('calculate_treated correctly samples treated and updates the drug stat variables <- list( state = list(queue_update = mockery::mock()), infectivity = list(queue_update = mockery::mock()), + recovery_rates = list(queue_update = mockery::mock()), drug = list(queue_update = mockery::mock()), - drug_time = list(queue_update = mockery::mock()), - dt = list(queue_update = mockery::mock()) + drug_time = list(queue_update = mockery::mock()) ) # Create a Bitset of individuals seeking treatment individuals: @@ -629,15 +629,15 @@ test_that('calculate_treated correctly samples treated and updates the drug stat # Check that update queued for dt for the non-slow parasite clearance individuals is correct: expect_bitset_update( - variables$dt$queue_update, - parameters$dt, + variables$recovery_rates$queue_update, + 1/parameters$dt, c(2, 3, 4, 5, 6, 7), 1) # Check that update queued for dt for the slow parasite clearance individuals is correct: expect_bitset_update( - variables$dt$queue_update, - unlist(parameters$dt_slow_parasite_clearance), + variables$recovery_rates$queue_update, + 1/unlist(parameters$dt_slow_parasite_clearance), c(1), 2) @@ -665,8 +665,8 @@ test_that('schedule_infections correctly schedules new infections', { parameters, 42 ) - - actual_infected <- mockery::mock_args(infection_mock)[[1]][[5]]$to_vector() + + actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector() actual_asymp_infected <- mockery::mock_args(asymp_mock)[[1]][[4]]$to_vector() expect_equal( @@ -705,7 +705,7 @@ test_that('prophylaxis is considered for medicated humans', { infection_outcome <- CompetingOutcome$new( targeted_process = function(timestep, target){ - infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) + infection_outcome_process(timestep, target, variables, renderer, parameters) }, size = parameters$human_population ) @@ -1053,9 +1053,9 @@ test_that('calculate_treated() successfully adds additional resistance columns t variables <- list( state = list(queue_update = mockery::mock()), infectivity = list(queue_update = mockery::mock()), + recovery_rates = list(queue_update = mockery::mock()), drug = list(queue_update = mockery::mock()), - drug_time = list(queue_update = mockery::mock()), - dt = list(queue_update = mockery::mock()) + drug_time = list(queue_update = mockery::mock()) ) renderer <- individual::Render$new(timesteps = 10) From c85ca933a0306065c0f41977ff4b067a7ae87069 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Mon, 24 Jun 2024 17:21:01 +0100 Subject: [PATCH 18/35] Re-added "age" to age outputs. --- tests/testthat/test-render.R | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R index ce3ebcac..8184c66f 100644 --- a/tests/testthat/test-render.R +++ b/tests/testthat/test-render.R @@ -170,7 +170,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 1, - 'n_0_729', + 'n_age_0_729', 2, timestep ) @@ -178,7 +178,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 2, - 'n_365_1824', + 'n_age_365_1824', 4, timestep ) @@ -186,7 +186,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 3, - 'n_1095_3649', + 'n_age_1095_3649', 7, timestep ) @@ -194,7 +194,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 4, - 'n_1825_2919', + 'n_age_1825_2919', 3, timestep ) @@ -202,7 +202,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 5, - 'n_1825_7299', + 'n_age_1825_7299', 15, timestep ) @@ -210,7 +210,7 @@ test_that('that age rendering works', { mockery::expect_args( renderer$render_mock(), 6, - 'n_3650_36499', + 'n_age_3650_36499', 90, timestep ) From 9142c585440e632397ad07d83e38ec089b351f8e Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Wed, 26 Jun 2024 10:13:32 +0100 Subject: [PATCH 19/35] Corrected age output names in Demography vignette to enable it to run (e.g., n_age... rather than n_...) --- vignettes/Demography.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/Demography.Rmd b/vignettes/Demography.Rmd index bb439836..9de15d9c 100644 --- a/vignettes/Demography.Rmd +++ b/vignettes/Demography.Rmd @@ -124,7 +124,7 @@ convert_to_long <- function(age_min, age_max, output) { data.frame( age_lower = age_min[[i]], age_upper = age_max[[i]], - n = output[,paste0('n_', age_min[[i]], '_',age_max[[i]])], + n = output[,paste0('n_age_', age_min[[i]], '_',age_max[[i]])], age_plot = age_min[[i]]/365, run = output$run, timestep = output$timestep) From b6787f7aa6f87d9d226dfb9e25a7aeff3e97f4bb Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Wed, 3 Jul 2024 15:04:53 +0100 Subject: [PATCH 20/35] Removed redundant commented lines. --- R/disease_progression.R | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/R/disease_progression.R b/R/disease_progression.R index f5cef8c5..11f63065 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -50,7 +50,6 @@ recovery_outcome_process <- function( variables$state$get_index_of("A")$and(target) ) - # Is there a reason we aren't doing this in one step...? These all recover to S... update_infection( variables$state, "S", @@ -61,15 +60,6 @@ recovery_outcome_process <- function( variables$state$get_index_of(c("U","Tr"))$and(target) ) - # update_infection( - # variables$state, - # "S", - # variables$infectivity, - # 0, - # variables$recovery_rates, - # 0, - # variables$state$get_index_of("Tr")$and(target) - # ) } #' @title Update the state of an individual as infection events occur From 5fec286f4556c1a4c8d0b56b606e60a43d13679c Mon Sep 17 00:00:00 2001 From: Giovanni Charles Date: Wed, 3 Jul 2024 15:12:20 +0100 Subject: [PATCH 21/35] Add print step to touchstone workflow --- .github/workflows/touchstone.yaml | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/.github/workflows/touchstone.yaml b/.github/workflows/touchstone.yaml index ff1124ad..bc36657f 100644 --- a/.github/workflows/touchstone.yaml +++ b/.github/workflows/touchstone.yaml @@ -20,6 +20,17 @@ env: R_GC_MEM_GROW: 3 jobs: + + print: + name: Print event + runs-on: ubuntu-latest + steps: + - name: Dump GitHub context + env: + GITHUB_CONTEXT: ${{ toJson(github) }} + run: | + echo "$GITHUB_CONTEXT" + prepare: # The other jobs all depend on this one succeeding. They'll implicitly get # skipped as well if this condition is not met. From 598ade0dfe86e027ee7c01564092896113c5ff55 Mon Sep 17 00:00:00 2001 From: Giovanni Charles Date: Wed, 3 Jul 2024 15:20:55 +0100 Subject: [PATCH 22/35] Revert "Add print step to touchstone workflow" This reverts commit 5fec286f4556c1a4c8d0b56b606e60a43d13679c. --- .github/workflows/touchstone.yaml | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/.github/workflows/touchstone.yaml b/.github/workflows/touchstone.yaml index bc36657f..ff1124ad 100644 --- a/.github/workflows/touchstone.yaml +++ b/.github/workflows/touchstone.yaml @@ -20,17 +20,6 @@ env: R_GC_MEM_GROW: 3 jobs: - - print: - name: Print event - runs-on: ubuntu-latest - steps: - - name: Dump GitHub context - env: - GITHUB_CONTEXT: ${{ toJson(github) }} - run: | - echo "$GITHUB_CONTEXT" - prepare: # The other jobs all depend on this one succeeding. They'll implicitly get # skipped as well if this condition is not met. From 09e6e8c166a9d5a6b07e88cf41dd7ec482d3a07c Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Thu, 4 Jul 2024 11:14:49 +0100 Subject: [PATCH 23/35] Improve cumsum performance --- R/competing_hazards.R | 13 ++++++------- R/disease_progression.R | 6 +++++- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index cdc6aea0..a4731f47 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -57,15 +57,14 @@ CompetingHazard <- R6::R6Class( occurs <- occur_rng < rate_to_prob(occur_rates) norm_probs <- event_rates / occur_rates norm_probs[is.na(norm_probs)] <- 0 - cum_probs <- apply(norm_probs, 1, cumsum) + + cumulative <- rep(0, private$size) event_rng <- private$rng(private$size) for(o in seq_along(private$outcomes)){ - if (o == 1) { - selected <- (event_rng <= cum_probs[o,]) - } else { - selected <- (event_rng > cum_probs[o - 1,]) & - (event_rng <= cum_probs[o,]) - } + next_cumulative <- cumulative + norm_probs[,o] + selected <- (event_rng > cumulative) & (event_rng <= next_cumulative) + cumulative <- next_cumulative + target <- individual::Bitset$new(private$size)$insert( which(selected & occurs) ) diff --git a/R/disease_progression.R b/R/disease_progression.R index 11f63065..22121d8a 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -5,12 +5,16 @@ #' @param variables the available human variables #' @param recovery_outcome competing hazards object for recovery rates #' @noRd +recovery_rates_process <- function(variables, recovery_outcome) { + recovery_outcome$set_rates(variables$recovery_rates$get_values()) +} + create_recovery_rates_process <- function( variables, recovery_outcome ) { function(timestep){ - recovery_outcome$set_rates(variables$recovery_rates$get_values()) + recovery_rates_process(variables, recovery_outcome) } } From 0f21bbcf876a5794ce490ea5e4755eaad823a824 Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Thu, 4 Jul 2024 11:18:29 +0100 Subject: [PATCH 24/35] Remove outlining of function --- R/disease_progression.R | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/R/disease_progression.R b/R/disease_progression.R index 22121d8a..11f63065 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -5,16 +5,12 @@ #' @param variables the available human variables #' @param recovery_outcome competing hazards object for recovery rates #' @noRd -recovery_rates_process <- function(variables, recovery_outcome) { - recovery_outcome$set_rates(variables$recovery_rates$get_values()) -} - create_recovery_rates_process <- function( variables, recovery_outcome ) { function(timestep){ - recovery_rates_process(variables, recovery_outcome) + recovery_outcome$set_rates(variables$recovery_rates$get_values()) } } From a58e064294f823ede85cfd6c542baf90abd10645 Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Thu, 4 Jul 2024 14:05:53 +0100 Subject: [PATCH 25/35] Optimize competing hazards. We only need to perform the selection for individuals that passed the occurs check, not the entire population. --- R/competing_hazards.R | 16 ++++++++-------- tests/testthat/test-competing-hazards.R | 2 +- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index a4731f47..d555bd9c 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -54,20 +54,20 @@ CompetingHazard <- R6::R6Class( ) occur_rates <- rowSums(event_rates) occur_rng <- private$rng(private$size) - occurs <- occur_rng < rate_to_prob(occur_rates) - norm_probs <- event_rates / occur_rates + occurs <- which(occur_rng < rate_to_prob(occur_rates)) + + norm_probs <- event_rates[occurs,] / occur_rates[occurs] norm_probs[is.na(norm_probs)] <- 0 - cumulative <- rep(0, private$size) - event_rng <- private$rng(private$size) - for(o in seq_along(private$outcomes)){ + cumulative <- rep(0, length(occurs)) + event_rng <- private$rng(length(occurs)) + + for (o in seq_along(private$outcomes)) { next_cumulative <- cumulative + norm_probs[,o] selected <- (event_rng > cumulative) & (event_rng <= next_cumulative) cumulative <- next_cumulative - target <- individual::Bitset$new(private$size)$insert( - which(selected & occurs) - ) + target <- individual::Bitset$new(private$size)$insert(occurs[selected]) if (target$size() > 0){ private$outcomes[[o]]$execute(t, target) } diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R index ef67e2d0..12203679 100644 --- a/tests/testthat/test-competing-hazards.R +++ b/tests/testthat/test-competing-hazards.R @@ -104,7 +104,7 @@ test_that('hazard resolves two outcomes when some events occur', { outcomes = list(outcome_1, outcome_2), rng = mockery::mock( c(0, 0, 1, 1), # some events occur - c(.05, .3, .2, .5), # event_rng + c(.05, .3), # event_rng ) ) From 1f7310b2f47fbd02b9d206c1467e526a5ea94a87 Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Thu, 4 Jul 2024 14:39:45 +0100 Subject: [PATCH 26/35] Fix empty occurs case. --- R/competing_hazards.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index d555bd9c..8acfc387 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -56,7 +56,7 @@ CompetingHazard <- R6::R6Class( occur_rng <- private$rng(private$size) occurs <- which(occur_rng < rate_to_prob(occur_rates)) - norm_probs <- event_rates[occurs,] / occur_rates[occurs] + norm_probs <- event_rates[occurs, , drop=FALSE] / occur_rates[occurs] norm_probs[is.na(norm_probs)] <- 0 cumulative <- rep(0, length(occurs)) From e2242b26ae26cbca33486392a00cda3df5be30bb Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Wed, 10 Jul 2024 13:44:59 +0100 Subject: [PATCH 27/35] Simplify competing hazards as a single random draw. --- R/competing_hazards.R | 18 ++- tests/testthat/test-competing-hazards.R | 183 ++++-------------------- 2 files changed, 33 insertions(+), 168 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index 8acfc387..b40ff175 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -52,22 +52,20 @@ CompetingHazard <- R6::R6Class( 'cbind', lapply(private$outcomes, function(x) x$rates) ) - occur_rates <- rowSums(event_rates) - occur_rng <- private$rng(private$size) - occurs <- which(occur_rng < rate_to_prob(occur_rates)) - norm_probs <- event_rates[occurs, , drop=FALSE] / occur_rates[occurs] - norm_probs[is.na(norm_probs)] <- 0 + total_rates <- rowSums(event_rates) + probs <- rate_to_prob(total_rates) * (event_rates / total_rates) + probs[is.na(probs)] <- 0 - cumulative <- rep(0, length(occurs)) - event_rng <- private$rng(length(occurs)) + rng <- private$rng(private$size) + cumulative <- rep(0, private$size) for (o in seq_along(private$outcomes)) { - next_cumulative <- cumulative + norm_probs[,o] - selected <- (event_rng > cumulative) & (event_rng <= next_cumulative) + next_cumulative <- cumulative + probs[,o] + selected <- which((rng > cumulative) & (rng <= next_cumulative)) cumulative <- next_cumulative - target <- individual::Bitset$new(private$size)$insert(occurs[selected]) + target <- individual::Bitset$new(private$size)$insert(selected) if (target$size() > 0){ private$outcomes[[o]]$execute(t, target) } diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R index 12203679..e1240d27 100644 --- a/tests/testthat/test-competing-hazards.R +++ b/tests/testthat/test-competing-hazards.R @@ -2,77 +2,32 @@ # Test the CompetingHazard class # ============================== -test_that('hazard resolves two normalised outcomes when all events occur', { - - size <- 4 - outcome_1_process <- mockery::mock() - outcome_1 <- CompetingOutcome$new( - targeted_process = outcome_1_process, - size = size - ) - outcome_2_process <- mockery::mock() - outcome_2 <- CompetingOutcome$new( - targeted_process = outcome_2_process, - size = size - ) - - hazard <- CompetingHazard$new( - outcomes = list(outcome_1, outcome_2), - rng = mockery::mock( - c(0, 0, 0, 0), # all events occur - c(.05, .3, .2, .5) # event_rng - ) - ) - - outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4)) - outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6)) - - hazard$resolve(0) - - mockery::expect_args( - outcome_1_process, - 1, - 0, - individual::Bitset$new(size)$insert(c(1, 3)) - ) - mockery::expect_args( - outcome_2_process, - 1, - 0, - individual::Bitset$new(size)$insert(c(2, 4)) - ) -}) -test_that('hazard resolves two unnormalised outcomes when all events occur', { - +test_that("hazard resolves two disjoint outcomes", { size <- 4 outcome_1_process <- mockery::mock() outcome_1 <- CompetingOutcome$new( targeted_process = outcome_1_process, size = size ) - outcome_2_process <- mockery::mock() outcome_2 <- CompetingOutcome$new( targeted_process = outcome_2_process, size = size ) - + hazard <- CompetingHazard$new( outcomes = list(outcome_1, outcome_2), - rng = mockery::mock( - c(0, 0, 0, 0), # all events occur - c(.05, .3, .2, .5) # event_rng - ) + rng = mockery::mock(c(.05, .3, .2, .5)) ) - outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) - outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + outcome_1$set_rates(c(10, 0, 10, 0)) + outcome_2$set_rates(c(0, 10, 0, 10)) hazard$resolve(0) - + mockery::expect_args( - outcome_1_process, + outcome_1_process, 1, 0, individual::Bitset$new(size)$insert(c(1, 3)) @@ -85,164 +40,76 @@ test_that('hazard resolves two unnormalised outcomes when all events occur', { ) }) -test_that('hazard resolves two outcomes when some events occur', { - +test_that("hazard resolves two competing outcomes", { size <- 4 outcome_1_process <- mockery::mock() outcome_1 <- CompetingOutcome$new( targeted_process = outcome_1_process, size = size ) - outcome_2_process <- mockery::mock() outcome_2 <- CompetingOutcome$new( targeted_process = outcome_2_process, size = size ) - + hazard <- CompetingHazard$new( outcomes = list(outcome_1, outcome_2), - rng = mockery::mock( - c(0, 0, 1, 1), # some events occur - c(.05, .3), # event_rng - ) + rng = mockery::mock(c(.7, .3, .2, .6)) ) - outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) - outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) + outcome_1$set_rates(c(5, 5, 5, 5)) + outcome_2$set_rates(c(5, 5, 5, 5)) hazard$resolve(0) - + mockery::expect_args( - outcome_1_process, + outcome_1_process, 1, 0, - individual::Bitset$new(size)$insert(1) + individual::Bitset$new(size)$insert(c(2, 3)) ) mockery::expect_args( outcome_2_process, 1, 0, - individual::Bitset$new(size)$insert(2) - ) -}) - - -test_that('hazard resolves when no rates are set', { - - size <- 4 - outcome_1_process <- mockery::mock() - outcome_1 <- CompetingOutcome$new( - targeted_process = outcome_1_process, - size = size - ) - - outcome_2_process <- mockery::mock() - outcome_2 <- CompetingOutcome$new( - targeted_process = outcome_2_process, - size = size - ) - - hazard <- CompetingHazard$new( - outcomes = list(outcome_1, outcome_2) - ) - - hazard$resolve(0) - - mockery::expect_called( - outcome_1_process, - 0 - ) - mockery::expect_called( - outcome_2_process, - 0 + individual::Bitset$new(size)$insert(c(1, 4)) ) }) -test_that('hazard resolves when rates are set to one', { - +test_that("hazard resolves partial outcomes", { size <- 4 outcome_1_process <- mockery::mock() outcome_1 <- CompetingOutcome$new( targeted_process = outcome_1_process, size = size ) - outcome_2_process <- mockery::mock() outcome_2 <- CompetingOutcome$new( targeted_process = outcome_2_process, size = size ) - - hazard <- CompetingHazard$new( - outcomes = list(outcome_1, outcome_2) - ) - - outcome_1$set_rates(rep(100, size)) - hazard$resolve(0) - mockery::expect_args( - outcome_1_process, - 1, - 0, - individual::Bitset$new(size)$insert(c(1, 2, 3, 4)) - ) - mockery::expect_called( - outcome_2_process, - 0 - ) -}) -test_that('hazard resolves three outcomes', { - - size <- 4 - outcome_1_process <- mockery::mock() - outcome_1 <- CompetingOutcome$new( - targeted_process = outcome_1_process, - size = size - ) - - outcome_2_process <- mockery::mock() - outcome_2 <- CompetingOutcome$new( - targeted_process = outcome_2_process, - size = size - ) - - outcome_3_process <- mockery::mock() - outcome_3 <- CompetingOutcome$new( - targeted_process = outcome_3_process, - size = size - ) - hazard <- CompetingHazard$new( - outcomes = list(outcome_1, outcome_2, outcome_3), - rng = mockery::mock( - c(0, 0, 0, 0), # all events occur - c(.04, .3, .14, .6), # event_rng - ) + outcomes = list(outcome_1, outcome_2), + rng = mockery::mock(c(.8, .4, .2, .6)) ) - outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2) - outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2) - outcome_3$set_rates(c(0.5, 0.5, 0.5, 0.5)) + outcome_1$set_rates(prob_to_rate(rep(0.5, size))) + outcome_2$set_rates(prob_to_rate(rep(0.5, size))) hazard$resolve(0) - + mockery::expect_args( - outcome_1_process, + outcome_1_process, 1, 0, - individual::Bitset$new(size)$insert(c(1, 3)) + individual::Bitset$new(size)$insert(c(3)) ) mockery::expect_args( outcome_2_process, 1, 0, - individual::Bitset$new(size)$insert(2) - ) - mockery::expect_args( - outcome_3_process, - 1, - 0, - individual::Bitset$new(size)$insert(4) + individual::Bitset$new(size)$insert(c(2, 4)) ) }) From 4c985591c86bafa4dec2e762984ee662b4a9b18b Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Wed, 10 Jul 2024 14:32:19 +0100 Subject: [PATCH 28/35] Try to optimize competing hazards --- R/competing_hazards.R | 56 ++++++++++++++------- R/disease_progression.R | 5 +- R/human_infection.R | 6 +-- R/processes.R | 3 +- tests/testthat/test-competing-hazards.R | 21 +++++--- tests/testthat/test-infection-integration.R | 2 +- tests/testthat/test-pev.R | 2 +- 7 files changed, 65 insertions(+), 30 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index b40ff175..d01bfc42 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -13,15 +13,25 @@ CompetingOutcome <- R6::R6Class( stop("size must be positive integer") } private$targeted_process <- targeted_process - self$rates <- rep(0, size) + + self$target <- individual::Bitset$new(size) + self$rates <- NULL }, - set_rates = function(rates){ + set_rates = function(target, rates){ + stopifnot(target$size() == length(rates)) + + # TODO: add an assign method to Bitset + self$target$or(target) self$rates <- rates }, execute = function(t, target){ private$targeted_process(t, target) - self$rates <- rep(0, length(self$rates)) }, + reset = function() { + self$target$clear() + self$rates <- NULL + }, + target = NULL, rates = NULL ) ) @@ -29,46 +39,58 @@ CompetingOutcome <- R6::R6Class( CompetingHazard <- R6::R6Class( "CompetingHazard", private = list( - outcomes = list(), size = NULL, + outcomes = list(), # RNG is passed in because mockery is not able to stub runif # TODO: change when fixed rng = NULL ), public = list( - initialize = function(outcomes, rng = runif){ + initialize = function(size, outcomes, rng = runif){ if (length(outcomes) == 0){ stop("At least one outcome must be provided") } if (!all(sapply(outcomes, function(x) inherits(x, "CompetingOutcome")))){ stop("All outcomes must be of class CompetingOutcome") } + private$size <- size private$outcomes <- outcomes - private$size <- length(outcomes[[1]]$rates) private$rng <- rng }, resolve = function(t){ - event_rates <- do.call( - 'cbind', - lapply(private$outcomes, function(x) x$rates) - ) + candidates <- individual::Bitset$new(private$size) + for (o in private$outcomes) { + candidates$or(o$target) + } + + rates <- matrix(ncol = length(private$outcomes), nrow = candidates$size(), 0) + for (i in seq_along(private$outcomes)) { + idx <- bitset_index( + candidates, + private$outcomes[[i]]$target) - total_rates <- rowSums(event_rates) - probs <- rate_to_prob(total_rates) * (event_rates / total_rates) + rates[idx, i] <- private$outcomes[[i]]$rates + } + + total_rates <- rowSums(rates) + probs <- rate_to_prob(total_rates) * (rates / total_rates) probs[is.na(probs)] <- 0 - rng <- private$rng(private$size) + rng <- private$rng(candidates$size()) + + cumulative <- rep(0, candidates$size()) - cumulative <- rep(0, private$size) for (o in seq_along(private$outcomes)) { next_cumulative <- cumulative + probs[,o] - selected <- which((rng > cumulative) & (rng <= next_cumulative)) + selected <- (rng > cumulative) & (rng <= next_cumulative) cumulative <- next_cumulative - target <- individual::Bitset$new(private$size)$insert(selected) - if (target$size() > 0){ + # TODO: change bitset_at to accept logical array + target <- bitset_at(candidates, which(selected)) + if (target$size() > 0) { private$outcomes[[o]]$execute(t, target) } + private$outcomes[[o]]$reset() } } ) diff --git a/R/disease_progression.R b/R/disease_progression.R index 11f63065..3269b6c5 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -10,7 +10,10 @@ create_recovery_rates_process <- function( recovery_outcome ) { function(timestep){ - recovery_outcome$set_rates(variables$recovery_rates$get_values()) + target <- variables$state$get_index_of(c("U", "Tr")) + recovery_outcome$set_rates( + target, + variables$recovery_rates$get_values(target)) } } diff --git a/R/human_infection.R b/R/human_infection.R index d8a39d90..6f16a7bb 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -121,9 +121,9 @@ 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) + infection_outcome$set_rates( + source_humans, + prob_to_rate(prob)) } #' @title Assigns infections to appropriate human states diff --git a/R/processes.R b/R/processes.R index c416a334..22793b7a 100644 --- a/R/processes.R +++ b/R/processes.R @@ -118,7 +118,8 @@ create_processes <- function( # Resolve competing hazards of infection with disease progression CompetingHazard$new( - outcomes = list(infection_outcome, recovery_outcome) + outcomes = list(infection_outcome, recovery_outcome), + size = parameters$human_population )$resolve ) diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R index e1240d27..8d12a176 100644 --- a/tests/testthat/test-competing-hazards.R +++ b/tests/testthat/test-competing-hazards.R @@ -5,6 +5,8 @@ test_that("hazard resolves two disjoint outcomes", { size <- 4 + population <- individual::Bitset$new(size)$not() + outcome_1_process <- mockery::mock() outcome_1 <- CompetingOutcome$new( targeted_process = outcome_1_process, @@ -17,12 +19,13 @@ test_that("hazard resolves two disjoint outcomes", { ) hazard <- CompetingHazard$new( + size = size, outcomes = list(outcome_1, outcome_2), rng = mockery::mock(c(.05, .3, .2, .5)) ) - outcome_1$set_rates(c(10, 0, 10, 0)) - outcome_2$set_rates(c(0, 10, 0, 10)) + outcome_1$set_rates(population, c(10, 0, 10, 0)) + outcome_2$set_rates(population, c(0, 10, 0, 10)) hazard$resolve(0) @@ -42,6 +45,8 @@ test_that("hazard resolves two disjoint outcomes", { test_that("hazard resolves two competing outcomes", { size <- 4 + population <- individual::Bitset$new(size)$not() + outcome_1_process <- mockery::mock() outcome_1 <- CompetingOutcome$new( targeted_process = outcome_1_process, @@ -54,12 +59,13 @@ test_that("hazard resolves two competing outcomes", { ) hazard <- CompetingHazard$new( + size = size, outcomes = list(outcome_1, outcome_2), rng = mockery::mock(c(.7, .3, .2, .6)) ) - outcome_1$set_rates(c(5, 5, 5, 5)) - outcome_2$set_rates(c(5, 5, 5, 5)) + outcome_1$set_rates(population, c(5, 5, 5, 5)) + outcome_2$set_rates(population, c(5, 5, 5, 5)) hazard$resolve(0) @@ -79,6 +85,8 @@ test_that("hazard resolves two competing outcomes", { test_that("hazard resolves partial outcomes", { size <- 4 + population <- individual::Bitset$new(size)$not() + outcome_1_process <- mockery::mock() outcome_1 <- CompetingOutcome$new( targeted_process = outcome_1_process, @@ -91,12 +99,13 @@ test_that("hazard resolves partial outcomes", { ) hazard <- CompetingHazard$new( + size = size, outcomes = list(outcome_1, outcome_2), rng = mockery::mock(c(.8, .4, .2, .6)) ) - outcome_1$set_rates(prob_to_rate(rep(0.5, size))) - outcome_2$set_rates(prob_to_rate(rep(0.5, size))) + outcome_1$set_rates(population, prob_to_rate(rep(0.5, size))) + outcome_2$set_rates(population, prob_to_rate(rep(0.5, size))) hazard$resolve(0) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index d399a527..d73b6657 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -707,7 +707,7 @@ test_that('prophylaxis is considered for medicated humans', { targeted_process = function(timestep, target){ infection_outcome_process(timestep, target, variables, renderer, parameters) }, - size = parameters$human_population + size = 4 ) infection_rates <- calculate_infections( diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R index 9a5f1175..f3c3d793 100644 --- a/tests/testthat/test-pev.R +++ b/tests/testthat/test-pev.R @@ -154,7 +154,7 @@ test_that('Infection considers pev efficacy', { targeted_process = function(timestep, target){ infection_process_resolved_hazard(timestep, target, variables, renderer, parameters) }, - size = parameters$human_population + size = 4 ) # remove randomness from infection sampling From 0bd7c10f83e14286f89f1ae5f92acd1cf63d578c Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Wed, 10 Jul 2024 15:07:25 +0100 Subject: [PATCH 29/35] All all states but S to create_recovery_rates_process --- R/disease_progression.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/disease_progression.R b/R/disease_progression.R index 3269b6c5..dde88bb2 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -10,7 +10,7 @@ create_recovery_rates_process <- function( recovery_outcome ) { function(timestep){ - target <- variables$state$get_index_of(c("U", "Tr")) + target <- variables$state$get_index_of("S")$not() recovery_outcome$set_rates( target, variables$recovery_rates$get_values(target)) From ded6876cfcd976fb7a3fd116e16b30df6a6957f0 Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Wed, 10 Jul 2024 15:32:38 +0100 Subject: [PATCH 30/35] Add a bitset_at implementation that works with a logical index. Eventually this should go into individual. --- R/RcppExports.R | 4 ++++ R/competing_hazards.R | 3 +-- R/utils.R | 2 ++ src/RcppExports.cpp | 13 +++++++++++++ src/utils.cpp | 22 ++++++++++++++++++++++ 5 files changed, 42 insertions(+), 2 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index ad658f2c..dcb14c70 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -101,6 +101,10 @@ bitset_index_cpp <- function(a, b) { .Call(`_malariasimulation_bitset_index_cpp`, a, b) } +bitset_at_logical_cpp <- function(source, filter) { + .Call(`_malariasimulation_bitset_at_logical_cpp`, source, filter) +} + fast_weighted_sample <- function(size, probs) { .Call(`_malariasimulation_fast_weighted_sample`, size, probs) } diff --git a/R/competing_hazards.R b/R/competing_hazards.R index d01bfc42..87bde320 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -85,8 +85,7 @@ CompetingHazard <- R6::R6Class( selected <- (rng > cumulative) & (rng <= next_cumulative) cumulative <- next_cumulative - # TODO: change bitset_at to accept logical array - target <- bitset_at(candidates, which(selected)) + target <- bitset_at_logical(candidates, selected) if (target$size() > 0) { private$outcomes[[o]]$execute(t, target) } diff --git a/R/utils.R b/R/utils.R index 0c98b099..1e0ab2b4 100644 --- a/R/utils.R +++ b/R/utils.R @@ -24,6 +24,8 @@ bernoulli_multi_p <- function(p) bernoulli_multi_p_cpp(p) #' @noRd bitset_index <- function(a, b) bitset_index_cpp(a$.bitset, b$.bitset) +bitset_at_logical <- function(a, b) individual::Bitset$new(from = bitset_at_logical_cpp(a$.bitset, b)) + #' @importFrom stats runif log_uniform <- function(size, rate) -rate * log(runif(size)) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index e2a6b9d6..c43084d6 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -336,6 +336,18 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// bitset_at_logical_cpp +Rcpp::XPtr bitset_at_logical_cpp(Rcpp::XPtr source, Rcpp::LogicalVector filter); +RcppExport SEXP _malariasimulation_bitset_at_logical_cpp(SEXP sourceSEXP, SEXP filterSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::XPtr >::type source(sourceSEXP); + Rcpp::traits::input_parameter< Rcpp::LogicalVector >::type filter(filterSEXP); + rcpp_result_gen = Rcpp::wrap(bitset_at_logical_cpp(source, filter)); + return rcpp_result_gen; +END_RCPP +} // fast_weighted_sample Rcpp::IntegerVector fast_weighted_sample(size_t size, std::vector probs); RcppExport SEXP _malariasimulation_fast_weighted_sample(SEXP sizeSEXP, SEXP probsSEXP) { @@ -376,6 +388,7 @@ static const R_CallMethodDef CallEntries[] = { {"_malariasimulation_random_restore_state", (DL_FUNC) &_malariasimulation_random_restore_state, 1}, {"_malariasimulation_bernoulli_multi_p_cpp", (DL_FUNC) &_malariasimulation_bernoulli_multi_p_cpp, 1}, {"_malariasimulation_bitset_index_cpp", (DL_FUNC) &_malariasimulation_bitset_index_cpp, 2}, + {"_malariasimulation_bitset_at_logical_cpp", (DL_FUNC) &_malariasimulation_bitset_at_logical_cpp, 2}, {"_malariasimulation_fast_weighted_sample", (DL_FUNC) &_malariasimulation_fast_weighted_sample, 2}, {"run_testthat_tests", (DL_FUNC) &run_testthat_tests, 0}, {NULL, NULL, 0} diff --git a/src/utils.cpp b/src/utils.cpp index d2b36043..ed5faf5f 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -43,6 +43,28 @@ std::vector bitset_index_cpp( return values; } +//[[Rcpp::export]] +Rcpp::XPtr bitset_at_logical_cpp( + Rcpp::XPtr source, + Rcpp::LogicalVector filter + ) { + if (source->size() != filter.size()) { + Rcpp::stop("vector of logicals must equal the size of the bitset"); + } + + individual_index_t result(source->max_size()); + + auto source_it = source->begin(); + auto filter_it = filter.begin(); + for (; source_it != source->end() && filter_it != filter.end(); ++source_it, ++filter_it) { + if (*filter_it) { + result.insert(*source_it); + } + } + + return Rcpp::XPtr(new individual_index_t(std::move(result)), true); +} + //[[Rcpp::export(rng = false)]] Rcpp::IntegerVector fast_weighted_sample( size_t size, From 8f9bcfc09b6eb03e4285a7f9cf6a7fddd6953a00 Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Wed, 10 Jul 2024 17:04:49 +0100 Subject: [PATCH 31/35] Fix tests --- tests/testthat/test-infection-integration.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index d73b6657..17972b6a 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -222,7 +222,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati targeted_process = function(timestep, target){ infection_outcome_process(timestep, target, variables, renderer, parameters) }, - size = parameters$human_population + size = 4 ) infections <- calculate_infections( From f363b94f135f2ac1c95bc401249a9bc616f6c86e Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Thu, 11 Jul 2024 14:01:03 +0100 Subject: [PATCH 32/35] Fix more tests --- tests/testthat/test-biting-integration.R | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R index 01138919..fbfd1fd4 100644 --- a/tests/testthat/test-biting-integration.R +++ b/tests/testthat/test-biting-integration.R @@ -30,9 +30,7 @@ test_that('biting_process integrates mosquito effects and human infection', { parameters, lagged_foim, lagged_eir, - 1, - 1, - infection_outcome + infection_outcome=infection_outcome ) bitten <- individual::Bitset$new(parameters$human_population) From b339b7a66b1ffc472511e9ad0381d0651a0d250c Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Fri, 12 Jul 2024 13:45:58 +0100 Subject: [PATCH 33/35] Add tests --- tests/testthat/test-competing-hazards.R | 124 +++++++++++++++++------- 1 file changed, 88 insertions(+), 36 deletions(-) diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R index 8d12a176..f7ada3fc 100644 --- a/tests/testthat/test-competing-hazards.R +++ b/tests/testthat/test-competing-hazards.R @@ -29,18 +29,10 @@ test_that("hazard resolves two disjoint outcomes", { hazard$resolve(0) - mockery::expect_args( - outcome_1_process, - 1, - 0, - individual::Bitset$new(size)$insert(c(1, 3)) - ) - mockery::expect_args( - outcome_2_process, - 1, - 0, - individual::Bitset$new(size)$insert(c(2, 4)) - ) + mockery::expect_args(outcome_1_process, 1, 0, + individual::Bitset$new(size)$insert(c(1, 3))) + mockery::expect_args(outcome_2_process, 1, 0, + individual::Bitset$new(size)$insert(c(2, 4))) }) test_that("hazard resolves two competing outcomes", { @@ -69,21 +61,13 @@ 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)) - ) - mockery::expect_args( - outcome_2_process, - 1, - 0, - individual::Bitset$new(size)$insert(c(1, 4)) - ) + mockery::expect_args(outcome_1_process, 1, 0, + individual::Bitset$new(size)$insert(c(2, 3))) + mockery::expect_args(outcome_2_process, 1, 0, + individual::Bitset$new(size)$insert(c(1, 4))) }) -test_that("hazard resolves partial outcomes", { +test_that("hazard may resolve to neither outcome", { size <- 4 population <- individual::Bitset$new(size)$not() @@ -109,16 +93,84 @@ test_that("hazard resolves partial outcomes", { hazard$resolve(0) - mockery::expect_args( - outcome_1_process, - 1, - 0, - individual::Bitset$new(size)$insert(c(3)) - ) - mockery::expect_args( - outcome_2_process, - 1, - 0, - individual::Bitset$new(size)$insert(c(2, 4)) + mockery::expect_args(outcome_1_process, 1, 0, + individual::Bitset$new(size)$insert(c(3))) + mockery::expect_args(outcome_2_process, 1, 0, + individual::Bitset$new(size)$insert(c(2, 4))) +}) + +test_that("outcomes can define a partial set of rates", { + size <- 4 + population <- individual::Bitset$new(size)$not() + + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + + hazard <- CompetingHazard$new( + size = size, + outcomes = list(outcome_1, outcome_2), + + # Only individuals 1, 3 and 4 get sampled + rng = mockery::mock(c(.2, .3, .6)) + ) + + outcome_1$set_rates(individual::Bitset$new(size)$insert(c(1,3)), c(5, 5)) + outcome_2$set_rates(individual::Bitset$new(size)$insert(c(1,4)), c(5, 5)) + + hazard$resolve(0) + + mockery::expect_args(outcome_1_process, 1, 0, + individual::Bitset$new(size)$insert(c(1, 3))) + mockery::expect_args(outcome_2_process, 1, 0, + individual::Bitset$new(size)$insert(c(4)) ) }) + +test_that("hazard resolves three competing outcomes", { + size <- 4 + population <- individual::Bitset$new(size)$not() + + outcome_1_process <- mockery::mock() + outcome_1 <- CompetingOutcome$new( + targeted_process = outcome_1_process, + size = size + ) + outcome_2_process <- mockery::mock() + outcome_2 <- CompetingOutcome$new( + targeted_process = outcome_2_process, + size = size + ) + outcome_3_process <- mockery::mock() + outcome_3 <- CompetingOutcome$new( + targeted_process = outcome_3_process, + size = size + ) + + hazard <- CompetingHazard$new( + size = size, + outcomes = list(outcome_1, outcome_2, outcome_3), + rng = mockery::mock(c(.1, .5, .8, .8)) + ) + + outcome_1$set_rates(population, c(5, 5, 5, 5)) + outcome_2$set_rates(population, c(5, 5, 5, 5)) + outcome_3$set_rates(population, c(5, 5, 5, 5)) + + hazard$resolve(0) + + mockery::expect_args(outcome_1_process, 1, 0, + individual::Bitset$new(size)$insert(c(1))) + mockery::expect_args(outcome_2_process, 1, 0, + individual::Bitset$new(size)$insert(c(2))) + mockery::expect_args(outcome_3_process, 1, 0, + individual::Bitset$new(size)$insert(c(3, 4))) +}) + From f2cacfba7d8651cf90d8a168076fa8500587fa5e Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Sat, 27 Jul 2024 09:49:57 +0100 Subject: [PATCH 34/35] Use new individual features --- DESCRIPTION | 2 +- R/RcppExports.R | 4 ---- R/competing_hazards.R | 5 ++--- R/utils.R | 2 -- src/RcppExports.cpp | 13 ------------- src/utils.cpp | 22 ---------------------- 6 files changed, 3 insertions(+), 45 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index ab68ccb0..c6d25b2d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -70,7 +70,7 @@ Remotes: mrc-ide/malariaEquilibrium, mrc-ide/individual Imports: - individual (>= 0.1.16), + individual (>= 0.1.17), malariaEquilibrium (>= 1.0.1), Rcpp, statmod, diff --git a/R/RcppExports.R b/R/RcppExports.R index dcb14c70..ad658f2c 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -101,10 +101,6 @@ bitset_index_cpp <- function(a, b) { .Call(`_malariasimulation_bitset_index_cpp`, a, b) } -bitset_at_logical_cpp <- function(source, filter) { - .Call(`_malariasimulation_bitset_at_logical_cpp`, source, filter) -} - fast_weighted_sample <- function(size, probs) { .Call(`_malariasimulation_fast_weighted_sample`, size, probs) } diff --git a/R/competing_hazards.R b/R/competing_hazards.R index 87bde320..13c522ad 100644 --- a/R/competing_hazards.R +++ b/R/competing_hazards.R @@ -20,8 +20,7 @@ CompetingOutcome <- R6::R6Class( set_rates = function(target, rates){ stopifnot(target$size() == length(rates)) - # TODO: add an assign method to Bitset - self$target$or(target) + self$target$copy_from(target) self$rates <- rates }, execute = function(t, target){ @@ -85,7 +84,7 @@ CompetingHazard <- R6::R6Class( selected <- (rng > cumulative) & (rng <= next_cumulative) cumulative <- next_cumulative - target <- bitset_at_logical(candidates, selected) + target <- bitset_at(candidates, selected) if (target$size() > 0) { private$outcomes[[o]]$execute(t, target) } diff --git a/R/utils.R b/R/utils.R index 1e0ab2b4..0c98b099 100644 --- a/R/utils.R +++ b/R/utils.R @@ -24,8 +24,6 @@ bernoulli_multi_p <- function(p) bernoulli_multi_p_cpp(p) #' @noRd bitset_index <- function(a, b) bitset_index_cpp(a$.bitset, b$.bitset) -bitset_at_logical <- function(a, b) individual::Bitset$new(from = bitset_at_logical_cpp(a$.bitset, b)) - #' @importFrom stats runif log_uniform <- function(size, rate) -rate * log(runif(size)) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index c43084d6..e2a6b9d6 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -336,18 +336,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// bitset_at_logical_cpp -Rcpp::XPtr bitset_at_logical_cpp(Rcpp::XPtr source, Rcpp::LogicalVector filter); -RcppExport SEXP _malariasimulation_bitset_at_logical_cpp(SEXP sourceSEXP, SEXP filterSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::XPtr >::type source(sourceSEXP); - Rcpp::traits::input_parameter< Rcpp::LogicalVector >::type filter(filterSEXP); - rcpp_result_gen = Rcpp::wrap(bitset_at_logical_cpp(source, filter)); - return rcpp_result_gen; -END_RCPP -} // fast_weighted_sample Rcpp::IntegerVector fast_weighted_sample(size_t size, std::vector probs); RcppExport SEXP _malariasimulation_fast_weighted_sample(SEXP sizeSEXP, SEXP probsSEXP) { @@ -388,7 +376,6 @@ static const R_CallMethodDef CallEntries[] = { {"_malariasimulation_random_restore_state", (DL_FUNC) &_malariasimulation_random_restore_state, 1}, {"_malariasimulation_bernoulli_multi_p_cpp", (DL_FUNC) &_malariasimulation_bernoulli_multi_p_cpp, 1}, {"_malariasimulation_bitset_index_cpp", (DL_FUNC) &_malariasimulation_bitset_index_cpp, 2}, - {"_malariasimulation_bitset_at_logical_cpp", (DL_FUNC) &_malariasimulation_bitset_at_logical_cpp, 2}, {"_malariasimulation_fast_weighted_sample", (DL_FUNC) &_malariasimulation_fast_weighted_sample, 2}, {"run_testthat_tests", (DL_FUNC) &run_testthat_tests, 0}, {NULL, NULL, 0} diff --git a/src/utils.cpp b/src/utils.cpp index ed5faf5f..d2b36043 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -43,28 +43,6 @@ std::vector bitset_index_cpp( return values; } -//[[Rcpp::export]] -Rcpp::XPtr bitset_at_logical_cpp( - Rcpp::XPtr source, - Rcpp::LogicalVector filter - ) { - if (source->size() != filter.size()) { - Rcpp::stop("vector of logicals must equal the size of the bitset"); - } - - individual_index_t result(source->max_size()); - - auto source_it = source->begin(); - auto filter_it = filter.begin(); - for (; source_it != source->end() && filter_it != filter.end(); ++source_it, ++filter_it) { - if (*filter_it) { - result.insert(*source_it); - } - } - - return Rcpp::XPtr(new individual_index_t(std::move(result)), true); -} - //[[Rcpp::export(rng = false)]] Rcpp::IntegerVector fast_weighted_sample( size_t size, From babc34b33a0758c04bd67067cc4be45ce8c2ad48 Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Sat, 27 Jul 2024 11:06:49 +0100 Subject: [PATCH 35/35] Fix output test --- tests/testthat/test-output.R | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/tests/testthat/test-output.R b/tests/testthat/test-output.R index 91d427e4..5cef37b9 100644 --- a/tests/testthat/test-output.R +++ b/tests/testthat/test-output.R @@ -17,15 +17,17 @@ test_that('Test age parameter function works', { sim <- run_simulation(timesteps = 1, parameters) - expect_true( - all( - paste0(rep(c("n_age", - "n", "n_inc", "p_inc", - "n","n_inc_clinical","p_inc_clinical", - "n","n_inc_severe","p_inc_severe", - "n","n_detect_lm","p_detect_lm","n_detect_pcr", - "ica_mean", "icm_mean","id_mean","ib_mean","iva_mean","ivm_mean"), each = 3),"_", - age_limits[-4]+rep(c(0,rep(c(1:3), each = 3), rep(4, 4), 5:10), each = 3),"_",age_limits[-1]-1+rep(c(0,rep(c(1:3), each = 3), rep(4, 4), 5:10), each = 3) - ) %in% - names(sim))) + prefixes <- c("n_age", + "n_inc", "p_inc", + "n_inc_clinical","p_inc_clinical", + "n_inc_severe","p_inc_severe", + "n_detect_lm","p_detect_lm","n_detect_pcr", + "ica_mean", "icm_mean","id_mean","ib_mean","iva_mean","ivm_mean") + offsets <- c(0, rep(1, 2), rep(2, 2), rep(3, 2), rep(4, 3), 5:10) + expect_equal(length(prefixes), length(offsets)) + + expected <- paste0(rep(prefixes, each = 3), + "_", age_limits[-4]+rep(offsets, each = 3), + "_", age_limits[-1]-1+rep(offsets, each = 3)) + expect_in(expected, names(sim)) })