Merge pull request #273 from InstituteforDiseaseModeling/charles/issu…

…e119

Rate limits caps for fertility, mortality, and prevalence rates

pacehrh/R/pace_exp_control.R

@@ -20,6 +20,7 @@
 #' pacehrh::InitializeScenarios()
 #' pacehrh::InitializeStochasticParameters()
 #' pacehrh::InitializeSeasonality()
+#' pacehrh::InitializeCadreRoles()
 #'
 #' scenario <- "ScenarioName"
 #'
@@ -148,7 +149,9 @@ SaveBaseSettings <- function(scenarioName = ""){
                    expression(BVE$stochasticParams),
                    expression(BVE$initialPopulation),
                    expression(BVE$populationLabels),
-                   expression(BVE$populationRangesTable)
+                   expression(BVE$populationRangesTable),
+                   expression(BVE$cadreRoles)
+#                   expression(BVE$changeRateLimits)
                    )
 
   checks <- sapply(varsToCheck, function(v){
@@ -302,7 +305,7 @@ ConfigureExperimentValues <- function(){
 #' \dontrun{
 #' initPop <- pacehrh:::loadInitialPopulation(sheetName = "Flat_Population")
 #' pcr <- pacehrh:::loadPopulationChangeRates(sheetName = "Flat_Rates")
-#' pars <- pacehrh:::loadStochasticParameters(sheetName = "Flat_StochasticParms")
+#' pars <- pacehrh:::loadStochasticParameters(stochasticParametersSheetName = "Flat_StochasticParms")
 #' years <- 2020:2040
 #' pcr <- pacehrh:::addRatesMatricesToPopulationChangeRates(pcr, years, NULL)
 #' }

pacehrh/R/pace_initialize.R

@@ -61,7 +61,7 @@ PaceInitialize <- function(globalConfigFile = NULL,
   if (is.null(stochasticParametersSheet)) {
     InitializeStochasticParameters()
   } else {
-    InitializeStochasticParameters(sheetName = stochasticParametersSheet)
+    InitializeStochasticParameters(stochasticParametersSheetName = stochasticParametersSheet)
   }
 
   # Local seasonality curves sheet and seasonality offsets sheets

pacehrh/R/pace_package_env.R

@@ -60,6 +60,7 @@ NULL
 .defaultStochasticParametersSheet <- "StochasticParameters"
 .defaultCadreRolesSheet <- "CadreRoles"
 .defaultTaskCadresSheet <- "TaskCadres"
+.defaultChangeRateLimitsSheet <- "ChangeRateLimits"
 
 # GLOBAL VARIABLES
 

pacehrh/R/pace_prevalence_matrix.R

@@ -5,12 +5,16 @@
 
 #' Generate A Matrix Of Stochastically Varied Prevalence/Incidence Rates
 #'
+#' @param debugEnv Caller-create environment into which the
+#'   `generatePrevalenceRatesMatrix()` call can insert some debugging
+#'   information to return to the caller. Intended for unit testing.
+#'
 #' @return Rates matrix
 #'
 #' @noRd
-generatePrevalenceRatesMatrix <- function(){
+generatePrevalenceRatesMatrix <- function(debugEnv = NULL){
   if (!GPE$stochasticity){
-    return(.generateNonStochPrm())
+    return(.generateNonStochPrm(debugEnv))
   }
 
   # Gather stuff we're going to need
@@ -33,24 +37,23 @@ generatePrevalenceRatesMatrix <- function(){
   # Mask for tasks affected by prevalence stochasticity
   sMask <- tasks$applyStochasticity
 
-  # Determine which rows refer to Malaria/HIV/TB
-  mhivtbMask <-
-    (
-      tasks$ServiceCat == "Malaria" |
-        tasks$ServiceCat == "Tuberculosis" | tasks$ServiceCat == "HIV"
-    )
-
-  # Determine which rows refer to childhood diseases
-  childDisMask <- (tasks$RelevantPop == "1-4")
-
   # Grab the initial prevalence/incidence values, and apply a stochastic tweak
   p = pars[pars$Value == "Incidence rates", ]$p
 
   initRates <- rep(1.0, nRows)
-
   initRates[sMask] <-
     .varyInitialIncidenceRates(tasks[sMask,]$StartingRateInPop, p)
 
+  # Calculate min/max constraints on rates
+  limits <- .getRatesLimits("Incidence")
+  rateLimits <- .getMinMaxRates(initRates, limits)
+  maskedRateLimits <- lapply(rateLimits, function(x){x[sMask]})
+
+  # Write some information back to the user-provided debugging environment
+  if (!is.null(debugEnv)) {
+    assign("rateLimits", rateLimits, envir = debugEnv)
+  }
+
   # Initialize a rates matrix
   m <-
     matrix(
@@ -66,17 +69,7 @@ generatePrevalenceRatesMatrix <- function(){
   q = pars[pars$Value == "Annual delta incidence rates", ]$q
   lims = p * q * c(-1, 1)
 
-  # Grab the annual prevalence deltas
-  deltaRatios <- tasks$AnnualDeltaRatio
-
-  # Set annual prevalence deltas to unity if the scenario requires it
-  if (BVE$scenario$o_MHIVTB_decr == FALSE){
-    deltaRatios[mhivtbMask] <- 1
-  }
-
-  if (BVE$scenario$o_ChildDis_decr == FALSE){
-    deltaRatios[childDisMask] <- 1
-  }
+  deltaRatios <- .computeAdjustedDeltaRatios(tasks)
 
   for (j in 2:nCols){
     e <- truncnorm::rtruncnorm(
@@ -87,13 +80,38 @@ generatePrevalenceRatesMatrix <- function(){
       b = lims[2]
     )
     deltas <- deltaRatios[sMask] * (1 + e)
-    m[sMask, j] <- m[sMask, j-1] * deltas
+    m[sMask, j] <- .applyRateLimits(m[sMask, j-1] * deltas, maskedRateLimits)
   }
 
   return(m)
 }
 
-.generateNonStochPrm <- function(){
+.computeAdjustedDeltaRatios <- function(tasks) {
+  deltaRatios <- tasks$AnnualDeltaRatio
+
+  # Set annual prevalence deltas to unity if the scenario requires it
+  if (BVE$scenario$o_MHIVTB_decr == FALSE){
+    # Determine which rows refer to Malaria/HIV/TB
+    mhivtbMask <-
+      (
+        tasks$ServiceCat == "Malaria" |
+          tasks$ServiceCat == "Tuberculosis" | tasks$ServiceCat == "HIV"
+      )
+
+    deltaRatios[mhivtbMask] <- 1
+  }
+
+  if (BVE$scenario$o_ChildDis_decr == FALSE){
+    # Determine which rows refer to childhood diseases
+    childDisMask <- (tasks$RelevantPop == "1-4")
+
+    deltaRatios[childDisMask] <- 1
+  }
+
+  return(deltaRatios)
+}
+
+.generateNonStochPrm <- function(debugEnv = NULL){
   years <- BVE$years
 
   tasks <-
@@ -112,19 +130,18 @@ generatePrevalenceRatesMatrix <- function(){
   # Mask for tasks affected by prevalence stochasticity
   sMask <- tasks$applyStochasticity
 
-  # Determine which rows refer to Malaria/HIV/TB
-  mhivtbMask <-
-    (
-      tasks$ServiceCat == "Malaria" |
-        tasks$ServiceCat == "Tuberculosis" | tasks$ServiceCat == "HIV"
-    )
-
-  # Determine which rows refer to childhood diseases
-  childDisMask <- (tasks$RelevantPop == "1-4")
-
   initRates <- rep(1.0, nRows)
   initRates[sMask] <- tasks[sMask,]$StartingRateInPop
 
+  # Calculate min/max constraints on rates
+  limits <- .getRatesLimits("Incidence")
+  rateLimits <- .getMinMaxRates(initRates, limits)
+  maskedRateLimits <- lapply(rateLimits, function(x){x[sMask]})
+
+  if (!is.null(debugEnv)) {
+    assign("rateLimits", rateLimits, envir = debugEnv)
+  }
+
   # Initialize a rates matrix
   m <-
     matrix(
@@ -135,22 +152,12 @@ generatePrevalenceRatesMatrix <- function(){
     )
   m[, 1] <- initRates
 
-  # Grab the annual prevalence deltas
-  deltaRatios <- tasks$AnnualDeltaRatio
-
-  # Set annual prevalence deltas to unity if the scenario requires it
-  if (BVE$scenario$o_MHIVTB_decr == FALSE){
-    deltaRatios[mhivtbMask] <- 1
-  }
-
-  if (BVE$scenario$o_ChildDis_decr == FALSE){
-    deltaRatios[childDisMask] <- 1
-  }
+  deltaRatios <- .computeAdjustedDeltaRatios(tasks)
 
   # Derive the rest of the matrix
   deltas <- deltaRatios[sMask]
   for (j in 2:nCols){
-    m[sMask, j] <- m[sMask, j-1] * deltas
+    m[sMask, j] <- .applyRateLimits(m[sMask, j-1] * deltas, maskedRateLimits)
   }
 
   return(m)