Merge pull request #151 from mrc-ide/dev

1.2.0

.Rbuildignore

@@ -9,4 +9,3 @@ docker
 LICENSE.md
 codecov.yml
 .github
-vignettes/*

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: malariasimulation
 Title: An individual based model for malaria
-Version: 1.1.0
+Version: 1.2.0
 Authors@R: c(
   person(
     given = "Giovanni",
@@ -19,9 +19,9 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 Remotes:
-  mrc-ide/malariaEquilibrium
+  mrc-ide/malariaEquilibrium,
 Imports:
-  individual,
+  individual (>= 0.1.7),
   malariaEquilibrium,
   Rcpp,
   statmod,
@@ -39,8 +39,9 @@ Suggests:
   DiagrammeR,
   cowplot,
   ggplot2,
-  covr
-RoxygenNote: 7.1.1
+  covr,
+  mgcv 
+RoxygenNote: 7.1.2
 Roxygen: list(markdown = TRUE)
 LinkingTo:
   Rcpp,

NEWS.md

@@ -0,0 +1,7 @@
+# malariasimulation 1.2.0
+
+  * added a `news.md` file to track changes to the package.
+  * n_inc_clinical includes treated humans
+  * disaggregate EIR and state counts by mosquito species
+  * remove redundant Total_M and and global EIR outputs
+  * new vignette on how to calibrate EIR to prevalence

R/RcppExports.R

@@ -65,6 +65,10 @@ bernoulli_multi_p_cpp <- function(p) {
     .Call(`_malariasimulation_bernoulli_multi_p_cpp`, p)
 }
 
+bitset_index_cpp <- function(a, b) {
+    .Call(`_malariasimulation_bitset_index_cpp`, a, b)
+}
+
 fast_weighted_sample <- function(size, probs) {
     .Call(`_malariasimulation_fast_weighted_sample`, size, probs)
 }

R/biting_process.R

@@ -85,12 +85,13 @@ simulate_bites <- function(
   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()
+    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]])
     p_bitten <- prob_bitten(timestep, variables, s_i, parameters)
     Q0 <- parameters$Q0[[s_i]]
@@ -118,20 +119,24 @@ simulate_bites <- function(
 
     lagged_eir[[s_i]]$save(n_infectious * a, timestep)
     species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
+    renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
-    n_bites <- rpois(1, species_eir * mean(psi))
-    if (n_bites > 0) {
-      bitten_humans$insert(
-        fast_weighted_sample(n_bites, lambda)
-      )
+    expected_bites <- species_eir * mean(psi)
+    if (expected_bites > 0) {
+      n_bites <- rpois(1, expected_bites)
+      if (n_bites > 0) {
+        bitten_humans$insert(
+          fast_weighted_sample(n_bites, lambda)
+        )
+      }
     }
 
     infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)
     lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
     foim <- calculate_foim(a, infectivity)
-    renderer$render(paste0('FOIM_', s_i), foim, timestep)
+    renderer$render(paste0('FOIM_', species_name), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
-    renderer$render(paste0('mu_', s_i), mu, timestep)
+    renderer$render(paste0('mu_', species_name), mu, timestep)
 
     if (parameters$individual_mosquitoes) {
       # update the ODE with stats for ovoposition calculations
@@ -167,7 +172,6 @@ simulate_bites <- function(
     }
   }
 
-  renderer$render('EIR', EIR, timestep)
   renderer$render('n_bitten', bitten_humans$size(), timestep)
   bitten_humans
 }
@@ -191,7 +195,7 @@ effective_biting_rates <- function(a, .pi, p_bitten) {
 
 calculate_infectious <- function(species, solvers, variables, parameters) {
   if (parameters$individual_mosquitoes) {
-    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not()
+    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not(TRUE)
     return(
       calculate_infectious_individual(
         species,
@@ -221,7 +225,7 @@ calculate_infectious_individual <- function(
 }
 
 calculate_infectious_compartmental <- function(solver_states) {
-  solver_states[[ADULT_ODE_INDICES['Im']]]
+  max(solver_states[[ADULT_ODE_INDICES['Im']]], 0)
 }
 
 intervention_coefficient <- function(p_bitten) {
@@ -234,7 +238,7 @@ human_pi <- function(zeta, psi) {
 
 blood_meal_rate <- function(v, z, parameters) {
   gonotrophic_cycle <- get_gonotrophic_cycle(v, parameters)
-  interrupted_foraging_time <- parameters$foraging_time / (1 - z)
+  interrupted_foraging_time <- parameters$foraging_time[[v]] / (1 - z)
   1 / (interrupted_foraging_time + gonotrophic_cycle)
 }
 

R/compartmental.R

@@ -77,7 +77,7 @@ parameterise_solvers <- function(models, parameters) {
   )
 }
 
-create_ode_rendering_process <- function(renderer, solvers, parameters) {
+create_compartmental_rendering_process <- function(renderer, solvers, parameters) {
   if (parameters$individual_mosquitoes) {
     indices <- ODE_INDICES
   } else {
@@ -86,18 +86,17 @@ create_ode_rendering_process <- function(renderer, solvers, parameters) {
 
   function(timestep) {
     counts <- rep(0, length(indices))
-    for (i in seq_along(solvers)) {
-      row <- solver_get_states(solvers[[i]])
+    for (s_i in seq_along(solvers)) {
+      row <- solver_get_states(solvers[[s_i]])
+      for (i in seq_along(indices)) {
+        renderer$render(
+          paste0(names(indices)[[i]], '_', parameters$species[[s_i]], '_count'),
+          row[[i]],
+          timestep
+        )
+      }
       counts <- counts + row
     }
-
-    for (i in seq_along(indices)) {
-      renderer$render(
-        paste0(names(indices)[[i]], '_count'),
-        counts[[i]],
-        timestep
-      )
-    }
   }
 }
 

R/events.R

@@ -9,9 +9,6 @@ create_events <- function(parameters) {
     clinical_infection = individual::TargetedEvent$new(parameters$human_population),
     asymptomatic_infection = individual::TargetedEvent$new(parameters$human_population),
 
-    # whether the infection is detected
-    detection = individual::TargetedEvent$new(parameters$human_population), 
-
     # MDA events
     mda_administer = individual::Event$new(),
     smc_administer = individual::Event$new(),
@@ -135,6 +132,11 @@ 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(
@@ -165,24 +167,6 @@ attach_event_listeners <- function(
   # Progression
   # ===========
   # When infection events fire, schedule the next stages of infection
-
-  events$clinical_infection$add_listener(
-    create_clinical_incidence_renderer(
-      variables$birth,
-      parameters,
-      renderer
-    )
-  )
-
-  events$detection$add_listener(
-    create_incidence_renderer(
-      variables$birth,
-      variables$is_severe,
-      parameters,
-      renderer
-    )
-  )
-
   if (parameters$individual_mosquitoes) {
     events$mosquito_infection$add_listener(
       individual::update_category_listener(variables$mosquito_state, 'Im')