diff --git a/NAMESPACE b/NAMESPACE index c31de7a4..46953fd3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,7 +1,7 @@ # Generated by roxygen2: do not edit by hand export(AL_params) -export(DHC_PQP_params) +export(DHA_PQP_params) export(SP_AQ_params) export(arab_params) export(fun_params) @@ -34,4 +34,5 @@ importFrom(stats,rexp) importFrom(stats,rnorm) importFrom(stats,rpois) importFrom(stats,runif) +importFrom(stats,weighted.mean) useDynLib(malariasimulation, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index 509518df..658e9057 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,37 +1,61 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -carrying_capacity <- function(timestep, model_seasonality, days_per_timestep, g0, g, h, K0, R_bar) { - .Call(`_malariasimulation_carrying_capacity`, timestep, model_seasonality, days_per_timestep, g0, g, h, K0, R_bar) +create_adult_mosquito_model <- function(growth_model, mu, tau, susceptible) { + .Call(`_malariasimulation_create_adult_mosquito_model`, growth_model, mu, tau, susceptible) } -rainfall <- function(t, days_per_timestep, g0, g, h) { - .Call(`_malariasimulation_rainfall`, t, days_per_timestep, g0, g, h) +adult_mosquito_model_update <- function(model, mu, foim, susceptible, f) { + invisible(.Call(`_malariasimulation_adult_mosquito_model_update`, model, mu, foim, susceptible, f)) +} + +create_adult_solver <- function(model, init) { + .Call(`_malariasimulation_create_adult_solver`, model, init) +} + +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) +} + +eggs_laid <- function(beta, mu, f) { + .Call(`_malariasimulation_eggs_laid`, beta, mu, f) +} + +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 odes a list of odes for each species of mosquito +#' @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(odes, state, species, species_names, dpl) { - .Call(`_malariasimulation_create_mosquito_emergence_process_cpp`, odes, state, species, species_names, dpl) +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) +} + +mosquito_model_update <- function(model, total_M, f, mum) { + invisible(.Call(`_malariasimulation_mosquito_model_update`, model, total_M, f, mum)) } -create_mosquito_model <- function(init, beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, days_per_timestep, g0, g, h, R_bar) { - .Call(`_malariasimulation_create_mosquito_model`, init, beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, days_per_timestep, g0, g, h, R_bar) +create_solver <- function(model, init) { + .Call(`_malariasimulation_create_solver`, model, init) } -mosquito_model_get_states <- function(model) { - .Call(`_malariasimulation_mosquito_model_get_states`, model) +solver_get_states <- function(solver) { + .Call(`_malariasimulation_solver_get_states`, solver) } -mosquito_model_step <- function(model, total_M) { - invisible(.Call(`_malariasimulation_mosquito_model_step`, model, total_M)) +solver_step <- function(solver) { + invisible(.Call(`_malariasimulation_solver_step`, solver)) } bernoulli_multi_p_cpp <- function(p) { diff --git a/R/biting_process.R b/R/biting_process.R index 737253c0..49200fb4 100644 --- a/R/biting_process.R +++ b/R/biting_process.R @@ -3,17 +3,28 @@ #' This is the biting process. It results in human and mosquito infection and #' mosquito death. #' @param renderer the model renderer object +#' @param solvers mosquito ode solvers +#' @param models mosquito ode models #' @param variables a list of all of the model variables #' @param events a list of all of the model events #' @param parameters model pararmeters #' @noRd -create_biting_process <- function(renderer, variables, events, parameters) { +create_biting_process <- function( + renderer, + solvers, + models, + variables, + events, + parameters + ) { function(timestep) { # Calculate combined EIR age <- get_age(variables$birth$get_values(), timestep) bitten_humans <- simulate_bites( renderer, + solvers, + models, variables, events, age, @@ -27,13 +38,23 @@ create_biting_process <- function(renderer, variables, events, parameters) { bitten_humans, age, parameters, - timestep + timestep, + renderer ) } } #' @importFrom stats rpois -simulate_bites <- function(renderer, variables, events, age, parameters, timestep) { +simulate_bites <- function( + renderer, + solvers, + models, + variables, + events, + age, + parameters, + timestep + ) { bitten_humans <- individual::Bitset$new(parameters$human_population) human_infectivity <- variables$infectivity$get_values() @@ -50,27 +71,47 @@ simulate_bites <- function(renderer, variables, events, age, parameters, timeste # Calculate pi (the relative biting rate for each human) psi <- unique_biting_rate(age, parameters) .pi <- human_pi(variables$zeta$get_values(), psi) - 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() + + # 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() + } + + EIR <- 0 for (s_i in seq_along(parameters$species)) { - # Calculate the probabilities of each human being bitten (given - # interventions) - species_index <- variables$species$get_index_of( - parameters$species[[s_i]] - ) + 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']]] + } + 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) + renderer$render( + paste0('p_repelled_', parameters$species[[s_i]]), + Z, + timestep + ) + renderer$render( + paste0('p_feed_survives_', parameters$species[[s_i]]), + W, + timestep + ) f <- blood_meal_rate(s_i, Z, parameters) - n_infectious <- infectious_index$copy()$and(species_index)$size() - lambda <- effective_biting_rate( .pi, age, @@ -81,115 +122,64 @@ simulate_bites <- function(renderer, variables, events, age, parameters, timeste parameters ) - renderer$render(paste0('lambda_', parameters$species[[s_i]]), mean(lambda), timestep) + renderer$render( + paste0('lambda_', parameters$species[[s_i]]), + mean(lambda), + timestep + ) n_bites <- rpois(1, n_infectious * sum(lambda)) + EIR <- EIR + n_bites if (n_bites > 0) { bitten_humans$insert( sample.int( parameters$human_population, n_bites, replace = TRUE, - prob=lambda + prob = lambda ) ) } - susceptible_species_index <- susceptible_index$copy()$and(species_index) - - calculate_mosquito_effects( - variables, - human_infectivity, - lambda, - events, - s_i, - susceptible_species_index, - species_index$and(adult_index), - W, - Z, - f, - parameters, - renderer, - timestep - ) + foim <- calculate_foim(human_infectivity, lambda) + 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) + + if (parameters$individual_mosquitoes) { + # update the ODE with stats for ovoposition calculations + mosquito_model_update(models[[s_i]], species_index$size(), f, mu) + + # update the individual mosquitoes + susceptible_species_index <- susceptible_index$copy()$and(species_index) + + biting_effects_individual( + variables, + foim, + events, + s_i, + susceptible_species_index, + species_index, + mu, + parameters, + timestep + ) + } else { + adult_mosquito_model_update( + models[[s_i]], + mu, + foim, + solver_states[[ADULT_ODE_INDICES['Sm']]], + f + ) + } } - renderer$render('EIR', bitten_humans$size(), timestep) + renderer$render('EIR', EIR, timestep) + renderer$render('n_bitten', bitten_humans$size(), timestep) bitten_humans } -#' @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 -#' @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 -#' @noRd -simulate_infection <- function( - variables, - events, - bitten_humans, - age, - parameters, - timestep - ) { - - if (bitten_humans$size() > 0) { - boost_immunity( - variables$ib, - bitten_humans, - variables$last_boosted_ib, - timestep, - parameters$ub - ) - } - - # Calculate Infected - infected_humans <- calculate_infections( - variables, - bitten_humans, - parameters, - timestep - ) - - clinical_infections <- calculate_clinical_infections( - variables, - infected_humans, - parameters, - timestep - ) - - if (parameters$severe_enabled) { - update_severe_disease( - timestep, - clinical_infections, - variables, - infected_humans, - parameters - ) - } - - treated <- calculate_treated( - variables, - clinical_infections, - events$recovery, - parameters, - timestep - ) - - schedule_infections( - events, - clinical_infections, - treated, - infected_humans, - parameters - ) -} - # ================= # Utility functions # ================= @@ -233,3 +223,17 @@ average_p_repelled <- function(p_repelled, .pi, Q0) { average_p_successful <- function(prob_bitten_survives, .pi, Q0) { (1 - Q0) + Q0 * sum(.pi * prob_bitten_survives) } + +# Unique biting rate (psi) for a human of a given age +unique_biting_rate <- function(age, parameters) { + 1 - parameters$rho * exp(- age / parameters$a0) +} + +#' @title Calculate the force of infection towards mosquitoes +#' +#' @param human_infectivity a vector of infectivities for each human +#' @param lambda a vector of biting rates for each human +#' @noRd +calculate_foim <- function(human_infectivity, lambda) { + sum(human_infectivity * lambda) +} diff --git a/R/drug_parameters.R b/R/drug_parameters.R index 9a22aed8..dde2b913 100644 --- a/R/drug_parameters.R +++ b/R/drug_parameters.R @@ -1,7 +1,7 @@ -#' @title Preset parameters for the DHC-PQP drug +#' @title Preset parameters for the DHA-PQP drug #' @description From SI of Commun. 5:5606 doi: 10.1038/ncomms6606 (2014) #' @export -DHC_PQP_params <- c(.95, 0.09434, 4.4, 28.1) +DHA_PQP_params <- c(.95, 0.09434, 4.4, 28.1) #' @title Preset parameters for the AL drug #' @description From SI of Commun. 5:5606 doi: 10.1038/ncomms6606 (2014) @@ -35,19 +35,44 @@ set_drugs <- function(parameters, drugs) { #' @title Parameterise clinical treatment #' #' @param parameters the model parameters -#' @param ft probability of seeking treatment -#' @param drugs vector of drugs indecies that are available -#' @param coverages vector of coverages for those drugs +#' @param drug the index of the drug to set as a treatment +#' @param timesteps vector of timesteps for each coverage change +#' @param coverages vector of coverages for this drug #' @export -set_clinical_treatment <- function(parameters, ft, drugs, coverages) { - parameters$ft <- ft - if (any(drugs < 1 | drugs > length(parameters$drug_efficacy))) { - stop('Drug indices are invalid') +set_clinical_treatment <- function(parameters, drug, timesteps, coverages) { + n_drugs <- length(parameters$drug_efficacy) + if (drug < 1 | drug > n_drugs) { + stop('Drug index is invalid, please set drugs using set_drugs') } - parameters$clinical_treatment_drugs <- drugs - if (!approx_sum(coverages, 1)) { - stop('Drug coverages do not sum to 1') + drug_index <- which(parameters$clinical_treatment_drugs == drug) + if (length(drug_index) == 0) { + drug_index <- length(parameters$clinical_treatment_drugs) + 1 + } + parameters$clinical_treatment_drugs[[drug_index]] <- drug + parameters$clinical_treatment_timesteps[[drug_index]] <- timesteps + parameters$clinical_treatment_coverages[[drug_index]] <- coverages + last_timestep <- max(unlist(parameters$clinical_treatment_timesteps)) + + for (t in seq(last_timestep)) { + if (sum(get_treatment_coverages(parameters, t)) > 1) { + stop('The sum of drug coverages cannot be greater than 1 at any timestep') + } } - parameters$clinical_treatment_coverages <- coverages parameters } + +get_treatment_coverages <- function(parameters, timestep) { + vnapply( + seq_along(parameters$clinical_treatment_drugs), + function(drug_index) { + previous <- which( + parameters$clinical_treatment_timesteps[[drug_index]] <= timestep + ) + if (length(previous) == 0) { + return(0) + } + last_set <- max(previous) + parameters$clinical_treatment_coverages[[drug_index]][[last_set]] + } + ) +} diff --git a/R/events.R b/R/events.R index 047f3408..f1d9790a 100644 --- a/R/events.R +++ b/R/events.R @@ -1,10 +1,10 @@ create_events <- function(parameters) { - list( + events <- list( # Human infection events clinical_infection = individual::TargetedEvent$new(parameters$human_population), asymptomatic_infection = individual::TargetedEvent$new(parameters$human_population), - # either clinical or asym infection - infection = individual::TargetedEvent$new(parameters$human_population), + # whether the infection is detected + detection = individual::TargetedEvent$new(parameters$human_population), subpatent_infection = individual::TargetedEvent$new(parameters$human_population), recovery = individual::TargetedEvent$new(parameters$human_population), @@ -18,12 +18,18 @@ create_events <- function(parameters) { smc_administer = individual::Event$new(), # Bednet events - throw_away_net = individual::TargetedEvent$new(parameters$human_population), - - # Mosquito events - mosquito_infection = individual::TargetedEvent$new(parameters$mosquito_limit), - mosquito_death = individual::TargetedEvent$new(parameters$mosquito_limit) + throw_away_net = individual::TargetedEvent$new(parameters$human_population) ) + if (parameters$individual_mosquitoes) { + events <- c( + events, + # Mosquito events + mosquito_infection = individual::TargetedEvent$new(parameters$mosquito_limit), + mosquito_death = individual::TargetedEvent$new(parameters$mosquito_limit) + ) + } + + events } initialise_events <- function(events, variables, parameters) { @@ -45,7 +51,6 @@ initialise_events <- function(events, variables, parameters) { 'U', parameters$du ) - initialise_progression( events$recovery, variables$state, @@ -53,11 +58,13 @@ initialise_events <- function(events, variables, parameters) { parameters$dt ) - incubating <- variables$mosquito_state$get_index_of('Pm') - events$mosquito_infection$schedule( - incubating, - log_uniform(incubating$size(), parameters$dem) - ) + if (parameters$individual_mosquitoes) { + incubating <- variables$mosquito_state$get_index_of('Pm') + events$mosquito_infection$schedule( + incubating, + log_uniform(incubating$size(), parameters$dem) + ) + } # Initialise interventions if (parameters$rtss) { @@ -142,6 +149,14 @@ attach_event_listeners <- function( parameters$dd ) ) + events$clinical_infection$add_listener( + create_clinical_incidence_renderer( + variables$birth, + parameters, + renderer + ) + ) + events$asymptomatic_infection$add_listener( create_progression_listener( events$subpatent_infection, @@ -155,7 +170,7 @@ attach_event_listeners <- function( ) ) - events$infection$add_listener( + events$detection$add_listener( create_incidence_renderer( variables$birth, variables$is_severe, @@ -164,17 +179,19 @@ attach_event_listeners <- function( ) ) - events$mosquito_infection$add_listener( - individual::update_category_listener(variables$mosquito_state, 'Im') - ) + if (parameters$individual_mosquitoes) { + events$mosquito_infection$add_listener( + individual::update_category_listener(variables$mosquito_state, 'Im') + ) - events$mosquito_death$add_listener( - function(timestep, target) { - variables$mosquito_state$queue_update('NonExistent', target) - events$mosquito_infection$clear_schedule(target) - renderer$render('mosquito_deaths', target$size(), timestep) - } - ) + events$mosquito_death$add_listener( + function(timestep, target) { + variables$mosquito_state$queue_update('NonExistent', target) + events$mosquito_infection$clear_schedule(target) + renderer$render('mosquito_deaths', target$size(), timestep) + } + ) + } if (parameters$bednets == 1) { events$throw_away_net$add_listener( diff --git a/R/human_infection.R b/R/human_infection.R index ca17a57c..db50e299 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -1,3 +1,96 @@ +#' @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 +#' @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 +#' @noRd +simulate_infection <- function( + variables, + events, + bitten_humans, + age, + parameters, + timestep, + renderer + ) { + if (bitten_humans$size() > 0) { + boost_immunity( + variables$ib, + bitten_humans, + variables$last_boosted_ib, + timestep, + parameters$ub + ) + } + + # Calculate Infected + infected_humans <- calculate_infections( + variables, + bitten_humans, + parameters, + 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 + ) + + if (parameters$severe_enabled) { + update_severe_disease( + timestep, + clinical_infections, + variables, + infected_humans, + parameters + ) + } + + + treated <- calculate_treated( + variables, + clinical_infections, + events$recovery, + events$detection, + parameters, + timestep, + renderer + ) + + renderer$render('n_treated', treated$size(), timestep) + + schedule_infections( + events, + clinical_infections, + treated, + infected_humans, + parameters + ) +} + #' @title Calculate overall infections for bitten humans #' @description #' Sample infected humans given prophylaxis and vaccination @@ -66,29 +159,10 @@ calculate_infections <- function( #' @param variables a list of all of the model variables #' @param infections bitset of infected humans #' @param parameters model parameters -#' @param timestep current timestep #' @noRd -calculate_clinical_infections <- function(variables, infections, parameters, timestep) { +calculate_clinical_infections <- function(variables, infections, parameters) { ica <- variables$ica$get_values(infections) icm <- variables$icm$get_values(infections) - - if (infections$size() > 0) { - boost_immunity( - variables$ica, - infections, - variables$last_boosted_ica, - timestep, - parameters$uc - ) - boost_immunity( - variables$id, - infections, - variables$last_boosted_id, - timestep, - parameters$ud - ) - } - phi <- clinical_immunity(ica, icm, parameters) bitset_at(infections, bernoulli_multi_p(phi)) } @@ -141,42 +215,48 @@ update_severe_disease <- function( #' @param recovery the recovery event #' @param parameters model parameters #' @param timestep the current timestep +#' @param renderer simulation renderer #' @noRd calculate_treated <- function( variables, clinical_infections, recovery, + detection, parameters, - timestep + timestep, + renderer ) { - if (length(parameters$clinical_treatment_coverages) == 0) { + treatment_coverages <- get_treatment_coverages(parameters, timestep) + ft <- sum(treatment_coverages) + + if (ft == 0) { return(individual::Bitset$new(parameters$human_population)) } - seek_treatment <- sample_bitset(clinical_infections, parameters$ft) + renderer$render('ft', ft, timestep) + seek_treatment <- sample_bitset(clinical_infections, ft) n_treat <- seek_treatment$size() - drugs <- parameters$clinical_treatment_drugs[ + drugs <- as.numeric(parameters$clinical_treatment_drugs[ sample.int( length(parameters$clinical_treatment_drugs), n_treat, - prob = parameters$clinical_treatment_coverages, + prob = treatment_coverages, replace = TRUE ) - ] + ]) successful <- bernoulli_multi_p(parameters$drug_efficacy[drugs]) treated_index <- bitset_at(seek_treatment, successful) # Update those who have been treated if (treated_index$size() > 0) { - successful_vector <- successful$to_vector() variables$state$queue_update('Tr', treated_index) variables$infectivity$queue_update( - parameters$cd * parameters$drug_rel_c[drugs[successful_vector]], + parameters$cd * parameters$drug_rel_c[drugs[successful]], treated_index ) variables$drug$queue_update( - drugs[successful_vector], + drugs[successful], treated_index ) variables$drug_time$queue_update( @@ -187,6 +267,7 @@ calculate_treated <- function( treated_index, log_uniform(treated_index$size(), parameters$dt) ) + detection$schedule(treated_index, 0) } treated_index } @@ -207,7 +288,9 @@ schedule_infections <- function( infections, parameters ) { - included <- events$infection$get_scheduled()$or(treated)$not() + included <- events$clinical_infection$get_scheduled()$or( + events$asymptomatic_infection$get_scheduled() + )$or(treated)$not() to_infect <- clinical_infections$and(included) to_infect_asym <- clinical_infections$not()$and(infections)$and(included) @@ -219,7 +302,7 @@ schedule_infections <- function( to_infect, infection_times ) - events$infection$schedule(to_infect, infection_times) + events$detection$schedule(to_infect, infection_times) events$subpatent_infection$clear_schedule(to_infect) events$recovery$clear_schedule(to_infect) } @@ -230,7 +313,7 @@ schedule_infections <- function( to_infect_asym, infection_times ) - events$infection$schedule(to_infect_asym, infection_times) + events$detection$schedule(to_infect_asym, infection_times) events$subpatent_infection$clear_schedule(to_infect_asym) events$recovery$clear_schedule(to_infect_asym) } @@ -310,11 +393,6 @@ asymptomatic_infectivity <- function(age, immunity, parameters) { parameters$cu + (parameters$cd - parameters$cu) * q ** parameters$gamma1 } -# Unique biting rate (psi) for a human of a given age -unique_biting_rate <- function(age, parameters) { - 1 - parameters$rho * exp(- age / parameters$a0) -} - # Implemented from Winskill 2017 - Supplementary Information page 4 blood_immunity <- function(ib, parameters) { parameters$b0 * ( diff --git a/R/model.R b/R/model.R index 08b863fa..c92ee386 100644 --- a/R/model.R +++ b/R/model.R @@ -29,14 +29,16 @@ run_simulation <- function(timesteps, parameters = NULL, correlations = NULL) { correlations, renderer ) - odes <- parameterise_ode(parameters) + vector_models <- parameterise_mosquito_models(parameters) + solvers <- parameterise_solvers(vector_models, parameters) individual::simulation_loop( processes = create_processes( renderer, variables, events, parameters, - odes, + vector_models, + solvers, correlations ), variables = variables, diff --git a/R/mortality_processes.R b/R/mortality_processes.R index 962fadea..5ec16c01 100644 --- a/R/mortality_processes.R +++ b/R/mortality_processes.R @@ -62,7 +62,7 @@ create_mortality_process <- function(variables, events, renderer, parameters) { } } - events$infection$clear_schedule(died) + events$detection$clear_schedule(died) events$clinical_infection$clear_schedule(died) events$asymptomatic_infection$clear_schedule(died) events$subpatent_infection$clear_schedule(died) diff --git a/R/mosquito_biology.R b/R/mosquito_biology.R index 7ba2a245..8235e388 100644 --- a/R/mosquito_biology.R +++ b/R/mosquito_biology.R @@ -3,33 +3,32 @@ #' "Modelling the impact of vector control interventions on Anopheles gambiae #' population dynamics" #' @param parameters model parameters +#' @param species the index of the species to find the equilibrium for #' @param foim equilibrium foim -#' @param m (optional) the total number of female adult mosquitos +#' @param m the total number of female adult mosquitos #' @noRd -initial_mosquito_counts <- function(parameters, foim = 0, m = NULL) { - if (is.null(m)) { - m <- parameters$total_M - } - omega <- calculate_omega(parameters) - n_E <- 2 * omega * parameters$mum * parameters$dl * ( +initial_mosquito_counts <- function(parameters, species, foim, m) { + omega <- calculate_omega(parameters, species) + mum <- parameters$mum[[species]] + n_E <- 2 * omega * mum * parameters$dl * ( 1. + parameters$dpl * parameters$mup ) * m - n_L <- 2 * parameters$mum * parameters$dl * ( + n_L <- 2 * mum * parameters$dl * ( 1. + parameters$dpl * parameters$mup ) * m - n_P <- 2 * parameters$dpl * parameters$mum * m + n_P <- 2 * parameters$dpl * mum * m - n_Sm <- m * parameters$mum / (foim + parameters$mum) + n_Sm <- m * parameters$mum / (foim + mum) - incubation_survival <- exp(-parameters$mum * parameters$dem) + incubation_survival <- exp(-mum * parameters$dem) - n_Pm <- m * foim / (foim + parameters$mum) * ( + n_Pm <- m * foim / (foim + mum) * ( 1. - incubation_survival ) - n_Im <- m * foim / (foim + parameters$mum) * incubation_survival + n_Im <- m * foim / (foim + mum) * incubation_survival c(n_E, n_L, n_P, n_Sm, n_Pm, n_Im) } @@ -40,18 +39,27 @@ initial_mosquito_counts <- function(parameters, foim = 0, m = NULL) { #' "Modelling the impact of vector control interventions on Anopheles gambiae #' population dynamics" #' @param parameters model parameters +#' @param species the index of the species to calculate for #' @noRd -calculate_omega <- function(parameters) { +calculate_omega <- function(parameters, species) { sub_omega <- parameters$gamma * parameters$ml / parameters$me - ( parameters$del / parameters$dl ) + ( (parameters$gamma - 1) * parameters$ml * parameters$del ) + mum <- parameters$mum[[species]] + + beta <- eggs_laid( + parameters$beta, + mum, + parameters$blood_meal_rates[[species]] + ) + -.5 * sub_omega + sqrt( .25 * sub_omega**2 + - .5 * parameters$gamma * parameters$beta * parameters$ml * parameters$del / - (parameters$me * parameters$mum * parameters$dl * ( + .5 * parameters$gamma * beta * parameters$ml * parameters$del / + (parameters$me * mum * parameters$dl * ( 1. + parameters$dpl * parameters$mup )) ) @@ -62,12 +70,13 @@ calculate_omega <- function(parameters) { #' "Modelling the impact of vector control interventions on Anopheles gambiae #' population dynamics" #' @param parameters model parameters +#' @param m number of adult mosquitoes +#' @param species index of the species to calculate for #' @noRd -calculate_carrying_capacity <- function(parameters) { - m <- parameters$total_M - omega <- calculate_omega(parameters) +calculate_carrying_capacity <- function(parameters, m, species) { + omega <- calculate_omega(parameters, species) - m * 2 * parameters$dl * parameters$mum * ( + m * 2 * parameters$dl * parameters$mum[[species]] * ( 1. + parameters$dpl * parameters$mup ) * parameters$gamma * (omega + 1) / ( omega / (parameters$ml * parameters$del) - ( @@ -82,7 +91,6 @@ calculate_carrying_capacity <- function(parameters) { calculate_R_bar <- function(parameters) { mean(vnapply(1:365, function(t) rainfall( t, - parameters$days_per_timestep, parameters$g0, parameters$g, parameters$h @@ -93,6 +101,7 @@ calculate_R_bar <- function(parameters) { #' #' @param parameters to work from #' @param EIR equilibrium to use, bites per person per year +#' @importFrom stats weighted.mean #' @noRd equilibrium_total_M <- function(parameters, EIR) { if (EIR == 0) { @@ -101,9 +110,10 @@ equilibrium_total_M <- function(parameters, EIR) { if (parameters$init_foim == 0) { stop('init_foim must be > 0 to calculate a non-zero equilibrium total_M') } + mum <- weighted.mean(parameters$mum, parameters$species_proportions) total_daily_eir <- EIR * parameters$human_population / 365 - lifetime <- parameters$init_foim * exp(-parameters$mum * parameters$dem) / ( - parameters$init_foim + parameters$mum + lifetime <- parameters$init_foim * exp(-mum * parameters$dem) / ( + parameters$init_foim + mum ) total_daily_eir / sum( parameters$species_proportions * parameters$blood_meal_rates * parameters$Q0 * lifetime @@ -116,59 +126,66 @@ equilibrium_total_M <- function(parameters, EIR) { #' @param parameters to work from #' @export peak_season_offset <- function(parameters) { - K0 <- calculate_carrying_capacity(parameters) - R_bar <- calculate_R_bar(parameters) + if (!parameters$model_seasonality) { + return(0) + } which.max(vnapply(seq(365), function(t) { - carrying_capacity( + rainfall( t, - parameters$model_seasonality, - parameters$days_per_timestep, parameters$g0, parameters$g, - parameters$h, - K0, - R_bar + parameters$h ) }))[[1]] } -#' @title Calculate the effects of biting on mosquito individuals +#' @title Calculate the death rate of mosquitoes given interventions +#' +#' @param f the feeding rate for this species of mosquito +#' @param W the mean probability that a mosquito feeds and survives +#' @param Z the mean probability that a mosquito is repelled +#' @param Z the mean probability that a mosquito is repelled +#' @noRd +death_rate <- function(f, W, Z, species, parameters) { + mum <- parameters$mum[[species]] + p1_0 <- exp(-mum * parameters$foraging_time) + gonotrophic_cycle <- get_gonotrophic_cycle(species, parameters) + p2 <- exp(-mum * gonotrophic_cycle) + p1 <- p1_0 * W / (1 - Z * p1_0) + -f * log(p1 * p2) +} + +get_gonotrophic_cycle <- function(v, parameters) { + f <- parameters$blood_meal_rates[[v]] + gonotrophic_cycle <- 1 / f - parameters$foraging_time +} + +#' @title Update the individual mosquito model after biting #' #' @param variables a list of variables in this simulation -#' @param human_infectivity the infectivity for each human -#' @param lambda the effective biting rate for this species on each human -#' @param mosquito_infection an event for mosquito infection +#' @param foim force of infection towards mosquitoes +#' @param events events in the simulation #' @param species the index of the species to calculate for #' @param susceptible_species the indices of susceptible mosquitos of the #' species -#' @param adult_species the indices of adult mosquitos of the -#' species -#' @param W the mean probability that a mosquito feeds and survives -#' @param Z the mean probability that a mosquito is repelled -#' @param f the feeding rate for this species of mosquito -#' @param renderer the model renderer object -#' @param timestep the current timestep +#' @param adult_species the indices of adult mosquitos of the species +#' @param mu the death rate of the current species #' @param parameters the model parameters +#' @param timestep the current timestep #' @noRd -calculate_mosquito_effects <- function( +biting_effects_individual <- function( variables, - human_infectivity, - lambda, + foim, events, species, susceptible_species, adult_species, - W, - Z, - f, + mu, parameters, - renderer, timestep ) { # deal with mosquito infections - lambda <- sum(human_infectivity * lambda) - renderer$render(paste0('FOIM_', species), lambda, timestep) - target <- sample_bitset(susceptible_species, lambda) + target <- sample_bitset(susceptible_species, foim) variables$mosquito_state$queue_update('Pm', target) events$mosquito_infection$schedule( target, @@ -176,18 +193,7 @@ calculate_mosquito_effects <- function( ) # deal with mosquito deaths - p1_0 <- exp(-parameters$mum * parameters$foraging_time) - gonotrophic_cycle <- get_gonotrophic_cycle(species, parameters) - p2 <- exp(-parameters$mum * gonotrophic_cycle) - p1 <- p1_0 * W / (1 - Z * p1_0) - mu <- -f * log(p1 * p2) died <- sample_bitset(adult_species, mu) - renderer$render(paste0('mu_', species), mu, timestep) events$mosquito_death$schedule(died, 0) } - -get_gonotrophic_cycle <- function(v, parameters) { - f <- parameters$blood_meal_rates[[v]] - gonotrophic_cycle <- 1 / f - parameters$foraging_time -} diff --git a/R/ode.R b/R/ode.R index ce229ff4..4711d064 100644 --- a/R/ode.R +++ b/R/ode.R @@ -1,16 +1,17 @@ ODE_INDICES <- c(E = 1, L = 2, P = 3) +ADULT_ODE_INDICES <- c(Sm = 4, Pm = 5, Im = 6) -parameterise_ode <- function(parameters) { +parameterise_mosquito_models <- function(parameters) { lapply( - parameters$species_proportions, - function(p) { + seq_along(parameters$species), + function(i) { + p <- parameters$species_proportions[[i]] m <- p * parameters$total_M - create_mosquito_model( - initial_mosquito_counts(parameters, 0, m)[ODE_INDICES], + growth_model <- create_mosquito_model( parameters$beta, parameters$del, parameters$me, - p * calculate_carrying_capacity(parameters), + calculate_carrying_capacity(parameters, m, i), parameters$gamma, parameters$dl, parameters$ml, @@ -18,26 +19,66 @@ parameterise_ode <- function(parameters) { parameters$mup, m, parameters$model_seasonality, - parameters$days_per_timestep, parameters$g0, parameters$g, parameters$h, calculate_R_bar(parameters) ) + + if (!parameters$individual_mosquitoes) { + susceptible <- initial_mosquito_counts( + parameters, + i, + parameters$init_foim, + m + )[ADULT_ODE_INDICES['Sm']] + return( + create_adult_mosquito_model( + growth_model, + parameters$mum[[i]], + parameters$dem, + susceptible * parameters$init_foim + ) + ) + } + growth_model } ) } -create_ode_rendering_process <- function(renderer, odes) { +parameterise_solvers <- function(models, parameters) { + lapply( + seq_along(models), + function(i) { + m <- parameters$species_proportions[[i]] * parameters$total_M + init <- initial_mosquito_counts(parameters, i, parameters$init_foim, m) + if (!parameters$individual_mosquitoes) { + return( + create_adult_solver(models[[i]], init) + ) + } + create_solver(models[[i]], init[ODE_INDICES]) + } + ) +} + +create_ode_rendering_process <- function(renderer, solvers, parameters) { + if (parameters$individual_mosquitoes) { + indices <- ODE_INDICES + } else { + indices <- c(ODE_INDICES, ADULT_ODE_INDICES) + } + function(timestep) { - counts <- rep(0, length(ODE_INDICES)) - for (ode in odes) { - row <- mosquito_model_get_states(ode) + counts <- rep(0, length(indices)) + for (i in seq_along(solvers)) { + row <- solver_get_states(solvers[[i]]) counts <- counts + row } - for (i in seq_along(ODE_INDICES)) { + + for (i in seq_along(indices)) { renderer$render( - paste0('mosquito_', names(ODE_INDICES)[[i]], '_count'), + paste0(names(indices)[[i]], '_count'), counts[[i]], timestep ) @@ -45,23 +86,15 @@ create_ode_rendering_process <- function(renderer, odes) { } } -#' @title Step mosquito ODE +#' @title Step mosquito solver #' @description calculates total_M per species and updates the vector ode #' -#' @param odes the models to step, one for each species -#' @param state the mosquito state variable -#' @param species the mosquito species variable -#' @param species_names the names of the mosquito species -#' @param renderer the model renderer +#' @param solvers for each species #' @noRd -create_ode_stepping_process <- function(odes, state, species, species_names, renderer) { +create_solver_stepping_process <- function(solvers) { function(timestep) { - adult <- state$get_index_of("NonExistent")$not() - for (s_i in seq_along(species_names)) { - total_M <- species$get_index_of(species_names[[s_i]])$and(adult)$size() - renderer$render(paste0('total_M_', s_i), total_M, timestep) - mosquito_model_step(odes[[s_i]], total_M) + for (solver in solvers) { + solver_step(solver) } - renderer$render('total_M', adult$size(), timestep) } } diff --git a/R/parameters.R b/R/parameters.R index aa641ac3..191cc0dd 100644 --- a/R/parameters.R +++ b/R/parameters.R @@ -61,8 +61,6 @@ #' * ad - scale parameter relating age to immunity #' * gammad - shape parameter relating age to immunity #' * d1 - minimum probability due to immunity -#' * dmin - minimum probability due to immunity NOTE: there appears to be a -#' mistake here! #' * id0 - scale parameter #' * kd - shape parameter #' @@ -232,11 +230,11 @@ #' * human_population - the number of humans to model #' * mosquito_limit - the maximum number of mosquitos to allow for in the #' simulation -#' * days_per_timestep - the number of days to model per timestep +#' * individual_mosquitoes - boolean whether adult mosquitoes are modelled +#' individually or compartmentaly #' #' @export get_parameters <- function(overrides = list()) { - days_per_timestep <- 1 parameters <- list( dd = 5, dt = 5, @@ -246,7 +244,7 @@ get_parameters <- function(overrides = list()) { dl = 3.72, dpl = .643, mup = .249, - mum = .249, #NOTE: set from sitefile + mum = .1253333, sigma_squared = 1.67, n_heterogeneity_groups = 5, # immunity decay rates @@ -275,8 +273,8 @@ get_parameters <- function(overrides = list()) { a0 = 8 * 365, rho = .85, # clinical immunity parameters - phi0 = .0749886, - phi1 = .0001191, + phi0 = .792, + phi1 = .00074, ic0 = 18.02366, kc = 2.36949, # severe disease immunity parameters @@ -300,7 +298,7 @@ get_parameters <- function(overrides = list()) { id0 = 1.577533, kd = .476614, # mortality parameters - average_age = 7663 / days_per_timestep, + average_age = 7663, v = .065, # NOTE: there are two definitions of this in the literature: one on line 124 and one in the parameters table pcm = .774368, pvm = .195768, @@ -354,9 +352,9 @@ get_parameters <- function(overrides = list()) { drug_rel_c = numeric(0), drug_prophilaxis_shape = numeric(0), drug_prophilaxis_scale = numeric(0), - clinical_treatment_drugs = numeric(0), - clinical_treatment_coverages = numeric(0), - ft = 0, + clinical_treatment_drugs = list(), + clinical_treatment_timesteps = list(), + clinical_treatment_coverages = list(), # rts,s rtss = FALSE, rtss_vmax = .93, @@ -407,6 +405,8 @@ get_parameters <- function(overrides = list()) { prevalence_rendering_max_ages = 10 * 365, incidence_rendering_min_ages = numeric(0), incidence_rendering_max_ages = numeric(0), + clinical_incidence_rendering_min_ages = numeric(0), + 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), @@ -414,7 +414,7 @@ get_parameters <- function(overrides = list()) { # misc human_population = 100, mosquito_limit = 100 * 1000, - days_per_timestep = days_per_timestep + individual_mosquitoes = TRUE ) # Override parameters with any client specified ones @@ -492,34 +492,42 @@ parameterise_mosquito_equilibrium <- function(parameters, EIR) { #' @export parameterise_total_M <- function(parameters, total_M) { parameters$total_M <- total_M - K0 <- calculate_carrying_capacity(parameters) - R_bar <- calculate_R_bar(parameters) - max_K <- max(vnapply(seq(365), function(t) { - carrying_capacity( - t, - parameters$model_seasonality, - parameters$days_per_timestep, - parameters$g0, - parameters$g, - parameters$h, - K0, - R_bar - ) - })) - omega <- calculate_omega(parameters) - max_total_M <- max_K * ( - 1 / ( - 2 * parameters$dl * parameters$mum * ( - 1 + parameters$dpl * parameters$mup + if (!parameters$individual_mosquitoes) { + return(parameters) + } + max_total_M <- 0 + for (i in seq_along(parameters$species)) { + species_M <- total_M * parameters$species_proportions[[i]] + K0 <- calculate_carrying_capacity(parameters, species_M, i) + R_bar <- calculate_R_bar(parameters) + max_K <- max(vnapply(seq(365), function(t) { + carrying_capacity( + t, + parameters$model_seasonality, + parameters$g0, + parameters$g, + parameters$h, + K0, + R_bar ) + })) + omega <- calculate_omega(parameters, i) + max_total_M <- max_total_M + max_K * ( + 1 / ( + 2 * parameters$dl * parameters$mum * ( + 1 + parameters$dpl * parameters$mup + ) + ) + ) * ( + 1 / ( + parameters$gamma * (omega + 1) + ) + ) * ( + omega / (parameters$ml * parameters$del) - ( + 1 / (parameters$ml * parameters$dl) + ) - 1 ) - ) * ( - 1 / ( - parameters$gamma * (omega + 1) - ) - ) * ( - omega / (parameters$ml * parameters$del) - (1 / (parameters$ml * parameters$dl)) - 1 - ) + } parameters$mosquito_limit <- ceiling(max_total_M * 5) #Allow for random fluctuations parameters } diff --git a/R/processes.R b/R/processes.R index a7ec7bc9..2e94da6c 100644 --- a/R/processes.R +++ b/R/processes.R @@ -7,7 +7,8 @@ #' @param variables a list of variables in the model #' @param events a list of events in the model #' @param parameters a list of model parameters -#' @param odes a list of vector ode models for each species +#' @param models a list of vector models, one for each species +#' @param solvers a list of ode solvers, one for each species #' @param correlations the intervention correlations object #' @noRd create_processes <- function( @@ -15,14 +16,14 @@ create_processes <- function( variables, events, parameters, - odes, + models, + solvers, correlations ) { + # ======== + # Immunity + # ======== processes <- list( - # ======== - # Immunity - # ======== - # Maternal immunity create_exponential_decay_process(variables$icm, parameters$rm), create_exponential_decay_process(variables$ivm, parameters$rvm), @@ -31,38 +32,55 @@ create_processes <- function( # Acquired immunity create_exponential_decay_process(variables$ica, parameters$rc), create_exponential_decay_process(variables$iva, parameters$rva), - create_exponential_decay_process(variables$id, parameters$rid), + create_exponential_decay_process(variables$id, parameters$rid) + ) - create_mosquito_emergence_process_cpp( - odes, - variables$mosquito_state$.variable, - variables$species$.variable, - parameters$species, - parameters$dpl - ), + if (parameters$individual_mosquitoes) { + processes <- c( + processes, + create_mosquito_emergence_process_cpp( + solvers, + variables$mosquito_state$.variable, + variables$species$.variable, + parameters$species, + parameters$dpl + ) + ) + } - # ============== - # Biting process - # ============== - # schedule infections for humans and set last_boosted_* - # move mosquitoes into incubating state - # kill mosquitoes caught in vector control - create_biting_process(renderer, variables, events, parameters), + # ============================== + # Biting and mortality processes + # ============================== + # schedule infections for humans and set last_boosted_* + # move mosquitoes into incubating state + # kill mosquitoes caught in vector control + processes <- c( + processes, + create_biting_process( + renderer, + solvers, + models, + variables, + events, + parameters + ), - create_mortality_process(variables, events, renderer, parameters), + create_mortality_process(variables, events, renderer, parameters) + ) - # =============== - # ODE integration - # =============== - create_ode_stepping_process( - odes, - variables$mosquito_state, - variables$species, - parameters$species, - renderer - ), + # =============== + # ODE integration + # =============== + processes <- c( + processes, + create_solver_stepping_process(solvers) + ) - # Rendering processes + # ========= + # Rendering + # ========= + processes <- c( + processes, individual::categorical_count_renderer_process( renderer, variables$state, @@ -77,18 +95,42 @@ create_processes <- function( variables$state, variables$birth, variables$is_severe, + variables$id, parameters, renderer ), - individual::categorical_count_renderer_process( - renderer, - variables$mosquito_state, - c('Sm', 'Pm', 'Im') - ), - - create_ode_rendering_process(renderer, odes) + create_ode_rendering_process(renderer, solvers, parameters) ) + if (parameters$individual_mosquitoes) { + processes <- c( + processes, + individual::categorical_count_renderer_process( + renderer, + variables$mosquito_state, + c('Sm', 'Pm', 'Im') + ), + create_total_M_renderer_individual( + variables$mosquito_state, + variables$species, + renderer, + parameters + ) + ) + } else { + processes <- c( + processes, + create_total_M_renderer_compartmental( + renderer, + solvers + ) + ) + } + + # ====================== + # Intervention processes + # ====================== + if (parameters$bednets) { processes <- c( processes, diff --git a/R/render.R b/R/render.R index 8357f0d8..17ab4759 100644 --- a/R/render.R +++ b/R/render.R @@ -14,18 +14,32 @@ create_prevelance_renderer <- function( state, birth, is_severe, + immunity, parameters, renderer ) { function(timestep) { - infected <- state$get_index_of(c('D', 'A')) - severe <- is_severe$get_index_of('yes')$and(infected) age <- get_age(birth$get_values(), timestep) + asymptomatic <- state$get_index_of('A') + asymptomatic_vector <- asymptomatic$to_vector() + asymptomatic_detected <- bernoulli_multi_p( + probability_of_detection( + age[asymptomatic_vector], + immunity$get_values(asymptomatic), + parameters + ) + ) + + detected <- state$get_index_of(c('Tr', 'D'))$or( + bitset_at(asymptomatic, asymptomatic_detected) + ) + + severe <- is_severe$get_index_of('yes')$and(detected) for (i in seq_along(parameters$prevalence_rendering_min_ages)) { lower <- parameters$prevalence_rendering_min_ages[[i]] upper <- parameters$prevalence_rendering_max_ages[[i]] p <- epi_from_subpopulation( - infected, + detected, age, lower, upper @@ -65,7 +79,7 @@ create_incidence_renderer <- function(birth, is_severe, parameters, renderer) { lower <- parameters$severe_incidence_rendering_min_ages[[i]] upper <- parameters$severe_incidence_rendering_max_ages[[i]] p <- epi_from_subpopulation( - severe, + severe$copy()$and(target), age, lower, upper @@ -75,6 +89,23 @@ create_incidence_renderer <- function(birth, is_severe, parameters, renderer) { } } +create_clinical_incidence_renderer <- function(birth, parameters, renderer) { + function(timestep, target) { + age <- get_age(birth$get_values(), timestep) + for (i in seq_along(parameters$clinical_incidence_rendering_min_ages)) { + lower <- parameters$clinical_incidence_rendering_min_ages[[i]] + upper <- parameters$clinical_incidence_rendering_max_ages[[i]] + p <- epi_from_subpopulation( + target, + age, + lower, + upper + ) + renderer$render(paste0('clin_inc_', lower, '_', upper), p, timestep) + } + } +} + create_variable_mean_renderer_process <- function( renderer, names, @@ -90,3 +121,37 @@ create_variable_mean_renderer_process <- function( } } } + +create_total_M_renderer_individual <- function( + mosquito_state, + species, + renderer, + parameters + ) { + function(timestep) { + adult <- mosquito_state$get_index_of('NonExistent')$not() + for (i in seq_along(parameters$species)) { + renderer$render( + paste0('total_M_', i), + species$get_index_of(parameters$species[[i]])$and(adult)$size(), + timestep + ) + } + + renderer$render('total_M', adult$size(), timestep) + } +} + +create_total_M_renderer_compartmental <- function(renderer, solvers) { + function(timestep) { + total_M <- 0 + for (i in seq_along(solvers)) { + row <- solver_get_states(solvers[[i]]) + species_M <- sum(row[ADULT_ODE_INDICES]) + total_M <- total_M + species_M + renderer$render(paste0('total_M_', i), species_M, timestep) + } + + renderer$render('total_M', total_M, timestep) + } +} diff --git a/R/utils.R b/R/utils.R index 7dcaa62f..15a6c865 100644 --- a/R/utils.R +++ b/R/utils.R @@ -12,10 +12,7 @@ bitset_at <- function(b, i) { individual::filter_bitset(b, i) } -#' @importFrom stats runif -bernoulli_multi_p <- function(p) { - individual::Bitset$new(from = bernoulli_multi_p_cpp(p)) -} +bernoulli_multi_p <- function(p) bernoulli_multi_p_cpp(p) #' @importFrom stats runif log_uniform <- function(size, rate) -rate * log(runif(size)) diff --git a/R/variables.R b/R/variables.R index 110a4cd0..e7e39fb5 100644 --- a/R/variables.R +++ b/R/variables.R @@ -190,37 +190,56 @@ create_variables <- function(parameters) { spray_time = spray_time ) - species <- individual::CategoricalVariable$new( - parameters$species, - sample( - parameters$species, - parameters$mosquito_limit, - prob = parameters$species_proportions, - replace = TRUE - ) - ) + # Add variables for individual mosquitoes + if (parameters$individual_mosquitoes) { + species_values <- rep(NA, parameters$mosquito_limit) + state_values <- rep(NA, parameters$mosquito_limit) + n_initialised <- 0 + for (i in seq_along(parameters$species)) { + mosquito_counts <- floor( + initial_mosquito_counts( + parameters, + i, + parameters$init_foim, + parameters$total_M * parameters$species_proportions[[i]] + ) + ) - mosquito_counts <- floor( - initial_mosquito_counts(parameters, parameters$init_foim) - ) + species_M <- sum(mosquito_counts[ADULT_ODE_INDICES]) - spare <- parameters$mosquito_limit - sum(mosquito_counts[c(4, 5, 6)]) + if (n_initialised + species_M > parameters$mosquito_limit) { + stop('Mosquito limit not high enough') + } - if (spare < 0) { - stop(paste('Mosquito limit not high enough. Short by', -spare, sep=' ')) - } + species_range <- seq(n_initialised, n_initialised + species_M) + species_values[species_range] <- parameters$species[[i]] + state_values[species_range] <- rep( + c('Sm', 'Pm', 'Im'), + times = mosquito_counts[ADULT_ODE_INDICES] + ) - mosquito_states <- c('Sm', 'Pm', 'Im', 'NonExistent') - mosquito_state <- individual::CategoricalVariable$new( - mosquito_states, - rep(mosquito_states, times = c(mosquito_counts[c(4, 5, 6)], spare)) - ) + n_initialised <- n_initialised + species_M + } - variables <- c( - variables, - species = species, - mosquito_state = mosquito_state - ) + # fill excess mosquitoes + species_values[is.na(species_values)] <- parameters$species[[1]] + state_values[is.na(state_values)] <- 'NonExistent' + + # initialise variables + species <- individual::CategoricalVariable$new( + parameters$species, + species_values + ) + mosquito_state <- individual::CategoricalVariable$new( + c('Sm', 'Pm', 'Im', 'NonExistent'), + state_values + ) + variables <- c( + variables, + species = species, + mosquito_state = mosquito_state + ) + } variables } diff --git a/R/vector_parameters.R b/R/vector_parameters.R index 32befa91..40a5026d 100644 --- a/R/vector_parameters.R +++ b/R/vector_parameters.R @@ -10,7 +10,8 @@ gamb_params <- list( dn0 = .533, phi_bednets = .89, rs = .2, - phi_indoors = .97 + phi_indoors = .97, + mum = .132 ) #' @title Preset parameters for the An. arabiensis vector @@ -25,7 +26,8 @@ arab_params <- list( dn0 = .533, phi_bednets = .9, rs = .2, - phi_indoors = .96 + phi_indoors = .96, + mum = .132 ) #' @title Preset parameters for the An. funestus vector @@ -40,7 +42,8 @@ fun_params <- list( dn0 = .533, phi_bednets = .9, rs = .2, - phi_indoors = .98 + phi_indoors = .98, + mum = .112 ) #' @title Parameterise the mosquito species to use in the model @@ -66,7 +69,8 @@ set_species <- function(parameters, species, proportions) { 'dn0', 'phi_bednets', 'rs', - 'phi_indoors' + 'phi_indoors', + 'mum' ) for (key in keys) { parameters[[key]] <- rep(NA, length(species)) diff --git a/man/DHC_PQP_params.Rd b/man/DHA_PQP_params.Rd similarity index 71% rename from man/DHC_PQP_params.Rd rename to man/DHA_PQP_params.Rd index 43f1d376..678dc95d 100644 --- a/man/DHC_PQP_params.Rd +++ b/man/DHA_PQP_params.Rd @@ -1,14 +1,14 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/drug_parameters.R \docType{data} -\name{DHC_PQP_params} -\alias{DHC_PQP_params} -\title{Preset parameters for the DHC-PQP drug} +\name{DHA_PQP_params} +\alias{DHA_PQP_params} +\title{Preset parameters for the DHA-PQP drug} \format{ An object of class \code{numeric} of length 4. } \usage{ -DHC_PQP_params +DHA_PQP_params } \description{ From SI of Commun. 5:5606 doi: 10.1038/ncomms6606 (2014) diff --git a/man/arab_params.Rd b/man/arab_params.Rd index a7642bf7..7a7dc1ab 100644 --- a/man/arab_params.Rd +++ b/man/arab_params.Rd @@ -5,7 +5,7 @@ \alias{arab_params} \title{Preset parameters for the An. arabiensis vector} \format{ -An object of class \code{list} of length 10. +An object of class \code{list} of length 11. } \usage{ arab_params diff --git a/man/create_mosquito_emergence_process_cpp.Rd b/man/create_mosquito_emergence_process_cpp.Rd index f2be9fff..71d7a606 100644 --- a/man/create_mosquito_emergence_process_cpp.Rd +++ b/man/create_mosquito_emergence_process_cpp.Rd @@ -4,10 +4,16 @@ \alias{create_mosquito_emergence_process_cpp} \title{Mosquito emergence process} \usage{ -create_mosquito_emergence_process_cpp(odes, state, species, species_names, dpl) +create_mosquito_emergence_process_cpp( + solvers, + state, + species, + species_names, + dpl +) } \arguments{ -\item{odes}{a list of odes for each species of mosquito} +\item{solvers}{a list of solver objects for each species of mosquito} \item{state}{the variable for the mosquito state} diff --git a/man/fun_params.Rd b/man/fun_params.Rd index 258d7851..1988e617 100644 --- a/man/fun_params.Rd +++ b/man/fun_params.Rd @@ -5,7 +5,7 @@ \alias{fun_params} \title{Preset parameters for the An. funestus vector} \format{ -An object of class \code{list} of length 10. +An object of class \code{list} of length 11. } \usage{ fun_params diff --git a/man/gamb_params.Rd b/man/gamb_params.Rd index a25408cc..bdac3df1 100644 --- a/man/gamb_params.Rd +++ b/man/gamb_params.Rd @@ -5,7 +5,7 @@ \alias{gamb_params} \title{Preset parameters for the An. gambiae s.s vector} \format{ -An object of class \code{list} of length 10. +An object of class \code{list} of length 11. } \usage{ gamb_params diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd index a660983d..6c040849 100644 --- a/man/get_parameters.Rd +++ b/man/get_parameters.Rd @@ -68,8 +68,6 @@ immunity reducing probability of detection: \item ad - scale parameter relating age to immunity \item gammad - shape parameter relating age to immunity \item d1 - minimum probability due to immunity -\item dmin - minimum probability due to immunity NOTE: there appears to be a -mistake here! \item id0 - scale parameter \item kd - shape parameter } @@ -254,7 +252,8 @@ miscellaneous: \item human_population - the number of humans to model \item mosquito_limit - the maximum number of mosquitos to allow for in the simulation -\item days_per_timestep - the number of days to model per timestep +\item individual_mosquitoes - boolean whether adult mosquitoes are modelled +individually or compartmentaly }} } \description{ diff --git a/man/set_clinical_treatment.Rd b/man/set_clinical_treatment.Rd index 0b35ff45..cfb469cf 100644 --- a/man/set_clinical_treatment.Rd +++ b/man/set_clinical_treatment.Rd @@ -4,16 +4,16 @@ \alias{set_clinical_treatment} \title{Parameterise clinical treatment} \usage{ -set_clinical_treatment(parameters, ft, drugs, coverages) +set_clinical_treatment(parameters, drug, timesteps, coverages) } \arguments{ \item{parameters}{the model parameters} -\item{ft}{probability of seeking treatment} +\item{drug}{the index of the drug to set as a treatment} -\item{drugs}{vector of drugs indecies that are available} +\item{timesteps}{vector of timesteps for each coverage change} -\item{coverages}{vector of coverages for those drugs} +\item{coverages}{vector of coverages for this drug} } \description{ Parameterise clinical treatment diff --git a/src/Random.cpp b/src/Random.cpp index fcebf9cd..183138ad 100644 --- a/src/Random.cpp +++ b/src/Random.cpp @@ -15,12 +15,12 @@ std::vector Random::bernoulli(size_t size, double p) { return Rcpp::as>(indices); } -individual_index_t Random::bernoulli_multi_p(const std::vector p) { - auto successes = individual_index_t(p.size()); +std::vector Random::bernoulli_multi_p(const std::vector p) { + auto successes = std::vector(); Rcpp::NumericVector r = Rcpp::runif(p.size()); for (auto i = 0u; i < p.size(); ++i) { if (r[i] < p[i]) { - successes.insert(i); + successes.push_back(i); } } return successes; diff --git a/src/Random.h b/src/Random.h index 3de791e5..771140d4 100644 --- a/src/Random.h +++ b/src/Random.h @@ -10,12 +10,11 @@ #define SRC_RANDOM_H_ #include -#include class RandomInterface { public: virtual std::vector bernoulli(size_t, double) = 0; - virtual individual_index_t bernoulli_multi_p(const std::vector) = 0; + virtual std::vector bernoulli_multi_p(const std::vector) = 0; virtual std::vector sample(size_t, size_t, bool) = 0; virtual ~RandomInterface() = default; }; @@ -27,7 +26,7 @@ class Random : public RandomInterface { return instance; } virtual std::vector bernoulli(size_t, double); - virtual individual_index_t bernoulli_multi_p(const std::vector); + virtual std::vector bernoulli_multi_p(const std::vector); virtual std::vector sample(size_t, size_t, bool); virtual ~Random() = default; Random(const Random &other) = delete; diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index f2ee52c1..0cdfd5d8 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -6,61 +6,111 @@ using namespace Rcpp; +// create_adult_mosquito_model +Rcpp::XPtr create_adult_mosquito_model(Rcpp::XPtr growth_model, double mu, double tau, double susceptible); +RcppExport SEXP _malariasimulation_create_adult_mosquito_model(SEXP growth_modelSEXP, SEXP muSEXP, SEXP tauSEXP, SEXP susceptibleSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::XPtr >::type growth_model(growth_modelSEXP); + Rcpp::traits::input_parameter< double >::type mu(muSEXP); + Rcpp::traits::input_parameter< double >::type tau(tauSEXP); + Rcpp::traits::input_parameter< double >::type susceptible(susceptibleSEXP); + rcpp_result_gen = Rcpp::wrap(create_adult_mosquito_model(growth_model, mu, tau, susceptible)); + return rcpp_result_gen; +END_RCPP +} +// adult_mosquito_model_update +void adult_mosquito_model_update(Rcpp::XPtr model, double mu, double foim, double susceptible, double f); +RcppExport SEXP _malariasimulation_adult_mosquito_model_update(SEXP modelSEXP, SEXP muSEXP, SEXP foimSEXP, SEXP susceptibleSEXP, SEXP fSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::XPtr >::type model(modelSEXP); + Rcpp::traits::input_parameter< double >::type mu(muSEXP); + Rcpp::traits::input_parameter< double >::type foim(foimSEXP); + Rcpp::traits::input_parameter< double >::type susceptible(susceptibleSEXP); + Rcpp::traits::input_parameter< double >::type f(fSEXP); + adult_mosquito_model_update(model, mu, foim, susceptible, f); + return R_NilValue; +END_RCPP +} +// create_adult_solver +Rcpp::XPtr create_adult_solver(Rcpp::XPtr model, std::vector init); +RcppExport SEXP _malariasimulation_create_adult_solver(SEXP modelSEXP, SEXP initSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::XPtr >::type model(modelSEXP); + Rcpp::traits::input_parameter< std::vector >::type init(initSEXP); + rcpp_result_gen = Rcpp::wrap(create_adult_solver(model, init)); + return rcpp_result_gen; +END_RCPP +} // carrying_capacity -double carrying_capacity(const size_t timestep, const bool model_seasonality, const double days_per_timestep, const double g0, const std::vector& g, const std::vector& h, const double K0, const double R_bar); -RcppExport SEXP _malariasimulation_carrying_capacity(SEXP timestepSEXP, SEXP model_seasonalitySEXP, SEXP days_per_timestepSEXP, SEXP g0SEXP, SEXP gSEXP, SEXP hSEXP, SEXP K0SEXP, SEXP R_barSEXP) { +double carrying_capacity(const size_t timestep, const bool model_seasonality, const double g0, const std::vector& g, const std::vector& h, const double K0, const double R_bar); +RcppExport SEXP _malariasimulation_carrying_capacity(SEXP timestepSEXP, SEXP model_seasonalitySEXP, SEXP g0SEXP, SEXP gSEXP, SEXP hSEXP, SEXP K0SEXP, SEXP R_barSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const size_t >::type timestep(timestepSEXP); Rcpp::traits::input_parameter< const bool >::type model_seasonality(model_seasonalitySEXP); - Rcpp::traits::input_parameter< const double >::type days_per_timestep(days_per_timestepSEXP); Rcpp::traits::input_parameter< const double >::type g0(g0SEXP); Rcpp::traits::input_parameter< const std::vector& >::type g(gSEXP); Rcpp::traits::input_parameter< const std::vector& >::type h(hSEXP); Rcpp::traits::input_parameter< const double >::type K0(K0SEXP); Rcpp::traits::input_parameter< const double >::type R_bar(R_barSEXP); - rcpp_result_gen = Rcpp::wrap(carrying_capacity(timestep, model_seasonality, days_per_timestep, g0, g, h, K0, R_bar)); + rcpp_result_gen = Rcpp::wrap(carrying_capacity(timestep, model_seasonality, g0, g, h, K0, R_bar)); + return rcpp_result_gen; +END_RCPP +} +// eggs_laid +double eggs_laid(double beta, double mu, double f); +RcppExport SEXP _malariasimulation_eggs_laid(SEXP betaSEXP, SEXP muSEXP, SEXP fSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< double >::type beta(betaSEXP); + Rcpp::traits::input_parameter< double >::type mu(muSEXP); + Rcpp::traits::input_parameter< double >::type f(fSEXP); + rcpp_result_gen = Rcpp::wrap(eggs_laid(beta, mu, f)); return rcpp_result_gen; END_RCPP } // rainfall -double rainfall(const size_t t, const double days_per_timestep, const double g0, const std::vector& g, const std::vector& h); -RcppExport SEXP _malariasimulation_rainfall(SEXP tSEXP, SEXP days_per_timestepSEXP, SEXP g0SEXP, SEXP gSEXP, SEXP hSEXP) { +double rainfall(const size_t t, const double g0, const std::vector& g, const std::vector& h); +RcppExport SEXP _malariasimulation_rainfall(SEXP tSEXP, SEXP g0SEXP, SEXP gSEXP, SEXP hSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const size_t >::type t(tSEXP); - Rcpp::traits::input_parameter< const double >::type days_per_timestep(days_per_timestepSEXP); Rcpp::traits::input_parameter< const double >::type g0(g0SEXP); Rcpp::traits::input_parameter< const std::vector& >::type g(gSEXP); Rcpp::traits::input_parameter< const std::vector& >::type h(hSEXP); - rcpp_result_gen = Rcpp::wrap(rainfall(t, days_per_timestep, g0, g, h)); + rcpp_result_gen = Rcpp::wrap(rainfall(t, g0, g, h)); return rcpp_result_gen; END_RCPP } // create_mosquito_emergence_process_cpp -Rcpp::XPtr create_mosquito_emergence_process_cpp(Rcpp::List odes, Rcpp::XPtr state, Rcpp::XPtr species, std::vector species_names, double dpl); -RcppExport SEXP _malariasimulation_create_mosquito_emergence_process_cpp(SEXP odesSEXP, SEXP stateSEXP, SEXP speciesSEXP, SEXP species_namesSEXP, SEXP dplSEXP) { +Rcpp::XPtr create_mosquito_emergence_process_cpp(Rcpp::List solvers, Rcpp::XPtr state, Rcpp::XPtr species, std::vector species_names, double dpl); +RcppExport SEXP _malariasimulation_create_mosquito_emergence_process_cpp(SEXP solversSEXP, SEXP stateSEXP, SEXP speciesSEXP, SEXP species_namesSEXP, SEXP dplSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::List >::type odes(odesSEXP); + Rcpp::traits::input_parameter< Rcpp::List >::type solvers(solversSEXP); Rcpp::traits::input_parameter< Rcpp::XPtr >::type state(stateSEXP); Rcpp::traits::input_parameter< Rcpp::XPtr >::type species(speciesSEXP); Rcpp::traits::input_parameter< std::vector >::type species_names(species_namesSEXP); Rcpp::traits::input_parameter< double >::type dpl(dplSEXP); - rcpp_result_gen = Rcpp::wrap(create_mosquito_emergence_process_cpp(odes, state, species, species_names, dpl)); + rcpp_result_gen = Rcpp::wrap(create_mosquito_emergence_process_cpp(solvers, state, species, species_names, dpl)); return rcpp_result_gen; END_RCPP } // create_mosquito_model -Rcpp::XPtr create_mosquito_model(std::vector init, double beta, double de, double mue, double K0, double gamma, double dl, double mul, double dp, double mup, size_t total_M, bool model_seasonality, double days_per_timestep, double g0, std::vector g, std::vector h, double R_bar); -RcppExport SEXP _malariasimulation_create_mosquito_model(SEXP initSEXP, SEXP betaSEXP, SEXP deSEXP, SEXP mueSEXP, SEXP K0SEXP, SEXP gammaSEXP, SEXP dlSEXP, SEXP mulSEXP, SEXP dpSEXP, SEXP mupSEXP, SEXP total_MSEXP, SEXP model_seasonalitySEXP, SEXP days_per_timestepSEXP, SEXP g0SEXP, SEXP gSEXP, SEXP hSEXP, SEXP R_barSEXP) { +Rcpp::XPtr create_mosquito_model(double beta, double de, double mue, double K0, double gamma, double dl, double mul, double dp, double mup, size_t total_M, bool model_seasonality, double g0, std::vector g, std::vector h, double R_bar); +RcppExport SEXP _malariasimulation_create_mosquito_model(SEXP betaSEXP, SEXP deSEXP, SEXP mueSEXP, SEXP K0SEXP, SEXP gammaSEXP, SEXP dlSEXP, SEXP mulSEXP, SEXP dpSEXP, SEXP mupSEXP, SEXP total_MSEXP, SEXP model_seasonalitySEXP, SEXP g0SEXP, SEXP gSEXP, SEXP hSEXP, SEXP R_barSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< std::vector >::type init(initSEXP); Rcpp::traits::input_parameter< double >::type beta(betaSEXP); Rcpp::traits::input_parameter< double >::type de(deSEXP); Rcpp::traits::input_parameter< double >::type mue(mueSEXP); @@ -72,39 +122,62 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double >::type mup(mupSEXP); Rcpp::traits::input_parameter< size_t >::type total_M(total_MSEXP); Rcpp::traits::input_parameter< bool >::type model_seasonality(model_seasonalitySEXP); - Rcpp::traits::input_parameter< double >::type days_per_timestep(days_per_timestepSEXP); Rcpp::traits::input_parameter< double >::type g0(g0SEXP); Rcpp::traits::input_parameter< std::vector >::type g(gSEXP); Rcpp::traits::input_parameter< std::vector >::type h(hSEXP); Rcpp::traits::input_parameter< double >::type R_bar(R_barSEXP); - rcpp_result_gen = Rcpp::wrap(create_mosquito_model(init, beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, days_per_timestep, g0, g, h, R_bar)); + rcpp_result_gen = Rcpp::wrap(create_mosquito_model(beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, g0, g, h, R_bar)); return rcpp_result_gen; END_RCPP } -// mosquito_model_get_states -std::vector mosquito_model_get_states(Rcpp::XPtr model); -RcppExport SEXP _malariasimulation_mosquito_model_get_states(SEXP modelSEXP) { +// mosquito_model_update +void mosquito_model_update(Rcpp::XPtr model, size_t total_M, double f, double mum); +RcppExport SEXP _malariasimulation_mosquito_model_update(SEXP modelSEXP, SEXP total_MSEXP, SEXP fSEXP, SEXP mumSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::XPtr >::type model(modelSEXP); + Rcpp::traits::input_parameter< size_t >::type total_M(total_MSEXP); + Rcpp::traits::input_parameter< double >::type f(fSEXP); + Rcpp::traits::input_parameter< double >::type mum(mumSEXP); + mosquito_model_update(model, total_M, f, mum); + return R_NilValue; +END_RCPP +} +// create_solver +Rcpp::XPtr create_solver(Rcpp::XPtr model, std::vector init); +RcppExport SEXP _malariasimulation_create_solver(SEXP modelSEXP, SEXP initSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< Rcpp::XPtr >::type model(modelSEXP); - rcpp_result_gen = Rcpp::wrap(mosquito_model_get_states(model)); + Rcpp::traits::input_parameter< std::vector >::type init(initSEXP); + rcpp_result_gen = Rcpp::wrap(create_solver(model, init)); return rcpp_result_gen; END_RCPP } -// mosquito_model_step -void mosquito_model_step(Rcpp::XPtr model, size_t total_M); -RcppExport SEXP _malariasimulation_mosquito_model_step(SEXP modelSEXP, SEXP total_MSEXP) { +// solver_get_states +std::vector solver_get_states(Rcpp::XPtr solver); +RcppExport SEXP _malariasimulation_solver_get_states(SEXP solverSEXP) { BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::XPtr >::type model(modelSEXP); - Rcpp::traits::input_parameter< size_t >::type total_M(total_MSEXP); - mosquito_model_step(model, total_M); + Rcpp::traits::input_parameter< Rcpp::XPtr >::type solver(solverSEXP); + rcpp_result_gen = Rcpp::wrap(solver_get_states(solver)); + return rcpp_result_gen; +END_RCPP +} +// solver_step +void solver_step(Rcpp::XPtr solver); +RcppExport SEXP _malariasimulation_solver_step(SEXP solverSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::XPtr >::type solver(solverSEXP); + solver_step(solver); return R_NilValue; END_RCPP } // bernoulli_multi_p_cpp -Rcpp::XPtr bernoulli_multi_p_cpp(const std::vector p); +std::vector bernoulli_multi_p_cpp(const std::vector p); RcppExport SEXP _malariasimulation_bernoulli_multi_p_cpp(SEXP pSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -118,12 +191,18 @@ END_RCPP RcppExport SEXP run_testthat_tests(); static const R_CallMethodDef CallEntries[] = { - {"_malariasimulation_carrying_capacity", (DL_FUNC) &_malariasimulation_carrying_capacity, 8}, - {"_malariasimulation_rainfall", (DL_FUNC) &_malariasimulation_rainfall, 5}, + {"_malariasimulation_create_adult_mosquito_model", (DL_FUNC) &_malariasimulation_create_adult_mosquito_model, 4}, + {"_malariasimulation_adult_mosquito_model_update", (DL_FUNC) &_malariasimulation_adult_mosquito_model_update, 5}, + {"_malariasimulation_create_adult_solver", (DL_FUNC) &_malariasimulation_create_adult_solver, 2}, + {"_malariasimulation_carrying_capacity", (DL_FUNC) &_malariasimulation_carrying_capacity, 7}, + {"_malariasimulation_eggs_laid", (DL_FUNC) &_malariasimulation_eggs_laid, 3}, + {"_malariasimulation_rainfall", (DL_FUNC) &_malariasimulation_rainfall, 4}, {"_malariasimulation_create_mosquito_emergence_process_cpp", (DL_FUNC) &_malariasimulation_create_mosquito_emergence_process_cpp, 5}, - {"_malariasimulation_create_mosquito_model", (DL_FUNC) &_malariasimulation_create_mosquito_model, 17}, - {"_malariasimulation_mosquito_model_get_states", (DL_FUNC) &_malariasimulation_mosquito_model_get_states, 1}, - {"_malariasimulation_mosquito_model_step", (DL_FUNC) &_malariasimulation_mosquito_model_step, 2}, + {"_malariasimulation_create_mosquito_model", (DL_FUNC) &_malariasimulation_create_mosquito_model, 15}, + {"_malariasimulation_mosquito_model_update", (DL_FUNC) &_malariasimulation_mosquito_model_update, 4}, + {"_malariasimulation_create_solver", (DL_FUNC) &_malariasimulation_create_solver, 2}, + {"_malariasimulation_solver_get_states", (DL_FUNC) &_malariasimulation_solver_get_states, 1}, + {"_malariasimulation_solver_step", (DL_FUNC) &_malariasimulation_solver_step, 1}, {"_malariasimulation_bernoulli_multi_p_cpp", (DL_FUNC) &_malariasimulation_bernoulli_multi_p_cpp, 1}, {"run_testthat_tests", (DL_FUNC) &run_testthat_tests, 0}, {NULL, NULL, 0} diff --git a/src/adult_mosquito_ode.cpp b/src/adult_mosquito_ode.cpp new file mode 100644 index 00000000..7fffe913 --- /dev/null +++ b/src/adult_mosquito_ode.cpp @@ -0,0 +1,91 @@ +/* + * mosquito_ode.cpp + * + * Created on: 11 Jun 2020 + * Author: gc1610 + */ + +#include +#include "adult_mosquito_ode.h" + +AdultMosquitoModel::AdultMosquitoModel( + MosquitoModel growth_model, + double mu, + double tau, + double incubating + ) : growth_model(growth_model), mu(mu), tau(tau) +{ + for (auto i = 0u; i < tau; ++i) { + lagged_incubating.push(incubating); + } +} + +integration_function_t create_ode(AdultMosquitoModel& model) { + //create original ode + auto growth_ode = create_ode(model.growth_model); + return [&model, growth_ode](const state_t& x, state_t& dxdt, double t) { + //set submodel total_M + model.growth_model.total_M = + x[get_idx(AdultODEState::S)] + + x[get_idx(AdultODEState::E)] + + x[get_idx(AdultODEState::I)]; + + //run the growth ode + growth_ode(x, dxdt, t); + + //run the adult ode + auto incubation_survival = exp(-model.mu * model.tau); + + dxdt[get_idx(AdultODEState::S)] = + .5 * x[get_idx(ODEState::P)] / model.growth_model.dp //growth to adult female + - x[get_idx(AdultODEState::S)] * model.foim //infections + - x[get_idx(AdultODEState::S)] * model.mu; //deaths + + dxdt[get_idx(AdultODEState::E)] = + x[get_idx(AdultODEState::S)] * model.foim //infections + - model.lagged_incubating.front() * incubation_survival //survived incubation period + - x[get_idx(AdultODEState::E)] * model.mu; // deaths + + dxdt[get_idx(AdultODEState::I)] = model.lagged_incubating.front() * incubation_survival //survived incubation period + - x[get_idx(AdultODEState::I)] * model.mu; // deaths + }; +} + +//[[Rcpp::export]] +Rcpp::XPtr create_adult_mosquito_model( + Rcpp::XPtr growth_model, + double mu, + double tau, + double susceptible + ) { + auto model = new AdultMosquitoModel(*growth_model, mu, tau, susceptible); + return Rcpp::XPtr(model, true); +} + +//[[Rcpp::export]] +void adult_mosquito_model_update( + Rcpp::XPtr model, + double mu, + double foim, + double susceptible, + double f + ) { + model->mu = mu; + model->foim = foim; + model->growth_model.f = f; + model->growth_model.mum = mu; + model->lagged_incubating.push(susceptible * foim); + if (model->lagged_incubating.size() > 0) { + model->lagged_incubating.pop(); + } +} + +//[[Rcpp::export]] +Rcpp::XPtr create_adult_solver( + Rcpp::XPtr model, + std::vector init) { + return Rcpp::XPtr( + new Solver(init, create_ode(*model)), + true + ); +} diff --git a/src/adult_mosquito_ode.h b/src/adult_mosquito_ode.h new file mode 100644 index 00000000..b97d0eb3 --- /dev/null +++ b/src/adult_mosquito_ode.h @@ -0,0 +1,32 @@ +/* + * adult_mosquito_ode.h + * + * Created on: 17 Mar 2021 + * Author: gc1610 + */ + +#ifndef SRC_ADULT_MOSQUITO_ODE_H_ +#define SRC_ADULT_MOSQUITO_ODE_H_ + +#include "mosquito_ode.h" + +/* + * The adult states are: + * 3 - S - Susceptible + * 4 - E - Incubating + * 5 - I - Infectious + */ +enum class AdultODEState : size_t {S = 3, E = 4, I = 5}; + +struct AdultMosquitoModel { + MosquitoModel growth_model; + std::queue lagged_incubating; //last tau values for incubating mosquitos + double mu; //death rate for adult female mosquitoes + const double tau; //extrinsic incubation period + double foim; //force of infection towards mosquitoes + AdultMosquitoModel(MosquitoModel, double, double, double); +}; + +integration_function_t create_ode(AdultMosquitoModel& model); + +#endif /* SRC_ADULT_MOSQUITO_ODE_H_ */ diff --git a/src/malariasimulation_types.h b/src/malariasimulation_types.h index a5e9db75..ca5c8f94 100644 --- a/src/malariasimulation_types.h +++ b/src/malariasimulation_types.h @@ -10,5 +10,7 @@ #include #include "mosquito_ode.h" +#include "adult_mosquito_ode.h" +#include "solver.h" #endif /* MALARIASIMULATION_TYPES_H_ */ diff --git a/src/mosquito_biology.cpp b/src/mosquito_biology.cpp index 1477306f..ef533f21 100644 --- a/src/mosquito_biology.cpp +++ b/src/mosquito_biology.cpp @@ -12,7 +12,6 @@ double carrying_capacity( const size_t timestep, const bool model_seasonality, - const double days_per_timestep, const double g0, const std::vector& g, const std::vector& h, @@ -20,16 +19,25 @@ double carrying_capacity( const double R_bar ) { if (model_seasonality) { - double r = rainfall(timestep, days_per_timestep, g0, g, h); - return std::max(K0 * r / R_bar, .01); + double r = rainfall(timestep, g0, g, h); + return std::max(K0 * r, .01) / R_bar; } return std::max(K0, .01); } +//[[Rcpp::export]] +double eggs_laid( + double beta, + double mu, + double f +) { + auto eov = beta / mu * (exp(mu / f) - 1); + return eov * mu * exp(-mu / f) / (1 - exp(-mu / f)); +} + //[[Rcpp::export]] double rainfall( const size_t t, - const double days_per_timestep, const double g0, const std::vector& g, const std::vector& h @@ -37,8 +45,8 @@ double rainfall( double result = g0; for (auto i = 0u; i < g.size(); ++i) { result += - g[i] * cos(2 * M_PI * t * days_per_timestep / 365 * (i + 1)) + - h[i] * sin(2 * M_PI * t * days_per_timestep / 365 * (i + 1)); + g[i] * cos(2 * M_PI * t * (i + 1) / 365) + + h[i] * sin(2 * M_PI * t * (i + 1) / 365); } - return result; + return std::max(result, 0.); } diff --git a/src/mosquito_biology.h b/src/mosquito_biology.h index 99cfef34..0fe33c61 100644 --- a/src/mosquito_biology.h +++ b/src/mosquito_biology.h @@ -14,7 +14,6 @@ double carrying_capacity( const size_t, const bool, const double, - const double, const std::vector&, const std::vector&, const double, @@ -24,9 +23,10 @@ double carrying_capacity( double rainfall( const size_t, const double, - const double, const std::vector&, const std::vector& ); +double eggs_laid(double, double, double); + #endif /* SRC_MOSQUITO_BIOLOGY_H_ */ diff --git a/src/mosquito_emergence.cpp b/src/mosquito_emergence.cpp index 6b46afbe..937fad47 100644 --- a/src/mosquito_emergence.cpp +++ b/src/mosquito_emergence.cpp @@ -13,14 +13,14 @@ //' @description Move mosquitos from NonExistent to Sm in line with the number of //' pupals in the ODE models //' -//' @param odes a list of odes for each species of mosquito +//' @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) //[[Rcpp::export]] Rcpp::XPtr create_mosquito_emergence_process_cpp( - Rcpp::List odes, + Rcpp::List solvers, Rcpp::XPtr state, Rcpp::XPtr species, std::vector species_names, @@ -30,8 +30,8 @@ Rcpp::XPtr create_mosquito_emergence_process_cpp( return Rcpp::XPtr( new process_t([=] (size_t t) { auto n = 0u; - for (Rcpp::XPtr ode : odes) { - n += ode->get_state()[get_idx(ODEState::P)] * rate; + for (Rcpp::XPtr solver: solvers) { + n += solver->state[get_idx(ODEState::P)] * rate; } auto source = state->get_index_of(std::vector{"NonExistent"}); if (source.size() < n) { @@ -46,8 +46,10 @@ Rcpp::XPtr create_mosquito_emergence_process_cpp( if (n > 0) { auto species_i = 0u; auto sourceit = source.begin(); - for (Rcpp::XPtr ode : odes) { - auto to_set = static_cast(ode->get_state()[get_idx(ODEState::P)] * rate); + for (Rcpp::XPtr solver : solvers) { + auto to_set = static_cast( + solver->state[get_idx(ODEState::P)] * rate + ); auto target = individual_index_t(state->size); for (auto i = 0u; i < to_set; ++i) { target.insert(*sourceit); diff --git a/src/mosquito_ode.cpp b/src/mosquito_ode.cpp index d9c217af..e710a6af 100644 --- a/src/mosquito_ode.cpp +++ b/src/mosquito_ode.cpp @@ -6,36 +6,37 @@ */ #include -#include #include "mosquito_ode.h" -#include integration_function_t create_ode(MosquitoModel& model) { - return [&model](const state_t& x , state_t& dxdt , double t) { + return [&model](const state_t& x, state_t& dxdt, double t) { auto K = carrying_capacity( t, model.model_seasonality, - model.days_per_timestep, model.g0, model.g, model.h, model.K0, model.R_bar ); - dxdt[get_idx(ODEState::E)] = model.beta * (model.total_M) //new eggs + auto beta = eggs_laid(model.beta, model.mum, model.f); + auto n_larvae = x[get_idx(ODEState::E)] + x[get_idx(ODEState::L)]; + + dxdt[get_idx(ODEState::E)] = beta * (model.total_M) //new eggs - x[get_idx(ODEState::E)] / model.de //growth to late larval stage - - x[get_idx(ODEState::E)] * model.mue * (1 + (x[get_idx(ODEState::E)] + x[get_idx(ODEState::L)]) / K); //early larval deaths - dxdt[1] = x[get_idx(ODEState::E)] / model.de //growth from early larval + - x[get_idx(ODEState::E)] * model.mue * (1 + n_larvae / K); //early larval deaths + + dxdt[get_idx(ODEState::L)] = x[get_idx(ODEState::E)] / model.de //growth from early larval - x[get_idx(ODEState::L)] / model.dl //growth to pupal - - x[get_idx(ODEState::L)] * model.mul * (1 + model.gamma * (x[get_idx(ODEState::E)] + x[get_idx(ODEState::L)]) / K); //late larval deaths - dxdt[2] = x[get_idx(ODEState::L)] / model.dl //growth to pupae + - x[get_idx(ODEState::L)] * model.mul * (1 + model.gamma * n_larvae / K); //late larval deaths + + dxdt[get_idx(ODEState::P)] = x[get_idx(ODEState::L)] / model.dl //growth to pupae - x[get_idx(ODEState::P)] / model.dp //growth to adult - x[get_idx(ODEState::P)] * model.mup; // death of pupae }; } MosquitoModel::MosquitoModel( - std::vector init, double beta, double de, double mue, @@ -47,7 +48,6 @@ MosquitoModel::MosquitoModel( double mup, size_t total_M, bool model_seasonality, - double days_per_timestep, double g0, std::vector g, std::vector h, @@ -64,36 +64,16 @@ MosquitoModel::MosquitoModel( mup(mup), total_M(total_M), model_seasonality(model_seasonality), - days_per_timestep(days_per_timestep), g0(g0), g(g), h(h), R_bar(R_bar) - { - for (auto i = 0u; i < state.size(); ++i) { - state[i] = init[i]; - } - ode = create_ode(*this); - rk = boost::numeric::odeint::make_dense_output( - a_tolerance, - r_tolerance, - boost::numeric::odeint::runge_kutta_dopri5() - ); -} + {} -void MosquitoModel::step(size_t new_total_M) { - total_M = new_total_M; - boost::numeric::odeint::integrate_adaptive(rk, ode, state, t, t + dt, dt); - ++t; -} -state_t MosquitoModel::get_state() { - return state; -} //[[Rcpp::export]] Rcpp::XPtr create_mosquito_model( - std::vector init, double beta, double de, double mue, @@ -105,14 +85,12 @@ Rcpp::XPtr create_mosquito_model( double mup, size_t total_M, bool model_seasonality, - double days_per_timestep, double g0, std::vector g, std::vector h, double R_bar ) { auto model = new MosquitoModel( - init, beta, de, mue, @@ -124,7 +102,6 @@ Rcpp::XPtr create_mosquito_model( mup, total_M, model_seasonality, - days_per_timestep, g0, g, h, @@ -134,12 +111,23 @@ Rcpp::XPtr create_mosquito_model( } //[[Rcpp::export]] -std::vector mosquito_model_get_states(Rcpp::XPtr model) { - auto state = model->get_state(); - return std::vector(state.cbegin(), state.cend()); +void mosquito_model_update( + Rcpp::XPtr model, + size_t total_M, + double f, + double mum + ) { + model->total_M = total_M; + model->f = f; + model->mum = mum; } //[[Rcpp::export]] -void mosquito_model_step(Rcpp::XPtr model, size_t total_M) { - model->step(total_M); +Rcpp::XPtr create_solver( + Rcpp::XPtr model, + std::vector init) { + return Rcpp::XPtr( + new Solver(init, create_ode(*model)), + true + ); } diff --git a/src/mosquito_ode.h b/src/mosquito_ode.h index 1867440f..9b4b3f20 100644 --- a/src/mosquito_ode.h +++ b/src/mosquito_ode.h @@ -8,17 +8,11 @@ #ifndef SRC_MOSQUITO_ODE_H_ #define SRC_MOSQUITO_ODE_H_ -#pragma GCC diagnostic ignored "-Wdeprecated-declarations" -// [[Rcpp::depends(BH)]] -#include -#pragma GCC diagnostic pop - -#include -#include +#include #include -#include #include #include "mosquito_biology.h" +#include "solver.h" /* * The states are: @@ -35,9 +29,6 @@ constexpr auto get_idx(T value) return static_cast>(value); } -using state_t = std::array; -using integration_function_t = std::function; - struct MosquitoModel { /* Parameters */ const double beta; //egg laying rate @@ -50,16 +41,15 @@ struct MosquitoModel { const double dp; //delay for for pupal growth const double mup; //death rate for pupae size_t total_M; //the number of adult female mosquitos in the model + double f; //biting rate + double mum; //adult mortality rate const bool model_seasonality; //whether to model seasonality - const double days_per_timestep; //scale of the fourier model for seasonality const double g0; //fourier shape parameter const std::vector g; //fourier shape parameters const std::vector h; //fourier shape parameters const double R_bar; //average rainfall - std::queue lagged_incubating; //last tau values for incubating mosquitos MosquitoModel( - std::vector init, double beta, double de, double mue, @@ -71,29 +61,12 @@ struct MosquitoModel { double mup, size_t total_M, bool model_seasonality, - double days_per_timestep, double g0, std::vector g, std::vector h, double R_bar ); - virtual void step(size_t); - virtual state_t get_state(); virtual ~MosquitoModel() {}; -private: - //solver fields - boost::numeric::odeint::dense_output_runge_kutta< - boost::numeric::odeint::controlled_runge_kutta< - boost::numeric::odeint::runge_kutta_dopri5 - > - >rk; - const double r_tolerance = 1.0e-6; - const double a_tolerance = 1.0e-6; - integration_function_t ode; - - double t = 0.; - const double dt = 1.; - state_t state; }; integration_function_t create_ode(MosquitoModel& model); diff --git a/src/solver.cpp b/src/solver.cpp new file mode 100644 index 00000000..1f0c3250 --- /dev/null +++ b/src/solver.cpp @@ -0,0 +1,43 @@ +/* + * mosquito_ode.cpp + * + * Created on: 11 Jun 2020 + * Author: gc1610 + */ + +#include +#include "solver.h" + +Solver::Solver( + const std::vector& init, + const integration_function_t& ode + ) : state(init), ode(ode) +{ + rk = boost::numeric::odeint::make_dense_output( + a_tolerance, + r_tolerance, + boost::numeric::odeint::runge_kutta_dopri5() + ); +} + +void Solver::step() { + boost::numeric::odeint::integrate_adaptive( + rk, + ode, + state, + t, + t + dt, + dt + ); + ++t; +} + +//[[Rcpp::export]] +std::vector solver_get_states(Rcpp::XPtr solver) { + return solver->state; +} + +//[[Rcpp::export]] +void solver_step(Rcpp::XPtr solver) { + solver->step(); +} diff --git a/src/solver.h b/src/solver.h new file mode 100644 index 00000000..112c2534 --- /dev/null +++ b/src/solver.h @@ -0,0 +1,42 @@ +/* + * solver.h + * + * Created on: 18 Mar 2020 + * Author: gc1610 + */ + +#ifndef SRC_SOLVER_H_ +#define SRC_SOLVER_H_ + +#pragma GCC diagnostic ignored "-Wdeprecated-declarations" +// [[Rcpp::depends(BH)]] +#include +#pragma GCC diagnostic pop + +#include +#include + +using state_t = std::vector; +using integration_function_t = std::function; + +struct Solver { + Solver( + const std::vector& init, + const integration_function_t& ode + ); + //solver fields + boost::numeric::odeint::dense_output_runge_kutta< + boost::numeric::odeint::controlled_runge_kutta< + boost::numeric::odeint::runge_kutta_dopri5 + > + >rk; + const double r_tolerance = 1.0e-6; + const double a_tolerance = 1.0e-6; + void step(); + double t = 0.; + const double dt = 1.; + state_t state; + integration_function_t ode; +}; + +#endif /* SRC_SOLVER_H_ */ diff --git a/src/test-emergence.cpp b/src/test-emergence.cpp index 78bbc2ef..8d5b1366 100644 --- a/src/test-emergence.cpp +++ b/src/test-emergence.cpp @@ -1,7 +1,7 @@ // All test files should include the // header file. -#include "mosquito_ode.h" +#include "solver.h" #include #include "test-mock.h" #include "mosquito_emergence.h" @@ -9,7 +9,9 @@ context("Emergence works") { test_that("emergence process fails when there are not enough individuals") { - MockODE ode; + //mock a solver + Solver solver(state_t{1000, 500, 100}, mock_integration); + //mock my variables auto state_values = std::vector(2000); for (auto n = 0; n < 2000; ++n) { @@ -25,20 +27,19 @@ context("Emergence works") { ); auto process = create_mosquito_emergence_process_cpp( - Rcpp::List::create(Rcpp::XPtr(&ode, false)), + Rcpp::List::create(Rcpp::XPtr(&solver, false)), Rcpp::XPtr(&state, false), Rcpp::XPtr(&species, false), {"gamb"}, 2 ); - REQUIRE_CALL(ode, get_state()) - .RETURN(state_t{1000, 500, 100}) - .TIMES(AT_LEAST(1)); CATCH_REQUIRE_THROWS((*process)(1)); } test_that("emergence process fails at the boundary") { - MockODE ode; + //mock a solver + Solver solver(state_t{1000, 500, 100}, mock_integration); + //mock my variables auto state_values = std::vector(2000); for (auto n = 0; n < 2000; ++n) { @@ -57,20 +58,18 @@ context("Emergence works") { ); auto process = create_mosquito_emergence_process_cpp( - Rcpp::List::create(Rcpp::XPtr(&ode, false)), + Rcpp::List::create(Rcpp::XPtr(&solver, false)), Rcpp::XPtr(&state, false), Rcpp::XPtr(&species, false), {"gamb"}, 2 ); - REQUIRE_CALL(ode, get_state()) - .RETURN(state_t{1000, 500, 100}) - .TIMES(AT_LEAST(1)); CATCH_REQUIRE_THROWS((*process)(1)); } test_that("emergence process creates the correct number of susceptibles") { - MockODE ode; + //mock a solver + Solver solver(state_t{1000, 500, 100}, mock_integration); //mock my variables auto state_values = std::vector(2000); @@ -90,7 +89,7 @@ context("Emergence works") { ); auto process = create_mosquito_emergence_process_cpp( - Rcpp::List::create(Rcpp::XPtr(&ode, false)), + Rcpp::List::create(Rcpp::XPtr(&solver, false)), Rcpp::XPtr(&state, false), Rcpp::XPtr(&species, false), {"gamb"}, @@ -103,16 +102,15 @@ context("Emergence works") { target.insert(i + 1000); } - REQUIRE_CALL(ode, get_state()) - .RETURN(state_t{1000, 500, 100}) - .TIMES(AT_LEAST(1)); REQUIRE_CALL(state, queue_update("Sm", target)); REQUIRE_CALL(species, queue_update("gamb", target)); (*process)(1); } test_that("emergence process works at the boundary") { - MockODE ode; + //mock a solver + Solver solver(state_t{1000, 500, 100}, mock_integration); + //mock my variables auto state_values = std::vector(2000); for (auto n = 0; n < 2000; ++n) { @@ -132,7 +130,7 @@ context("Emergence works") { ); auto process = create_mosquito_emergence_process_cpp( - Rcpp::List::create(Rcpp::XPtr(&ode, false)), + Rcpp::List::create(Rcpp::XPtr(&solver, false)), Rcpp::XPtr(&state, false), Rcpp::XPtr(&species, false), {"gamb"}, @@ -143,9 +141,6 @@ context("Emergence works") { target.insert(i + 1000); } - REQUIRE_CALL(ode, get_state()) - .RETURN(state_t{1000, 500, 100}) - .TIMES(AT_LEAST(1)); REQUIRE_CALL(state, queue_update("Sm", target)); REQUIRE_CALL(species, queue_update("gamb", target)); (*process)(1); diff --git a/src/test-mock.h b/src/test-mock.h index 99cf2fbd..dcc5ffad 100644 --- a/src/test-mock.h +++ b/src/test-mock.h @@ -9,7 +9,7 @@ #ifndef SRC_TEST_MOCK_H_ #define SRC_TEST_MOCK_H_ -#include "mosquito_ode.h" +#include "solver.h" #include #include #include "Random.h" @@ -22,34 +22,17 @@ struct MockCategory : public CategoricalVariable { MAKE_MOCK2(queue_update, void(const std::string, const individual_index_t&), override); }; -class MockODE : public MosquitoModel { -public: - MockODE() : MosquitoModel( - std::vector{0, 0, 0}, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - false, - 1, - 0, - std::vector{}, - std::vector{}, - 0 - ) {}; - MAKE_MOCK0(get_state, state_t(), override); -}; - class MockRandom : public RandomInterface { public: MAKE_MOCK2(bernoulli, std::vector(size_t, double), override); MAKE_MOCK3(sample, std::vector(size_t, size_t, bool), override); }; + +integration_function_t mock_integration = []( + const state_t&, + state_t&, + double t +) {}; + #endif /* SRC_TEST_MOCK_H_ */ diff --git a/src/utils.cpp b/src/utils.cpp index 794269e6..33cf634e 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -3,9 +3,10 @@ #include "Random.h" //[[Rcpp::export]] -Rcpp::XPtr bernoulli_multi_p_cpp(const std::vector p) { - return Rcpp::XPtr( - new individual_index_t(Random::get_instance().bernoulli_multi_p(p)), - true - ); +std::vector bernoulli_multi_p_cpp(const std::vector p) { + auto values = Random::get_instance().bernoulli_multi_p(p); + for (auto i = 0u; i < values.size(); ++i) { + values[i]++; + } + return values; } diff --git a/tests/performance/profile.R b/tests/performance/profile.R index 78896fc0..5d203f54 100644 --- a/tests/performance/profile.R +++ b/tests/performance/profile.R @@ -22,7 +22,7 @@ simparams <- get_parameters(c( eq <- human_equilibrium(EIR = eir, ft = ft, p = jamie_params, age = 0:99) simparams <- parameterise_human_equilibrium(simparams, eq) simparams <- parameterise_mosquito_equilibrium(simparams, EIR=eir) -simparams <- set_drugs(simparams, list(AL_params, DHC_PQP_params)) +simparams <- set_drugs(simparams, list(AL_params, DHA_PQP_params)) simparams <- set_clinical_treatment(simparams, ft, c(1, 2), c(.5, .5)) #profvis({output <- run_simulation(sim_length, simparams)}) diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R index edd54148..d57aa1ff 100644 --- a/tests/testthat/test-biology.R +++ b/tests/testthat/test-biology.R @@ -19,7 +19,7 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR', { ) )) parameters <- parameterise_mosquito_equilibrium(parameters, EIR) - m_eq <- initial_mosquito_counts(parameters, foim) + m_eq <- initial_mosquito_counts(parameters, 1, foim, parameters$total_M) #set up arguments for EIR calculation variables <- create_variables(parameters) @@ -73,7 +73,7 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR (with h ) )) parameters <- parameterise_mosquito_equilibrium(parameters, EIR) - m_eq <- initial_mosquito_counts(parameters, foim) + m_eq <- initial_mosquito_counts(parameters, 1, foim, parameters$total_M) #set up arguments for EIR calculation variables <- create_variables(parameters) @@ -141,32 +141,25 @@ test_that('mosquito_effects correctly samples mortalities and infections without parameters <- get_parameters() events <- create_events(parameters) variables <- create_variables(parameters) - infectivity <- c(.6, 0, .2, .3) - lambda <- c(.1, .2, .3, .4) - f <- parameters$blood_meal_rates[[1]] infected <- individual::Bitset$new(100)$insert(1:25) dead <- individual::Bitset$new(100)$insert(c(1:15, 51:60, 76:85)) bernoulli_mock = mockery::mock(infected, dead) - mockery::stub(calculate_mosquito_effects, 'sample_bitset', bernoulli_mock) + mockery::stub(biting_effects_individual, 'sample_bitset', bernoulli_mock) events$mosquito_infection <- mock_event(events$mosquito_infection) events$mosquito_death <- mock_event(events$mosquito_death) variables$mosquito_state <- mock_category( c('Sm', 'Pm', 'Im', 'NonExistent'), rep('Sm', 100) ) - calculate_mosquito_effects( + biting_effects_individual( variables, - infectivity, - lambda, + .5, events, 1, individual::Bitset$new(100)$insert(1:50), individual::Bitset$new(100)$insert(1:100), - 1, - 0, - f, + parameters$mum, parameters, - individual::Render$new(1), timestep = 1 ) diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R index c40e5c21..585056be 100644 --- a/tests/testthat/test-biting-integration.R +++ b/tests/testthat/test-biting-integration.R @@ -12,6 +12,8 @@ test_that('biting_process integrates mosquito effects and human infection', { biting_process <- create_biting_process( renderer, + list(), #solvers + list(), #models variables, events, parameters @@ -29,6 +31,8 @@ test_that('biting_process integrates mosquito effects and human infection', { bites_mock, 1, renderer, + list(), + list(), variables, events, age, @@ -44,7 +48,8 @@ test_that('biting_process integrates mosquito effects and human infection', { bitten, age, parameters, - timestep + timestep, + renderer ) }) @@ -74,14 +79,26 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects', lambda_mock <- mockery::mock(.5, .5, .5) mosquito_effects_mock <- mockery::mock() + ode_update <- mockery::mock() mockery::stub(simulate_bites, 'effective_biting_rate', lambda_mock) - mockery::stub(simulate_bites, 'calculate_mosquito_effects', mosquito_effects_mock) + mockery::stub(simulate_bites, 'biting_effects_individual', mosquito_effects_mock) mockery::stub(simulate_bites, 'rpois', mockery::mock(2)) mockery::stub(simulate_bites, 'sample.int', mockery::mock(c(2, 3))) .pi <- rep(1 / population, population) mockery::stub(simulate_bites, 'human_pi', mockery::mock(.pi)) - bitten <- simulate_bites(renderer, variables, events, age, parameters, timestep) + mockery::stub(simulate_bites, 'mosquito_model_update', ode_update) + ode_model <- mockery::mock() + bitten <- simulate_bites( + renderer, + list(), #solvers + list(ode_model), #models + variables, + events, + age, + parameters, + timestep + ) expect_equal(bitten$to_vector(), c(2, 3)) @@ -106,15 +123,14 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects', effects_args <- mockery::mock_args(mosquito_effects_mock) expect_equal(effects_args[[1]][[1]], variables) - expect_equal(effects_args[[1]][[2]], infectivity) - expect_equal(effects_args[[1]][[3]], .5) - expect_equal(effects_args[[1]][[4]], events) - expect_equal(effects_args[[1]][[5]], 1) - expect_equal(effects_args[[1]][[6]]$to_vector(), 11:25) - expect_equal(effects_args[[1]][[7]]$to_vector(), c(1:10, 11:25)) - expect_equal(effects_args[[1]][[8]], 1) - expect_equal(effects_args[[1]][[9]], 0) - expect_equal(effects_args[[1]][[10]], f) - expect_equal(effects_args[[1]][[11]], parameters) - expect_equal(effects_args[[1]][[12]], renderer) + expect_equal(effects_args[[1]][[2]], .55) + expect_equal(effects_args[[1]][[3]], events) + expect_equal(effects_args[[1]][[4]], 1) + expect_equal(effects_args[[1]][[5]]$to_vector(), 11:25) + expect_equal(effects_args[[1]][[6]]$to_vector(), c(1:10, 11:25)) + expect_equal(effects_args[[1]][[7]], parameters$mum) + expect_equal(effects_args[[1]][[8]], parameters) + expect_equal(effects_args[[1]][[9]], timestep) + + mockery::expect_args(ode_update, 1, ode_model, 25, f, parameters$mum) }) diff --git a/tests/testthat/test-emergence.R b/tests/testthat/test-emergence.R index dae540a0..a60221ca 100644 --- a/tests/testthat/test-emergence.R +++ b/tests/testthat/test-emergence.R @@ -3,7 +3,6 @@ test_that('carrying_capacity is calculated correctly', { carrying_capacity( 100, TRUE, - 1, 2, c(.3, .6, .9), c(.1, .4, .7), @@ -15,23 +14,6 @@ test_that('carrying_capacity is calculated correctly', { ) }) -test_that('carrying_capacity is takes into account the timescale', { - expect_equal( - carrying_capacity( - 100, - TRUE, - 5, - 2, - c(.3, .6, .9), - c(.1, .4, .7), - 10, - 2 - ), - 12.8, - tolerance = 1e-1 - ) -}) - test_that('carrying_capacity cycles every year', { time_points <- c(1, 30, 160, 240, 365) for (t in time_points) { @@ -40,7 +22,6 @@ test_that('carrying_capacity cycles every year', { carrying_capacity( t, TRUE, - 5, 2, c(.3, .6, .9), c(.1, .4, .7), @@ -50,7 +31,6 @@ test_that('carrying_capacity cycles every year', { carrying_capacity( t + 365 * y, TRUE, - 5, 2, c(.3, .6, .9), c(.1, .4, .7), @@ -72,7 +52,6 @@ test_that('carrying_capacity can avoid seasonality', { carrying_capacity( t + 365 * y, FALSE, - 5, 2, c(.3, .6, .9), c(.1, .4, .7), diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index 3a949678..4485f6ef 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -1,12 +1,20 @@ test_that('simulate_infection integrates different types of infection and scheduling', { population <- 8 timestep <- 5 - parameters <- get_parameters(list(human_population = population, severe_enabled = TRUE)) + parameters <- get_parameters(list( + human_population = population, + severe_enabled = TRUE + )) 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(c(.2, .3, .5, .9, .2, .3, .5, .9)) + ib = individual::DoubleVariable$new(immunity), + id = individual::DoubleVariable$new(immunity), + state = list(get_index_of = mockery::mock(asymptomatics)) ) total_eir <- 5 @@ -14,10 +22,14 @@ test_that('simulate_infection integrates different types of infection and schedu bitten <- individual::Bitset$new(population)$insert(c(1, 3, 5, 7)) boost_immunity_mock <- mockery::mock() - infection_mock <- mockery::mock(c(1, 3, 5)) - clinical_infection_mock <- mockery::mock(c(1, 3)) - severe_infection_mock <- mockery::mock(c(1)) - treated_mock <- mockery::mock(c(3)) + 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() mockery::stub(simulate_infection, 'boost_immunity', boost_immunity_mock) @@ -26,7 +38,15 @@ test_that('simulate_infection integrates different types of infection and schedu 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) - simulate_infection(variables, events, bitten, age, parameters, timestep) + simulate_infection( + variables, + events, + bitten, + age, + parameters, + timestep, + renderer + ) mockery::expect_args( boost_immunity_mock, @@ -51,18 +71,17 @@ test_that('simulate_infection integrates different types of infection and schedu clinical_infection_mock, 1, variables, - c(1, 3, 5), - parameters, - timestep + infected, + parameters ) mockery::expect_args( severe_infection_mock, 1, timestep, - c(1, 3), + clinical, variables, - c(1, 3, 5), + infected, parameters ) @@ -70,19 +89,21 @@ test_that('simulate_infection integrates different types of infection and schedu treated_mock, 1, variables, - c(1, 3), + clinical, events$recovery, + events$detection, parameters, - timestep + timestep, + renderer ) mockery::expect_args( schedule_mock, 1, events, - c(1, 3), - c(3), - c(1, 3, 5), + clinical, + treated, + infected, parameters ) }) @@ -90,8 +111,8 @@ test_that('simulate_infection integrates different types of infection and schedu test_that('calculate_infections works various combinations of drug and vaccination', { timestep <- 50 parameters <- get_parameters() - parameters <- set_drugs(parameters, list(AL_params, DHC_PQP_params)) - parameters <- set_clinical_treatment(parameters, .5, 2, 1) + parameters <- set_drugs(parameters, list(AL_params, DHA_PQP_params)) + parameters <- set_clinical_treatment(parameters, 2, 1, .5) variables <- list( state = individual::CategoricalVariable$new( @@ -113,7 +134,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati weibull_mock <- mockery::mock(.2) rtss_antibodies_mock <- mockery::mock(c(2, 3)) rtss_efficacy_mock <- mockery::mock(c(.2, .3)) - bernoulli_mock <- mockery::mock(individual::Bitset$new(4)$insert(2)) + bernoulli_mock <- mockery::mock(2) mockery::stub(calculate_infections, 'blood_immunity', immunity_mock) mockery::stub(calculate_infections, 'dweibull', weibull_mock) mockery::stub(calculate_infections, 'calculate_rtss_antibodies', rtss_antibodies_mock) @@ -181,7 +202,7 @@ test_that('calculate_clinical_infections correctly samples clinically infected', mockery::stub(calculate_clinical_infections, 'boost_immunity', boost_mock) mockery::stub(calculate_clinical_infections, 'clinical_immunity', immunity_mock) - bernoulli_mock <- mockery::mock(individual::Bitset$new(4)$insert(c(1, 3))) + bernoulli_mock <- mockery::mock(c(1, 3)) mockery::stub(calculate_clinical_infections, 'bernoulli_multi_p', bernoulli_mock) infections <- individual::Bitset$new(4)$insert(c(2, 3, 4)) @@ -189,32 +210,11 @@ test_that('calculate_clinical_infections correctly samples clinically infected', clinical_infections <- calculate_clinical_infections( variables, infections, - parameters, - timestep + parameters ) expect_equal(clinical_infections$to_vector(), c(2, 4)) - mockery::expect_args( - boost_mock, - 1, - variables$ica, - infections, - variables$last_boosted_ica, - 5, - parameters$uc - ) - - mockery::expect_args( - boost_mock, - 2, - variables$id, - infections, - variables$last_boosted_id, - 5, - parameters$ud - ) - mockery::expect_args( immunity_mock, 1, @@ -232,8 +232,9 @@ test_that('calculate_clinical_infections correctly samples clinically infected', test_that('calculate_treated correctly samples treated and updates the drug state', { parameters <- get_parameters() - parameters <- set_drugs(parameters, list(AL_params, DHC_PQP_params)) - parameters <- set_clinical_treatment(parameters, .5, c(1, 2), c(.5, .5)) + parameters <- set_drugs(parameters, list(AL_params, DHA_PQP_params)) + parameters <- set_clinical_treatment(parameters, 1, 1, .25) + parameters <- set_clinical_treatment(parameters, 2, 1, .25) timestep <- 5 events <- create_events(parameters) variables <- list( @@ -243,8 +244,10 @@ test_that('calculate_treated correctly samples treated and updates the drug stat drug_time = list(queue_update = mockery::mock()) ) - schedule_mock <- mockery::mock() - mockery::stub(calculate_treated, 'recovery$schedule', schedule_mock) + recovery_mock <- mockery::mock() + detection_mock <- mockery::mock() + mockery::stub(calculate_treated, 'recovery$schedule', recovery_mock) + mockery::stub(calculate_treated, 'detection$schedule', detection_mock) seek_treatment <- individual::Bitset$new(4)$insert(c(1, 2, 4)) mockery::stub( @@ -254,7 +257,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat ) sample_mock <- mockery::mock(c(2, 1, 1, 1)) mockery::stub(calculate_treated, 'sample.int', sample_mock) - bernoulli_mock <- mockery::mock(individual::Bitset$new(4)$insert(c(1, 3))) + bernoulli_mock <- mockery::mock(c(1, 3)) mockery::stub(calculate_treated, 'bernoulli_multi_p', bernoulli_mock) mockery::stub(calculate_treated, 'log_uniform', mockery::mock(c(3, 4))) @@ -264,8 +267,10 @@ test_that('calculate_treated correctly samples treated and updates the drug stat variables, clinical_infections, events$recovery, + events$detection, parameters, - timestep + timestep, + mock_render(timestep) ) mockery::expect_args( @@ -273,7 +278,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat 1, 2, 3, - c(.5, .5), + c(.25, .25), TRUE ) @@ -293,28 +298,35 @@ test_that('calculate_treated correctly samples treated and updates the drug stat expect_bitset_update(variables$drug_time$queue_update, 5, c(1, 4)) expect_bitset_schedule( - schedule_mock, + recovery_mock, c(1, 4), c(3, 4) ) + expect_bitset_schedule( + detection_mock, + c(1, 4), + 0 + ) }) test_that('schedule_infections correctly schedules new infections', { parameters <- get_parameters() - scheduled <- individual::Bitset$new(20) - scheduled$insert(c(1, 3, 7, 15)) + scheduled <- individual::Bitset$new(20)$insert(c(1, 3, 7, 15)) clinical_mock <- mockery::mock() asym_mock <- mockery::mock() - all_mock <- mockery::mock(cycle = TRUE) + detection_mock <- mockery::mock() events <- list( - infection = list( - schedule = all_mock, - get_scheduled = mockery::mock(scheduled) + detection = list(schedule = detection_mock), + clinical_infection = list( + schedule = clinical_mock, + get_scheduled = mockery::mock(individual::Bitset$new(20)$insert(c(1, 3))) + ), + asymptomatic_infection = list( + schedule = asym_mock, + get_scheduled = mockery::mock(individual::Bitset$new(20)$insert(c(7, 15))) ), - clinical_infection = list(schedule = clinical_mock), - asymptomatic_infection = list(schedule = asym_mock), subpatent_infection = list(clear_schedule = mockery::mock()), recovery = list(clear_schedule = mockery::mock()) ) @@ -328,12 +340,15 @@ test_that('schedule_infections correctly schedules new infections', { ) ) - infections <- individual::Bitset$new(20) - infections$insert(1:20) - clinical_infections <- individual::Bitset$new(20) - clinical_infections$insert(5:15) - treated <- individual::Bitset$new(20) - treated$insert(7:12) + mockery::stub( + schedule_infections, + 'bernoulli_multi_p', + mockery::mock(c(1, 3, 6)) + ) + + infections <- individual::Bitset$new(20)$insert(1:20) + clinical_infections <- individual::Bitset$new(20)$insert(5:15) + treated <- individual::Bitset$new(20)$insert(7:12) schedule_infections( events, @@ -350,7 +365,7 @@ test_that('schedule_infections correctly schedules new infections', { ) expect_bitset_schedule( - all_mock, + detection_mock, c(5, 6, 13, 14), c(5, 6, 13, 14) ) @@ -362,7 +377,7 @@ test_that('schedule_infections correctly schedules new infections', { ) expect_bitset_schedule( - all_mock, + detection_mock, c(2, 4, 16, 17, 18, 19, 20), c(2, 4, 16, 17, 18, 19, 20), call = 2 @@ -371,7 +386,7 @@ test_that('schedule_infections correctly schedules new infections', { test_that('prophylaxis is considered for medicated humans', { parameters <- get_parameters() - parameters <- set_drugs(parameters, list(AL_params, DHC_PQP_params)) + parameters <- set_drugs(parameters, list(AL_params, DHA_PQP_params)) events <- create_events(parameters) timestep <- 50 @@ -388,7 +403,7 @@ test_that('prophylaxis is considered for medicated humans', { ) bitten <- individual::Bitset$new(4)$insert(seq(4)) - m <- mockery::mock(individual::Bitset$new(4)$insert(seq(3))) + m <- mockery::mock(seq(3)) mockery::stub(calculate_infections, 'bernoulli_multi_p', m) calculate_infections(variables, bitten, parameters, timestep) diff --git a/tests/testthat/test-ode.R b/tests/testthat/test-ode.R index 4fdc1a94..f1f12754 100644 --- a/tests/testthat/test-ode.R +++ b/tests/testthat/test-ode.R @@ -1,20 +1,68 @@ test_that('ODE stays at equilibrium with a constant total_M', { + parameters <- get_parameters() + total_M <- 1000 + f <- parameters$blood_meal_rates + models <- parameterise_mosquito_models(parameters) + solvers <- parameterise_solvers(models, parameters) + timesteps <- 365 * 10 + + counts <- c() + + for (t in seq(timesteps)) { + counts <- rbind(counts, c(t, solver_get_states(solvers[[1]]))) + mosquito_model_update(models[[1]], total_M, f, parameters$mum) + solver_step(solvers[[1]]) + } + + expected <- c() + equilibrium <- initial_mosquito_counts( + parameters, + 1, + parameters$init_foim, + parameters$total_M + )[ODE_INDICES] + for (t in seq(timesteps)) { + expected <- rbind(expected, c(t, equilibrium)) + } + + expect_equal(counts, expected, tolerance=1e-4) +}) + +test_that('Adult ODE stays at equilibrium with a constant foim and mu', { + foim = 0.5 parameters <- get_parameters(list( - species_proportions = 1 + individual_mosquitoes = FALSE, + init_foim = foim )) total_M <- 1000 - model <- parameterise_ode(parameters)[[1]] + f <- parameters$blood_meal_rates + models <- parameterise_mosquito_models(parameters) + solvers <- parameterise_solvers(models, parameters) timesteps <- 365 * 10 counts <- c() for (t in seq(timesteps)) { - counts <- rbind(counts, c(t, mosquito_model_get_states(model))) - mosquito_model_step(model, total_M) + states <- solver_get_states(solvers[[1]]) + counts <- rbind(counts, c(t, states)) + adult_mosquito_model_update( + models[[1]], + parameters$mum, + foim, + states[ADULT_ODE_INDICES['Sm']], + f + ) + solver_step(solvers[[1]]) } expected <- c() - equilibrium <- initial_mosquito_counts(parameters)[seq(3)] + equilibrium <- initial_mosquito_counts( + parameters, + 1, + parameters$init_foim, + total_M + ) + for (t in seq(timesteps)) { expected <- rbind(expected, c(t, equilibrium)) } @@ -25,22 +73,29 @@ test_that('ODE stays at equilibrium with a constant total_M', { test_that('ODE stays at equilibrium with low total_M', { total_M <- 10 parameters <- get_parameters(list( - total_M = total_M, - species = 'all', - species_proportions = 1 + total_M = total_M )) - model <- parameterise_ode(parameters)[[1]] + f <- parameters$blood_meal_rates + models <- parameterise_mosquito_models(parameters) + solvers <- parameterise_solvers(models, parameters) timesteps <- 365 * 10 counts <- c() for (t in seq(timesteps)) { - counts <- rbind(counts, c(t, mosquito_model_get_states(model))) - mosquito_model_step(model, total_M) + counts <- rbind(counts, c(t, solver_get_states(solvers[[1]]))) + mosquito_model_update(models[[1]], total_M, f, parameters$mum) + solver_step(solvers[[1]]) } expected <- c() - equilibrium <- initial_mosquito_counts(parameters)[seq(3)] + equilibrium <- initial_mosquito_counts( + parameters, + 1, + parameters$init_foim, + parameters$total_M + )[ODE_INDICES] + for (t in seq(timesteps)) { expected <- rbind(expected, c(t, equilibrium)) } @@ -53,11 +108,11 @@ test_that('Changing total_M stabilises', { total_M_0 <- 500 total_M_1 <- 400 parameters <- get_parameters(list( - total_M = total_M_0, - species = 'all', - species_proportions = 1 + total_M = total_M_0 )) - model <- parameterise_ode(parameters)[[1]] + f <- parameters$blood_meal_rates + models <- parameterise_mosquito_models(parameters) + solvers <- parameterise_solvers(models, parameters) timesteps <- 365 * 10 change <- 50 burn_in <- 365 * 5 @@ -65,17 +120,23 @@ test_that('Changing total_M stabilises', { counts <- c() for (t in seq(timesteps)) { - counts <- rbind(counts, c(t, mosquito_model_get_states(model))) + counts <- rbind(counts, c(t, solver_get_states(solvers[[1]]))) if (t < change) { total_M <- total_M_0 } else { total_M <- total_M_1 } - mosquito_model_step(model, total_M) + mosquito_model_update(models[[1]], total_M, f, parameters$mum) + solver_step(solvers[[1]]) } - initial_eq <- initial_mosquito_counts(parameters)[seq(3)] - final_eq <- counts[burn_in, seq(3) + 1] + initial_eq <- initial_mosquito_counts( + parameters, + 1, + parameters$init_foim, + parameters$total_M + )[ODE_INDICES] + final_eq <- counts[burn_in, ODE_INDICES + 1] expect_equal(counts[1,], c(1, initial_eq), tolerance=1e-4) expect_equal(counts[timesteps,], c(timesteps, final_eq), tolerance=1e-4) @@ -86,6 +147,7 @@ test_that('Changing total_M stabilises', { test_that('You can parameterise Total_M = 0 🤯', { total_M <- 1000 parameters <- get_parameters() + f <- parameters$blood_meal_rates parameters <- set_species( parameters, species = list(arab_params, fun_params), @@ -95,14 +157,21 @@ test_that('You can parameterise Total_M = 0 🤯', { parameters, total_M ) - model <- parameterise_ode(parameters)[[2]] + models <- parameterise_mosquito_models(parameters) + solvers <- parameterise_solvers(models, parameters) timesteps <- 365 * 10 counts <- c() for (t in seq(timesteps)) { - counts <- rbind(counts, c(t, mosquito_model_get_states(model))) - mosquito_model_step(model, 0) + counts <- rbind(counts, c(t, solver_get_states(solvers[[2]]))) + mosquito_model_update( + models[[2]], + total_M * parameters$species_proportions[[2]], + f, + parameters$mum[[2]] + ) + solver_step(solvers[[2]]) } expected <- cbind( diff --git a/tests/testthat/test-process-integration.R b/tests/testthat/test-process-integration.R index 1413127f..9e8fc901 100644 --- a/tests/testthat/test-process-integration.R +++ b/tests/testthat/test-process-integration.R @@ -3,14 +3,16 @@ test_that('create_processes makes valid process functions', { parameters <- get_parameters() events <- create_events(parameters) variables <- create_variables(parameters) - odes <- parameterise_ode(parameters) + vector_models <- parameterise_mosquito_models(parameters) + solvers <- parameterise_solvers(vector_models, parameters) renderer <- individual::Render$new(1) processes <- create_processes( renderer, variables, events, parameters, - odes + vector_models, + solvers ) for (process in processes) { expect(is.function(process) || inherits(process, 'externalptr'), 'Process is not a function') diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R index 6a4a1932..6723c862 100644 --- a/tests/testthat/test-render.R +++ b/tests/testthat/test-render.R @@ -2,26 +2,30 @@ test_that('that default rendering works', { timestep <- 0 parameters <- get_parameters() state <- individual::CategoricalVariable$new( - c('U', 'A', 'D', 'S'), + c('U', 'A', 'D', 'S', 'Tr'), c('U', 'A', 'D', 'S') ) birth <- individual::DoubleVariable$new( -c(2, 5, 10, 11) * 365 ) + immunity <- individual::DoubleVariable$new(rep(1, 4)) is_severe <- individual::CategoricalVariable$new( c('yes', 'no'), rep('no', 4) ) + renderer <- mock_render(1) process <- create_prevelance_renderer( state, birth, is_severe, + immunity, parameters, renderer ) + mockery::stub(process, 'bernoulli_multi_p', mockery::mock(1)) process(timestep) mockery::expect_args( @@ -37,12 +41,13 @@ test_that('that default rendering works when no one is in the age range', { timestep <- 0 parameters <- get_parameters() state <- individual::CategoricalVariable$new( - c('U', 'A', 'D', 'S'), + c('U', 'A', 'D', 'S', 'Tr'), rep('S', 4) ) birth <- individual::DoubleVariable$new( -c(1, 11, 21, 11) * 365 ) + immunity <- individual::DoubleVariable$new(rep(1, 4)) is_severe <- individual::CategoricalVariable$new( c('yes', 'no'), rep('no', 4) @@ -53,10 +58,12 @@ test_that('that default rendering works when no one is in the age range', { state, birth, is_severe, + immunity, parameters, renderer ) process(timestep) + mockery::stub(process, 'bernoulli_multi_p', mockery::mock(1)) mockery::expect_args( renderer$render, @@ -77,12 +84,13 @@ test_that('that severe rendering works', { prevalence_rendering_max_ages = NULL )) state <- individual::CategoricalVariable$new( - c('U', 'A', 'D', 'S'), + c('U', 'A', 'D', 'S', 'Tr'), c('U', 'D', 'D', 'S') ) birth <- individual::DoubleVariable$new( -c(2, 5, 10, 11) * 365 ) + immunity <- individual::DoubleVariable$new(rep(1, 4)) is_severe <- individual::CategoricalVariable$new( c('yes', 'no'), c('no', 'yes', 'no', 'no') @@ -92,6 +100,7 @@ test_that('that severe rendering works', { state, birth, is_severe, + immunity, parameters, renderer ) @@ -114,6 +123,43 @@ test_that('that severe rendering works', { ) }) +test_that('that clinical incidence rendering works', { + timestep <- 0 + year <- 365 + parameters <- get_parameters(list( + clinical_incidence_rendering_min_ages = c(0, 2) * year, + clinical_incidence_rendering_max_ages = c(5, 10) * year, + prevalence_rendering_min_ages = NULL, + prevalence_rendering_max_ages = NULL + )) + + birth <- individual::DoubleVariable$new( + -c(2, 5, 10, 11) * year + ) + + renderer <- mock_render(1) + process <- create_clinical_incidence_renderer(birth, parameters, renderer) + + process(timestep, individual::Bitset$new(4)$insert(c(1, 2, 4))) + + mockery::expect_args( + renderer$render, + 1, + 'clin_inc_0_1825', + 1, + timestep + ) + + mockery::expect_args( + renderer$render, + 2, + 'clin_inc_730_3650', + 2/3, + timestep + ) +}) + + test_that('that incidence rendering works', { timestep <- 0 year <- 365 diff --git a/tests/testthat/test-rtss.R b/tests/testthat/test-rtss.R index 1292393f..78b9e187 100644 --- a/tests/testthat/test-rtss.R +++ b/tests/testthat/test-rtss.R @@ -84,7 +84,7 @@ test_that('Infection considers vaccine efficacy', { rep(.2, 4) ) - bernoulli_mock <- mockery::mock(individual::Bitset$new(4)$insert(c(1, 2))) + bernoulli_mock <- mockery::mock(c(1, 2)) mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock) calculate_infections( variables, diff --git a/tests/testthat/test-seasonality.R b/tests/testthat/test-seasonality.R index 79fc2ad7..1b90639a 100644 --- a/tests/testthat/test-seasonality.R +++ b/tests/testthat/test-seasonality.R @@ -7,14 +7,21 @@ test_that('Seasonality correctly affects P', { species_proportions = 1 )) total_M <- 1000 - model <- parameterise_ode(parameters)[[1]] + models <- parameterise_mosquito_models(parameters) + solvers <- parameterise_solvers(models, parameters) timesteps <- 365 * 40 counts <- c() for (t in seq(timesteps)) { - counts <- rbind(counts, c(t, mosquito_model_get_states(model))) - mosquito_model_step(model, total_M) + counts <- rbind(counts, c(t, solver_get_states(solvers[[1]]))) + mosquito_model_update( + models[[1]], + total_M, + parameters$blood_meal_rates, + parameters$mum + ) + solver_step(solvers[[1]]) } burn_in <- 20 diff --git a/tests/testthat/test-treatment.R b/tests/testthat/test-treatment.R new file mode 100644 index 00000000..9c0a95be --- /dev/null +++ b/tests/testthat/test-treatment.R @@ -0,0 +1,61 @@ + +test_that('You can set time varying clinical treatment', { + parameters <- get_parameters() + parameters <- set_drugs(parameters, list(DHA_PQP_params, AL_params)) + parameters <- set_clinical_treatment( + parameters, + drug = 1, + timesteps = c(50, 100), + coverages = c(.2, .5) + ) + parameters <- set_clinical_treatment( + parameters, + drug = 2, + timesteps = c(25, 75), + coverages = c(.3, .4) + ) + + expect_equal( + get_treatment_coverages(parameters, 10), + c(0, 0) + ) + + expect_equal( + get_treatment_coverages(parameters, 25), + c(0, .3) + ) + + expect_equal( + get_treatment_coverages(parameters, 30), + c(0, .3) + ) + + expect_equal( + get_treatment_coverages(parameters, 50), + c(.2, .3) + ) + + expect_equal( + get_treatment_coverages(parameters, 80), + c(.2, .4) + ) +}) + +test_that('You cannot set invalid coverages', { + parameters <- get_parameters() + parameters <- set_drugs(parameters, list(DHA_PQP_params, AL_params)) + parameters <- set_clinical_treatment( + parameters, + drug = 1, + timesteps = c(50, 100), + coverages = c(.2, .5) + ) + expect_error( + set_clinical_treatment( + parameters, + drug = 2, + timesteps = c(25, 75), + coverages = c(.3, .6) + ) + ) +}) diff --git a/vignettes/Treatment.Rmd b/vignettes/Treatment.Rmd index 294f1c44..a7ca9ff7 100644 --- a/vignettes/Treatment.Rmd +++ b/vignettes/Treatment.Rmd @@ -35,7 +35,7 @@ simparams <- get_parameters(c( human_population = human_population ) )) -simparams <- set_drugs(simparams, list(AL_params, DHC_PQP_params)) +simparams <- set_drugs(simparams, list(AL_params, DHA_PQP_params)) # NOTE: this equilibrium is calculated using a different prophylaxis model # it is just used as a starting point @@ -55,11 +55,9 @@ for (ft in fts) { eq <- get_equilibrium(starting_EIR, ft) simparams <- parameterise_human_equilibrium(simparams, eq) simparams <- parameterise_mosquito_equilibrium(simparams, starting_EIR) - simparams <- set_clinical_treatment(simparams, ft, c(1, 2), c(.5, .5)) - output <- run_simulation( - sim_length, - simparams - ) + simparams <- set_clinical_treatment(simparams, 1, 1, ft / 2) + simparams <- set_clinical_treatment(simparams, 2, 1, ft / 2) + output <- run_simulation(sim_length, simparams) outputs[[length(outputs) + 1]] <- output } ```