Bug/spraying

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: malariasimulation
 Title: An individual based model for malaria
-Version: 1.6.0
+Version: 1.6.1
 Authors@R: c(
   person(
     given = "Giovanni",
@@ -70,12 +70,12 @@ Remotes:
   mrc-ide/malariaEquilibrium,
   mrc-ide/individual
 Imports:
-    individual (>= 0.1.7),
+    individual (>= 0.1.16),
     malariaEquilibrium (>= 1.0.1),
     Rcpp,
     statmod,
     MASS,
-    dqrng,
+    dqrng (>= 0.4),
     sitmo,
     BH,
     R6,
@@ -91,7 +91,7 @@ Suggests:
   ggplot2,
   covr,
   mgcv 
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 Roxygen: list(markdown = TRUE)
 LinkingTo:
   Rcpp,

NAMESPACE

@@ -18,6 +18,7 @@ export(rtss_profile)
 export(run_metapop_simulation)
 export(run_simulation)
 export(run_simulation_with_repetitions)
+export(set_antimalarial_resistance)
 export(set_bednets)
 export(set_carrying_capacity)
 export(set_clinical_treatment)

NEWS.md

@@ -1,3 +1,10 @@
+# malariasimulation 1.6.1 (wip)
+
+ * Fix bug in competing hazards between mass and EPI vaccines. Where individuals
+   can be enrolled onto both strategies if applied on the same timestep.
+ * Fix bug with min_wait. Min wait was working off of the final primary dose. It
+   now works of of the first dose.
+
 # malariasimulation 1.6.0
 
   * Fix MDA bug where undetectable asymptomatics are treated

R/RcppExports.R

@@ -9,6 +9,14 @@ adult_mosquito_model_update <- function(model, mu, foim, susceptible, f) {
     invisible(.Call(`_malariasimulation_adult_mosquito_model_update`, model, mu, foim, susceptible, f))
 }
 
+adult_mosquito_model_save_state <- function(model) {
+    .Call(`_malariasimulation_adult_mosquito_model_save_state`, model)
+}
+
+adult_mosquito_model_restore_state <- function(model, state) {
+    invisible(.Call(`_malariasimulation_adult_mosquito_model_restore_state`, model, state))
+}
+
 create_adult_solver <- function(model, init, r_tol, a_tol, max_steps) {
     .Call(`_malariasimulation_create_adult_solver`, model, init, r_tol, a_tol, max_steps)
 }
@@ -37,10 +45,18 @@ rainfall <- function(t, g0, g, h, floor) {
     .Call(`_malariasimulation_rainfall`, t, g0, g, h, floor)
 }
 
+exponential_process_cpp <- function(variable, rate) {
+    .Call(`_malariasimulation_exponential_process_cpp`, variable, rate)
+}
+
 solver_get_states <- function(solver) {
     .Call(`_malariasimulation_solver_get_states`, solver)
 }
 
+solver_set_states <- function(solver, t, state) {
+    invisible(.Call(`_malariasimulation_solver_set_states`, solver, t, state))
+}
+
 solver_step <- function(solver) {
     invisible(.Call(`_malariasimulation_solver_step`, solver))
 }
@@ -57,10 +73,26 @@ timeseries_push <- function(timeseries, value, timestep) {
     invisible(.Call(`_malariasimulation_timeseries_push`, timeseries, value, timestep))
 }
 
+timeseries_save_state <- function(timeseries) {
+    .Call(`_malariasimulation_timeseries_save_state`, timeseries)
+}
+
+timeseries_restore_state <- function(timeseries, state) {
+    invisible(.Call(`_malariasimulation_timeseries_restore_state`, timeseries, state))
+}
+
 random_seed <- function(seed) {
     invisible(.Call(`_malariasimulation_random_seed`, seed))
 }
 
+random_save_state <- function() {
+    .Call(`_malariasimulation_random_save_state`)
+}
+
+random_restore_state <- function(state) {
+    invisible(.Call(`_malariasimulation_random_restore_state`, state))
+}
+
 bernoulli_multi_p_cpp <- function(p) {
     .Call(`_malariasimulation_bernoulli_multi_p_cpp`, p)
 }

R/antimalarial_resistance.R

@@ -0,0 +1,146 @@
+#' @title Parameterise antimalarial resistance
+#' @description
+#' Parameterise antimalarial resistance
+#' 
+#' @param parameters the model parameters
+#' @param drug the index of the drug which resistance is being set, as set by the set_drugs() function, in the parameter list
+#' @param timesteps vector of time steps for each update to resistance proportion and resistance outcome probability
+#' @param artemisinin_resistance_proportion vector of updates to the proportions of infections that are artemisinin resistant at time t
+#' @param partner_drug_resistance_proportion vector of updates to the proportions of infections that are partner-drug resistant at time t
+#' @param slow_parasite_clearance_probability vector of updates to the proportion of artemisinin-resistant infections that result in early treatment failure
+#' @param early_treatment_failure_probability vector of updates to the proportion of artemisinin-resistant infections that result in slow parasite clearance
+#' @param late_clinical_failure_probability vector of updates to the proportion of partner-drug-resistant infections that result in late clinical failure
+#' @param late_parasitological_failure_probability vector of updates to the proportion of partner-drug-resistant infections that result in late parasitological failure
+#' @param reinfection_during_prophylaxis_probability vector of updates to the proportion of partner-drug-resistant infections that result in reinfection during prophylaxis
+#' @param slow_parasite_clearance_time single value representing the mean time individual's experiencing slow parasite clearance reside in the treated state
+#' @export
+set_antimalarial_resistance <- function(parameters,
+                                        drug,
+                                        timesteps,
+                                        artemisinin_resistance_proportion,
+                                        partner_drug_resistance_proportion,
+                                        slow_parasite_clearance_probability,
+                                        early_treatment_failure_probability,
+                                        late_clinical_failure_probability,
+                                        late_parasitological_failure_probability,
+                                        reinfection_during_prophylaxis_probability,
+                                        slow_parasite_clearance_time) {
+
+  if(any(partner_drug_resistance_proportion > 0,
+         late_clinical_failure_probability > 0,
+         late_parasitological_failure_probability > 0,
+         reinfection_during_prophylaxis_probability > 0)) {
+    stop("Parameters set for unimplemented feature - late clinical failure, late parasitological failure, or reinfection during prophylaxis")
+  }
+
+  if(any(c(length(artemisinin_resistance_proportion),
+           length(partner_drug_resistance_proportion), 
+           length(slow_parasite_clearance_probability), 
+           length(early_treatment_failure_probability), 
+           length(late_clinical_failure_probability), 
+           length(late_parasitological_failure_probability), 
+           length(reinfection_during_prophylaxis_probability)) != length(timesteps))) {
+    stop("Length of one or more resistance parameter vectors does not match time steps specified for update")
+  }
+
+  if(any(artemisinin_resistance_proportion < 0 | artemisinin_resistance_proportion > 1 |
+         partner_drug_resistance_proportion < 0 | partner_drug_resistance_proportion > 1)) {
+    stop("Artemisinin and partner-drug resistance proportions must fall between 0 and 1")
+  }
+
+  if(any(slow_parasite_clearance_probability < 0 | slow_parasite_clearance_probability > 1 |
+         early_treatment_failure_probability < 0 | early_treatment_failure_probability > 1 |
+         late_clinical_failure_probability < 0 | late_clinical_failure_probability > 1 |
+         late_parasitological_failure_probability < 0 | late_parasitological_failure_probability > 1 |
+         reinfection_during_prophylaxis_probability < 0 | reinfection_during_prophylaxis_probability > 1)) {
+    stop("Resistance outcome probabilities must fall between 0 and 1")
+  }
+
+  if(length(slow_parasite_clearance_time) != 1) {
+    stop("Error: length of slow_parasite_clearance_time not equal to 1")
+  }
+
+  if(slow_parasite_clearance_time <= 0) {
+    stop("Error: slow_parasite_clearance_time is non-positive")
+  }
+
+  parameters$antimalarial_resistance <- TRUE
+
+  n_drugs <- length(parameters$drug_efficacy)
+
+  if (drug < 1 | drug > n_drugs) {
+    stop('Drug index is invalid, please set drugs using set_drugs')
+  }
+
+  drug_index <- which(parameters$antimalarial_resistance_drug == drug)
+
+  if (length(drug_index) == 0) {
+    drug_index <- length(parameters$antimalarial_resistance_drug) + 1
+  }
+
+  parameters$antimalarial_resistance_drug[[drug_index]] <- drug
+  parameters$antimalarial_resistance_timesteps[[drug_index]] <- timesteps
+  parameters$artemisinin_resistance_proportion[[drug_index]] <- artemisinin_resistance_proportion
+  parameters$partner_drug_resistance_proportion[[drug_index]] <- partner_drug_resistance_proportion
+  parameters$slow_parasite_clearance_probability[[drug_index]] <- slow_parasite_clearance_probability
+  parameters$early_treatment_failure_probability[[drug_index]] <- early_treatment_failure_probability
+  parameters$late_clinical_failure_probability[[drug_index]] <- late_clinical_failure_probability
+  parameters$late_parasitological_failure_probability[[drug_index]] <- late_parasitological_failure_probability
+  parameters$reinfection_during_prophylaxis_probability[[drug_index]] <- reinfection_during_prophylaxis_probability
+  parameters$dt_slow_parasite_clearance[[drug_index]] <- slow_parasite_clearance_time
+
+  return(parameters)
+
+}
+
+#' @title Retrieve resistance parameters
+#' @description
+#' Retrieve the resistance parameters associated with the drug each individual receiving clinical 
+#' treatment has been administered in the current time step.
+#' 
+#' @param parameters the model parameters
+#' @param drugs vector of integers representing the drugs administered to each individual receiving treatment
+#' @param timestep the current time step
+get_antimalarial_resistance_parameters <- function(parameters, drugs, timestep) {
+
+  if(!parameters$antimalarial_resistance) {
+    stop("Error: Antimalarial resistance has not been parameterised; antimalarial_resistance = FALSE")
+  }
+
+  blank_vector <- numeric(length = length(drugs))
+  artemisinin_resistance_proportion <- blank_vector
+  partner_drug_resistance_proportion <- blank_vector
+  slow_parasite_clearance_probability <- blank_vector
+  early_treatment_failure_probability <- blank_vector
+  late_clinical_failure_probability <- blank_vector
+  late_parasitological_failure_probability <- blank_vector
+  reinfection_during_prophylaxis_probability <- blank_vector
+  dt_slow_parasite_clearance <- rep(parameters$dt, length = length(drugs))
+
+  for(i in seq_along(parameters$antimalarial_resistance_drug)) {
+    drug <- parameters$antimalarial_resistance_drug[[i]]
+    treated_with_drug <- which(drugs == drug)
+    resistance_timestep <- match_timestep(ts = parameters$antimalarial_resistance_timesteps[[i]], t = timestep)
+    artemisinin_resistance_proportion[treated_with_drug] <- parameters$artemisinin_resistance_proportion[[i]][resistance_timestep]
+    partner_drug_resistance_proportion[treated_with_drug] <- parameters$partner_drug_resistance_proportion[[i]][resistance_timestep]
+    slow_parasite_clearance_probability[treated_with_drug] <- parameters$slow_parasite_clearance_probability[[i]][resistance_timestep]
+    early_treatment_failure_probability[treated_with_drug] <- parameters$early_treatment_failure_probability[[i]][resistance_timestep]
+    late_clinical_failure_probability[treated_with_drug] <- parameters$late_clinical_failure_probability[[i]][resistance_timestep]
+    late_parasitological_failure_probability[treated_with_drug] <- parameters$late_parasitological_failure_probability[[i]][resistance_timestep]
+    reinfection_during_prophylaxis_probability[treated_with_drug] <- parameters$reinfection_during_prophylaxis_probability[[i]][resistance_timestep]
+    dt_slow_parasite_clearance[treated_with_drug] <- parameters$dt_slow_parasite_clearance[[i]]
+  }
+
+  resistance_parameters <- list()
+  resistance_parameters$artemisinin_resistance_proportion <- artemisinin_resistance_proportion
+  resistance_parameters$partner_drug_resistance_proportion <- partner_drug_resistance_proportion
+  resistance_parameters$slow_parasite_clearance_probability <- slow_parasite_clearance_probability
+  resistance_parameters$early_treatment_failure_probability <- early_treatment_failure_probability
+  resistance_parameters$late_clinical_failure_probability <- late_clinical_failure_probability
+  resistance_parameters$late_parasitological_failure_probability <- late_parasitological_failure_probability
+  resistance_parameters$reinfection_during_prophylaxis_probability <- reinfection_during_prophylaxis_probability
+  resistance_parameters$dt_slow_parasite_clearance <- dt_slow_parasite_clearance
+
+  return(resistance_parameters)
+
+}

R/biting_process.R

@@ -103,7 +103,7 @@ simulate_bites <- function(
 
   for (s_i in seq_along(parameters$species)) {
     species_name <- parameters$species[[s_i]]
-    solver_states <- solver_get_states(solvers[[s_i]])
+    solver_states <- solvers[[s_i]]$get_states()
     p_bitten <- prob_bitten(timestep, variables, s_i, parameters)
     Q0 <- parameters$Q0[[s_i]]
     W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, Q0)
@@ -167,7 +167,7 @@ simulate_bites <- function(
     if (parameters$individual_mosquitoes) {
       # update the ODE with stats for ovoposition calculations
       aquatic_mosquito_model_update(
-        models[[s_i]],
+        models[[s_i]]$.model,
         species_index$size(),
         f,
         mu
@@ -189,7 +189,7 @@ simulate_bites <- function(
       )
     } else {
       adult_mosquito_model_update(
-        models[[s_i]],
+        models[[s_i]]$.model,
         mu,
         foim,
         solver_states[[ADULT_ODE_INDICES['Sm']]],
@@ -210,9 +210,7 @@ simulate_bites <- function(
 calculate_eir <- function(species, solvers, variables, parameters, timestep) {
   a <- human_blood_meal_rate(species, variables, parameters, timestep)
   infectious <- calculate_infectious(species, solvers, variables, parameters)
-  age <- get_age(variables$birth$get_values(), timestep)
-  psi <- unique_biting_rate(age, parameters)
-  infectious * a * mean(psi)
+  infectious * a
 }
 
 effective_biting_rates <- function(a, .pi, p_bitten) {
@@ -235,7 +233,7 @@ calculate_infectious <- function(species, solvers, variables, parameters) {
       )
     )
   }
-  calculate_infectious_compartmental(solver_get_states(solvers[[species]]))
+  calculate_infectious_compartmental(solvers[[species]]$get_states())
 }
 
 calculate_infectious_individual <- function(