Use rounding for treatments

[wenzeslaus]

Replaced flooring by rounding in all treatments.

• PR in rpops with this change Use rounding for treatments rpops#203
• PR in r.pops.spread with this change Use rounding for treatments r.pops.spread#79

[wenzeslaus]

This causes some, but not much changes in r.pops.spread test results:

======================================================================
FAIL: test_outputs_mortality_pesticide_treatment (__main__.TestSpread.test_outputs_mortality_pesticide_treatment)
Check mortality together with pesticide treatment (all_infected_in_cell)
----------------------------------------------------------------------
...
AssertionError: r.univar map=average percentile=90.0 nprocs=1 separator== format=plain -g difference:
mismatch values (key, reference, actual): [('max', 5.6, 5.8), ('mean', 0.526, 0.552422813036034)]
command: r.univar map=average percentile=90.0 nprocs=1 separator== format=plain -g {'map': 'average', 'separator': '=', 'flags': 'g'}

======================================================================
FAIL: test_outputs_mortality_treatment (__main__.TestSpread.test_outputs_mortality_treatment)
Check mortality together with treatment
----------------------------------------------------------------------
...
AssertionError: r.univar map=dead_2020_12_31 percentile=90.0 nprocs=1 separator== format=plain -g difference:
mismatch values (key, reference, actual): [('sum', 2483, 2505)]
command: r.univar map=dead_2020_12_31 percentile=90.0 nprocs=1 separator== format=plain -g {'map': 'dead_2020_12_31', 'separator': '=', 'flags': 'g'}

----------------------------------------------------------------------

[wenzeslaus]

Two tests with mortality in r.pops.spread fail:

def test_outputs_mortality_treatment(self):
    """Check mortality together with treatment"""
# ...
def test_outputs_mortality_pesticide_treatment(self):
    """Check mortality together with pesticide treatment (all_infected_in_cell)"""
# ...

The other two treatment tests without mortality pass even with this change:

def test_sei_treatments_removal_ratio_to_all(self):
    """Check outputs with SEI and treatment_length == 0"""
# ...
def test_sei_treatments_pesticide_ratio_to_all(self):
    """Check outputs with SEI and treatment_length != 0 (pesticide)"""
# ...

[wenzeslaus]

I merged the explicit use_mortality code here. Unfortunately, I didn't notice I took different approach in making the conditions, so the diff in the PR is now larger than it was before, i.e., there is a lot of lines which are different because of indent, but the code is the same otherwise. This is a host_pool.hpp diff between the main branch before the explicit use_mortality branch and this PR. While it shows the explicit use_mortality changes from #231, it shows the changes for rounding more clearly.

diff --git a/include/pops/host_pool.hpp b/include/pops/host_pool.hpp
index 6fc3b79..f22cebe 100644
--- a/include/pops/host_pool.hpp
+++ b/include/pops/host_pool.hpp
@@ -85,10 +85,16 @@ public:
      * host infection over time. Expectation is that mortality tracker is of
      * length (1/mortality_rate + mortality_time_lag).
      *
+     * If *use_mortality* is set to false, the *mortality_tracker_vector* is not
+     * used or otherwise accessed, and any attempts to retrieve mortality information
+     * will provide infected hosts as if there is one cohort or empty results
+     * (but will not fail).
+     *
      * Host is added to the environment by the constructor. Afterwards, the environment
      * is not modified.
      *
      * @param model_type Type of the model (SI or SEI)
+     * @param use_mortality true if morality should be used
      * @param susceptible Raster of susceptible hosts
      * @param exposed Raster of exposed or infected hosts
      * @param latency_period Length of the latency period in steps
@@ -109,6 +115,7 @@ public:
      */
     HostPool(
         ModelType model_type,
+        bool use_mortality,
         IntegerRaster& susceptible,
         std::vector<IntegerRaster>& exposed,
         unsigned latency_period,
@@ -137,6 +144,7 @@ public:
           total_hosts_(total_hosts),
           environment_(environment),
           model_type_(model_type),
+          use_mortality_(use_mortality),
           dispersers_stochasticity_(dispersers_stochasticity),
           reproductive_rate_(reproductive_rate),
           establishment_stochasticity_(establishment_stochasticity),
@@ -148,6 +156,42 @@ public:
         environment.add_host(this);
     }
 
+    /** Constructor overload which takes Config instead of the individual values. */
+    HostPool(
+        const Config& config,
+        IntegerRaster& susceptible,
+        std::vector<IntegerRaster>& exposed,
+        IntegerRaster& infected,
+        IntegerRaster& total_exposed,
+        IntegerRaster& resistant,
+        std::vector<IntegerRaster>& mortality_tracker_vector,
+        IntegerRaster& died,
+        IntegerRaster& total_hosts,
+        Environment& environment,
+        std::vector<std::vector<int>>& suitable_cells)
+        : susceptible_(susceptible),
+          infected_(infected),
+          exposed_(exposed),
+          latency_period_(config.latency_period_steps),
+          total_exposed_(total_exposed),
+          resistant_(resistant),
+          mortality_tracker_vector_(mortality_tracker_vector),
+          died_(died),
+          total_hosts_(total_hosts),
+          environment_(environment),
+          model_type_(config.model_type_as_enum()),
+          use_mortality_(config.use_mortality),
+          dispersers_stochasticity_(config.generate_stochasticity),
+          reproductive_rate_(config.reproductive_rate),
+          establishment_stochasticity_(config.establishment_stochasticity),
+          deterministic_establishment_probability_(config.establishment_probability),
+          rows_(config.rows),
+          cols_(config.cols),
+          suitable_cells_(suitable_cells)
+    {
+        environment.add_host(this);
+    }
+
     /**
      * @brief Set pest-host table
      *
@@ -266,7 +310,8 @@ public:
         susceptible_(row, col) -= 1;
         if (model_type_ == ModelType::SusceptibleInfected) {
             infected_(row, col) += 1;
-            mortality_tracker_vector_.back()(row, col) += 1;
+            if (use_mortality_)
+                mortality_tracker_vector_.back()(row, col) += 1;
         }
         else if (model_type_ == ModelType::SusceptibleExposedInfected) {
             exposed_.back()(row, col) += 1;
@@ -472,7 +517,7 @@ public:
                 index += 1;
             }
         }
-        if (infected_moved > 0) {
+        if (use_mortality_ && infected_moved > 0) {
             std::vector<int> mortality_draw = draw_n_from_cohorts(
                 mortality_tracker_vector_,
                 infected_moved,
@@ -517,6 +562,7 @@ public:
      * @brief Completely remove any hosts
      *
      * Removes hosts completely (as opposed to moving them to another pool).
+     * If mortality is not active, the *mortality* parameter is ignored.
      *
      * @param row Row index of the cell
      * @param col Column index of the cell
@@ -525,9 +571,6 @@ public:
      * @param infected Number of infected hosts to remove.
      * @param mortality Number of infected hosts in each mortality cohort.
      *
-     * @note Counts are doubles, so that handling of floating point values is managed
-     * here in the same way as in the original treatment code.
-     *
      * @note This does not remove resistant just like the original implementation in
      * treatments.
      */
@@ -558,6 +601,11 @@ public:
         // Possibly reuse in the I->S removal.
         if (infected <= 0)
             return;
+        if (!use_mortality_) {
+            infected_(row, col) -= infected;
+            reset_total_host(row, col);
+            return;
+        }
         if (mortality_tracker_vector_.size() != mortality.size()) {
             throw std::invalid_argument(
                 "mortality is not the same size as the internal mortality tracker ("
@@ -565,9 +613,8 @@ public:
                 + " != " + std::to_string(mortality.size()) + ") for cell ("
                 + std::to_string(row) + ", " + std::to_string(col) + ")");
         }
-
-        double mortality_total = 0;
-        for (size_t i = 0; i < mortality.size(); ++i) {
+        int mortality_total = 0;
+        for (size_t i = 0; i < mortality_tracker_vector_.size(); ++i) {
             if (mortality_tracker_vector_[i](row, col) < mortality[i]) {
                 throw std::invalid_argument(
                     "Mortality value [" + std::to_string(i) + "] is too high ("
@@ -580,24 +627,20 @@ public:
                 mortality_tracker_vector_[i](row, col) - mortality[i];
             mortality_total += mortality[i];
         }
-        // These two values will only match if we actually compute one from another
-        // and once we don't need to keep the exact same double to int results for
-        // tests. First condition always fails the tests. The second one may potentially
-        // fail.
-        if (false && infected != mortality_total) {
+        if (infected != mortality_total) {
             throw std::invalid_argument(
                 "Total of removed mortality values differs from removed infected "
                 "count ("
                 + std::to_string(mortality_total) + " != " + std::to_string(infected)
-                + " for cell (" + std::to_string(row) + ", " + std::to_string(col)
+                + ") for cell (" + std::to_string(row) + ", " + std::to_string(col)
                 + ")");
         }
-        if (false && infected_(row, col) < mortality_total) {
+        if (infected_(row, col) < mortality_total) {
             throw std::invalid_argument(
                 "Total of removed mortality values is higher than current number "
-                "of infected hosts for cell ("
-                + std::to_string(row) + ", " + std::to_string(col) + ") is too high ("
+                "of infected hosts ("
                 + std::to_string(mortality_total) + " > " + std::to_string(infected)
+                + ") for cell (" + std::to_string(row) + ", " + std::to_string(col)
                 + ")");
         }
         infected_(row, col) -= infected;
@@ -623,13 +666,15 @@ public:
         // remove percentage of infestation/infection in the infected class
         infected_(row, col) -= count;
         // remove the removed infected from mortality cohorts
-        if (count > 0) {
-            std::vector<int> mortality_draw = draw_n_from_cohorts(
-                mortality_tracker_vector_, count, row, col, generator);
-            int index = 0;
-            for (auto& raster : mortality_tracker_vector_) {
-                raster(row, col) -= mortality_draw[index];
-                index += 1;
+        if (use_mortality_) {
+            if (count > 0) {
+                std::vector<int> mortality_draw = draw_n_from_cohorts(
+                    mortality_tracker_vector_, count, row, col, generator);
+                int index = 0;
+                for (auto& raster : mortality_tracker_vector_) {
+                    raster(row, col) -= mortality_draw[index];
+                    index += 1;
+                }
             }
         }
         // move infested/infected host back to susceptible pool
@@ -702,6 +747,8 @@ public:
     /**
      * @brief Make hosts resistant in a given cell
      *
+     * If mortality is not active, the *mortality* parameter is ignored.
+     *
      * @param row Row index of the cell
      * @param col Column index of the cell
      * @param susceptible Number of susceptible hosts to make resistant
@@ -748,6 +795,12 @@ public:
             total_resistant += exposed[i];
         }
         infected_(row, col) -= infected;
+        total_resistant += infected;
+        resistant_(row, col) += total_resistant;
+        if (!use_mortality_) {
+            reset_total_host(row, col);
+            return;
+        }
         if (mortality_tracker_vector_.size() != mortality.size()) {
             throw std::invalid_argument(
                 "mortality is not the same size as the internal mortality tracker ("
@@ -757,7 +810,7 @@ public:
         }
         int mortality_total = 0;
         // no simple zip in C++, falling back to indices
-        for (size_t i = 0; i < mortality.size(); ++i) {
+        for (size_t i = 0; i < mortality_tracker_vector_.size(); ++i) {
             mortality_tracker_vector_[i](row, col) -= mortality[i];
             mortality_total += mortality[i];
         }
@@ -773,8 +826,7 @@ public:
                 + " for cell (" + std::to_string(row) + ", " + std::to_string(col)
                 + "))");
         }
-        total_resistant += infected;
-        resistant_(row, col) += total_resistant;
+        reset_total_host(row, col);
     }
 
     /**
@@ -814,7 +866,7 @@ public:
     void apply_mortality_at(
         RasterIndex row, RasterIndex col, double mortality_rate, int mortality_time_lag)
     {
-        if (mortality_rate <= 0)
+        if (mortality_rate <= 0 || !use_mortality_)
             return;
         int max_index = mortality_tracker_vector_.size() - mortality_time_lag - 1;
         for (int index = 0; index <= max_index; index++) {
@@ -890,6 +942,8 @@ public:
      */
     void step_forward_mortality()
     {
+        if (!use_mortality_)
+            return;
         rotate_left_by_one(mortality_tracker_vector_);
     }
 
@@ -977,6 +1031,9 @@ public:
     /**
      * @brief Get infected hosts in each mortality cohort at a given cell
      *
+     * If mortality is not active, it returns number of all infected individuals
+     * as if there would be only one cohort.
+     *
      * @param row Row index of the cell
      * @param col Column index of the cell
      *
@@ -985,6 +1042,12 @@ public:
     std::vector<int> mortality_by_group_at(RasterIndex row, RasterIndex col) const
     {
         std::vector<int> all;
+
+        if (!use_mortality_) {
+            all.push_back(infected_at(row, col));
+            return all;
+        }
+
         all.reserve(mortality_tracker_vector_.size());
         for (const auto& raster : mortality_tracker_vector_)
             all.push_back(raster(row, col));
@@ -1077,7 +1140,8 @@ public:
                 auto& oldest = exposed_.front();
                 // Move hosts to infected raster
                 infected_ += oldest;
-                mortality_tracker_vector_.back() += oldest;
+                if (use_mortality_)
+                    mortality_tracker_vector_.back() += oldest;
                 total_exposed_ += (oldest * (-1));
                 // Reset the raster
                 // (hosts moved from the raster)
@@ -1160,6 +1224,7 @@ private:
     const Environment& environment_;
 
     ModelType model_type_;
+    bool use_mortality_{false};
 
     bool dispersers_stochasticity_{false};
     double reproductive_rate_{0};

[ChrisJones687]

This all works as intended in rpops. The rounding changes results as expected based on tests here.