Use rounding for treatments

include/pops/host_pool.hpp

@@ -87,7 +87,8 @@ class HostPool : public HostPoolInterface<RasterIndex>
      *
      * 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 empty results (but will not fail).
+     * 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.
@@ -561,6 +562,7 @@ class HostPool : public HostPoolInterface<RasterIndex>
      * @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
@@ -569,9 +571,6 @@ class HostPool : public HostPoolInterface<RasterIndex>
      * @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.
      */
@@ -602,49 +601,47 @@ class HostPool : public HostPoolInterface<RasterIndex>
         // Possibly reuse in the I->S removal.
         if (infected <= 0)
             return;
-        if (use_mortality_) {
-            if (mortality_tracker_vector_.size() != mortality.size()) {
-                throw std::invalid_argument(
-                    "mortality is not the same size as the internal mortality tracker ("
-                    + std::to_string(mortality_tracker_vector_.size())
-                    + " != " + 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) {
-                if (mortality_tracker_vector_[i](row, col) < mortality[i]) {
-                    throw std::invalid_argument(
-                        "Mortality value [" + std::to_string(i) + "] is too high ("
-                        + std::to_string(mortality[i]) + " > "
-                        + std::to_string(mortality_tracker_vector_[i](row, col))
-                        + ") for cell (" + std::to_string(row) + ", "
-                        + std::to_string(col) + ")");
-                }
-                mortality_tracker_vector_[i](row, col) =
-                    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) {
-                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) + ")");
-            }
-            if (false && infected_(row, col) < mortality_total) {
+        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 ("
+                + std::to_string(mortality_tracker_vector_.size())
+                + " != " + std::to_string(mortality.size()) + ") for cell ("
+                + std::to_string(row) + ", " + std::to_string(col) + ")");
+        }
+        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(
-                    "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 (" + std::to_string(mortality_total) + " > "
-                    + std::to_string(infected) + ")");
+                    "Mortality value [" + std::to_string(i) + "] is too high ("
+                    + std::to_string(mortality[i]) + " > "
+                    + std::to_string(mortality_tracker_vector_[i](row, col))
+                    + ") for cell (" + std::to_string(row) + ", " + std::to_string(col)
+                    + ")");
             }
+            mortality_tracker_vector_[i](row, col) =
+                mortality_tracker_vector_[i](row, col) - mortality[i];
+            mortality_total += mortality[i];
+        }
+        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)
+                + ")");
+        }
+        if (infected_(row, col) < mortality_total) {
+            throw std::invalid_argument(
+                "Total of removed mortality values is higher than current number "
+                "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;
         reset_total_host(row, col);
@@ -750,6 +747,8 @@ class HostPool : public HostPoolInterface<RasterIndex>
     /**
      * @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
@@ -796,35 +795,38 @@ class HostPool : public HostPoolInterface<RasterIndex>
             total_resistant += exposed[i];
         }
         infected_(row, col) -= infected;
-        if (use_mortality_) {
-            if (mortality_tracker_vector_.size() != mortality.size()) {
-                throw std::invalid_argument(
-                    "mortality is not the same size as the internal mortality tracker ("
-                    + std::to_string(mortality_tracker_vector_.size())
-                    + " != " + std::to_string(mortality.size()) + ") for cell ("
-                    + std::to_string(row) + ", " + std::to_string(col) + ")");
-            }
-            int mortality_total = 0;
-            // no simple zip in C++, falling back to indices
-            for (size_t i = 0; i < mortality.size(); ++i) {
-                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) {
-                throw std::invalid_argument(
-                    "Total of mortality values differs from formerly infected, now resistant "
-                    "count ("
-                    + std::to_string(mortality_total)
-                    + " != " + std::to_string(infected) + " for cell ("
-                    + std::to_string(row) + ", " + std::to_string(col) + "))");
-            }
-        }
         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 ("
+                + std::to_string(mortality_tracker_vector_.size())
+                + " != " + std::to_string(mortality.size()) + ") for cell ("
+                + std::to_string(row) + ", " + std::to_string(col) + ")");
+        }
+        int mortality_total = 0;
+        // no simple zip in C++, falling back to indices
+        for (size_t i = 0; i < mortality_tracker_vector_.size(); ++i) {
+            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) {
+            throw std::invalid_argument(
+                "Total of mortality values differs from formerly infected, now resistant "
+                "count ("
+                + std::to_string(mortality_total) + " != " + std::to_string(infected)
+                + " for cell (" + std::to_string(row) + ", " + std::to_string(col)
+                + "))");
+        }
+        reset_total_host(row, col);
     }
 
     /**
@@ -1029,6 +1031,9 @@ class HostPool : public HostPoolInterface<RasterIndex>
     /**
      * @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
      *
@@ -1037,11 +1042,15 @@ class HostPool : public HostPoolInterface<RasterIndex>
     std::vector<int> mortality_by_group_at(RasterIndex row, RasterIndex col) const
     {
         std::vector<int> all;
-        if (use_mortality_) {
-            all.reserve(mortality_tracker_vector_.size());
-            for (const auto& raster : mortality_tracker_vector_)
-                all.push_back(raster(row, col));
+
+        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));
         return all;
     }
 

include/pops/treatments.hpp

@@ -125,15 +125,15 @@ class BaseTreatment : public AbstractTreatment<HostPool, FloatRaster>
     }
 
     // returning double allows identical results with the previous version
-    double get_treated(int i, int j, int count)
+    int get_treated(int i, int j, int count)
     {
         return get_treated(i, j, count, this->application_);
     }
 
-    double get_treated(int i, int j, int count, TreatmentApplication application)
+    int get_treated(int i, int j, int count, TreatmentApplication application)
     {
         if (application == TreatmentApplication::Ratio) {
-            return count * this->map_(i, j);
+            return std::lround(count * this->map_(i, j));
         }
         else if (application == TreatmentApplication::AllInfectedInCell) {
             return static_cast<bool>(this->map_(i, j)) ? count : 0;
@@ -173,20 +173,21 @@ class SimpleTreatment : public BaseTreatment<HostPool, FloatRaster>
         for (const auto& indices : host_pool.suitable_cells()) {
             int i = indices[0];
             int j = indices[1];
-            int remove_susceptible = static_cast<int>(std::ceil(this->get_treated(
-                i, j, host_pool.susceptible_at(i, j), TreatmentApplication::Ratio)));
-            int remove_infected = static_cast<int>(
-                std::ceil(this->get_treated(i, j, host_pool.infected_at(i, j))));
+            int remove_susceptible = this->get_treated(
+                i, j, host_pool.susceptible_at(i, j), TreatmentApplication::Ratio);
+            // Treated infected are computed as a sum of treated in mortality groups.
+            int remove_infected = 0;
             std::vector<int> remove_mortality;
             for (int count : host_pool.mortality_by_group_at(i, j)) {
-                remove_mortality.push_back(
-                    static_cast<int>(std::ceil(this->get_treated(i, j, count))));
+                int remove = this->get_treated(i, j, count);
+                remove_mortality.push_back(remove);
+                remove_infected += remove;
             }
+            // Will need to use infected directly if not mortality.
 
             std::vector<int> remove_exposed;
             for (int count : host_pool.exposed_by_group_at(i, j)) {
-                remove_exposed.push_back(
-                    static_cast<int>(std::ceil(this->get_treated(i, j, count))));
+                remove_exposed.push_back(this->get_treated(i, j, count));
             }
             host_pool.completely_remove_hosts_at(
                 i,
@@ -240,26 +241,19 @@ class PesticideTreatment : public BaseTreatment<HostPool, FloatRaster>
         for (const auto& indices : host_pool.suitable_cells()) {
             int i = indices[0];
             int j = indices[1];
-            // Given how the original code was written (everything was first converted
-            // to ints and subtractions happened only afterwards), this needs ints,
-            // not doubles to pass the r.pops.spread test (unlike the other code which
-            // did substractions before converting to ints), so the conversion to ints
-            // happened only later. Now get_treated returns double and floor or ceil is
-            // applied to the result to get the same results as before.
-            int susceptible_resistant = static_cast<int>(std::floor(this->get_treated(
-                i, j, host_pool.susceptible_at(i, j), TreatmentApplication::Ratio)));
+            int susceptible_resistant = this->get_treated(
+                i, j, host_pool.susceptible_at(i, j), TreatmentApplication::Ratio);
             std::vector<int> resistant_exposed_list;
             for (const auto& number : host_pool.exposed_by_group_at(i, j)) {
-                resistant_exposed_list.push_back(
-                    static_cast<int>(std::floor(this->get_treated(i, j, number))));
+                resistant_exposed_list.push_back(this->get_treated(i, j, number));
             }
+            int infected = 0;
             std::vector<int> resistant_mortality_list;
             for (const auto& number : host_pool.mortality_by_group_at(i, j)) {
-                resistant_mortality_list.push_back(
-                    static_cast<int>(std::floor(this->get_treated(i, j, number))));
+                int remove = this->get_treated(i, j, number);
+                resistant_mortality_list.push_back(remove);
+                infected += remove;
             }
-            int infected = static_cast<int>(
-                std::floor(this->get_treated(i, j, host_pool.infected_at(i, j))));
             host_pool.make_resistant_at(
                 i,
                 j,

tests/test_model.cpp

@@ -625,7 +625,7 @@ int test_model_sei_deterministic_with_treatments()
     config.set_date_end(2020, 12, 31);
     config.set_step_unit(StepUnit::Month);
     config.set_step_num_units(1);
-    config.use_mortality = false;
+    config.use_mortality = true;
     config.mortality_frequency = "year";
     config.mortality_frequency_n = 1;
     config.use_treatments = true;
@@ -667,6 +667,13 @@ int test_model_sei_deterministic_with_treatments()
         suitable_cells);
     std::vector<TestModel::StandardSingleHostPool*> host_pools = {&host_pool};
     TestModel::StandardMultiHostPool multi_host_pool(host_pools, config);
+    PestHostTable<TestModel::StandardSingleHostPool> pest_host_table(
+        model.environment());
+    pest_host_table.add_host_info(
+        config.establishment_probability,  // using as host susceptibility
+        config.mortality_rate,
+        config.mortality_time_lag);
+    multi_host_pool.set_pest_host_table(pest_host_table);
     TestModel::StandardPestPool pest_pool{
         dispersers, established_dispersers, outside_dispersers};
     SpreadRateAction<TestModel::StandardMultiHostPool, int> spread_rate(
@@ -698,10 +705,12 @@ int test_model_sei_deterministic_with_treatments()
     for (int row = 0; row < expected_infected.rows(); ++row)
         for (int col = 0; col < expected_infected.rows(); ++col)
             if (pesticide_treatment(row, col) > 0)
-                expected_infected(row, col) = static_cast<int>(
-                    std::floor(2 * pesticide_treatment(row, col) * infected(row, col)));
-    // Valus is based on the result which is considered correct.
-    Raster<int> expected_dispersers = {{0, 0, 0}, {0, 5, 0}, {0, 0, 2}};
+                expected_infected(row, col) = std::lround(
+                    2 * pesticide_treatment(row, col) * expected_infected(row, col));
+    expected_infected(0, 0) += 5;  // based on what is considered a correct result
+    expected_infected(1, 1) -= 5;  // based on what is considered a correct result
+    // Values are based on the result which is considered correct.
+    Raster<int> expected_dispersers = {{5, 0, 0}, {0, 10, 0}, {0, 0, 2}};
 
     for (unsigned int step = 0; step < config.scheduler().get_num_steps(); ++step) {
         model.run_step(