Testing foi

R/events.R

@@ -18,12 +18,6 @@ create_events <- function(variables_list, parameters_list) {
 
   )
 
-  ## Add RS_event if endemic pathogen is required
-  if (parameters_list$endemic_or_epidemic == "endemic") {
-    RS_event <- TargetedEvent$new(population_size = parameters_list$human_population)
-    events_list <- c(events_list, list(RS_event = RS_event))
-  }
-
   # Add listener to the EI event:
   events_list$EI_event$add_listener(
     function(t, target) {
@@ -38,13 +32,18 @@ create_events <- function(variables_list, parameters_list) {
     }
   )
 
-  ## Add listener to the RS event if endemic pathogen is required
+  # Add RS_event and listener if endemic pathogen is required (i.e. individuals going R->S)
   if (parameters_list$endemic_or_epidemic == "endemic") {
+
+    RS_event <- TargetedEvent$new(population_size = parameters_list$human_population)
+    events_list <- c(events_list, list(RS_event = RS_event))
+
     events_list$RS_event$add_listener(
       function(t, target) {
         variables_list$disease_state$queue_update("S", target)
       }
     )
+
   }
 
   # Return the list of model events:

R/model.R

@@ -8,8 +8,8 @@ run_simulation <- function(parameters_list) {
 
   # Generate the model variables:
   variables_list <- create_variables(parameters_list)
-  parameters_list <- variables_list$parameters_list
-  variables_list <- variables_list$variables_list
+  parameters_list <- variables_list$parameters_list # note: this could be written more nicely and in a way
+  variables_list <- variables_list$variables_list   #       that doesn't require recursive modification
 
   # Generate the model events:
   events_list <- create_events(variables_list = variables_list, parameters_list = parameters_list)

R/parameters.R

@@ -145,6 +145,8 @@ get_parameters <- function(overrides = list()) {
     stop("duration_immune must be specified if endemic_or_epidemic is set to endemic")
   }
 
+  ## ADD MORE CHECKS IN HERE FOR PARAMETERS ##
+
   # Return the list of parameters
   parameters
 

R/processes.R

@@ -115,9 +115,6 @@ create_SE_process <- function(variables_list, events_list, parameters_list){
   function(t) {
 
     ## Timestep printing
-    # if (t %% 10 == 0) {
-    #   print(paste0("timestep = ", (t * parameters_list$dt)))
-    # }
     print(t)
 
     ## Bitset for all infectious individuals
@@ -135,7 +132,7 @@ create_SE_process <- function(variables_list, events_list, parameters_list){
       spec_household <- household_bitset_list[[i]]
 
       # Retrieve the bitset of all infectious individuals in the i-th household
-      spec_household_I <- I$and(spec_household)
+      spec_household_I <- I$copy()$and(spec_household)
 
       #  Calculate the FOI for the i-th household - with and without farUVC installed
       if (parameters_list$far_uvc_household) {
@@ -165,7 +162,7 @@ create_SE_process <- function(variables_list, events_list, parameters_list){
       spec_workplace <- workplace_bitset_list[[i]]
 
       # Get the indices of infectious individuals in the i-th workplace:
-      spec_workplace_I <- I$and(spec_workplace)
+      spec_workplace_I <- I$copy()$and(spec_workplace)
 
       # Calculate the workplace-specific FOI of the i-th workplace - with and without farUVC installed
       if (parameters_list$far_uvc_workplace) {
@@ -195,7 +192,7 @@ create_SE_process <- function(variables_list, events_list, parameters_list){
       spec_school <- school_bitset_list[[i]]
 
       # Retrieve the indices of all infectious individuals in the i-th school:
-      spec_school_I <- I$and(spec_school)
+      spec_school_I <- I$copy()$and(spec_school)
 
       # Calculate the school-specific FOI for the i-th school - with and without farUVC installed
       if (parameters_list$far_uvc_school) {
@@ -256,7 +253,7 @@ create_SE_process <- function(variables_list, events_list, parameters_list){
         spec_leisure <- variables_list$specific_leisure$get_index_of(as.character(spec_leisure_setting))
 
         # Retrieve the indices of all infectious individuals in the particular leisure setting being considered
-        spec_leisure_I <- I$and(spec_leisure)
+        spec_leisure_I <- I$copy()$and(spec_leisure)
 
         # Calculate the leisure-specific FOI for the i-th leisure setting - with and without farUVC installed
         if (parameters_list$far_uvc_leisure) {
@@ -287,7 +284,7 @@ create_SE_process <- function(variables_list, events_list, parameters_list){
     # Sum the household, workplace, school leisure, and community FOIs to get the total FOI for each
     # individual:
     total_FOI <- household_FOI + workplace_FOI + school_FOI + leisure_FOI + community_FOI
-    print(c(max(household_FOI), max(workplace_FOI), max(school_FOI), max(leisure_FOI), max(community_FOI), max(total_FOI)))
+    print(c(max(household_FOI), max(workplace_FOI), max(school_FOI), max(leisure_FOI), max(community_FOI)))
 
     # Calculate the probability of getting infected in the current interval for each individual:
     p_inf <- 1 - exp(-total_FOI * parameters_list$dt)

R/variables.R

@@ -11,8 +11,9 @@ create_variables <- function(parameters_list) {
   disease_state_variable <- individual::CategoricalVariable$new(categories = disease_states,
                                                                 initial_values = initial_disease_states)
 
-  # If user wants to use empirical distribution of households from ONS - this generates
-  # both the household variable AND the age class variable
+  # Initialise and populate the age and household variables
+
+  # If user wants to use empirical distribution of households and ages from ONS sample
   if (parameters_list$household_distribution_generation == "empirical") {
 
     # Bootstrap sampling of households from ONS 2011 Census reference panel of household sizes and age composition
@@ -27,7 +28,7 @@ create_variables <- function(parameters_list) {
     household_variable <- individual::CategoricalVariable$new(categories = as.character(1:max(household_age_list$individual_households)),
                                                               initial_values = as.character(household_age_list$individual_households))
 
-    ## If user wants to specify age-class proportions manually
+  # If user wants to specify age-class proportions and associated households manually
   } else {
 
     # Initialise and populate the age class variable:
@@ -54,7 +55,7 @@ create_variables <- function(parameters_list) {
   school_variable <- CategoricalVariable$new(categories = as.character(0:num_schools),
                                              initial_values = initial_school_settings)
 
-  # Initialise and populate the leisure setting variable that stores all the leisure locations individual COULD go to
+  # Initialise and populate the leisure setting variable that stores all the leisure locations an individual COULD go to
   # Generating the number and sizes of each leisure setting
   leisure_setting_sizes <- sample_negbinom(N = parameters_list$human_population,
                                            prop_max = parameters_list$leisure_prop_max,
@@ -63,7 +64,7 @@ create_variables <- function(parameters_list) {
   initial_leisure_settings <- generate_initial_leisure(parameters_list = parameters_list, leisure_setting_sizes = leisure_setting_sizes) # returns list to initialise RaggedInteger
   leisure_variable <- RaggedInteger$new(initial_values = initial_leisure_settings)
 
-  # Initialise and populate the leisure setting variable for where individuals go to on a particular day
+  # Initialise and populate the leisure setting variable for where individuals go to on a particular day (this is dynamically updated in the model)
   possible_leisure_settings <- unique(unlist(initial_leisure_settings))
   possible_leisure_settings <- possible_leisure_settings[order(possible_leisure_settings)]
   specific_day_leisure_variable <- CategoricalVariable$new(categories = as.character(possible_leisure_settings),
@@ -80,8 +81,7 @@ create_variables <- function(parameters_list) {
     specific_leisure = specific_day_leisure_variable
   )
 
-  # If any setting has UVC parameterised, retrieve the setting sizes:
-  # Decide whether we want to do this step regardless of UVC and append sizes to variables_list
+  # If any setting has UVC installed, retrieve the sizes of all of the settings:
   if(any(parameters_list$far_uvc_workplace,
          parameters_list$far_uvc_school,
          parameters_list$far_uvc_leisure,

inst/farUVC_testing.R

@@ -11,68 +11,71 @@ source("R/model.R")
 source("R/utils.R")
 
 ## Running the model with no farUVC
-nofarUVC_parameters_list <- get_parameters(overrides = list(beta_household = 1,
-                                                            beta_workplace = 1,
-                                                            beta_school = 1,
-                                                            beta_community = 1,
-                                                            beta_leisure = 1,
-                                                            simulation_time = 35))
+nofarUVC_parameters_list <- get_parameters(overrides = list(beta_household = 0.2,
+                                                            beta_workplace = 0.2,
+                                                            beta_school = 0.2,
+                                                            beta_leisure = 0.2,
+                                                            beta_community = 0.05,
+                                                            simulation_time = 80,
+                                                            number_initially_exposed = 50))
 nofarUVC_output <- run_simulation(nofarUVC_parameters_list)
 health_cols <-  c("royalblue3","firebrick3","darkorchid3", "orange2")
 matplot(
   x = nofarUVC_output[[1]]*nofarUVC_parameters_list$dt, y = nofarUVC_output[-1],
   type="l",lwd=2,lty = 1,col = adjustcolor(col = health_cols, alpha.f = 0.85),
   xlab = "Time",ylab = "Count"
 )
-
-
-
+# betas = 0.1 ~ AR = 38% (finished by 200)
+# betas = 0.2 ~ R0 =
 
 ## Running the model with farUVC
-farUVC_parameters_list <- get_parameters(overrides = list(beta_household = 1,
-                                                          beta_workplace = 1,
-                                                          beta_school = 1,
-                                                          beta_community = 0,
-                                                          beta_leisure = 1,
-                                                          simulation_time = 35,
+farUVC_parameters_list <- get_parameters(overrides = list(beta_household = 0.1,
+                                                          beta_workplace = 0.1,
+                                                          beta_school = 0.1,
+                                                          beta_leisure = 0.1,
+                                                          beta_community = 0.05,
+                                                          simulation_time = 400,
                                                           far_uvc_workplace = TRUE,
                                                           far_uvc_workplace_coverage_type = "random",
-                                                          far_uvc_workplace_coverage = 1,
-                                                          far_uvc_workplace_efficacy = 1,
+                                                          far_uvc_workplace_coverage = 0.25,
+                                                          far_uvc_workplace_efficacy = 0.6,
                                                           far_uvc_workplace_timestep = 0,
                                                           far_uvc_school = TRUE,
                                                           far_uvc_school_coverage_type = "random",
-                                                          far_uvc_school_coverage = 1,
-                                                          far_uvc_school_efficacy = 1,
+                                                          far_uvc_school_coverage = 0.25,
+                                                          far_uvc_school_efficacy = 0.6,
                                                           far_uvc_school_timestep = 0,
                                                           far_uvc_leisure = TRUE,
                                                           far_uvc_leisure_coverage_type = "random",
-                                                          far_uvc_leisure_coverage = 1,
-                                                          far_uvc_leisure_efficacy = 1,
+                                                          far_uvc_leisure_coverage = 0.25,
+                                                          far_uvc_leisure_efficacy = 0.6,
                                                           far_uvc_leisure_timestep = 0,
                                                           far_uvc_household = TRUE,
                                                           far_uvc_household_coverage_type = "random",
-                                                          far_uvc_household_coverage = 1,
-                                                          far_uvc_household_efficacy = 1,
+                                                          far_uvc_household_coverage = 0.25,
+                                                          far_uvc_household_efficacy = 0.6,
                                                           far_uvc_household_timestep = 0))
 farUVC_output <- run_simulation(farUVC_parameters_list)
-# parameters_list <- farUVC_parameters_list
-
 
 par(mfrow = c(1, 2))
 health_cols <-  c("royalblue3","firebrick3","darkorchid3", "orange2")
 matplot(
-  x = farUVC_output[[1]]*farUVC_parameters_list$dt, y = farUVC_output[-1],
+  x = nofarUVC_output[[1]]*nofarUVC_parameters_list$dt, y = nofarUVC_output[-c(1, 2, 5)],
   type="l",lwd=2,lty = 1,col = adjustcolor(col = health_cols, alpha.f = 0.85),
   xlab = "Time",ylab = "Count"
 )
 matplot(
-  x = nofarUVC_output[[1]]*nofarUVC_parameters_list$dt, y = nofarUVC_output[-1],
+  x = farUVC_output[[1]]*farUVC_parameters_list$dt, y = farUVC_output[-c(1, 2, 5)],
   type="l",lwd=2,lty = 1,col = adjustcolor(col = health_cols, alpha.f = 0.85),
   xlab = "Time",ylab = "Count"
 )
 
+plot(nofarUVC_output[[1]]*nofarUVC_parameters_list$dt, nofarUVC_output[, 3], type = "l")
+lines(farUVC_output[[1]]*farUVC_parameters_list$dt, farUVC_output[, 3], col = "red")
 
+par(mfrow = c(1, 1))
+plot(nofarUVC_output[[1]]*nofarUVC_parameters_list$dt, nofarUVC_output[, 5], type = "l", xlab = "Time (Days)", ylab = "Total Infected")
+lines(farUVC_output[[1]]*farUVC_parameters_list$dt, farUVC_output[, 5], col = "red")
 
 
 parameters_list$seed <- 10

tests/testthat/test-interventions.R

@@ -42,12 +42,12 @@ test_that("set_uvc() errors when setting input not from allowed list of Far UVC
   # Check that set_uvc() errors when the setting input not either "workplace",
   # "school", or "leisure"
   expect_error(object = set_uvc(parameters_list = parameters,
-                                setting = "hosptial",
+                                setting = "hospital",
                                 coverage = c(0.8),
                                 coverage_type = "random",
                                 efficacy = c(0.8),
                                 timestep = c(1)),
-               regexp = "Error: Input setting invalid - far UvC only deployable in workplace, school, leisure, or household settings")
+               regexp = "Error: Input setting invalid - far UVC only deployable in workplace, school, leisure, or household settings")
 })
 
 test_that("set_uvc() errors when coverage_type input not from allowed list of Far UVC settings", {