Skip to content

Commit

Permalink
Merge pull request #101 from mrc-ide/dev
Browse files Browse the repository at this point in the history 
Version 1.0.0
  • Loading branch information
@giovannic
giovannic authored May 12, 2021
2 parents 9bb47af + 009e5c3 commit 78743fb
Show file tree
Hide file tree
Showing 41 changed files with 1,150 additions and 772 deletions.
4 changes: 1 addition & 3 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,23 +8,21 @@ export(fun_params)
export(gamb_params)
export(get_correlation_parameters)
export(get_parameters)
export(parameterise_human_equilibrium)
export(parameterise_mosquito_equilibrium)
export(parameterise_total_M)
export(peak_season_offset)
export(remove_unused_jamie)
export(run_simulation)
export(run_simulation_with_repetitions)
export(set_bednets)
export(set_clinical_treatment)
export(set_drugs)
export(set_equilibrium)
export(set_mda)
export(set_rtss)
export(set_smc)
export(set_species)
export(set_spraying)
export(set_tbv)
export(translate_jamie)
importFrom(MASS,mvrnorm)
importFrom(Rcpp,sourceCpp)
importFrom(stats,dweibull)
Expand Down
25 changes: 12 additions & 13 deletions R/RcppExports.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,18 @@ create_adult_solver <- function(model, init) {
.Call(`_malariasimulation_create_adult_solver`, model, init)
}

create_history <- function(size, default_value) {
.Call(`_malariasimulation_create_history`, size, default_value)
}

history_at <- function(history, timestep) {
.Call(`_malariasimulation_history_at`, history, timestep)
}

history_push <- function(history, value, timestep) {
invisible(.Call(`_malariasimulation_history_push`, history, value, timestep))
}

carrying_capacity <- function(timestep, model_seasonality, g0, g, h, K0, R_bar) {
.Call(`_malariasimulation_carrying_capacity`, timestep, model_seasonality, g0, g, h, K0, R_bar)
}
Expand All @@ -25,19 +37,6 @@ rainfall <- function(t, g0, g, h) {
.Call(`_malariasimulation_rainfall`, t, g0, g, h)
}

#' @title Mosquito emergence process
#' @description Move mosquitos from NonExistent to Sm in line with the number of
#' pupals in the ODE models
#'
#' @param solvers a list of solver objects for each species of mosquito
#' @param state the variable for the mosquito state
#' @param species the variable for the mosquito species
#' @param species_names a vector of category names for the species variable
#' @param dpl the delay for pupal growth (in timesteps)
create_mosquito_emergence_process_cpp <- function(solvers, state, species, species_names, dpl) {
.Call(`_malariasimulation_create_mosquito_emergence_process_cpp`, solvers, state, species, species_names, dpl)
}

create_mosquito_model <- function(beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, g0, g, h, R_bar) {
.Call(`_malariasimulation_create_mosquito_model`, beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, g0, g, h, R_bar)
}
Expand Down
164 changes: 113 additions & 51 deletions R/biting_process.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,9 @@ create_biting_process <- function(
models,
variables,
events,
parameters
parameters,
lagged_foim,
lagged_eir
) {
function(timestep) {
# Calculate combined EIR
Expand All @@ -29,7 +31,9 @@ create_biting_process <- function(
events,
age,
parameters,
timestep
timestep,
lagged_foim,
lagged_eir
)

simulate_infection(
Expand All @@ -53,7 +57,9 @@ simulate_bites <- function(
events,
age,
parameters,
timestep
timestep,
lagged_foim,
lagged_eir
) {
bitten_humans <- individual::Bitset$new(parameters$human_population)

Expand All @@ -70,7 +76,8 @@ simulate_bites <- function(

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

# Get some indices for later
if (parameters$individual_mosquitoes) {
Expand All @@ -82,43 +89,22 @@ simulate_bites <- function(
EIR <- 0

for (s_i in seq_along(parameters$species)) {
if (!parameters$individual_mosquitoes) {
solver_states <- solver_get_states(solvers[[s_i]])
}

if (parameters$individual_mosquitoes) {
species_index <- variables$species$get_index_of(
parameters$species[[s_i]]
)$and(adult_index)
n_infectious <- infectious_index$copy()$and(species_index)$size()
} else {
n_infectious <- solver_states[[ADULT_ODE_INDICES['Im']]]
}

solver_states <- solver_get_states(solvers[[s_i]])
p_bitten <- prob_bitten(timestep, variables, s_i, parameters)

Q0 <- parameters$Q0[[s_i]]
Z <- average_p_repelled(p_bitten$prob_repelled, .pi, Q0)
W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, Q0)
renderer$render(
paste0('p_repelled_', parameters$species[[s_i]]),
Z,
timestep
)
renderer$render(
paste0('p_feed_survives_', parameters$species[[s_i]]),
W,
timestep
)
Z <- average_p_repelled(p_bitten$prob_repelled, .pi, Q0)
f <- blood_meal_rate(s_i, Z, parameters)

lambda <- effective_biting_rate(
.pi,
age,
lambda <- .effective_biting_rates(
s_i,
psi,
.pi,
zeta,
p_bitten,
f,
W,
Z,
f,
parameters
)

Expand All @@ -128,8 +114,26 @@ simulate_bites <- function(
timestep
)

n_bites <- rpois(1, n_infectious * sum(lambda))
EIR <- EIR + n_bites
if (parameters$individual_mosquitoes) {
species_index <- variables$species$get_index_of(
parameters$species[[s_i]]
)$and(adult_index)
n_infectious <- calculate_infectious_individual(
s_i,
variables,
infectious_index,
adult_index,
species_index,
parameters
)
} else {
n_infectious <- calculate_infectious_compartmental(solver_states)
}

lagged_eir[[s_i]]$save(n_infectious * sum(lambda), timestep)
species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
EIR <- EIR + species_eir / sum(psi)
n_bites <- rpois(1, species_eir)
if (n_bites > 0) {
bitten_humans$insert(
sample.int(
Expand All @@ -142,6 +146,7 @@ simulate_bites <- function(
}

foim <- calculate_foim(human_infectivity, lambda)
lagged_foim$save(foim, timestep)
renderer$render(paste0('FOIM_', s_i), foim, timestep)
mu <- death_rate(f, W, Z, s_i, parameters)
renderer$render(paste0('mu_', s_i), mu, timestep)
Expand All @@ -155,7 +160,7 @@ simulate_bites <- function(

biting_effects_individual(
variables,
foim,
lagged_foim$get(timestep - parameters$delay_gam),
events,
s_i,
susceptible_species_index,
Expand All @@ -168,7 +173,7 @@ simulate_bites <- function(
adult_mosquito_model_update(
models[[s_i]],
mu,
foim,
lagged_foim$get(timestep - parameters$delay_gam),
solver_states[[ADULT_ODE_INDICES['Sm']]],
f
)
Expand All @@ -180,25 +185,82 @@ simulate_bites <- function(
bitten_humans
}


# =================
# Utility functions
# =================

#' @title Calculate the effective biting rate for a species on each human given
#' vector control interventions
#' @description
#' Implemented from Griffin et al 2010 S2 page 6
#' @param .pi relative biting rate for each human
#' @param age of each human (timesteps)
#' @param species to model
#' @param p_bitten the probabilities of feeding given vector controls
#' @param f blood meal rate
#' @param W average probability of a successful bite
#' @param parameters of the model
#' @noRd
effective_biting_rate <- function(.pi, age, species, p_bitten, f, W, parameters) {
calculate_eir <- function(species, solvers, variables, parameters, timestep) {
lambda <- effective_biting_rates(species, variables, parameters, timestep)
infectious <- calculate_infectious(species, solvers, variables, parameters)
infectious * sum(lambda)
}

effective_biting_rates <- function(species, variables, parameters, timestep) {
age <- get_age(variables$birth$get_values(), timestep)
psi <- unique_biting_rate(age, parameters)
zeta <- variables$zeta$get_values()
p_bitten <- prob_bitten(timestep, variables, species, parameters)
.pi <- human_pi(zeta, psi)
Q0 <- parameters$Q0[[species]]
W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, Q0)
Z <- average_p_repelled(p_bitten$prob_repelled, .pi, Q0)
f <- blood_meal_rate(species, Z, parameters)
.effective_biting_rates(species, psi, .pi, zeta, p_bitten, W, Z, f, parameters)
}

.effective_biting_rates <- function(
species,
psi,
.pi,
zeta,
p_bitten,
W,
Z,
f,
parameters
) {
a <- human_blood_meal_rate(f, species, W, parameters)
a * .pi * p_bitten$prob_bitten / sum(.pi * p_bitten$prob_bitten_survives)
a * intervention_coefficient(p_bitten) * zeta * psi
}

calculate_infectious <- function(species, solvers, variables, parameters) {
if (parameters$individual_mosquitoes) {
adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not()
return(
calculate_infectious_individual(
species,
variables,
variables$mosquito_state$get_index_of('Im'),
adult_index,
variables$species$get_index_of(
parameters$species[[species]]
)$and(adult_index),
parameters
)
)
}
calculate_infectious_compartmental(solver_get_states(solvers[[species]]))
}

calculate_infectious_individual <- function(
species,
variables,
infectious_index,
adult_index,
species_index,
parameters
) {

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

calculate_infectious_compartmental <- function(solver_states) {
solver_states[[ADULT_ODE_INDICES['Im']]]
}

intervention_coefficient <- function(p_bitten) {
p_bitten$prob_bitten / sum(p_bitten$prob_bitten_survives)
}

human_pi <- function(zeta, psi) {
Expand Down
Loading

0 comments on commit 78743fb

Please sign in to comment.