Merge pull request #300 from mrc-ide/competing_hazards_recovery_infec…

…tion_resolution

Competing hazards recovery infection resolution

DESCRIPTION

@@ -70,7 +70,7 @@ Remotes:
   mrc-ide/malariaEquilibrium,
   mrc-ide/individual
 Imports:
-    individual (>= 0.1.16),
+    individual (>= 0.1.17),
     malariaEquilibrium (>= 1.0.1),
     Rcpp,
     statmod,

R/biting_process.R

@@ -15,6 +15,8 @@
 #' on the other populations
 #' @param mixing_index an index for this population's position in the
 #' lagged_infectivity list (default: 1)
+#' @param infection_outcome competing hazards object for infection rates
+#' @param timestep the current timestep
 #' @noRd
 create_biting_process <- function(
   renderer,
@@ -26,12 +28,12 @@ create_biting_process <- function(
   lagged_infectivity,
   lagged_eir,
   mixing_fn = NULL,
-  mixing_index = 1
+  mixing_index = 1,
+  infection_outcome
   ) {
   function(timestep) {
     # Calculate combined EIR
     age <- get_age(variables$birth$get_values(), timestep)
-
     bitten_humans <- simulate_bites(
       renderer,
       solvers,
@@ -54,7 +56,8 @@ create_biting_process <- function(
       age,
       parameters,
       timestep,
-      renderer
+      renderer,
+      infection_outcome
     )
   }
 }

R/competing_hazards.R

@@ -0,0 +1,95 @@
+## Define classes to resolve competing hazards
+CompetingOutcome <- R6::R6Class(
+  "CompetingOutcome",
+  private = list(
+    targeted_process = NULL
+  ),
+  public = list(
+    initialize = function(targeted_process, size){
+      if (!is.function(targeted_process)){
+        stop("targeted_process must be a function")
+      }
+      if (!is.numeric(size) || size <= 0){
+        stop("size must be positive integer")
+      }
+      private$targeted_process <- targeted_process
+
+      self$target <- individual::Bitset$new(size)
+      self$rates <- NULL
+    },
+    set_rates = function(target, rates){
+      stopifnot(target$size() == length(rates))
+
+      self$target$copy_from(target)
+      self$rates <- rates
+    },
+    execute = function(t, target){
+      private$targeted_process(t, target)
+    },
+    reset = function() {
+      self$target$clear()
+      self$rates <- NULL
+    },
+    target = NULL,
+    rates = NULL
+  )
+)
+
+CompetingHazard <- R6::R6Class(
+  "CompetingHazard",
+  private = list(
+    size = NULL,
+    outcomes = list(),
+    # RNG is passed in because mockery is not able to stub runif
+    # TODO: change when fixed
+    rng = NULL
+  ),
+  public = list(
+    initialize = function(size, outcomes, rng = runif){
+      if (length(outcomes) == 0){
+        stop("At least one outcome must be provided")
+      }
+      if (!all(sapply(outcomes, function(x) inherits(x, "CompetingOutcome")))){
+        stop("All outcomes must be of class CompetingOutcome")
+      }
+      private$size <- size
+      private$outcomes <- outcomes
+      private$rng <- rng
+    },
+    resolve = function(t){
+      candidates <- individual::Bitset$new(private$size)
+      for (o in private$outcomes) {
+        candidates$or(o$target)
+      }
+
+      rates <- matrix(ncol = length(private$outcomes), nrow = candidates$size(), 0)
+      for (i in seq_along(private$outcomes)) {
+        idx <- bitset_index(
+          candidates,
+          private$outcomes[[i]]$target)
+
+        rates[idx, i] <- private$outcomes[[i]]$rates
+      }
+
+      total_rates <- rowSums(rates)
+      probs <- rate_to_prob(total_rates) * (rates / total_rates)
+      probs[is.na(probs)] <- 0
+
+      rng <- private$rng(candidates$size())
+
+      cumulative <- rep(0, candidates$size())
+
+      for (o in seq_along(private$outcomes)) {
+        next_cumulative <- cumulative + probs[,o]
+        selected <- (rng > cumulative) & (rng <= next_cumulative)
+        cumulative <- next_cumulative
+
+        target <- bitset_at(candidates, selected)
+        if (target$size() > 0) {
+          private$outcomes[[o]]$execute(t, target)
+        }
+        private$outcomes[[o]]$reset()
+      }
+    }
+  )
+)

R/disease_progression.R

@@ -1,3 +1,70 @@
+#' @title Calculate recovery rates
+#' @description Calculates recovery rates for each individual in the population
+#' for storage in competing hazards object and subsequent resolution
+#'
+#' @param variables the available human variables
+#' @param recovery_outcome competing hazards object for recovery rates
+#' @noRd
+create_recovery_rates_process <- function(
+  variables,
+  recovery_outcome
+) {
+  function(timestep){
+    target <- variables$state$get_index_of("S")$not()
+    recovery_outcome$set_rates(
+      target,
+      variables$recovery_rates$get_values(target))
+  }
+}
+
+
+#' @title Disease progression outcomes (recovery)
+#' @description Following resolution of competing hazards, update state and
+#' infectivity of sampled individuals
+#'
+#' @param timestep the current timestep
+#' @param target the sampled recovering individuals
+#' @param variables the available human variables
+#' @param parameters model parameters
+#' @param renderer competing hazards object for recovery rates
+#' @noRd
+recovery_outcome_process <- function(
+    timestep,
+    target,
+    variables,
+    parameters,
+    renderer
+){
+
+  update_to_asymptomatic_infection(
+    variables,
+    parameters,
+    timestep,
+    variables$state$get_index_of("D")$and(target)
+  )
+
+  update_infection(
+    variables$state,
+    "U",
+    variables$infectivity,
+    parameters$cu,
+    variables$recovery_rates,
+    1/parameters$du,
+    variables$state$get_index_of("A")$and(target)
+  )
+
+  update_infection(
+    variables$state,
+    "S",
+    variables$infectivity,
+    0,
+    variables$recovery_rates,
+    0,
+    variables$state$get_index_of(c("U","Tr"))$and(target)
+  )
+
+}
+
 #' @title Update the state of an individual as infection events occur
 #' @description Randomly moves individuals towards the later stages of disease
 #' and updates their infectivity
@@ -8,63 +75,17 @@
 #' @param new_infectivity the new infectivity of the progressed individuals
 #' @noRd
 update_infection <- function(
-  state,
-  to_state,
-  infectivity,
-  new_infectivity,
-  to_move
-  ) {
-  state$queue_update(to_state, to_move)
-  infectivity$queue_update(new_infectivity, to_move)
-}
-
-create_progression_process <- function(
     state,
-    from_state,
     to_state,
-    rate,
     infectivity,
-    new_infectivity
+    new_infectivity,
+    recovery_rates,
+    new_recovery_rate,
+    to_move
 ) {
-  function(timestep) {
-
-    # Retrieve the indices of all individuals in the to_move state:
-    index <- state$get_index_of(from_state)
-
-    # If the length of rate is greater than 1 (when it's a variable):
-    if (inherits(rate, "DoubleVariable")) {
-      rate <- rate$get_values(index)
-    }
-
-    # Sample the individuals to be moved into a new Bitset using the transition rate(s):
-    to_move <- index$sample(1/rate)
-
-    # Update the infection status of those individuals who are moving:
-    update_infection(
-      state,
-      to_state,
-      infectivity,
-      new_infectivity,
-      to_move
-    )
-  }
-}
-
-create_asymptomatic_progression_process <- function(
-  state,
-  rate,
-  variables,
-  parameters
-  ) {
-  function(timestep) {
-    to_move <- state$get_index_of('D')$sample(1/rate)
-    update_to_asymptomatic_infection(
-      variables,
-      parameters,
-      timestep,
-      to_move
-    )
-  }
+  state$queue_update(to_state, to_move)
+  infectivity$queue_update(new_infectivity, to_move)
+  recovery_rates$queue_update(new_recovery_rate, to_move)
 }
 
 #' @title Modelling the progression to asymptomatic disease
@@ -75,11 +96,11 @@ create_asymptomatic_progression_process <- function(
 #' @param parameters model parameters
 #' @noRd
 update_to_asymptomatic_infection <- function(
-  variables,
-  parameters,
-  timestep,
-  to_move
-  ) {
+    variables,
+    parameters,
+    timestep,
+    to_move
+) {
   if (to_move$size() > 0) {
     variables$state$queue_update('A', to_move)
     new_infectivity <- asymptomatic_infectivity(
@@ -94,5 +115,9 @@ update_to_asymptomatic_infection <- function(
       new_infectivity,
       to_move
     )
+    variables$recovery_rates$queue_update(
+      1/parameters$da,
+      to_move
+    )
   }
 }