Merge pull request #30 from mrc-ide/Make_asc_rate_params_variable

Make asc rate params variable

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: cowflu
 Title: Cows With Flu
-Version: 1.0.1
+Version: 1.0.2
 Authors@R: c(person("Thomas", "Rawson", role = c("aut", "cre"),
                     email = "t.rawson@imperial.ac.uk"),
              person("Rich", "FitzJohn", role = "aut"),

inst/dust/cows.cpp

@@ -42,12 +42,11 @@ void sum_over_regions(real_type *cows,
 
 struct outbreak_detection_parameters {
   bool proportion_only;
-  // These all fixed for now:
-  double proportion_scaling = 10;
-  double N_scaling = 0.7;
-  double strength_scaling = 0.95;
-  double I_scaling = 150;
-} outbreak_detection ;
+  double N_scaling;
+  double proportion_scaling;
+  double strength_scaling;
+  double I_scaling;
+};
 
 
 template <typename real_type, typename rng_state_type>
@@ -440,7 +439,12 @@ class cows {
     const auto likelihood_choice = read_likelihood_type(pars, "likelihood_choice");
 
     const bool outbreak_detection_proportion_only = dust2::r::read_bool(pars, "outbreak_detection_proportion_only", false);
-    const auto outbreak_detection_parameters{outbreak_detection_proportion_only};
+    const real_type N_scaling = dust2::r::read_real(pars, "N_scaling", 0.7);
+    const real_type proportion_scaling = dust2::r::read_real(pars, "proportion_scaling", 10);
+    const real_type strength_scaling = dust2::r::read_real(pars, "strength_scaling", 0.95);
+    const real_type I_scaling = dust2::r::read_real(pars, "I_scaling", 150);
+
+    const outbreak_detection_parameters outbreak_detection{outbreak_detection_proportion_only, N_scaling, proportion_scaling, strength_scaling, I_scaling};
 
     const size_t n_seed = dust2::r::read_size(pars, "n_seed");
     std::vector<size_t> seed_time(n_seed);
@@ -462,7 +466,7 @@ class cows {
     // dust2::r::read_real_vector(pars, n_seed, seed_herd.data(), "seed_herd", true);
     // dust2::r::read_real_vector(pars, n_seed, seed_amount.data(), "seed_amount", true);
 
-    return shared_state{n_seed, seed_time, seed_herd, seed_amount, n_herds, n_regions, gamma, sigma, beta, alpha, time_test, n_test, likelihood_choice, region_start, herd_to_region_lookup, p_region_export, p_cow_export, n_cows_per_herd, movement_matrix, start_count, start_herd, asc_rate, dispersion, condition_on_export, outbreak_detection_parameters};
+    return shared_state{n_seed, seed_time, seed_herd, seed_amount, n_herds, n_regions, gamma, sigma, beta, alpha, time_test, n_test, likelihood_choice, region_start, herd_to_region_lookup, p_region_export, p_cow_export, n_cows_per_herd, movement_matrix, start_count, start_herd, asc_rate, dispersion, condition_on_export, outbreak_detection};
   }
 
   static internal_state build_internal(const shared_state& shared) {
@@ -494,6 +498,12 @@ class cows {
     } else {
       dust2::r::read_real_vector(pars, shared.n_regions, shared.asc_rate.data(), "asc_rate", false);
     }
+
+    shared.outbreak_detection.N_scaling = dust2::r::read_real(pars, "N_scaling", shared.outbreak_detection.N_scaling);
+    shared.outbreak_detection.proportion_scaling = dust2::r::read_real(pars, "proportion_scaling", shared.outbreak_detection.proportion_scaling);
+    shared.outbreak_detection.strength_scaling = dust2::r::read_real(pars, "strength_scaling", shared.outbreak_detection.strength_scaling);
+    shared.outbreak_detection.I_scaling = dust2::r::read_real(pars, "I_scaling", shared.outbreak_detection.I_scaling);
+
   }
 
   // This is a reasonable default implementation in the no-internal

tests/testthat/test-ascertainment_parameters.R

@@ -0,0 +1,42 @@
+test_that("The outbreak detection parameters work as expected.", {
+  pars <- test_toy_inputs(time_test = 1)
+  n_particles <- 3
+  times <- 0:75
+
+  ## Make the probability of declaring an outbreak very high:
+  pars$N_scaling <- 1e-6
+  pars$strength_scaling <- 1e-6
+  pars$I_scaling <- 1e-6
+
+  sys <- dust2::dust_system_create(cows(), pars, n_particles = n_particles, dt = 0.25)
+  dust2::dust_system_set_state_initial(sys)
+  ## sys$packer_state$unpack(dust2::dust_system_state(sys))
+  ## We subset only the model states, and not additionally added exported variables:
+  end_of_core_states <- (pars$n_herds + pars$n_regions)*4
+  s <- dust2::dust_system_simulate(sys, times)
+  ## Extract number of outbreaks
+  baseline_outbreaks <- s[(end_of_core_states + pars$n_herds + 1):(end_of_core_states + pars$n_herds + pars$n_regions),,]
+  ## Sum over time
+  baseline_outbreaks <- apply(baseline_outbreaks, c(1,2), sum)
+
+  ## Now repeat with parameters that greatly reduce the probability of declaring an outbreak close to 0:
+  pars$N_scaling <- 1e6
+  pars$strength_scaling <- 1000
+  pars$I_scaling <- 1e6
+
+  sys <- dust2::dust_system_create(cows(), pars, n_particles = n_particles, dt = 0.25)
+  dust2::dust_system_set_state_initial(sys)
+  s <- dust2::dust_system_simulate(sys, times)
+
+  less_outbreaks <- s[(end_of_core_states + pars$n_herds + 1):(end_of_core_states + pars$n_herds + pars$n_regions),,]
+  ## Sum over time
+  less_outbreaks <- apply(less_outbreaks, c(1,2), sum)
+
+  for(i in 1:3){
+    for(j in 1:3){
+      expect_true(less_outbreaks[i,j] <= baseline_outbreaks[i,j],
+                  "More outbreaks declared despite lower probability.")
+    }
+  }
+
+})