Merge pull request #162 from mrc-ide/dev

Dev

.github/workflows/R-CMD-check.yaml

@@ -31,7 +31,7 @@ jobs:
     steps:
       - uses: actions/checkout@v2
 
-      - uses: r-lib/actions/setup-r@v1
+      - uses: r-lib/actions/setup-r@v2
         with:
           r-version: ${{ matrix.config.r }}
 
@@ -78,4 +78,4 @@ jobs:
         uses: actions/upload-artifact@main
         with:
           name: ${{ runner.os }}-r${{ matrix.config.r }}-results
-          path: check
+          path: check

DESCRIPTION

@@ -1,13 +1,37 @@
 Package: malariasimulation
 Title: An individual based model for malaria
-Version: 1.2.0
+Version: 1.3.0
 Authors@R: c(
   person(
     given = "Giovanni",
     family = "Charles",
     role = c('aut', 'cre'),
     email = '[email protected]'
   ),
+  person(
+    given = "Peter",
+    family = "Windskill",
+    role = c('aut'),
+    email = '[email protected]'
+  ),
+  person(
+    given = "Hillary",
+    family = "Topazian",
+    role = c('aut'),
+    email = '[email protected]'
+  ),
+  person(
+    given = "Joseph",
+    family = "Challenger",
+    role = c('aut'),
+    email = '[email protected]'
+  ),
+  person(
+    given = "Richard",
+    family = "Fitzjohn",
+    role = c('aut'),
+    email = '[email protected]'
+  ),
   person(
     given = "Imperial College of Science, Technology and Medicine",
     family = "",
@@ -20,9 +44,10 @@ Encoding: UTF-8
 LazyData: true
 Remotes:
   mrc-ide/malariaEquilibrium,
+  mrc-ide/individual
 Imports:
   individual (>= 0.1.7),
-  malariaEquilibrium,
+  malariaEquilibrium (>= 1.0.1),
   Rcpp,
   statmod,
   MASS,

NAMESPACE

@@ -4,6 +4,7 @@ export(AL_params)
 export(DHA_PQP_params)
 export(SP_AQ_params)
 export(arab_params)
+export(find_birthrates)
 export(fun_params)
 export(gamb_params)
 export(get_correlation_parameters)
@@ -15,6 +16,7 @@ export(run_simulation)
 export(run_simulation_with_repetitions)
 export(set_bednets)
 export(set_clinical_treatment)
+export(set_demography)
 export(set_drugs)
 export(set_equilibrium)
 export(set_mass_rtss)
@@ -32,5 +34,6 @@ importFrom(stats,rexp)
 importFrom(stats,rnorm)
 importFrom(stats,rpois)
 importFrom(stats,runif)
+importFrom(stats,uniroot)
 importFrom(stats,weighted.mean)
 useDynLib(malariasimulation, .registration = TRUE)

NEWS.md

@@ -1,3 +1,11 @@
+# malariasimulation 1.3.0
+
+  * Custom demography
+  * Rainfall floor
+  * Documentation fixes for vector controls
+  * Simple model for severe disease
+  * Alignment of immunity function with legacy model
+
 # malariasimulation 1.2.0
 
   * added a `news.md` file to track changes to the package.

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) {
-    .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)
+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)
 }
 
 aquatic_mosquito_model_update <- function(model, total_M, f, mum) {
@@ -37,16 +37,16 @@ history_push <- function(history, value, timestep) {
     invisible(.Call(`_malariasimulation_history_push`, history, value, timestep))
 }
 
-carrying_capacity <- function(timestep, model_seasonality, g0, g, h, K0, R_bar) {
-    .Call(`_malariasimulation_carrying_capacity`, timestep, model_seasonality, g0, g, h, K0, R_bar)
+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)
 }
 
 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)
+rainfall <- function(t, g0, g, h, floor) {
+    .Call(`_malariasimulation_rainfall`, t, g0, g, h, floor)
 }
 
 solver_get_states <- function(solver) {

R/compartmental.R

@@ -24,7 +24,8 @@ parameterise_mosquito_models <- function(parameters) {
         parameters$h,
         calculate_R_bar(parameters),
         parameters$mum[[i]],
-        parameters$blood_meal_rates[[i]]
+        parameters$blood_meal_rates[[i]],
+        parameters$rainfall_floor
       )
 
       if (!parameters$individual_mosquitoes) {
@@ -87,7 +88,11 @@ create_compartmental_rendering_process <- function(renderer, solvers, parameters
   function(timestep) {
     counts <- rep(0, length(indices))
     for (s_i in seq_along(solvers)) {
-      row <- solver_get_states(solvers[[s_i]])
+      if (parameters$species_proportions[[s_i]] > 0) {
+        row <- solver_get_states(solvers[[s_i]])
+      } else {
+        row <- rep(0, length(indices))
+      }
       for (i in seq_along(indices)) {
         renderer$render(
           paste0(names(indices)[[i]], '_', parameters$species[[s_i]], '_count'),
@@ -105,10 +110,12 @@ create_compartmental_rendering_process <- function(renderer, solvers, parameters
 #'
 #' @param solvers for each species
 #' @noRd
-create_solver_stepping_process <- function(solvers) {
+create_solver_stepping_process <- function(solvers, parameters) {
   function(timestep) {
-    for (solver in solvers) {
-      solver_step(solver)
+    for (i in seq_along(solvers)) {
+      if (parameters$species_proportions[[i]] > 0) {
+        solver_step(solvers[[i]])
+      }
     }
   }
 }

R/compatibility.R

@@ -51,7 +51,18 @@ translations = list(
   tau = 'dem',
   Q0 = 'Q0',
   f = 'blood_meal_rates',
-  tl = 'delay_gam'
+  tl = 'delay_gam',
+  uv = 'uv',
+  dv = 'rva',
+  PVM = 'pvm',
+  dvm = 'rvm',
+  fv0 = 'fv0',
+  av = 'av',
+  gammav = 'gammav',
+  theta0 = 'theta0',
+  theta1 = 'theta1',
+  IV0 = 'iv0',
+  kv = 'kv'
 )
 
 back_translations = list(
@@ -95,7 +106,18 @@ back_translations = list(
   dem = 'tau',
   Q0 = mean_param('Q0', 'Q0', 'species_proportions'),
   blood_meal_rates = mean_param('f', 'blood_meal_rates', 'species_proportions'),
-  delay_gam = 'tl'
+  delay_gam = 'tl',
+  uv = 'uv',
+  rva = 'dv',
+  pvm = 'PVM',
+  rvm = 'dvm',
+  fv0 = 'fv0',
+  av = 'av',
+  gammav = 'gammav',
+  theta0 = 'theta0',
+  theta1 = 'theta1',
+  IV0 = 'iv0',
+  kv = 'kv'
 )
 
 #' @description translate parameter keys from the malariaEquilibrium format
@@ -162,7 +184,8 @@ remove_unused_equilibrium <- function(params) {
 #' equilibrium parameters and set up the initial human and mosquito population
 #' to acheive init_EIR
 #' @param parameters model parameters to update
-#' @param init_EIR the desired initial EIR
+#' @param init_EIR the desired initial EIR (infectious bites per day over the entire human
+#' population)
 #' @param eq_params parameters from the malariaEquilibrium package, if null.
 #' The default malariaEquilibrium parameters will be used
 #' @export

R/events.R

@@ -132,11 +132,6 @@ attach_event_listeners <- function(
       parameters
     )
   )
-  events$asymptomatic_progression$add_listener(
-    function(timestep, target) {
-      variables$is_severe$queue_update('no', target)
-    }
-  )
 
   # Recovery events
   events$subpatent_progression$add_listener(
@@ -147,6 +142,7 @@ attach_event_listeners <- function(
       parameters$cu
     )
   )
+
   events$subpatent_progression$add_listener(
     create_rate_listener('A', 'U', renderer)
   )

R/human_infection.R

@@ -62,16 +62,13 @@ simulate_infection <- function(
     timestep
   )
 
-  if (parameters$severe_enabled) {
-    update_severe_disease(
-      timestep,
-      clinical_infections,
-      variables,
-      infected_humans,
-      parameters,
-      renderer
-    )
-  }
+  update_severe_disease(
+    timestep,
+    infected_humans,
+    variables,
+    parameters,
+    renderer
+  )
 
   treated <- calculate_treated(
     variables,
@@ -210,55 +207,49 @@ calculate_clinical_infections <- function(
 #' @description
 #' Sample severely infected humans from clinically infected
 #' @param timestep current timestep
-#' @param clinical_infections indices of clinically infected humans
+#' @param infections indices of all infected humans
 #' @param variables a list of all of the model variables
-#' @param infections indices of all infected humans (for immunity boosting)
 #' @param parameters model parameters
 #' @param renderer model outputs
 #' @noRd
 update_severe_disease <- function(
   timestep,
-  clinical_infections,
-  variables,
   infections,
+  variables,
   parameters,
   renderer
   ) {
-  if (clinical_infections$size() > 0) {
-    age <- get_age(variables$birth$get_values(clinical_infections), timestep)
-    iva <- variables$iva$get_values(clinical_infections)
+  if (infections$size() > 0) {
+    age <- get_age(variables$birth$get_values(infections), timestep)
+    iva <- variables$iva$get_values(infections)
+    ivm <- variables$ivm$get_values(infections)
     theta <- severe_immunity(
       age,
       iva,
-      variables$ivm$get_values(clinical_infections),
+      ivm,
       parameters
     )
     develop_severe <- bernoulli_multi_p(theta)
-    severe_infections <- bitset_at(clinical_infections, develop_severe)
-    variables$is_severe$queue_update('yes', severe_infections)
-    variables$is_severe$queue_update(
-      'no',
-      severe_infections$copy()$not(TRUE)$and(infections)
-    ) # Prevent any unwanted severe reinfections
+    severe_infections <- bitset_at(infections, develop_severe)
     incidence_renderer(
       variables$birth,
       renderer,
       severe_infections,
-      clinical_infections,
+      infections,
       theta,
       'inc_severe_',
       parameters$severe_incidence_rendering_min_ages,
       parameters$severe_incidence_rendering_max_ages,
       timestep
     )
-    boost_immunity(
-      variables$iva,
-      infections,
-      variables$last_boosted_iva,
-      timestep,
-      parameters$uv
-    )
   }
+  boost_immunity(
+    variables$iva,
+    infections,
+    variables$last_boosted_iva,
+    timestep,
+    parameters$uv
+  )
 }
 
 #' @title Calculate treated humans
@@ -386,6 +377,8 @@ boost_immunity <- function(
 
 # Implemented from Winskill 2017 - Supplementary Information page 4
 clinical_immunity <- function(acquired_immunity, maternal_immunity, parameters) {
+  acquired_immunity[acquired_immunity > 0] <- acquired_immunity[acquired_immunity > 0] + 0.5
+
   parameters$phi0 * (
     parameters$phi1 +
       (1 - parameters$phi1) /
@@ -399,6 +392,8 @@ clinical_immunity <- function(acquired_immunity, maternal_immunity, parameters)
 
 # Implemented from Winskill 2017 - Supplementary Information page 5
 severe_immunity <- function(age, acquired_immunity, maternal_immunity, parameters) {
+  acquired_immunity[acquired_immunity > 0] <- acquired_immunity[acquired_immunity > 0] + 0.5
+
   fv <- 1 - (1 - parameters$fv0) / (
     1 + (age / parameters$av) ** parameters$gammav
   )
@@ -431,6 +426,8 @@ asymptomatic_infectivity <- function(age, immunity, parameters) {
 
 # Implemented from Winskill 2017 - Supplementary Information page 4
 blood_immunity <- function(ib, parameters) {
+  ib[ib > 0] <- ib[ib > 0] + 0.5
+
   parameters$b0 * (
     parameters$b1 +
       (1 - parameters$b1) /

R/model.R

@@ -56,8 +56,6 @@
 #'  * p_inc_severe: the sum of probabilities of severe infection for humans between an inclusive age range at this timestep.
 #' severe incidence columns can be set with
 #' severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
-#'  * severe_deaths: number of deaths due to severe malaria. severe_enabled must be
-#' set to TRUE
 #'  * E_count: number of mosquitoes in the early larval stage (per species)
 #'  * L_count: number of mosquitoes in the late larval stage (per species)
 #'  * P_count: number of mosquitoes in the pupal stage (per species)