Use rounding for non-stochastic dispersers from cell

inst/include/actions.hpp

@@ -90,16 +90,14 @@ class SpreadAction
                     // From all the generated dispersers, some go to the soil in the
                     // same cell and don't participate in the kernel-driven dispersal.
                     auto dispersers_to_soil =
-                        std::round(to_soil_percentage_ * dispersers_from_cell);
+                        std::lround(to_soil_percentage_ * dispersers_from_cell);
                     soil_pool_->dispersers_to(dispersers_to_soil, i, j, generator);
                     dispersers_from_cell -= dispersers_to_soil;
                 }
-                pests.set_dispersers_at(i, j, dispersers_from_cell);
-                pests.set_established_dispersers_at(i, j, dispersers_from_cell);
+                pests.set_dispersers_at(i, j, dispersers_from_cell, 0);
             }
             else {
-                pests.set_dispersers_at(i, j, 0);
-                pests.set_established_dispersers_at(i, j, 0);
+                pests.set_dispersers_at(i, j, 0, 0);
             }
         }
     }
@@ -123,17 +121,15 @@ class SpreadAction
             if (pests.dispersers_at(i, j) > 0) {
                 for (int k = 0; k < pests.dispersers_at(i, j); k++) {
                     std::tie(row, col) = dispersal_kernel_(generator, i, j);
-                    // if (row < 0 || row >= rows_ || col < 0 || col >= cols_) {
                     if (host_pool.is_outside(row, col)) {
                         pests.add_outside_disperser_at(row, col);
-                        pests.remove_established_dispersers_at(i, j, 1);
                         continue;
                     }
                     // Put a disperser to the host pool.
                     auto dispersed =
                         host_pool.disperser_to(row, col, generator.establishment());
-                    if (!dispersed) {
-                        pests.remove_established_dispersers_at(i, j, 1);
+                    if (dispersed) {
+                        pests.add_established_dispersers_at(i, j, 1);
                     }
                 }
             }
@@ -370,7 +366,8 @@ class MoveOverpopulatedPests
                 // for leaving_percentage == 0.5
                 // 2 infected -> 1 leaving
                 // 3 infected -> 1 leaving
-                int leaving = original_count * leaving_percentage_;
+                int leaving =
+                    static_cast<int>(std::floor(original_count * leaving_percentage_));
                 leaving = hosts.pests_from(i, j, leaving, generator.overpopulation());
                 if (row < 0 || row >= rows_ || col < 0 || col >= cols_) {
                     pests.add_outside_dispersers_at(row, col, leaving);
@@ -510,7 +507,7 @@ class Mortality
     void action(Hosts& hosts)
     {
         for (auto indices : hosts.suitable_cells()) {
-            if (action_mortality_) {
+            if (static_cast<bool>(action_mortality_)) {
                 hosts.apply_mortality_at(
                     indices[0], indices[1], mortality_rate_, mortality_time_lag_);
             }

inst/include/config.hpp

@@ -569,12 +569,19 @@ class Config
         for (const auto& row : values) {
             if (row.size() < 3) {
                 throw std::invalid_argument(
-                    "3 values are required for each pest-host table row");
+                    "3 values are required for each pest-host table row "
+                    "(but row size is "
+                    + std::to_string(row.size()) + ")");
             }
             PestHostTableDataRow resulting_row;
             resulting_row.susceptibility = row[0];
             resulting_row.mortality_rate = row[1];
             resulting_row.mortality_time_lag = row[2];
+            if (resulting_row.susceptibility < 0 || resulting_row.susceptibility > 1) {
+                throw std::invalid_argument(
+                    "Susceptibility needs to be >=0 and <=1, not "
+                    + std::to_string(resulting_row.susceptibility));
+            }
             pest_host_table_data_.push_back(std::move(resulting_row));
         }
     }

inst/include/deterministic_kernel.hpp

@@ -172,8 +172,10 @@ class DeterministicDispersalKernel
             // The invalid state is checked later, in this case using the kernel type.
             return;
         }
-        number_of_columns = ceil(max_distance / east_west_resolution) * 2 + 1;
-        number_of_rows = ceil(max_distance / north_south_resolution) * 2 + 1;
+        number_of_columns =
+            static_cast<int>(ceil(max_distance / east_west_resolution)) * 2 + 1;
+        number_of_rows =
+            static_cast<int>(ceil(max_distance / north_south_resolution)) * 2 + 1;
         Raster<double> prob_size(number_of_rows, number_of_columns, 0);
         probability = prob_size;
         probability_copy = prob_size;

inst/include/host_pool.hpp

@@ -26,6 +26,7 @@
 #include "environment_interface.hpp"
 #include "competency_table.hpp"
 #include "pest_host_table.hpp"
+#include "utils.hpp"
 
 namespace pops {
 
@@ -306,7 +307,7 @@ class HostPool : public HostPoolInterface<RasterIndex>
             }
         }
         else {
-            dispersers_from_cell = lambda * infected_at(row, col);
+            dispersers_from_cell = std::lround(lambda * infected_at(row, col));
         }
         return dispersers_from_cell;
     }
@@ -326,7 +327,17 @@ class HostPool : public HostPoolInterface<RasterIndex>
         if (pest_host_table_) {
             suitability *= pest_host_table_->susceptibility(this);
         }
-        return environment_.influence_suitability_at(row, col, suitability);
+        suitability = environment_.influence_suitability_at(row, col, suitability);
+        if (suitability < 0 || suitability > 1) {
+            throw std::invalid_argument(
+                "Suitability should be >=0 and <=1, not " + std::to_string(suitability)
+                + " (susceptible: " + std::to_string(susceptible_(row, col))
+                + ", total population: "
+                + std::to_string(environment_.total_population_at(row, col))
+                + ", susceptibility: "
+                + std::to_string(pest_host_table_->susceptibility(this)) + ")");
+        }
+        return suitability;
     }
 
     /**
@@ -477,16 +488,9 @@ class HostPool : public HostPoolInterface<RasterIndex>
         // Since suitable cells originally comes from the total hosts, check first total
         // hosts and proceed only if there was no host.
         if (total_hosts_(row_to, col_to) == 0) {
-            for (auto indices : suitable_cells_) {
-                int i = indices[0];
-                int j = indices[1];
-                // TODO: This looks like a bug. Flag is needed for found and push back
-                // should happen only after the loop.
-                if ((i == row_to) && (j == col_to)) {
-                    std::vector<int> added_index = {row_to, col_to};
-                    suitable_cells_.push_back(added_index);
-                    break;
-                }
+            std::vector<int> new_index = {row_to, col_to};
+            if (!container_contains(suitable_cells_, new_index)) {
+                suitable_cells_.push_back(new_index);
             }
         }
 
@@ -528,10 +532,10 @@ class HostPool : public HostPoolInterface<RasterIndex>
     void completely_remove_hosts_at(
         RasterIndex row,
         RasterIndex col,
-        double susceptible,
-        std::vector<double> exposed,
-        double infected,
-        const std::vector<double>& mortality)
+        int susceptible,
+        std::vector<int> exposed,
+        int infected,
+        const std::vector<int>& mortality)
     {
         if (susceptible > 0)
             susceptible_(row, col) = susceptible_(row, col) - susceptible;
@@ -795,6 +799,9 @@ class HostPool : public HostPoolInterface<RasterIndex>
      * individuals is multiplied by the mortality rate to calculate the number of hosts
      * that die that time step.
      *
+     * If mortality rate is zero (<=0), no mortality is applied and mortality tracker
+     * vector stays as is, i.e., no hosts die.
+     *
      * To be used together with step_forward_mortality().
      *
      * @param row Row index of the cell
@@ -805,6 +812,8 @@ class HostPool : public HostPoolInterface<RasterIndex>
     void apply_mortality_at(
         RasterIndex row, RasterIndex col, double mortality_rate, int mortality_time_lag)
     {
+        if (mortality_rate <= 0)
+            return;
         int max_index = mortality_tracker_vector_.size() - mortality_time_lag - 1;
         for (int index = 0; index <= max_index; index++) {
             int mortality_in_index = 0;
@@ -815,8 +824,8 @@ class HostPool : public HostPoolInterface<RasterIndex>
                     mortality_in_index = mortality_tracker_vector_[index](row, col);
                 }
                 else {
-                    mortality_in_index =
-                        mortality_rate * mortality_tracker_vector_[index](row, col);
+                    mortality_in_index = static_cast<int>(std::floor(
+                        mortality_rate * mortality_tracker_vector_[index](row, col)));
                 }
                 mortality_tracker_vector_[index](row, col) -= mortality_in_index;
                 died_(row, col) += mortality_in_index;

inst/include/network.hpp

@@ -58,7 +58,7 @@ class EdgeGeometry : public std::vector<RasterCell>
     /** Get cost of the whole segment (edge). */
     double cost() const
     {
-        if (total_cost_)
+        if (static_cast<bool>(total_cost_))
             return total_cost_;
         // This is short for ((size - 2) + (2 * 1/2)) * cost per cell.
         return (this->size() - 1) * cost_per_cell_;
@@ -72,7 +72,7 @@ class EdgeGeometry : public std::vector<RasterCell>
     /** Get cost per cell for the segment (edge). */
     double cost_per_cell() const
     {
-        if (total_cost_)
+        if (static_cast<bool>(total_cost_))
             return total_cost_ / (this->size() - 1);
         return cost_per_cell_;
     }
@@ -879,14 +879,9 @@ class Network
             auto node_1_id = node_id_from_text(node_1_text);
             auto node_2_id = node_id_from_text(node_2_text);
             if (node_1_id < 1 || node_2_id < 1) {
-                throw std::runtime_error(std::string(
-                    "Node ID must be greater than zero: " + node_1_text + " "
-                    + node_2_text));
-            }
-            if (node_1_id == node_2_id) {
                 throw std::runtime_error(
-                    std::string("Edge cannot begin and end with the same node: ")
-                    + node_1_text + " " + node_2_text);
+                    std::string("Node ID must be greater than zero (node 1, node 2): ")
+                    + node_1_text + ", " + node_2_text + ", line: " + line);
             }
             Segment segment;
 

inst/include/pest_host_table.hpp

@@ -103,7 +103,7 @@ class PestHostTable
      * @param host Pointer to the host to get the information for
      * @return Mortality time lag value
      */
-    double mortality_time_lag(const HostPool* host) const
+    int mortality_time_lag(const HostPool* host) const
     {
         auto host_index = environment_.host_index(host);
         return mortality_time_lags_.at(host_index);

inst/include/pest_pool.hpp

@@ -53,13 +53,17 @@ class PestPool
     {}
     /**
      * @brief Set number of dispersers
+     *
      * @param row Row number
      * @param col Column number
-     * @param value The new value
+     * @param dispersers Number of dispersers
+     * @param established_dispersers Number of established dispersers
      */
-    void set_dispersers_at(RasterIndex row, RasterIndex col, int value)
+    void set_dispersers_at(
+        RasterIndex row, RasterIndex col, int dispersers, int established_dispersers)
     {
-        dispersers_(row, col) = value;
+        dispersers_(row, col) = dispersers;
+        established_dispersers_(row, col) = established_dispersers;
     }
     /**
      * @brief Return number of dispersers
@@ -82,32 +86,21 @@ class PestPool
     {
         return dispersers_;
     }
+
     /**
-     * @brief Set number of established dispersers
+     * @brief Add established dispersers
      *
-     * Established are dispersers which left cell (row, col) and established themselves
-     * elsewhere, i.e., origin of the established dispersers is tracked.
+     * Established dispersers are dispersers which left cell (row, col) and
+     * established themselves elsewhere, i.e., origin of the established dispersers
+     * is tracked.
      *
      * @param row Row number
      * @param col Column number
-     * @param value The new value
-     */
-    void set_established_dispersers_at(RasterIndex row, RasterIndex col, int value)
-    {
-        established_dispersers_(row, col) = value;
-    }
-    // TODO: The following function should not be necessary because pests can't
-    // un-establish. It exists just because it mirrors how the raster was handled in the
-    // original Simulation code.
-    /**
-     * @brief Remove established dispersers
-     * @param row Row number
-     * @param col Column number
-     * @param count How many dispers to remove
+     * @param count How many dispersers to add
      */
-    void remove_established_dispersers_at(RasterIndex row, RasterIndex col, int count)
+    void add_established_dispersers_at(RasterIndex row, RasterIndex col, int count)
     {
-        established_dispersers_(row, col) -= count;
+        established_dispersers_(row, col) += count;
     }
     /**
      * @brief Add a disperser which left the study area

inst/include/quarantine.hpp

@@ -147,20 +147,20 @@ class QuarantineEscapeAction
         DistDir closest;
         if (directions_.at(Direction::N)
             && (i - n) * north_south_resolution_ < mindist) {
-            mindist = (i - n) * north_south_resolution_;
+            mindist = static_cast<int>(std::floor((i - n) * north_south_resolution_));
             closest = std::make_tuple(mindist, Direction::N);
         }
         if (directions_.at(Direction::S)
             && (s - i) * north_south_resolution_ < mindist) {
-            mindist = (s - i) * north_south_resolution_;
+            mindist = static_cast<int>(std::floor((s - i) * north_south_resolution_));
             closest = std::make_tuple(mindist, Direction::S);
         }
         if (directions_.at(Direction::E) && (e - j) * west_east_resolution_ < mindist) {
-            mindist = (e - j) * west_east_resolution_;
+            mindist = static_cast<int>(std::floor((e - j) * west_east_resolution_));
             closest = std::make_tuple(mindist, Direction::E);
         }
         if (directions_.at(Direction::W) && (j - w) * west_east_resolution_ < mindist) {
-            mindist = (j - w) * west_east_resolution_;
+            mindist = static_cast<int>(std::floor((j - w) * west_east_resolution_));
             closest = std::make_tuple(mindist, Direction::W);
         }
         return closest;

inst/include/radial_kernel.hpp

@@ -198,8 +198,8 @@ class RadialDispersalKernel
         }
         theta = von_mises(generator);
 
-        row -= round(distance * cos(theta) / north_south_resolution);
-        col += round(distance * sin(theta) / east_west_resolution);
+        row -= lround(distance * cos(theta) / north_south_resolution);
+        col += lround(distance * sin(theta) / east_west_resolution);
 
         return std::make_tuple(row, col);
     }