Skip to content

Commit

Permalink
competing_hazards.R: Giovanni designed two objects: CompetingOutcome …
Browse files Browse the repository at this point in the history
…(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.
  • Loading branch information
RJSheppard committed May 8, 2024
1 parent a9b495c commit 903547a
Show file tree
Hide file tree
Showing 12 changed files with 952 additions and 387 deletions.
100 changes: 51 additions & 49 deletions R/biting_process.R
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand All @@ -74,8 +77,9 @@ simulate_bites <- function(
mixing = 1,
mixing_index = 1
) {

bitten_humans <- individual::Bitset$new(parameters$human_population)

human_infectivity <- variables$infectivity$get_values()
if (parameters$tbv) {
human_infectivity <- account_for_tbv(
Expand All @@ -86,21 +90,21 @@ simulate_bites <- function(
)
}
renderer$render('infectivity', mean(human_infectivity), timestep)

# Calculate pi (the relative biting rate for each human)
psi <- unique_biting_rate(age, parameters)
zeta <- variables$zeta$get_values()
.pi <- human_pi(zeta, psi)

# Get some indices for later
if (parameters$individual_mosquitoes) {
infectious_index <- variables$mosquito_state$get_index_of('Im')
susceptible_index <- variables$mosquito_state$get_index_of('Sm')
adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not(TRUE)
}

EIR <- 0

for (s_i in seq_along(parameters$species)) {
species_name <- parameters$species[[s_i]]
solver_states <- solvers[[s_i]]$get_states()
Expand All @@ -111,7 +115,7 @@ simulate_bites <- function(
f <- blood_meal_rate(s_i, Z, parameters)
a <- .human_blood_meal_rate(f, s_i, W, parameters)
lambda <- effective_biting_rates(a, .pi, p_bitten)

if (parameters$individual_mosquitoes) {
species_index <- variables$species$get_index_of(
parameters$species[[s_i]]
Expand All @@ -127,18 +131,18 @@ simulate_bites <- function(
} else {
n_infectious <- calculate_infectious_compartmental(solver_states)
}

# store the current population's EIR for later
lagged_eir[[mixing_index]][[s_i]]$save(n_infectious * a, timestep)

# calculated the EIR for this timestep after mixing
species_eir <- sum(
vnapply(
lagged_eir,
function(l) l[[s_i]]$get(timestep - parameters$de)
) * mixing
)

renderer$render(paste0('EIR_', species_name), species_eir, timestep)
EIR <- EIR + species_eir
expected_bites <- species_eir * mean(psi)
Expand All @@ -150,7 +154,7 @@ simulate_bites <- function(
)
}
}

infectivity <- vnapply(
lagged_infectivity,
function(l) l$get(timestep - parameters$delay_gam)
Expand All @@ -163,7 +167,7 @@ simulate_bites <- function(
renderer$render(paste0('FOIM_', species_name), foim, timestep)
mu <- death_rate(f, W, Z, s_i, parameters)
renderer$render(paste0('mu_', species_name), mu, timestep)

if (parameters$individual_mosquitoes) {
# update the ODE with stats for ovoposition calculations
aquatic_mosquito_model_update(
Expand All @@ -172,10 +176,10 @@ simulate_bites <- function(
f,
mu
)

# update the individual mosquitoes
susceptible_species_index <- susceptible_index$copy()$and(species_index)

biting_effects_individual(
variables,
foim,
Expand All @@ -197,7 +201,7 @@ simulate_bites <- function(
)
}
}

renderer$render('n_bitten', bitten_humans$size(), timestep)
bitten_humans
}
Expand All @@ -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) {
Expand Down Expand Up @@ -246,7 +248,7 @@ calculate_infectious_individual <- function(
species_index,
parameters
) {

infectious_index$copy()$and(species_index)$size()
}

Expand Down
78 changes: 78 additions & 0 deletions R/competing_hazards.R
Original file line number Diff line number Diff line change
@@ -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_probs <- do.call(
'cbind',
lapply(private$outcomes, function(x) x$rates)
)
occur_prob <- rowSums(event_probs)
occur_rng <- private$rng(private$size)
occurs <- occur_rng < occur_prob
norm_probs <- event_probs / occur_prob
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)
}
}
}
)
)
Loading

0 comments on commit 903547a

Please sign in to comment.