Add option to move percentage of exposed/infected hosts to susceptibles (#175)

Similarly to lethal temperature and mortality. Removes certain percentage of exposed/infected host and makes them susceptible.

The scheduling is done in the same way as lethal temperature, so it is scheduled yearly on the first of certain month.

Refactors relevant common code from mortality.

Closes #173.

Author: petrasovaa

include/pops/config.hpp

@@ -54,6 +54,10 @@ class Config
     int lethal_temperature_month{0};
     bool weather{false};
     double reproductive_rate{0};
+    // survival rate
+    bool use_survival_rate{false};
+    int survival_rate_month{0};
+    int survival_rate_day{0};
     // SI/SEI
     std::string model_type;
     int latency_period_steps;
@@ -110,6 +114,9 @@ class Config
         if (use_lethal_temperature)
             lethal_schedule_ =
                 scheduler_.schedule_action_yearly(lethal_temperature_month, 1);
+        if (use_survival_rate)
+            survival_rate_schedule_ = scheduler_.schedule_action_yearly(
+                survival_rate_month, survival_rate_day);
         if (use_spreadrates)
             spread_rate_schedule_ = schedule_from_string(
                 scheduler_, spreadrate_frequency, spreadrate_frequency_n);
@@ -154,6 +161,17 @@ class Config
         return lethal_schedule_;
     }
 
+    const std::vector<bool>& survival_rate_schedule() const
+    {
+        if (!use_survival_rate)
+            throw std::logic_error(
+                "survival_rate_schedule() not available when use_survival_rate is false");
+        if (!schedules_created_)
+            throw std::logic_error(
+                "Schedules were not created before calling survival_rate_schedule()");
+        return survival_rate_schedule_;
+    }
+
     const std::vector<bool>& spread_rate_schedule() const
     {
         if (!use_spreadrates)
@@ -203,6 +221,17 @@ class Config
         return get_number_of_scheduled_actions(lethal_schedule_);
     }
 
+    unsigned num_survival_rate()
+    {
+        if (!use_survival_rate)
+            throw std::logic_error(
+                "num_survival_rate() not available when use_survival_rate is false");
+        if (!schedules_created_)
+            throw std::logic_error(
+                "Schedules were not created before calling num_survival_rate()");
+        return get_number_of_scheduled_actions(survival_rate_schedule_);
+    }
+
     unsigned rate_num_steps()
     {
         if (!use_spreadrates)
@@ -302,6 +331,7 @@ class Config
     std::vector<bool> output_schedule_;
     std::vector<bool> mortality_schedule_;
     std::vector<bool> lethal_schedule_;
+    std::vector<bool> survival_rate_schedule_;
     std::vector<bool> spread_rate_schedule_;
     std::vector<bool> quarantine_schedule_;
 };

include/pops/model.hpp

@@ -189,6 +189,7 @@ class Model
         std::vector<IntegerRaster>& mortality_tracker,
         IntegerRaster& died,
         const std::vector<FloatRaster>& temperatures,
+        const std::vector<FloatRaster>& survival_rates,
         const FloatRaster& weather_coefficient,
         Treatments<IntegerRaster, FloatRaster>& treatments,
         IntegerRaster& resistant,
@@ -212,6 +213,19 @@ class Model
                 config_.lethal_temperature,
                 suitable_cells);
         }
+        // removal of percentage of dispersers
+        if (config_.use_survival_rate && config_.survival_rate_schedule()[step]) {
+            int survival_step =
+                simulation_step_to_action_step(config_.survival_rate_schedule(), step);
+            simulation_.remove_percentage(
+                infected,
+                susceptible,
+                mortality_tracker,
+                exposed,
+                total_exposed,
+                survival_rates[survival_step],
+                suitable_cells);
+        }
         // actual spread
         if (config_.spread_schedule()[step]) {
             simulation_.generate(

include/pops/simulation.hpp

@@ -55,6 +55,31 @@ std::vector<int> draw_n_from_v(std::vector<int> v, unsigned n, Generator& genera
     return v;
 }
 
+/** Draws n elements from a cohort of rasters. Expects n to be equal or less than
+ *  sum of cohorts at cell (i, j).
+ */
+template<typename Generator, typename IntegerRaster, typename RasterIndex = int>
+std::vector<int> draw_n_from_cohorts(
+    std::vector<IntegerRaster>& cohorts,
+    int n,
+    RasterIndex row,
+    RasterIndex col,
+    Generator& generator)
+{
+    std::vector<int> categories;
+    unsigned index = 0;
+    for (auto& raster : cohorts) {
+        categories.insert(categories.end(), raster(row, col), index);
+        index += 1;
+    }
+    std::vector<int> draw = draw_n_from_v(categories, n, generator);
+    std::vector<int> cohort_counts;
+    for (index = 0; index < cohorts.size(); index++) {
+        cohort_counts.push_back(std::count(draw.begin(), draw.end(), index));
+    }
+    return cohort_counts;
+}
+
 /** The type of a epidemiological model (SI or SEI)
  */
 enum class ModelType
@@ -207,6 +232,67 @@ class Simulation
         }
     }
 
+    /** Removes percentage of exposed and infected
+     *
+     * @param infected Currently infected hosts
+     * @param susceptible Currently susceptible hosts
+     * @param mortality_tracker_vector Hosts that are infected at a specific time step
+     * @param exposed Exposed hosts per cohort
+     * @param total_exposed Total exposed in all exposed cohorts
+     * @param survival_rate Raster between 0 and 1 representing pest survival rate
+     * @param suitable_cells used to run model only where host are known to occur
+     */
+    void remove_percentage(
+        IntegerRaster& infected,
+        IntegerRaster& susceptible,
+        std::vector<IntegerRaster>& mortality_tracker_vector,
+        std::vector<IntegerRaster>& exposed,
+        IntegerRaster& total_exposed,
+        const FloatRaster& survival_rate,
+        const std::vector<std::vector<int>>& suitable_cells)
+    {
+        for (auto indices : suitable_cells) {
+            int i = indices[0];
+            int j = indices[1];
+            if (survival_rate(i, j) < 1) {
+                int removed = 0;
+                // remove percentage of infestation/infection in the infected class
+                int removed_infected =
+                    infected(i, j) - std::lround(infected(i, j) * survival_rate(i, j));
+                infected(i, j) -= removed_infected;
+                removed += removed_infected;
+                // remove the removed infected from mortality cohorts
+                if (removed_infected > 0) {
+                    std::vector<int> mortality_draw = draw_n_from_cohorts(
+                        mortality_tracker_vector, removed_infected, i, j, generator_);
+                    int index = 0;
+                    for (auto& raster : mortality_tracker_vector) {
+                        raster(i, j) -= mortality_draw[index];
+                        index += 1;
+                    }
+                }
+                // remove the same percentage for total exposed and remove randomly from
+                // each cohort
+                int total_removed_exposed =
+                    total_exposed(i, j)
+                    - std::lround(total_exposed(i, j) * survival_rate(i, j));
+                total_exposed(i, j) -= total_removed_exposed;
+                removed += total_removed_exposed;
+                if (total_removed_exposed > 0) {
+                    std::vector<int> exposed_draw = draw_n_from_cohorts(
+                        exposed, total_removed_exposed, i, j, generator_);
+                    int index = 0;
+                    for (auto& raster : exposed) {
+                        raster(i, j) -= exposed_draw[index];
+                        index += 1;
+                    }
+                }
+                // move infested/infected host back to susceptible pool
+                susceptible(i, j) += removed;
+            }
+        }
+    }
+
     /** kills infected hosts based on mortality rate and timing. In the last year
      * of mortality tracking the first index all remaining tracked infected hosts
      * are removed. In indexes that are in the mortality_time_lag no mortality occurs.
@@ -347,43 +433,26 @@ class Simulation
             resistant_moved = std::count(draw.begin(), draw.end(), 4);
 
             if (exposed_moved > 0) {
+                std::vector<int> exposed_draw = draw_n_from_cohorts(
+                    exposed, exposed_moved, row_from, col_from, generator_);
                 int index = 0;
                 for (auto& raster : exposed) {
-                    auto exposed_count = raster(row_from, col_from);
-                    exposed_categories.insert(
-                        exposed_categories.end(), exposed_count, index);
-
-                    index += 1;
-                }
-                std::vector<int> exposed_draw =
-                    draw_n_from_v(exposed_categories, exposed_moved, generator_);
-                index = 0;
-                for (auto& raster : exposed) {
-                    auto exposed_moved_in_cohort =
-                        std::count(exposed_draw.begin(), exposed_draw.end(), index);
-                    raster(row_from, col_from) -= exposed_moved_in_cohort;
-                    raster(row_to, col_to) += exposed_moved_in_cohort;
+                    raster(row_from, col_from) -= exposed_draw[index];
+                    raster(row_to, col_to) += exposed_draw[index];
                     index += 1;
                 }
             }
-
             if (infected_moved > 0) {
-                std::vector<int> mortality_categories;
+                std::vector<int> mortality_draw = draw_n_from_cohorts(
+                    mortality_tracker_vector,
+                    infected_moved,
+                    row_from,
+                    col_from,
+                    generator_);
                 int index = 0;
                 for (auto& raster : mortality_tracker_vector) {
-                    auto mortality_count = raster(row_from, col_from);
-                    mortality_categories.insert(
-                        mortality_categories.end(), mortality_count, index);
-                    index += 1;
-                }
-                std::vector<int> mortality_draw =
-                    draw_n_from_v(mortality_categories, infected_moved, generator_);
-                index = 0;
-                for (auto& raster : mortality_tracker_vector) {
-                    auto mortality_moved_in_cohort =
-                        std::count(mortality_draw.begin(), mortality_draw.end(), index);
-                    raster(row_from, col_from) -= mortality_moved_in_cohort;
-                    raster(row_to, col_to) += mortality_moved_in_cohort;
+                    raster(row_from, col_from) -= mortality_draw[index];
+                    raster(row_to, col_to) += mortality_draw[index];
                     index += 1;
                 }
             }

tests/CMakeLists.txt

@@ -23,6 +23,7 @@ add_pops_test(test_deterministic)
 add_pops_test(test_model)
 add_pops_test(test_mortality)
 add_pops_test(test_movements)
+add_pops_test(test_survival_rate)
 add_pops_test(test_raster)
 add_pops_test(test_scheduling)
 add_pops_test(test_simulation)

tests/test_model.cpp

@@ -59,6 +59,7 @@ int test_with_reduced_stochasticity()
     config.model_type = "SI";
     config.latency_period_steps = 0;
     config.use_lethal_temperature = false;
+    config.use_survival_rate = false;
     config.use_quarantine = true;
     config.quarantine_frequency = "year";
     config.quarantine_frequency_n = 1;
@@ -125,6 +126,7 @@ int test_with_reduced_stochasticity()
         mortality_tracker,
         died,
         empty_float,
+        empty_float,
         empty_float[0],
         treatments,
         zeros,
@@ -212,6 +214,7 @@ int test_deterministic()
     config.model_type = "SI";
     config.latency_period_steps = 0;
     config.use_lethal_temperature = false;
+    config.use_survival_rate = false;
     config.use_quarantine = false;
     config.use_spreadrates = true;
     config.spreadrate_frequency = "year";
@@ -265,6 +268,7 @@ int test_deterministic()
         mortality_tracker,
         died,
         empty_float,
+        empty_float,
         empty_float[0],
         treatments,
         zeros,
@@ -351,6 +355,7 @@ int test_deterministic_exponential()
     config.model_type = "SI";
     config.latency_period_steps = 0;
     config.use_lethal_temperature = false;
+    config.use_survival_rate = false;
     config.use_quarantine = false;
     config.use_spreadrates = true;
     config.spreadrate_frequency = "year";
@@ -405,6 +410,7 @@ int test_deterministic_exponential()
         mortality_tracker,
         died,
         empty_float,
+        empty_float,
         empty_float[0],
         treatments,
         zeros,
@@ -487,6 +493,7 @@ int test_model_sei_deterministic()
     config.model_type = "SEI";
     config.latency_period_steps = 11;
     config.use_lethal_temperature = false;
+    config.use_survival_rate = false;
     config.use_quarantine = false;
     config.use_spreadrates = true;
     config.spreadrate_frequency = "year";
@@ -550,6 +557,7 @@ int test_model_sei_deterministic()
             mortality_tracker,
             died,
             empty_float,
+            empty_float,
             empty_float[0],
             treatments,
             zeros,
@@ -627,6 +635,7 @@ int test_model_sei_deterministic_with_treatments()
     config.model_type = "SEI";
     config.latency_period_steps = 11;
     config.use_lethal_temperature = false;
+    config.use_survival_rate = false;
     config.use_quarantine = false;
     config.use_spreadrates = true;
     config.spreadrate_frequency = "year";
@@ -708,6 +717,7 @@ int test_model_sei_deterministic_with_treatments()
             mortality_tracker,
             died,
             empty_float,
+            empty_float,
             empty_float[0],
             treatments,
             zeros,

tests/test_network_kernel.cpp

@@ -128,6 +128,7 @@ int test_model_with_network()
             empty_ints,
             zeros,
             empty_floats,
+            empty_floats,
             empty_floats[0],
             treatments,
             zeros,

tests/test_overpopulation_movements.cpp

@@ -148,6 +148,7 @@ int test_model()
         empty_ints,
         zeros,
         empty_floats,
+        empty_floats,
         empty_floats[0],
         treatments,
         zeros,

tests/test_simulation_kernels.cpp

@@ -248,6 +248,7 @@ int test_simulation_with_kernels_generic(
     config.model_type = "SI";
     config.latency_period_steps = 3;
     config.use_lethal_temperature = false;
+    config.use_survival_rate = false;
     config.use_quarantine = false;
     config.use_spreadrates = true;
     config.spreadrate_frequency = "year";
@@ -355,6 +356,7 @@ int test_model_with_kernels_generic(
     config.model_type = "SI";
     config.latency_period_steps = 3;
     config.use_lethal_temperature = false;
+    config.use_survival_rate = false;
     config.use_quarantine = false;
     config.use_spreadrates = true;
     config.spreadrate_frequency = "year";
@@ -432,6 +434,7 @@ int test_model_with_kernels_generic(
             mortality_tracker,
             died,
             empty_float,
+            empty_float,
             empty_float[0],
             treatments,
             zeros,