Merge pull request #235 from mrc-ide/feat/carrying_capacity

Feat/carrying capacity
mrc-ide · May 30, 2023 · c9aa9d7 · c9aa9d7
2 parents cdaeae8 + e016850
commit c9aa9d7
Show file tree

Hide file tree

Showing 32 changed files with 947 additions and 490 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -8,6 +8,7 @@ export(create_pev_profile)
 export(fun_params)
 export(gamb_params)
 export(get_correlation_parameters)
+export(get_init_carrying_capacity)
 export(get_parameters)
 export(parameterise_mosquito_equilibrium)
 export(parameterise_total_M)
@@ -18,6 +19,7 @@ export(run_metapop_simulation)
 export(run_simulation)
 export(run_simulation_with_repetitions)
 export(set_bednets)
+export(set_carrying_capacity)
 export(set_clinical_treatment)
 export(set_demography)
 export(set_drugs)
@@ -31,6 +33,7 @@ export(set_smc)
 export(set_species)
 export(set_spraying)
 export(set_tbv)
+export(steph_params)
 importFrom(MASS,mvrnorm)
 importFrom(Rcpp,sourceCpp)
 importFrom(stats,qnorm)

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -13,8 +13,8 @@ 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)
 }
 
-create_aquatic_mosquito_model <- function(beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, g0, g, h, R_bar, mum, f, rainfall_floor) {
-    .Call(`_malariasimulation_create_aquatic_mosquito_model`, beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, g0, g, h, R_bar, mum, f, rainfall_floor)
+create_aquatic_mosquito_model <- function(beta, de, mue, k_timeseries, gamma, dl, mul, dp, mup, total_M, model_seasonality, g0, g, h, R_bar, mum, f, rainfall_floor) {
+    .Call(`_malariasimulation_create_aquatic_mosquito_model`, beta, de, mue, k_timeseries, gamma, dl, mul, dp, mup, total_M, model_seasonality, g0, g, h, R_bar, mum, f, rainfall_floor)
 }
 
 aquatic_mosquito_model_update <- function(model, total_M, f, mum) {
@@ -25,18 +25,6 @@ create_aquatic_solver <- function(model, init, r_tol, a_tol, max_steps) {
     .Call(`_malariasimulation_create_aquatic_solver`, model, init, r_tol, a_tol, max_steps)
 }
 
-create_history <- function(size, default_value) {
-    .Call(`_malariasimulation_create_history`, size, default_value)
-}
-
-history_at <- function(history, timestep) {
-    .Call(`_malariasimulation_history_at`, history, timestep)
-}
-
-history_push <- function(history, value, timestep) {
-    invisible(.Call(`_malariasimulation_history_push`, history, value, timestep))
-}
-
 carrying_capacity <- function(timestep, model_seasonality, g0, g, h, K0, R_bar, rainfall_floor) {
     .Call(`_malariasimulation_carrying_capacity`, timestep, model_seasonality, g0, g, h, K0, R_bar, rainfall_floor)
 }
@@ -57,6 +45,18 @@ solver_step <- function(solver) {
     invisible(.Call(`_malariasimulation_solver_step`, solver))
 }
 
+create_timeseries <- function(size, default_value) {
+    .Call(`_malariasimulation_create_timeseries`, size, default_value)
+}
+
+timeseries_at <- function(timeseries, timestep, linear) {
+    .Call(`_malariasimulation_timeseries_at`, timeseries, timestep, linear)
+}
+
+timeseries_push <- function(timeseries, value, timestep) {
+    invisible(.Call(`_malariasimulation_timeseries_push`, timeseries, value, timestep))
+}
+
 random_seed <- function(seed) {
     invisible(.Call(`_malariasimulation_random_seed`, seed))
 }

diff --git a/R/biting_process.R b/R/biting_process.R
@@ -17,21 +17,21 @@
 #' lagged_infectivity list (default: 1)
 #' @noRd
 create_biting_process <- function(
-  renderer,
-  solvers,
-  models,
-  variables,
-  events,
-  parameters,
-  lagged_infectivity,
-  lagged_eir,
-  mixing = 1,
-  mixing_index = 1
-  ) {
+    renderer,
+    solvers,
+    models,
+    variables,
+    events,
+    parameters,
+    lagged_infectivity,
+    lagged_eir,
+    mixing = 1,
+    mixing_index = 1
+) {
   function(timestep) {
     # Calculate combined EIR
     age <- get_age(variables$birth$get_values(), timestep)
-
+    
     bitten_humans <- simulate_bites(
       renderer,
       solvers,
@@ -46,7 +46,7 @@ create_biting_process <- function(
       mixing,
       mixing_index
     )
-
+    
     simulate_infection(
       variables,
       events,
@@ -61,21 +61,21 @@ create_biting_process <- function(
 
 #' @importFrom stats rpois
 simulate_bites <- function(
-  renderer,
-  solvers,
-  models,
-  variables,
-  events,
-  age,
-  parameters,
-  timestep,
-  lagged_infectivity,
-  lagged_eir,
-  mixing = 1,
-  mixing_index = 1
-  ) {
+    renderer,
+    solvers,
+    models,
+    variables,
+    events,
+    age,
+    parameters,
+    timestep,
+    lagged_infectivity,
+    lagged_eir,
+    mixing = 1,
+    mixing_index = 1
+) {
   bitten_humans <- individual::Bitset$new(parameters$human_population)
-
+  
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
     human_infectivity <- account_for_tbv(
@@ -86,21 +86,21 @@ simulate_bites <- function(
     )
   }
   renderer$render('infectivity', mean(human_infectivity), timestep)
-
+  
   # Calculate pi (the relative biting rate for each human)
   psi <- unique_biting_rate(age, parameters)
   zeta <- variables$zeta$get_values()
   .pi <- human_pi(zeta, psi)
-
+  
   # Get some indices for later
   if (parameters$individual_mosquitoes) {
     infectious_index <- variables$mosquito_state$get_index_of('Im')
     susceptible_index <- variables$mosquito_state$get_index_of('Sm')
     adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not(TRUE)
   }
-
+  
   EIR <- 0
-
+  
   for (s_i in seq_along(parameters$species)) {
     species_name <- parameters$species[[s_i]]
     solver_states <- solver_get_states(solvers[[s_i]])
@@ -111,7 +111,7 @@ simulate_bites <- function(
     f <- blood_meal_rate(s_i, Z, parameters)
     a <- .human_blood_meal_rate(f, s_i, W, parameters)
     lambda <- effective_biting_rates(a, .pi, p_bitten)
-
+    
     if (parameters$individual_mosquitoes) {
       species_index <- variables$species$get_index_of(
         parameters$species[[s_i]]
@@ -127,18 +127,18 @@ simulate_bites <- function(
     } else {
       n_infectious <- calculate_infectious_compartmental(solver_states)
     }
-
+    
     # store the current population's EIR for later
     lagged_eir[[mixing_index]][[s_i]]$save(n_infectious * a, timestep)
-
+    
     # calculated the EIR for this timestep after mixing
     species_eir <- sum(
       vnapply(
         lagged_eir,
         function(l) l[[s_i]]$get(timestep - parameters$de)
       ) * mixing
     )
-
+    
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     expected_bites <- species_eir * mean(psi)
@@ -150,7 +150,7 @@ simulate_bites <- function(
         )
       }
     }
-
+    
     infectivity <- vnapply(
       lagged_infectivity,
       function(l) l$get(timestep - parameters$delay_gam)
@@ -163,7 +163,7 @@ simulate_bites <- function(
     renderer$render(paste0('FOIM_', species_name), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
     renderer$render(paste0('mu_', species_name), mu, timestep)
-
+    
     if (parameters$individual_mosquitoes) {
       # update the ODE with stats for ovoposition calculations
       aquatic_mosquito_model_update(
@@ -172,10 +172,10 @@ simulate_bites <- function(
         f,
         mu
       )
-
+      
       # update the individual mosquitoes
       susceptible_species_index <- susceptible_index$copy()$and(species_index)
-
+      
       biting_effects_individual(
         variables,
         foim,
@@ -197,7 +197,7 @@ simulate_bites <- function(
       )
     }
   }
-
+  
   renderer$render('n_bitten', bitten_humans$size(), timestep)
   bitten_humans
 }
@@ -239,13 +239,13 @@ calculate_infectious <- function(species, solvers, variables, parameters) {
 }
 
 calculate_infectious_individual <- function(
-  species,
-  variables,
-  infectious_index,
-  adult_index,
-  species_index,
-  parameters
-  ) {
+    species,
+    variables,
+    infectious_index,
+    adult_index,
+    species_index,
+    parameters
+) {
 
   infectious_index$copy()$and(species_index)$size()
 }
@@ -255,7 +255,7 @@ calculate_infectious_compartmental <- function(solver_states) {
 }
 
 intervention_coefficient <- function(p_bitten) {
-   p_bitten$prob_bitten / sum(p_bitten$prob_bitten_survives)
+  p_bitten$prob_bitten / sum(p_bitten$prob_bitten_survives)
 }
 
 human_pi <- function(zeta, psi) {

diff --git a/R/compartmental.R b/R/compartmental.R
@@ -1,17 +1,31 @@
 ODE_INDICES <- c(E = 1, L = 2, P = 3)
 ADULT_ODE_INDICES <- c(Sm = 4, Pm = 5, Im = 6)
 
-parameterise_mosquito_models <- function(parameters) {
+parameterise_mosquito_models <- function(parameters, timesteps) {
+
   lapply(
     seq_along(parameters$species),
     function(i) {
       p <- parameters$species_proportions[[i]]
       m <- p * parameters$total_M
+      # Baseline carrying capacity
+      k0 <- calculate_carrying_capacity(parameters, m, i)
+      # Create the carrying capacity object
+      k_timeseries <- create_timeseries(size = length(parameters$carrying_capacity_timesteps), k0)
+      if(parameters$carrying_capacity){
+        for(j in 1:length(parameters$carrying_capacity_timesteps)){
+          timeseries_push(
+            k_timeseries,
+            parameters$carrying_capacity_values[j,i],
+            parameters$carrying_capacity_timesteps[j]
+          )
+        }
+      }
       growth_model <- create_aquatic_mosquito_model(
         parameters$beta,
         parameters$del,
         parameters$me,
-        calculate_carrying_capacity(parameters, m, i),
+        k_timeseries,
         parameters$gamma,
         parameters$dl,
         parameters$ml,
@@ -27,7 +41,7 @@ parameterise_mosquito_models <- function(parameters) {
         parameters$blood_meal_rates[[i]],
         parameters$rainfall_floor
       )
-
+      
       if (!parameters$individual_mosquitoes) {
         susceptible <- initial_mosquito_counts(
           parameters,
@@ -84,7 +98,7 @@ create_compartmental_rendering_process <- function(renderer, solvers, parameters
   } else {
     indices <- c(ODE_INDICES, ADULT_ODE_INDICES)
   }
-
+  
   function(timestep) {
     counts <- rep(0, length(indices))
     for (s_i in seq_along(solvers)) {

diff --git a/R/lag.R b/R/lag.R
@@ -5,15 +5,15 @@ LaggedValue <- R6::R6Class(
   ),
   public = list(
     initialize = function(max_lag, default) {
-      private$history <- create_history(max_lag, default)
+      private$history <- create_timeseries(max_lag, default)
     },
 
     save = function(value, timestep) {
-      history_push(private$history, value, timestep)
+      timeseries_push(private$history, value, timestep)
     },
 
     get = function(timestep) {
-      history_at(private$history, timestep)
+      timeseries_at(private$history, timestep, TRUE)
     }
   )
 )