Merge pull request #88 from mrc-ide/feat/prevalence

Feat/prevalence

R/events.R

@@ -3,8 +3,8 @@ create_events <- function(parameters) {
     # 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),
 
@@ -149,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,
@@ -162,7 +170,7 @@ attach_event_listeners <- function(
     )
   )
 
-  events$infection$add_listener(
+  events$detection$add_listener(
     create_incidence_renderer(
       variables$birth,
       variables$is_severe,

R/human_infection.R

@@ -18,7 +18,6 @@ simulate_infection <- function(
   timestep,
   renderer
   ) {
-
   if (bitten_humans$size() > 0) {
     boost_immunity(
       variables$ib,
@@ -37,11 +36,27 @@ simulate_infection <- function(
     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,
-    timestep
+    parameters
   )
 
   if (parameters$severe_enabled) {
@@ -59,6 +74,7 @@ simulate_infection <- function(
     variables,
     clinical_infections,
     events$recovery,
+    events$detection,
     parameters,
     timestep,
     renderer
@@ -143,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))
 }
@@ -224,6 +221,7 @@ calculate_treated <- function(
   variables,
   clinical_infections,
   recovery,
+  detection,
   parameters,
   timestep,
   renderer
@@ -252,14 +250,13 @@ calculate_treated <- function(
 
   # 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(
@@ -270,6 +267,7 @@ calculate_treated <- function(
       treated_index,
       log_uniform(treated_index$size(), parameters$dt)
     )
+    detection$schedule(treated_index, 0)
   }
   treated_index
 }
@@ -290,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)
@@ -302,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)
   }
@@ -313,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)
   }

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)

R/parameters.R

@@ -405,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),

R/processes.R

@@ -95,6 +95,7 @@ create_processes <- function(
       variables$state,
       variables$birth,
       variables$is_severe,
+      variables$id,
       parameters,
       renderer
     ),

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,

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))

src/Random.cpp

@@ -15,12 +15,12 @@ std::vector<size_t> Random::bernoulli(size_t size, double p) {
     return Rcpp::as<std::vector<size_t>>(indices);
 }
 
-individual_index_t Random::bernoulli_multi_p(const std::vector<double> p) {
-    auto successes = individual_index_t(p.size());
+std::vector<size_t> Random::bernoulli_multi_p(const std::vector<double> p) {
+    auto successes = std::vector<size_t>();
     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;

src/Random.h

@@ -10,12 +10,11 @@
 #define SRC_RANDOM_H_
 
 #include <vector>
-#include <individual.h>
 
 class RandomInterface {
 public:
     virtual std::vector<size_t> bernoulli(size_t, double) = 0;
-    virtual individual_index_t bernoulli_multi_p(const std::vector<double>) = 0;
+    virtual std::vector<size_t> bernoulli_multi_p(const std::vector<double>) = 0;
     virtual std::vector<size_t> sample(size_t, size_t, bool) = 0;
     virtual ~RandomInterface() = default;
 };
@@ -27,7 +26,7 @@ class Random : public RandomInterface {
         return instance;
     }
     virtual std::vector<size_t> bernoulli(size_t, double);
-    virtual individual_index_t bernoulli_multi_p(const std::vector<double>);
+    virtual std::vector<size_t> bernoulli_multi_p(const std::vector<double>);
     virtual std::vector<size_t> sample(size_t, size_t, bool);
     virtual ~Random() = default;
     Random(const Random &other) = delete;

src/RcppExports.cpp

@@ -177,7 +177,7 @@ BEGIN_RCPP
 END_RCPP
 }
 // bernoulli_multi_p_cpp
-Rcpp::XPtr<individual_index_t> bernoulli_multi_p_cpp(const std::vector<double> p);
+std::vector<size_t> bernoulli_multi_p_cpp(const std::vector<double> p);
 RcppExport SEXP _malariasimulation_bernoulli_multi_p_cpp(SEXP pSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;

src/utils.cpp

@@ -3,9 +3,10 @@
 #include "Random.h"
 
 //[[Rcpp::export]]
-Rcpp::XPtr<individual_index_t> bernoulli_multi_p_cpp(const std::vector<double> p) {
-    return Rcpp::XPtr<individual_index_t>(
-        new individual_index_t(Random::get_instance().bernoulli_multi_p(p)),
-        true
-    );
+std::vector<size_t> bernoulli_multi_p_cpp(const std::vector<double> p) {
+    auto values = Random::get_instance().bernoulli_multi_p(p);
+    for (auto i = 0u; i < values.size(); ++i) {
+        values[i]++;
+    }
+    return values;
 }