From ec471b15773b37631ed0aadda3eda85e152c3056 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Thu, 14 Nov 2024 09:41:47 -0500 Subject: [PATCH] Use rounding for non-stochastic dispersers from cell Update to PoPS Core with flooring instead of rounding in computation of number of non-stochastic dispersers from cell. --- inst/cpp/pops-core | 2 +- inst/include/actions.hpp | 19 +++---- inst/include/config.hpp | 9 +++- inst/include/deterministic_kernel.hpp | 6 ++- inst/include/host_pool.hpp | 45 +++++++++------- inst/include/network.hpp | 13 ++--- inst/include/pest_host_table.hpp | 2 +- inst/include/pest_pool.hpp | 37 ++++++------- inst/include/quarantine.hpp | 8 +-- inst/include/radial_kernel.hpp | 4 +- inst/include/raster.hpp | 77 ++++++++++++++++++++++++--- inst/include/simulation.hpp | 9 ++-- inst/include/soils.hpp | 2 +- inst/include/treatments.hpp | 38 +++++++------ 14 files changed, 174 insertions(+), 97 deletions(-) diff --git a/inst/cpp/pops-core b/inst/cpp/pops-core index 19eab166..19b77bfd 160000 --- a/inst/cpp/pops-core +++ b/inst/cpp/pops-core @@ -1 +1 @@ -Subproject commit 19eab166224bb5f75a285ec2d96920ec8e2d3f4c +Subproject commit 19b77bfd3e67a6e3e2cd19801848a59ada697394 diff --git a/inst/include/actions.hpp b/inst/include/actions.hpp index e92bbc54..1b05188e 100644 --- a/inst/include/actions.hpp +++ b/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(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(action_mortality_)) { hosts.apply_mortality_at( indices[0], indices[1], mortality_rate_, mortality_time_lag_); } diff --git a/inst/include/config.hpp b/inst/include/config.hpp index 79ea58b1..a5382a9e 100644 --- a/inst/include/config.hpp +++ b/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)); } } diff --git a/inst/include/deterministic_kernel.hpp b/inst/include/deterministic_kernel.hpp index a3df69eb..edcd02a9 100644 --- a/inst/include/deterministic_kernel.hpp +++ b/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(ceil(max_distance / east_west_resolution)) * 2 + 1; + number_of_rows = + static_cast(ceil(max_distance / north_south_resolution)) * 2 + 1; Raster prob_size(number_of_rows, number_of_columns, 0); probability = prob_size; probability_copy = prob_size; diff --git a/inst/include/host_pool.hpp b/inst/include/host_pool.hpp index 2c600b31..ec7221d5 100644 --- a/inst/include/host_pool.hpp +++ b/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 } } 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 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 // 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 added_index = {row_to, col_to}; - suitable_cells_.push_back(added_index); - break; - } + std::vector 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 void completely_remove_hosts_at( RasterIndex row, RasterIndex col, - double susceptible, - std::vector exposed, - double infected, - const std::vector& mortality) + int susceptible, + std::vector exposed, + int infected, + const std::vector& mortality) { if (susceptible > 0) susceptible_(row, col) = susceptible_(row, col) - susceptible; @@ -795,6 +799,9 @@ class HostPool : public HostPoolInterface * 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 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 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(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; diff --git a/inst/include/network.hpp b/inst/include/network.hpp index 56caf671..6349cccb 100644 --- a/inst/include/network.hpp +++ b/inst/include/network.hpp @@ -58,7 +58,7 @@ class EdgeGeometry : public std::vector /** Get cost of the whole segment (edge). */ double cost() const { - if (total_cost_) + if (static_cast(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 /** Get cost per cell for the segment (edge). */ double cost_per_cell() const { - if (total_cost_) + if (static_cast(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; diff --git a/inst/include/pest_host_table.hpp b/inst/include/pest_host_table.hpp index 0d70ba3e..a0017e43 100644 --- a/inst/include/pest_host_table.hpp +++ b/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); diff --git a/inst/include/pest_pool.hpp b/inst/include/pest_pool.hpp index 7eb3e6bb..e548cf3f 100644 --- a/inst/include/pest_pool.hpp +++ b/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 diff --git a/inst/include/quarantine.hpp b/inst/include/quarantine.hpp index 71da68b1..6720ab5c 100644 --- a/inst/include/quarantine.hpp +++ b/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(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(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(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(std::floor((j - w) * west_east_resolution_)); closest = std::make_tuple(mindist, Direction::W); } return closest; diff --git a/inst/include/radial_kernel.hpp b/inst/include/radial_kernel.hpp index b6e4722b..ce265cb9 100644 --- a/inst/include/radial_kernel.hpp +++ b/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); } diff --git a/inst/include/raster.hpp b/inst/include/raster.hpp index 03ff60d8..54ccdc2c 100644 --- a/inst/include/raster.hpp +++ b/inst/include/raster.hpp @@ -275,7 +275,11 @@ class Raster } template - Raster& operator+=(OtherNumber value) + typename std::enable_if< + !(std::is_floating_point::value + && std::is_integral::value), + Raster&>::type + operator+=(OtherNumber value) { std::for_each( data_, data_ + (cols_ * rows_), [&value](Number& a) { a += value; }); @@ -283,7 +287,11 @@ class Raster } template - Raster& operator-=(OtherNumber value) + typename std::enable_if< + !(std::is_floating_point::value + && std::is_integral::value), + Raster&>::type + operator-=(OtherNumber value) { std::for_each( data_, data_ + (cols_ * rows_), [&value](Number& a) { a -= value; }); @@ -291,7 +299,11 @@ class Raster } template - Raster& operator*=(OtherNumber value) + typename std::enable_if< + !(std::is_floating_point::value + && std::is_integral::value), + Raster&>::type + operator*=(OtherNumber value) { std::for_each( data_, data_ + (cols_ * rows_), [&value](Number& a) { a *= value; }); @@ -299,13 +311,65 @@ class Raster } template - Raster& operator/=(OtherNumber value) + typename std::enable_if< + !(std::is_floating_point::value + && std::is_integral::value), + Raster&>::type + operator/=(OtherNumber value) { std::for_each( data_, data_ + (cols_ * rows_), [&value](Number& a) { a /= value; }); return *this; } + template + typename std::enable_if< + std::is_floating_point::value && std::is_integral::value, + Raster&>::type + operator+=(OtherNumber value) + { + std::for_each(data_, data_ + (cols_ * rows_), [&value](Number& a) { + a += static_cast(std::floor(value)); + }); + return *this; + } + + template + typename std::enable_if< + std::is_floating_point::value && std::is_integral::value, + Raster&>::type + operator-=(OtherNumber value) + { + std::for_each(data_, data_ + (cols_ * rows_), [&value](Number& a) { + a -= static_cast(std::floor(value)); + }); + return *this; + } + + template + typename std::enable_if< + std::is_floating_point::value && std::is_integral::value, + Raster&>::type + operator*=(OtherNumber value) + { + std::for_each(data_, data_ + (cols_ * rows_), [&value](Number& a) { + a *= static_cast(std::floor(value)); + }); + return *this; + } + + template + typename std::enable_if< + std::is_floating_point::value && std::is_integral::value, + Raster&>::type + operator/=(OtherNumber value) + { + std::for_each(data_, data_ + (cols_ * rows_), [&value](Number& a) { + a /= static_cast(std::floor(value)); + }); + return *this; + } + template typename std::enable_if< std::is_floating_point::value @@ -496,12 +560,13 @@ class Raster return out; } - friend inline Raster pow(Raster image, double value) + friend inline Raster pow(const Raster& image, double value) { image.for_each([value](Number& a) { a = std::pow(a, value); }); return image; } - friend inline Raster sqrt(Raster image) + + friend inline Raster sqrt(const Raster& image) { image.for_each([](Number& a) { a = std::sqrt(a); }); return image; diff --git a/inst/include/simulation.hpp b/inst/include/simulation.hpp index 582185bb..e4d5db8d 100644 --- a/inst/include/simulation.hpp +++ b/inst/include/simulation.hpp @@ -35,10 +35,11 @@ namespace pops { /*! A class to control the spread simulation. * - * \deprecated - * The class is deprecated in favor of individual action classes and a higher-level - * Model. The class corresponding to the original Simulation class before too much code - * accumulated in Simulation is SpreadAction. The class is now used only in tests. + * @note + * The class is deprecated for external use in favor of individual action classes and a + * higher-level Model. The class corresponding to the original Simulation class before + * too much code accumulated in Simulation is SpreadAction. The class is now used only + * in tests. * * The Simulation class handles the mechanics of the model, but the * timing of the events or steps should be handled outside of this diff --git a/inst/include/soils.hpp b/inst/include/soils.hpp index 8e23435d..3c4802d3 100644 --- a/inst/include/soils.hpp +++ b/inst/include/soils.hpp @@ -90,7 +90,7 @@ class SoilPool } } else { - dispersers = lambda * count; + dispersers = static_cast(std::floor(lambda * count)); } auto draw = draw_n_from_cohorts(*rasters_, dispersers, row, col, generator); size_t index = 0; diff --git a/inst/include/treatments.hpp b/inst/include/treatments.hpp index 39381a65..d6f15df3 100644 --- a/inst/include/treatments.hpp +++ b/inst/include/treatments.hpp @@ -136,7 +136,7 @@ class BaseTreatment : public AbstractTreatment return count * this->map_(i, j); } else if (application == TreatmentApplication::AllInfectedInCell) { - return this->map_(i, j) ? count : 0; + return static_cast(this->map_(i, j)) ? count : 0; } throw std::runtime_error( "BaseTreatment::get_treated: unknown TreatmentApplication"); @@ -173,18 +173,20 @@ class SimpleTreatment : public BaseTreatment for (auto indices : host_pool.suitable_cells()) { int i = indices[0]; int j = indices[1]; - double remove_susceptible = this->get_treated( - i, j, host_pool.susceptible_at(i, j), TreatmentApplication::Ratio); - double remove_infected = - this->get_treated(i, j, host_pool.infected_at(i, j)); - std::vector remove_mortality; + int remove_susceptible = static_cast(std::ceil(this->get_treated( + i, j, host_pool.susceptible_at(i, j), TreatmentApplication::Ratio))); + int remove_infected = static_cast( + std::ceil(this->get_treated(i, j, host_pool.infected_at(i, j)))); + std::vector remove_mortality; for (int count : host_pool.mortality_by_group_at(i, j)) { - remove_mortality.push_back(this->get_treated(i, j, count)); + remove_mortality.push_back( + static_cast(std::ceil(this->get_treated(i, j, count)))); } - std::vector remove_exposed; + std::vector remove_exposed; for (int count : host_pool.exposed_by_group_at(i, j)) { - remove_exposed.push_back(this->get_treated(i, j, count)); + remove_exposed.push_back( + static_cast(std::ceil(this->get_treated(i, j, count)))); } host_pool.completely_remove_hosts_at( i, @@ -241,23 +243,29 @@ class PesticideTreatment : public BaseTreatment // 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). - int susceptible_resistant = this->get_treated( - i, j, host_pool.susceptible_at(i, j), TreatmentApplication::Ratio); + // 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(std::floor(this->get_treated( + i, j, host_pool.susceptible_at(i, j), TreatmentApplication::Ratio))); std::vector resistant_exposed_list; for (const auto& number : host_pool.exposed_by_group_at(i, j)) { - resistant_exposed_list.push_back(this->get_treated(i, j, number)); + resistant_exposed_list.push_back( + static_cast(std::floor(this->get_treated(i, j, number)))); } std::vector resistant_mortality_list; for (const auto& number : host_pool.mortality_by_group_at(i, j)) { - resistant_mortality_list.push_back(this->get_treated(i, j, number)); + resistant_mortality_list.push_back( + static_cast(std::floor(this->get_treated(i, j, number)))); } + int infected = static_cast( + std::floor(this->get_treated(i, j, host_pool.infected_at(i, j)))); host_pool.make_resistant_at( i, j, susceptible_resistant, resistant_exposed_list, - this->get_treated(i, j, host_pool.infected_at(i, j)), + infected, resistant_mortality_list); } }