From 3bec6214cb3978e0442fdcb2d31b6cbbd4610790 Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Mon, 27 Nov 2023 15:37:42 -0500 Subject: [PATCH 01/12] Use multihost API with single host --- inst/cpp/pops-core | 2 +- inst/include/actions.hpp | 532 +++++++++++ inst/include/anthropogenic_kernel.hpp | 2 + inst/include/competency_table.hpp | 109 +++ inst/include/config.hpp | 83 ++ inst/include/date.hpp | 7 +- inst/include/environment.hpp | 107 ++- inst/include/environment_interface.hpp | 55 ++ inst/include/gamma_kernel.hpp | 1 - inst/include/generator_provider.hpp | 84 +- inst/include/host_pool.hpp | 1122 ++++++++++++++++++++++++ inst/include/host_pool_interface.hpp | 50 ++ inst/include/model.hpp | 315 ++++--- inst/include/model_type.hpp | 63 ++ inst/include/multi_host_pool.hpp | 461 ++++++++++ inst/include/natural_kernel.hpp | 2 + inst/include/network.hpp | 11 +- inst/include/pest_host_use_table.hpp | 81 ++ inst/include/pest_pool.hpp | 134 +++ inst/include/quarantine.hpp | 27 +- inst/include/simulation.hpp | 995 +++++++++------------ inst/include/soils.hpp | 14 +- inst/include/spread_rate.hpp | 180 ++-- inst/include/treatments.hpp | 247 ++---- inst/include/utils.hpp | 31 +- src/pops.cpp | 78 +- 26 files changed, 3751 insertions(+), 1042 deletions(-) create mode 100644 inst/include/actions.hpp create mode 100644 inst/include/competency_table.hpp create mode 100644 inst/include/environment_interface.hpp create mode 100644 inst/include/host_pool.hpp create mode 100644 inst/include/host_pool_interface.hpp create mode 100644 inst/include/model_type.hpp create mode 100644 inst/include/multi_host_pool.hpp create mode 100644 inst/include/pest_host_use_table.hpp create mode 100644 inst/include/pest_pool.hpp diff --git a/inst/cpp/pops-core b/inst/cpp/pops-core index e62d22c0..bf490a9f 160000 --- a/inst/cpp/pops-core +++ b/inst/cpp/pops-core @@ -1 +1 @@ -Subproject commit e62d22c0f7d3d9eed6436c3ce47e2c8b75511645 +Subproject commit bf490a9f41851ea0e980cd05f386d8a041cff892 diff --git a/inst/include/actions.hpp b/inst/include/actions.hpp new file mode 100644 index 00000000..e92bbc54 --- /dev/null +++ b/inst/include/actions.hpp @@ -0,0 +1,532 @@ +/* + * PoPS model - pest or pathogen spread simulation + * + * Copyright (C) 2023 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * Chris Jones (cjones1688 gmail com) + * Anna Petrasova (kratochanna gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_ACTIONS_HPP +#define POPS_ACTIONS_HPP + +#include +#include +#include +#include +#include +#include +#include + +#include "utils.hpp" +#include "model_type.hpp" +#include "soils.hpp" + +namespace pops { + +/** + * Spread of pest or pathogens + * + * The spread starts with generation of dispersers in the hosts. + * This is followed by dispersion and possible establishment. + */ +template< + typename Hosts, + typename Pests, + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex, + typename DispersalKernel, + typename Generator> +class SpreadAction +{ +public: + /** + * @brief Create the object with a given dispersal kernel. + * + * DispersalKernel is callable object or function with one parameter + * which is the random number engine (generator). The return value + * is row and column in the raster (or outside of it). The current + * position is passed as parameters. The return value is in the + * form of a tuple with row and column so that std::tie() is usable + * on the result, i.e. function returning + * `std::make_tuple(row, column)` fulfills this requirement. + * + * @param dispersal_kernel Dispersal kernel to move the dispersers + */ + SpreadAction(DispersalKernel& dispersal_kernel) + : dispersal_kernel_(dispersal_kernel) + {} + + /** Perform the generation and spread */ + void action(Hosts& host_pool, Pests& pests, Generator& generator) + { + this->generate(host_pool, pests, generator); + this->disperse(host_pool, pests, generator); + } + + /** Generates dispersers based on hosts + * + * @note Unlike other functions, this resets the state of dispersers in the pest + * pool. + */ + void generate(Hosts& host_pool, Pests& pests, Generator& generator) + { + for (auto indices : host_pool.suitable_cells()) { + int i = indices[0]; + int j = indices[1]; + int dispersers_from_cell = + host_pool.dispersers_from(i, j, generator.disperser_generation()); + if (dispersers_from_cell > 0) { + if (soil_pool_) { + // 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); + 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); + } + else { + pests.set_dispersers_at(i, j, 0); + pests.set_established_dispersers_at(i, j, 0); + } + } + } + + /** Moves dispersers (dispersing individuals) to dispersal locations + * + * The function sends dispersers previously generated by generate() to locations in + * host pool generated by the dispersal kernel. Any SI-SEI differences are resolved + * by the host. + */ + void disperse(Hosts& host_pool, Pests& pests, Generator& generator) + { + // The interaction does not happen over the member variables yet, use empty + // variables. This requires SI/SEI to be fully resolved in host and not in + // disperse_and_infect. + int row; + int col; + for (auto indices : host_pool.suitable_cells()) { + int i = indices[0]; + int j = indices[1]; + 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 (soil_pool_) { + // Get dispersers from the soil if there are any. + auto num_dispersers = soil_pool_->dispersers_from(i, j, generator); + // Put each disperser to the host pool. + for (int k = 0; k < num_dispersers; k++) { + host_pool.disperser_to(i, j, generator.establishment()); + } + } + } + } + + /** + * @brief Activate storage of dispersers in soil + * + * Calling this function activates the soils. By default, the soil pool is not used. + * The parameters are soil pool used to store the dispersers and + * a percentage (0-1 ratio) of dispersers which will be sent to the soil (and may + * establish or not depending on the soil pool object). + * + * Soil pool is optional and nothing is done when it is not set. + * This function needs to be called separately some time after the object is created + * to activate the soil part of the simulation. This avoids the need for many + * more constructors or for many optional parameters which need default values. + * + * @param soil_pool Soils pool object to use for storage + * @param dispersers_percentage Percentage of dispersers moving to the soil + */ + void activate_soils( + std::shared_ptr> + soil_pool, + double dispersers_percentage) + { + this->soil_pool_ = soil_pool; + this->to_soil_percentage_ = dispersers_percentage; + } + +private: + /** + * Dispersal kernel + */ + DispersalKernel& dispersal_kernel_; + /** + * Optional soil pool + */ + std::shared_ptr> + soil_pool_{nullptr}; + /** + * Percentage (0-1 ratio) of disperers to be send to soil + */ + double to_soil_percentage_{0}; +}; + +/** + * Pest or pathogen individuals survive only with a given survival rate + * + * This removes a given percentage of exposed and infected in host pool. + */ +template +class SurvivalRateAction +{ +public: + /** + * @brief Create object with a given survival rate. + * + * Survival rate is a spatially variable float raster with values between 0 and 1. + * From pest perspective this is the ratio of surviving individuals on hosts, + * while from the host perspective this is the ratio of hosts to keep infected or + * exposed. + * + * @param survival_rate Survival rate + */ + SurvivalRateAction(const FloatRaster& survival_rate) : survival_rate_(survival_rate) + {} + /** + * Reduce the infection based on the pest survival rate. + * + * Infected and total number of exposed hosts are removed directly, while mortality + * cohorts and exposed cohorts are left to be managed by the host pool. In other + * words, the details of how total infected and total exposed are distributed among + * the cohorts is managed by the host pools. + */ + template + void action(Hosts& hosts, Generator& generator) + { + for (auto indices : hosts.suitable_cells()) { + int i = indices[0]; + int j = indices[1]; + if (survival_rate_(i, j) < 1) { + hosts.remove_infection_by_ratio_at( + i, j, survival_rate_(i, j), generator.survival_rate()); + } + } + } + +private: + const FloatRaster& survival_rate_; +}; + +/** + * Remove infection or infestation by temperature threshold + */ +template< + typename Hosts, + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex, + typename Generator> +class RemoveByTemperature +{ +public: + /** + * @brief Create an object with reference to the environment and the temperature + * threshold + * + * All infection in a cell with temperature under the *lethal_temperature* value is + * removed from hosts. + * + * @param environment Environment which supplies the temperature per cell + * @param lethal_temperature Temperature at which lethal conditions occur + */ + RemoveByTemperature( + const EnvironmentInterface& + environment, + double lethal_temperature) + : environment_(environment), lethal_temperature_(lethal_temperature) + {} + /** Perform the removal of infection */ + void action(Hosts& hosts, Generator& generator) + { + for (auto indices : hosts.suitable_cells()) { + int i = indices[0]; + int j = indices[1]; + if (environment_.temperature_at(i, j) < lethal_temperature_) { + // now this includes also mortality, but it does not include exposed + hosts.remove_all_infected_at(i, j, generator.lethal_temperature()); + } + } + } + +private: + const EnvironmentInterface& + environment_; + const double lethal_temperature_; +}; + +/** + * Move overpopulated pests in one cell to a different cell. + * + * Overpopulation is measured by the ratio of infected hosts to all hosts. + */ +template< + typename Hosts, + typename Pests, + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex, + typename DispersalKernel> +class MoveOverpopulatedPests +{ +public: + /** + * @brief Create object with given dispersal kernel and overpopulation parameters + * + * @param dispersal_kernel Dispersal kernel used for pest individuals leaving a cell + * @param overpopulation_percentage Ratio of infected to all hosts considered as + * overpopulation + * @param leaving_percentage Ratio of pest individuals leaving the cell + * @param rows Number of rows in the whole area (for tracking outside dispersers) + * @param cols Number of columns in the whole area (for tracking outside dispersers) + */ + MoveOverpopulatedPests( + DispersalKernel& dispersal_kernel, + double overpopulation_percentage, + double leaving_percentage, + RasterIndex rows, + RasterIndex cols) + : dispersal_kernel_(dispersal_kernel), + overpopulation_percentage_(overpopulation_percentage), + leaving_percentage_(leaving_percentage), + rows_(rows), + cols_(cols) + {} + /** Move overflowing pest population to other hosts. + * + * When the number of pests (pest population) is too high, part of them moves + * to a different location. Number of infected/infested hosts is considered to be + * the number of pests (groups of pest) in a raster cell. + * + * The movement happens in two stages. First, all the leaving pests are identified + * and removed from the source cells. Second, the move to the target cells is + * performed. This means that even if the resulting number of pests in the target + * cell is considered too high, it is left as is and the move is performed the next + * time this function is called. + * + * If the pests (pest population) cannot be accommodated in the target cell due to + * the insufficient number of susceptible hosts, the excessive pests die. + * + * @note Exposed hosts do not count towards total number of pest, + * i.e., *total_host* is assumed to be S + E in SEI model. + * @note Mortality and exposure are not supported by this function, i.e., the + * mortality rasters are not modified while the infected are. This is left to be + * managed by the host pool, but the current host pool does not support that. + */ + template + void action(Hosts& hosts, Pests& pests, Generator& generator) + { + struct Move + { + RasterIndex row; + RasterIndex col; + int count; + }; + std::vector moves; + + // Identify the moves. Remove from source cells. + for (auto indices : hosts.suitable_cells()) { + int i = indices[0]; + int j = indices[1]; + int original_count = hosts.infected_at(i, j); + // No move with only one infected host (one unit). + if (original_count <= 1) + continue; + // r = I / (I + S) + // r = I / (I + S + E_1 + E_2 + ...) + double ratio = original_count / double(hosts.total_hosts_at(i, j)); + if (ratio >= overpopulation_percentage_) { + int row; + int col; + std::tie(row, col) = + dispersal_kernel_(generator.overpopulation(), i, j); + // for leaving_percentage == 0.5 + // 2 infected -> 1 leaving + // 3 infected -> 1 leaving + int leaving = 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); + continue; + } + // Doing the move here would create inconsistent results as some + // target cells would be evaluated after the moved pest arrived, + // possibly triggering another move. + // So, instead, we just collect them and apply later. (The pest is in + // the air when in the moves vector.) + moves.push_back(Move({row, col, leaving})); + } + } + // Perform the moves to target cells. + for (const auto& move : moves) { + // Pests which won't fit are ignored (disappear). This can happen if there + // is simply not enough S hosts to accommodate all the pests moving from the + // source or if multiple sources end up in the same target cell and there is + // not enough S hosts to accommodate all of them. The decision is made in + // the host pool. Here, we ignore the return value specifying the number of + // accepted pests. + hosts.pests_to(move.row, move.col, move.count, generator.overpopulation()); + } + } + +private: + DispersalKernel& dispersal_kernel_; + double overpopulation_percentage_; + double leaving_percentage_; + RasterIndex rows_; + RasterIndex cols_; +}; + +/** + * Moves hosts from one location to another + * + * @note Mortality and non-host individuals are not supported in movements. + */ +template< + typename Hosts, + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex> +class HostMovement +{ +public: + /** + * @brief Creates an object with state pointing to the movements to be used + * + * @param step the current step of the simulation + * @param last_index the last index to not be used from movements + * @param movements a vector of ints with row_from, col_from, row_to, col_to, and + * num_hosts + * @param movement_schedule a vector matching movements with the step at which the + * movement from movements are applied + * + * @note There is a mix of signed and unsigned ints for step and movement schedule. + */ + HostMovement( + unsigned step, + unsigned last_index, + const std::vector>& movements, + const std::vector& movement_schedule) + : step_(step), + last_index_(last_index), + movements_(movements), + movement_schedule_(movement_schedule) + {} + + /** Perform the movement + * + * @note This one returns a value, so to create a consistent interface for all + * actions, the tracking would have to happen in the action which would be better + * design anyway. + */ + template + unsigned action(Hosts& hosts, Generator& generator) + { + for (unsigned i = last_index_; i < movements_.size(); i++) { + const auto& moved = movements_[i]; + unsigned move_schedule = movement_schedule_[i]; + if (move_schedule != step_) { + return i; + } + int row_from = moved[0]; + int col_from = moved[1]; + int row_to = moved[2]; + int col_to = moved[3]; + hosts.move_hosts_from_to( + row_from, col_from, row_to, col_to, moved[4], generator.movement()); + } + return movements_.size(); + } + +private: + const unsigned step_; + const unsigned last_index_; + const std::vector>& movements_; + const std::vector& movement_schedule_; +}; + +/** + * Mortality of the hosts + * + * Kills infected hosts based on mortality rate and timing. The host pool implementation + * of mortality is used to perform the actual process while this class provides + * parameters and takes care of the two steps in the mortality process (dying and + * aging). + * + * The *mortality_tracker_vector* used by hosts must fit with the *mortality_time_lag* + * here. It needs to have a minimum size of mortality_time_lag + 1. See + * HostPool::apply_mortality_at() and HostPool::step_forward_mortality(). + */ +template +class Mortality +{ +public: + /** + * @brief Create object which will let host decide its mortality parameters. + */ + Mortality() : action_mortality_(false) {} + /** + * @brief Create object with fixed mortality rate and time lag. + * + * @param mortality_rate Percent of infected hosts that die each time period + * @param mortality_time_lag Time lag prior to mortality beginning + */ + Mortality(double mortality_rate, int mortality_time_lag) + : mortality_rate_(mortality_rate), + mortality_time_lag_(mortality_time_lag), + action_mortality_(true) + {} + /** + * Perform the action by applying mortality and moving the mortality tracker + * forward. + */ + void action(Hosts& hosts) + { + for (auto indices : hosts.suitable_cells()) { + if (action_mortality_) { + hosts.apply_mortality_at( + indices[0], indices[1], mortality_rate_, mortality_time_lag_); + } + else { + hosts.apply_mortality_at(indices[0], indices[1]); + } + } + hosts.step_forward_mortality(); + } + +private: + const double mortality_rate_ = 0; + const int mortality_time_lag_ = 0; + const double action_mortality_ = false; +}; + +} // namespace pops + +#endif // POPS_ACTIONS_HPP diff --git a/inst/include/anthropogenic_kernel.hpp b/inst/include/anthropogenic_kernel.hpp index 477b8bc8..16ff727b 100644 --- a/inst/include/anthropogenic_kernel.hpp +++ b/inst/include/anthropogenic_kernel.hpp @@ -35,8 +35,10 @@ namespace pops { * Same structure as the natural kernel, but the parameters are for anthropogenic * kernel when available. * + * @param config Configuration for the kernel * @param dispersers The disperser raster (reference, for deterministic kernel) * @param network Network (initialized or not) + * * @return Created kernel */ template diff --git a/inst/include/competency_table.hpp b/inst/include/competency_table.hpp new file mode 100644 index 00000000..90aa496b --- /dev/null +++ b/inst/include/competency_table.hpp @@ -0,0 +1,109 @@ +/* + * PoPS model - Competency table for hosts and pest + * + * Copyright (C) 2023 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_COMPETENCY_TABLE_HPP +#define POPS_COMPETENCY_TABLE_HPP + +#include +#include +#include + +#include "config.hpp" + +namespace pops { + +template +class CompetencyTable +{ +public: + using Environment = typename HostPool::Environment; + + CompetencyTable(const Environment& environment) : environment_(environment) {} + CompetencyTable(const Config& config, const Environment& environment) + : environment_(environment) + { + for (const auto& row : config.competency_table_data()) { + competency_table_.emplace_back(row.presence_absence, row.competency); + } + } + + void + add_host_competencies(const std::vector& presence_absence, double competency) + { + competency_table_.emplace_back(presence_absence, competency); + } + + double competency_at(RasterIndex row, RasterIndex col, const HostPool* host) const + { + auto presence_absence = environment_.host_presence_at(row, col); + auto host_index = environment_.host_index(host); + return find_competency(presence_absence, host_index); + } + + double + find_competency(const std::vector& presence_absence, size_t host_index) const + { + // Go over all the rows and find the highest competency which fulfilled the + // presence criteria. + double competency = 0; + for (const auto& table_row : competency_table_) { + // probably faster if we just wait for the iteration over the whole thing + if (!table_row.presence_absence[host_index]) + continue; + if (presence_absence.size() != table_row.presence_absence.size()) { + throw std::invalid_argument( + "Number of hosts in the environment is not the same as " + "the number of hosts in the competency table (" + + std::to_string(presence_absence.size()) + + " != " + std::to_string(table_row.presence_absence.size()) + ")"); + } + // Pick always the highest competency. Don't test row which has lower one. + if (table_row.competency <= competency) + continue; + // Row is considered only if all required hosts are present. + bool condition_fulfilled = true; + for (size_t i = 0; i < presence_absence.size(); ++i) { + if (!table_row.presence_absence[i]) { + continue; + } + if (!presence_absence[i]) { + condition_fulfilled = false; + break; + } + } + if (condition_fulfilled) + competency = table_row.competency; + } + return competency; + } + +private: + struct TableRow + { + std::vector presence_absence; + double competency; + + // Constructor can be removed with C++20, GCC 10, clang 16. + TableRow(const std::vector& presence_absence, double competency) + : presence_absence(presence_absence), competency(competency) + {} + }; + std::vector competency_table_; + const Environment& environment_; +}; + +} // namespace pops + +#endif // POPS_COMPETENCY_TABLE_HPP diff --git a/inst/include/config.hpp b/inst/include/config.hpp index 2d7f535f..14c4c770 100644 --- a/inst/include/config.hpp +++ b/inst/include/config.hpp @@ -118,6 +118,18 @@ std::map read_key_value_pairs( class Config { public: + struct PestHostUseTableDataRow + { + double susceptibility; + double mortality_rate; + double mortality_time_lag; + }; + struct CompetencyTableDataRow + { + std::vector presence_absence; + double competency; + }; + // Seed int random_seed{0}; bool multiple_random_seeds{false}; @@ -438,6 +450,20 @@ class Config season_end_month_ = std::stoi(end); } + void set_arrival_behavior(std::string value) + { + if (value != "infect" && value != "land") { + throw std::invalid_argument( + "arrival behavior can be 'infect' or 'land' but not: " + value); + } + arrival_behavior_ = value; + } + + const std::string& arrival_behavior() const + { + return arrival_behavior_; + } + /** * Read seeds from text. * @@ -468,6 +494,7 @@ class Config "anthropogenic_dispersal", "establishment", "weather", + "lethal_temperature", "movement", "overpopulation", "survival_rate", @@ -485,6 +512,57 @@ class Config this->multiple_random_seeds = true; } + const std::vector& pest_host_use_table_data() const + { + return pest_host_use_table_data_; + } + const std::vector& competency_table_data() const + { + return competency_table_data_; + } + + void read_pest_host_use_table(const std::vector>& values) + { + for (const auto& row : values) { + if (row.size() < 3) { + throw std::invalid_argument( + "3 values are required for each pest-host-use table row"); + } + PestHostUseTableDataRow resulting_row; + resulting_row.susceptibility = row[0]; + resulting_row.mortality_rate = row[1]; + resulting_row.mortality_time_lag = row[2]; + pest_host_use_table_data_.push_back(std::move(resulting_row)); + } + } + + void create_pest_host_use_table_from_parameters(int num_of_hosts) + { + for (int i = 0; i < num_of_hosts; ++i) { + PestHostUseTableDataRow resulting_row; + resulting_row.susceptibility = 1; + resulting_row.mortality_rate = this->mortality_rate; + resulting_row.mortality_time_lag = this->mortality_time_lag; + pest_host_use_table_data_.push_back(std::move(resulting_row)); + } + } + + void read_competency_table(const std::vector>& values) + { + for (const auto& row : values) { + if (row.size() < 2) { + throw std::invalid_argument( + "At least 2 values are required for each competency table row"); + } + CompetencyTableDataRow resulting_row; + for (auto it = row.begin(); it < std::prev(row.end()); ++it) { + resulting_row.presence_absence.push_back(bool(*it)); + } + resulting_row.competency = row.back(); + competency_table_data_.push_back(std::move(resulting_row)); + } + } + private: Date date_start_{"0-01-01"}; Date date_end_{"0-01-02"}; @@ -498,6 +576,8 @@ class Config Scheduler scheduler_{date_start_, date_end_, step_unit_, step_num_units_}; bool schedules_created_{false}; + std::string arrival_behavior_{"infect"}; + std::vector spread_schedule_; std::vector output_schedule_; std::vector mortality_schedule_; @@ -506,6 +586,9 @@ class Config std::vector spread_rate_schedule_; std::vector quarantine_schedule_; std::vector weather_table_; + + std::vector pest_host_use_table_data_; + std::vector competency_table_data_; }; } // namespace pops diff --git a/inst/include/date.hpp b/inst/include/date.hpp index 3b5b5412..32567279 100644 --- a/inst/include/date.hpp +++ b/inst/include/date.hpp @@ -446,9 +446,10 @@ int Date::weeks_from_date(Date start) /*! * Returns the current date as a string - * Will be replaced by C++20 format function - * \param Date, current date - * \return String, current date as string + * + * Will be replaced by C++20 format function. + * + * \return current date as a string */ std::string Date::to_string() { diff --git a/inst/include/environment.hpp b/inst/include/environment.hpp index 30f5aa18..1cdedf9b 100644 --- a/inst/include/environment.hpp +++ b/inst/include/environment.hpp @@ -24,8 +24,10 @@ #include #include +#include "environment_interface.hpp" #include "normal_distribution_with_uniform_fallback.hpp" #include "utils.hpp" +#include "host_pool_interface.hpp" #include "generator_provider.hpp" namespace pops { @@ -64,8 +66,13 @@ inline WeatherType weather_type_from_string(const std::string& text) * * Currently, only handles weather coefficient for soils. Holds only the current state. */ -template +template< + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex, + typename Generator> class Environment + : public EnvironmentInterface { public: Environment() {} @@ -75,9 +82,10 @@ class Environment * * @param raster Raster with the weather coefficient. */ - void update_weather_coefficient(const FloatRaster& raster) + void update_weather_coefficient(const FloatRaster& raster) override { current_weather_coefficient = &raster; + weather_ = true; } /** @@ -99,9 +107,10 @@ class Environment * @throw std::invalid_argument when dimensions of *mean* and *stddev* differ or * when mean is out of range */ - template void update_weather_from_distribution( - const FloatRaster& mean, const FloatRaster& stddev, Generator& generator) + const FloatRaster& mean, + const FloatRaster& stddev, + Generator& generator) override { if (mean.rows() != stddev.rows()) { throw std::invalid_argument( @@ -142,6 +151,7 @@ class Environment } } current_weather_coefficient = &stored_weather_coefficient; + weather_ = true; } /** @@ -153,7 +163,7 @@ class Environment * * @throw std::logic_error when coefficient is not set */ - double weather_coefficient_at(RasterIndex row, RasterIndex col) const + double weather_coefficient_at(RasterIndex row, RasterIndex col) const override { if (!current_weather_coefficient) { throw std::logic_error("Weather coefficient used, but not provided"); @@ -161,6 +171,71 @@ class Environment return current_weather_coefficient->operator()(row, col); } + double influence_reproductive_rate_at( + RasterIndex row, RasterIndex col, double value) const override + { + if (!weather_) + return value; + return value * weather_coefficient_at(row, col); + } + + double influence_probability_of_establishment_at( + RasterIndex row, RasterIndex col, double value) const override + { + if (!weather_) + return value; + return value * weather_coefficient_at(row, col); + } + + int total_population_at(RasterIndex row, RasterIndex col) const override + { + // If total population is used, use that instead of computing it. + if (total_population_) + return total_population_->operator()(row, col); + int sum = 0; + if (other_individuals_) + sum += other_individuals_->operator()(row, col); + for (const auto& host : hosts_) + sum += host->total_hosts_at(row, col); + return sum; + } + + std::vector host_presence_at(RasterIndex row, RasterIndex col) const override + { + std::vector presence; + presence.reserve(hosts_.size()); + for (const auto& host : hosts_) + presence.push_back(host->total_hosts_at(row, col)); + return presence; + } + + void set_other_individuals(const IntegerRaster* individuals) override + { + other_individuals_ = individuals; + } + + void set_total_population(const IntegerRaster* individuals) override + { + total_population_ = individuals; + } + + void add_host(const HostPoolInterface* host) override + { + // no-op if already there, may become an error in the future + if (container_contains(hosts_, host)) + return; + hosts_.push_back(host); + } + + size_t host_index(const HostPoolInterface* host) const override + { + auto it = std::find(hosts_.begin(), hosts_.end(), host); + if (it == hosts_.end()) + throw std::invalid_argument( + "Environment::host_index: Host is not in the environment"); + return std::distance(hosts_.begin(), it); + } + /** * @brief Get weather coefficient raster * @@ -168,7 +243,7 @@ class Environment * * @throw std::logic_error when coefficient is not set */ - const FloatRaster& weather_coefficient() const + const FloatRaster& weather_coefficient() const override { if (!current_weather_coefficient) { throw std::logic_error("Weather coefficient used, but not provided"); @@ -176,6 +251,19 @@ class Environment return *current_weather_coefficient; } + void update_temperature(const FloatRaster& raster) override + { + temperature_ = &raster; + } + + double temperature_at(RasterIndex row, RasterIndex col) const override + { + if (!temperature_) { + throw std::logic_error("Temperature used, but not provided"); + } + return temperature_->operator()(row, col); + } + protected: static constexpr double weather_coefficient_min = 0; static constexpr double weather_coefficient_max = 1; @@ -187,6 +275,13 @@ class Environment */ const FloatRaster* current_weather_coefficient{nullptr}; FloatRaster stored_weather_coefficient; + bool weather_{false}; + + std::vector*> hosts_; // host, non-owning + const IntegerRaster* other_individuals_{nullptr}; // non-hosts, non-owning + const IntegerRaster* total_population_{nullptr}; // non-hosts, non-owning + + const FloatRaster* temperature_{nullptr}; }; } // namespace pops diff --git a/inst/include/environment_interface.hpp b/inst/include/environment_interface.hpp new file mode 100644 index 00000000..62e3e0f2 --- /dev/null +++ b/inst/include/environment_interface.hpp @@ -0,0 +1,55 @@ +/* + * PoPS model - pest or pathogen spread simulation + * + * Copyright (C) 2023 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_ENVIRONMENT_INTERFACE_HPP +#define POPS_ENVIRONMENT_INTERFACE_HPP + +#include "host_pool_interface.hpp" + +namespace pops { + +template< + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex, + typename Generator> +class EnvironmentInterface +{ +public: + virtual ~EnvironmentInterface() {} + virtual void update_weather_coefficient(const FloatRaster& raster) = 0; + virtual void update_weather_from_distribution( + const FloatRaster& mean, const FloatRaster& stddev, Generator& generator) = 0; + virtual double weather_coefficient_at(RasterIndex row, RasterIndex col) const = 0; + virtual double influence_reproductive_rate_at( + RasterIndex row, RasterIndex col, double value) const = 0; + virtual double influence_probability_of_establishment_at( + RasterIndex row, RasterIndex col, double value) const = 0; + virtual int total_population_at(RasterIndex row, RasterIndex col) const = 0; + virtual std::vector + host_presence_at(RasterIndex row, RasterIndex col) const = 0; + // in general, we may have other individuals-non const + virtual void set_other_individuals(const IntegerRaster* individuals) = 0; + virtual void set_total_population(const IntegerRaster* individuals) = 0; + virtual void add_host(const HostPoolInterface* host) = 0; + virtual size_t host_index(const HostPoolInterface* host) const = 0; + virtual const FloatRaster& weather_coefficient() const = 0; + virtual void update_temperature(const FloatRaster& raster) = 0; + virtual double temperature_at(RasterIndex row, RasterIndex col) const = 0; +}; + +} // namespace pops + +#endif // POPS_ENVIRONMENT_INTERFACE_HPP diff --git a/inst/include/gamma_kernel.hpp b/inst/include/gamma_kernel.hpp index af57e1c2..2d9e31b9 100644 --- a/inst/include/gamma_kernel.hpp +++ b/inst/include/gamma_kernel.hpp @@ -139,7 +139,6 @@ class GammaKernel } } throw std::invalid_argument("unable to find solution to gamma icdf "); - return -1; } }; diff --git a/inst/include/generator_provider.hpp b/inst/include/generator_provider.hpp index 2ab56157..b751add3 100644 --- a/inst/include/generator_provider.hpp +++ b/inst/include/generator_provider.hpp @@ -51,6 +51,7 @@ class RandomNumberGeneratorProviderInterface virtual Generator& anthropogenic_dispersal() = 0; virtual Generator& establishment() = 0; virtual Generator& weather() = 0; + virtual Generator& lethal_temperature() = 0; virtual Generator& movement() = 0; virtual Generator& overpopulation() = 0; virtual Generator& survival_rate() = 0; @@ -66,10 +67,12 @@ class RandomNumberGeneratorProviderInterface * as this generator object can be used directly or, more importantly, * standard generator can be used in its place. */ -template -class SingleGeneratorProvider : public RandomNumberGeneratorProviderInterface +template +class SingleGeneratorProvider + : public RandomNumberGeneratorProviderInterface { public: + using Generator = GeneratorType; /** * @brief Seeds the underlying generator * @param seed for the underlying generator @@ -80,14 +83,14 @@ class SingleGeneratorProvider : public RandomNumberGeneratorProviderInterface& seeds) + void seed(const std::map& seeds) override { UNUSED(seeds); throw std::invalid_argument( @@ -98,7 +101,7 @@ class SingleGeneratorProvider : public RandomNumberGeneratorProviderInterface& seeds) + void seed(const std::map& seeds) override { this->set_seed_by_name( seeds, "disperser_generation", disperser_generation_generator_); @@ -254,6 +263,7 @@ class MultiRandomNumberGeneratorProvider seeds, "anthropogenic_dispersal", anthropogenic_dispersal_generator_); this->set_seed_by_name(seeds, "establishment", establishment_generator_); this->set_seed_by_name(seeds, "weather", weather_generator_); + this->set_seed_by_name(seeds, "lethal_temperature", lethal_temperature_); this->set_seed_by_name(seeds, "movement", movement_generator_); this->set_seed_by_name(seeds, "overpopulation", overpopulation_generator_); this->set_seed_by_name(seeds, "survival_rate", survival_rate_generator_); @@ -261,7 +271,7 @@ class MultiRandomNumberGeneratorProvider } /** Re-seed using named seeds, otherwise increment single seed */ - void seed(const Config& config) + void seed(const Config& config) override { if (!config.random_seeds.empty()) { this->seed(config.random_seeds); @@ -271,47 +281,52 @@ class MultiRandomNumberGeneratorProvider } } - Generator& disperser_generation() + Generator& disperser_generation() override { return disperser_generation_generator_; } - Generator& natural_dispersal() + Generator& natural_dispersal() override { return natural_dispersal_generator_; } - Generator& anthropogenic_dispersal() + Generator& anthropogenic_dispersal() override { return anthropogenic_dispersal_generator_; } - Generator& establishment() + Generator& establishment() override { return establishment_generator_; } - Generator& weather() + Generator& weather() override { return weather_generator_; } - Generator& movement() + Generator& lethal_temperature() override + { + return lethal_temperature_; + } + + Generator& movement() override { return movement_generator_; } - Generator& overpopulation() + Generator& overpopulation() override { return overpopulation_generator_; } - Generator& survival_rate() + Generator& survival_rate() override { return survival_rate_generator_; } - Generator& soil() + Generator& soil() override { return soil_generator_; } @@ -337,6 +352,7 @@ class MultiRandomNumberGeneratorProvider Generator anthropogenic_dispersal_generator_; Generator establishment_generator_; Generator weather_generator_; + Generator lethal_temperature_; Generator movement_generator_; Generator overpopulation_generator_; Generator survival_rate_generator_; @@ -359,10 +375,11 @@ class MultiRandomNumberGeneratorProvider * an exception if used directly as UniformRandomBitGenerator, but the * object was seeded with multiple seeds. */ -template +template class RandomNumberGeneratorProvider { public: + using Generator = GeneratorType; /** * Seeds first generator with the seed and then each subsequent generator with * seed += 1. *multi* decides if single generator is created or if multiple @@ -424,6 +441,11 @@ class RandomNumberGeneratorProvider return impl->weather(); } + Generator& lethal_temperature() + { + return impl->lethal_temperature(); + } + Generator& movement() { return impl->movement(); @@ -448,12 +470,12 @@ class RandomNumberGeneratorProvider using result_type = typename Generator::result_type; - static result_type min() + static constexpr result_type min() { return Generator::min(); } - static result_type max() + static constexpr result_type max() { return Generator::max(); } diff --git a/inst/include/host_pool.hpp b/inst/include/host_pool.hpp new file mode 100644 index 00000000..dafb8b20 --- /dev/null +++ b/inst/include/host_pool.hpp @@ -0,0 +1,1122 @@ +/* + * PoPS model - pest or pathogen spread simulation + * + * Copyright (C) 2023 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_HOST_POOL_HPP +#define POPS_HOST_POOL_HPP + +#include +#include +#include +#include + +#include "host_pool_interface.hpp" +#include "model_type.hpp" +#include "environment_interface.hpp" +#include "competency_table.hpp" +#include "pest_host_use_table.hpp" + +namespace pops { + +/** + * Host pool managing susceptible, exposed, infected, and resistant hosts. + * + * @tparam IntegerRaster Integer raster type + * @tparam FloatRaster Floating point raster type + * @tparam RasterIndex Type for indexing the rasters + * @tparam GeneratorProvider Provider of random number generators + * + * GeneratorProvider needs to provide Generator memeber which is the type of the + * underlying random number generators. + */ +template< + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex, + typename GeneratorProvider> +class HostPool : public HostPoolInterface +{ +public: + /** + * Type of environment object providing information about weather and other + * environmental properties. + */ + using Environment = EnvironmentInterface< + IntegerRaster, + FloatRaster, + RasterIndex, + GeneratorProvider>; + /** + * Standard random number generator to be passed directly to the methods. + */ + using Generator = typename GeneratorProvider::Generator; + + /** + * @brief Creates an object with stored references and host properties. + * + * The *exposed* vector is a list of hosts exposed in the previous steps. + * The length of the vector is the number of steps of the latency + * period plus one. See the step_forward() function for details on how + * the vector is used. + * + * The *total_populations* can be total number of hosts in the basic case + * or it can be the total size of population of all relevant species + * both host and non-host if dilution effect should be applied. + * + * If establishment stochasticity is disabled, + * *establishment_probability* is used to decide whether or not + * a disperser is established in a cell. Value 1 means that all + * dispersers will establish and value 0 means that no dispersers + * will establish. + * + * The *mortality_tracker_vector* is a vector of matrices for tracking infected + * host infection over time. Expectation is that mortality tracker is of + * length (1/mortality_rate + mortality_time_lag). + * + * @param model_type Type of the model (SI or SEI) + * @param susceptible Raster of susceptible hosts + * @param exposed Raster of exposed or infected hosts + * @param latency_period Length of the latency period in steps + * @param infected Infected hosts + * @param total_exposed Raster tracking all exposed hosts + * @param resistant Resistant hosts + * @param mortality_tracker_vector Raster tracking hosts for mortality + * @param died Raster tracking all hosts which died + * @param total_hosts Total number of hosts + * @param environment Environment which influences the processes + * @param dispersers_stochasticity Enable stochasticity in generating of dispersers + * @param reproductive_rate Reproductive rate in ideal conditions + * @param establishment_stochasticity Enable stochasticity in establishment step + * @param establishment_probability Fixed probability disperser establishment + * @param rows Number of rows for the study area + * @param cols Number of columns for the study area + * @param suitable_cells Cells where hosts are known to occur + */ + HostPool( + ModelType model_type, + IntegerRaster& susceptible, + std::vector& exposed, + unsigned latency_period, + IntegerRaster& infected, + IntegerRaster& total_exposed, + IntegerRaster& resistant, + std::vector& mortality_tracker_vector, + IntegerRaster& died, + IntegerRaster& total_hosts, + const Environment& environment, + bool dispersers_stochasticity, + double reproductive_rate, + bool establishment_stochasticity, + double establishment_probability, + RasterIndex rows, + RasterIndex cols, + std::vector>& suitable_cells) + : susceptible_(susceptible), + infected_(infected), + exposed_(exposed), + latency_period_(latency_period), + total_exposed_(total_exposed), + resistant_(resistant), + mortality_tracker_vector_(mortality_tracker_vector), + died_(died), + total_hosts_(total_hosts), + environment_(environment), + model_type_(model_type), + dispersers_stochasticity_(dispersers_stochasticity), + reproductive_rate_(reproductive_rate), + establishment_stochasticity_(establishment_stochasticity), + deterministic_establishment_probability_(establishment_probability), + rows_(rows), + cols_(cols), + suitable_cells_(suitable_cells) + {} + + void set_pest_host_use_table(const PestHostUseTable& pest_host_use_table) + { + this->pest_host_use_table_ = &pest_host_use_table; + } + + void + set_competency_table(const CompetencyTable& competency_table) + { + this->competency_table_ = &competency_table; + } + + /** + * @brief Move disperser to a cell in the host pool + * + * Processes event when a disperser lands in a cell potentially establishing on a + * host. The disperser may or may not establish a based on host availability, + * weather, establishment probability, and stochasticity. + * + * Any new dispersers targeting host in the host pool should be processed using this + * function. + * + * @param row Row number of the target cell + * @param col Column number of the target cell + * @param generator Random number generator + * + * @return true if disperser has established in the cell, false otherwise + * + * @throw std::runtime_error if model type is unsupported (i.e., not SI or SEI) + */ + int disperser_to(RasterIndex row, RasterIndex col, Generator& generator) + { + if (susceptible_(row, col) <= 0) + return 0; + double probability_of_establishment = establishment_probability_at(row, col); + bool establish = can_disperser_establish( + probability_of_establishment, + establishment_stochasticity_, + deterministic_establishment_probability_, + generator); + if (establish) + return add_disperser_at(row, col); + return 0; + } + + static bool can_disperser_establish( + double probability_of_establishment, + bool establishment_stochasticity, + double deterministic_establishment_probability, + Generator& generator) + { + std::uniform_real_distribution distribution_uniform(0.0, 1.0); + double establishment_tester = 1 - deterministic_establishment_probability; + if (establishment_stochasticity) + establishment_tester = distribution_uniform(generator); + if (establishment_tester < probability_of_establishment) + return true; + return false; + } + + /** + * @brief Add disperser to a cell + * + * Turns disperser into infection considering model type (SI, SEI). + * + * Unlike disperser_to(), this is transforming a disperser into infection right away + * without any further evaluation of establishment or stochasticity. + * + * @param row Row number of the target cell + * @param col Column number of the target cell + * + * @return 1 if disperser was added, 0 if there was not available host. + * + * @throw std::runtime_error if model type is unsupported (i.e., not SI or SEI) + * + * @note This may be merged with pests_to() in the future. + */ + int add_disperser_at(RasterIndex row, RasterIndex col) + { + if (susceptible_(row, col) <= 0) + return 0; + susceptible_(row, col) -= 1; + if (model_type_ == ModelType::SusceptibleInfected) { + infected_(row, col) += 1; + mortality_tracker_vector_.back()(row, col) += 1; + } + else if (model_type_ == ModelType::SusceptibleExposedInfected) { + exposed_.back()(row, col) += 1; + total_exposed_(row, col) += 1; + } + else { + throw std::runtime_error( + "Unknown ModelType value in HostPool::add_disperser_at()"); + } + return 1; + } + + /** + * @brief Get dispersers produced in a cell + * + * Each time the function is called it generates number of dispersers based on the + * current infection, host attributes, and the environment. + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param generator Random number generator + * @return Number of generated dispersers + */ + int dispersers_from(RasterIndex row, RasterIndex col, Generator& generator) const + { + if (infected_at(row, col) <= 0) + return 0; + double lambda = + environment_.influence_reproductive_rate_at(row, col, reproductive_rate_); + if (competency_table_) { + lambda *= competency_table_->competency_at(row, col, this); + } + int dispersers_from_cell = 0; + if (dispersers_stochasticity_) { + std::poisson_distribution distribution(lambda); + for (int k = 0; k < infected_at(row, col); k++) { + dispersers_from_cell += distribution(generator); + } + } + else { + dispersers_from_cell = lambda * infected_at(row, col); + } + return dispersers_from_cell; + } + + /** + * @brief Get establishment probability for a cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Establishment probability + */ + double establishment_probability_at(RasterIndex row, RasterIndex col) const + { + double probability_of_establishment = + (double)(susceptible_(row, col)) + / environment_.total_population_at(row, col); + if (pest_host_use_table_) { + probability_of_establishment *= pest_host_use_table_->susceptibility(this); + } + return environment_.influence_probability_of_establishment_at( + row, col, probability_of_establishment); + } + + /** + * @brief Move pests from a cell + * + * Moves pest from a cell reducing number of infected and increasing number of + * susceptible hosts in the cell. + * + * This is a request to move the number of pests given by *count* from a cell. + * The function may reduce the number based on the availability of pests in the cell + * or internal mechanics of the pool. However, the current implementation simply + * moves the given number of pests without any checks. + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param count Number of pests requested to move from the cell + * @param generator Random number generator (for compatibility with multi-host API) + * + * @return Number of pests actually moved from the cell + * + * @note For consitency with the previous implementation, this does not modify + * mortality cohorts nor touches the exposed cohorts. + */ + int pests_from(RasterIndex row, RasterIndex col, int count, Generator& generator) + { + UNUSED(generator); + susceptible_(row, col) += count; + infected_(row, col) -= count; + return count; + } + + /** + * @brief Move pests to a cell + * + * This directly modifies the number of infected hosts. Unlike disperser_to(), this + * is assuming the pests will establish right away without any further evaluation of + * establishment or stochasticity. + * + * When there is more pests than the target cell can accept based on the number of + * susceptible hosts, all susceptible hosts are infected. The number of actually + * infected hosts (accepted pests) is returned. + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param count Number of pests requested to move to the cell + * @param generator Random number generator (for compatibility with multi-host API) + * + * @return Number of accepted pests + * + * @note For consistency with the previous implementation, this does not make hosts + * exposed in the SEI model. Instead, the hosts are infected right away. This may + * become a feature in the future. + * + * @note For consistency with the previous implementation, this does not modify the + * mortality cohorts. This wil need to be fixed in the future. + * + * @note This may be merged with add_disperser_at() in the future. + */ + int pests_to(RasterIndex row, RasterIndex col, int count, Generator& generator) + { + UNUSED(generator); + // The target cell can accept all. + if (susceptible_(row, col) >= count) { + susceptible_(row, col) -= count; + infected_(row, col) += count; + } + // More pests than the target cell can accept. + else { + count = susceptible_(row, col); + susceptible_(row, col) -= count; + infected_(row, col) += count; + } + return count; + } + + // + /** + * @brief Move hosts from a cell to a cell + * + * Moves specified number of hosts from a cell to another cell. The distribution + * of the hosts among susceptible, exposed, and other pools is stochastic. + * + * @param row_from Row index of the source cell + * @param col_from Column index of the source cell + * @param row_to Row index of the target cell + * @param col_to Column index of the target cell + * @param count Number of pests to move + * @param generator Random number generator for distribution of host among pools + * + * @return Number of pests actually moved between cells + * + * @note Mortality is not supported. + */ + int move_hosts_from_to( + RasterIndex row_from, + RasterIndex col_from, + RasterIndex row_to, + RasterIndex col_to, + int count, + Generator& generator) + { + int total_hosts_moved{count}; + if (count > total_hosts_(row_from, col_from)) { + total_hosts_moved = total_hosts_(row_from, col_from); + } + int total_infecteds = infected_(row_from, col_from); + int suscepts = susceptible_(row_from, col_from); + int expose = total_exposed_(row_from, col_from); + int resist = resistant_(row_from, col_from); + // Set up vector of numeric categories (infected = 1, susceptible = 2, + // exposed = 3) for drawing number of moved in each category. + std::vector categories(total_infecteds, 1); + categories.insert(categories.end(), suscepts, 2); + categories.insert(categories.end(), expose, 3); + categories.insert(categories.end(), resist, 4); + + std::vector draw = draw_n_from_v(categories, total_hosts_moved, generator); + int infected_moved = std::count(draw.begin(), draw.end(), 1); + int susceptible_moved = std::count(draw.begin(), draw.end(), 2); + int exposed_moved = std::count(draw.begin(), draw.end(), 3); + int resistant_moved = std::count(draw.begin(), draw.end(), 4); + + if (exposed_moved > 0) { + std::vector exposed_draw = draw_n_from_cohorts( + exposed_, exposed_moved, row_from, col_from, generator); + int index = 0; + for (auto& raster : exposed_) { + raster(row_from, col_from) -= exposed_draw[index]; + raster(row_to, col_to) += exposed_draw[index]; + index += 1; + } + } + if (infected_moved > 0) { + std::vector 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_) { + raster(row_from, col_from) -= mortality_draw[index]; + raster(row_to, col_to) += mortality_draw[index]; + index += 1; + } + } + // Ensure that the target cell of host movement is in suitable cells. + // 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; + } + } + } + + infected_(row_from, col_from) -= infected_moved; + susceptible_(row_from, col_from) -= susceptible_moved; + total_hosts_(row_from, col_from) -= total_hosts_moved; + total_exposed_(row_from, col_from) -= exposed_moved; + resistant_(row_from, col_from) -= resistant_moved; + infected_(row_to, col_to) += infected_moved; + susceptible_(row_to, col_to) += susceptible_moved; + total_hosts_(row_to, col_to) += total_hosts_moved; + total_exposed_(row_to, col_to) += exposed_moved; + resistant_(row_to, col_to) += resistant_moved; + + // Returned total hosts actually moved is based only on the total host and no + // other checks are performed. This assumes that the counts are correct in the + // object (precodition). + return total_hosts_moved; + } + + /** + * @brief Completely remove any hosts + * + * Removes hosts completely (as opposed to moving them to another pool). + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param susceptible Number of susceptible hosts to remove. + * @param exposed Number of exposed hosts to remove by cohort. + * @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 threatment code. + * + * @note This does not remove resistant just like the original implementation in + * treatments. + */ + void completely_remove_hosts_at( + RasterIndex row, + RasterIndex col, + double susceptible, + std::vector exposed, + double infected, + const std::vector& mortality) + { + if (susceptible > 0) + susceptible_(row, col) = susceptible_(row, col) - susceptible; + + if (exposed.size() != exposed_.size()) { + throw std::invalid_argument( + "counts is not the same size as the internal list of exposed (" + + std::to_string(exposed.size()) + + " != " + std::to_string(exposed_.size()) + ") for cell (" + + std::to_string(row) + ", " + std::to_string(col) + ")"); + } + + // no simple zip in C++, falling back to indices + for (size_t i = 0; i < exposed.size(); ++i) { + exposed_[i](row, col) -= exposed[i]; + } + + // Possibly reuse in the I->S removal. + if (infected <= 0) + 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) + ")"); + } + + 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) { + 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) + + ")"); + } + infected_(row, col) -= infected; + reset_total_host(row, col); + } + + /** + * @brief Remove infected hosts and make the hosts susceptible + * + * Distribution of removed infected hosts among mortality groups is stochastic. + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param count Number of host to remove infection from + * @param generator Random number generator to provide stochasticity for mortality + * + * @note Using this method either assumes the SI model or it needs to be used + * together with remove_exposed_at() to handle SEI model. + */ + void remove_infected_at( + RasterIndex row, RasterIndex col, int count, Generator& generator) + { + // 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 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 + susceptible_(row, col) += count; + } + + /** + * @brief Make all infected hosts susceptible at the given cell + * @param row Row index of the cell + * @param col Column index of the cell + * @param generator Random number generator to provide stochasticity for mortality + */ + void remove_all_infected_at(RasterIndex row, RasterIndex col, Generator& generator) + { + auto count = this->infected_at(row, col); + this->remove_infected_at(row, col, count, generator); + } + + /** + * @brief Remove percentage of infestation/infection + * + * remove the same percentage for total exposed and remove randomly from each cohort + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param ratio Ratio of removed infection + * @param generator Random number generator to provide stochasticity for mortality + */ + void remove_infection_by_ratio_at( + RasterIndex row, RasterIndex col, double ratio, Generator& generator) + { + auto infected = this->infected_at(row, col); + int removed_infected = infected - std::lround(infected * ratio); + this->remove_infected_at(row, col, removed_infected, generator); + auto exposed = this->exposed_at(row, col); + int total_removed_exposed = exposed - std::lround(exposed * ratio); + this->remove_exposed_at(row, col, total_removed_exposed, generator); + } + + /** + * @brief Remove exposed hosts and make the hosts susceptible + * + * Distribution of removed hosts among exposed groups is stochastic. + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param count Number of host to remove infection from + * @param generator Random number generator to provide stochasticity for mortality + */ + void + remove_exposed_at(RasterIndex row, RasterIndex col, int count, Generator& generator) + { + // remove the same percentage for total exposed and remove randomly from + // each cohort + total_exposed_(row, col) -= count; + if (count > 0) { + std::vector exposed_draw = + draw_n_from_cohorts(exposed_, count, row, col, generator); + int index = 0; + for (auto& raster : exposed_) { + raster(row, col) -= exposed_draw[index]; + index += 1; + } + } + // move infested/infected host back to susceptible pool + susceptible_(row, col) += count; + } + + /** + * @brief Make hosts resistant in a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param susceptible Number of susceptible hosts to make resistant + * @param exposed Number of exposed hosts in each cohort to make resistant + * @param infected Number of infected hosts to make resistant + * @param mortality Number of infected hosts in each mortality cohort to make + * resistant + * + * @throws std::invalid_argument if counts exceed the maximum amounts possible given + * the current state (currently checked only for susceptible) + * @throws std::invalid_argument if sizes don't equal to number of cohorts + */ + void make_resistant_at( + RasterIndex row, + RasterIndex col, + int susceptible, + const std::vector& exposed, + int infected, + const std::vector& mortality) + { + int total_resistant = 0; + + if (susceptible_(row, col) < susceptible) { + throw std::invalid_argument( + "Total of newly resistant is higher than current number (" + + std::to_string(susceptible) + " > " + + std::to_string(susceptible_(row, col)) + ") for cell (" + + std::to_string(row) + ", " + std::to_string(col) + ")"); + } + + susceptible_(row, col) -= susceptible; + total_resistant += susceptible; + + if (exposed.size() != exposed_.size()) { + throw std::invalid_argument( + "exposed is not the same size as the internal list of exposed (" + + std::to_string(exposed.size()) + + " != " + std::to_string(exposed_.size()) + ") for cell (" + + std::to_string(row) + ", " + std::to_string(col) + ")"); + } + // no simple zip in C++, falling back to indices + for (size_t i = 0; i < exposed.size(); ++i) { + exposed_[i](row, col) -= exposed[i]; + total_resistant += exposed[i]; + } + infected_(row, col) -= infected; + 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 formely 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; + } + + /** + * @brief Remove resistance of host at a cell + * + * All resistant hosts are moved to the susceptible group. + * + * @param row Row index of the cell + * @param col Column index of the cell + */ + void remove_resistance_at(RasterIndex row, RasterIndex col) + { + susceptible_(row, col) += resistant_(row, col); + resistant_(row, col) = 0; + } + + /** + * @brief Apply mortality at a given cell + * + * Each cohort in exposed to mortality rate except those inside the mortality time + * lag. In indexes of the cohort vector that are in the mortality_time_lag, no + * mortality occurs. In the last year of mortality tracking, the all remaining + * tracked infected hosts are removed. In all other indexes the number of tracked + * individuals is multiplied by the mortality rate to calculate the number of hosts + * that die that time step. + * + * To be used together with step_forward_mortality(). + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param mortality_rate Percent of infected hosts that die each time period + * @param mortality_time_lag Time lag prior to mortality beginning + */ + void apply_mortality_at( + RasterIndex row, RasterIndex col, double mortality_rate, int mortality_time_lag) + { + int max_index = mortality_tracker_vector_.size() - mortality_time_lag - 1; + for (int index = 0; index <= max_index; index++) { + int mortality_in_index = 0; + if (mortality_tracker_vector_[index](row, col) > 0) { + // used to ensure that all infected hosts in the last year of + // tracking mortality + if (index == 0) { + mortality_in_index = mortality_tracker_vector_[index](row, col); + } + else { + mortality_in_index = + mortality_rate * mortality_tracker_vector_[index](row, col); + } + mortality_tracker_vector_[index](row, col) -= mortality_in_index; + died_(row, col) += mortality_in_index; + if (mortality_in_index > infected_(row, col)) { + throw std::runtime_error( + "Mortality[" + std::to_string(index) + + "] is higher than current number of infected hosts (" + + std::to_string(mortality_in_index) + " > " + + std::to_string(infected_(row, col)) + ") for cell (" + + std::to_string(row) + ", " + std::to_string(col) + ")"); + } + if (mortality_in_index > total_hosts_(row, col)) { + throw std::runtime_error( + "Mortality[" + std::to_string(index) + + "] is higher than current number of total hosts (" + + std::to_string(mortality_in_index) + " > " + + std::to_string(total_hosts_(row, col)) + ") for cell (" + + std::to_string(row) + ", " + std::to_string(col) + ")"); + } + if (infected_(row, col) > 0) { + infected_(row, col) -= mortality_in_index; + } + if (total_hosts_(row, col) > 0) { + total_hosts_(row, col) -= mortality_in_index; + } + } + } + } + + /** + * @brief Apply mortality at a given cell + * + * Uses pest-host-use table for mortality parameters. + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @see apply_mortality_at(RasterIndex, RasterIndex, double, int) + */ + void apply_mortality_at(RasterIndex row, RasterIndex col) + { + if (!pest_host_use_table_) { + throw std::invalid_argument( + "Set pest-host-use table before calling apply_mortality_at " + "or provide parameters in the function call"); + } + this->apply_mortality_at( + row, + col, + pest_host_use_table_->mortality_rate(this), + pest_host_use_table_->mortality_time_lag(this)); + } + + /** + * @brief Make a step forward in tracking mortality cohorts + * + * This makes the mortality cohorts age. + * + * To be used together with apply_mortality_at(). + */ + void step_forward_mortality() + { + rotate_left_by_one(mortality_tracker_vector_); + } + + /** + * @brief Get number of infected hosts at a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Number of infected hosts + */ + int infected_at(RasterIndex row, RasterIndex col) const + { + return infected_(row, col); + } + + /** + * @brief Get number of susceptible hosts at a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Number of susceptible hosts + */ + int susceptible_at(RasterIndex row, RasterIndex col) const + { + return susceptible_(row, col); + } + + /** + * @brief Get number of exposed hosts at a given cell + * + * The number is the total number of exposed hosts in all exposed cohorts. + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Total number of exposed hosts + * + * @note The number is taken from the internally managed total. This will be merged + * with computed_exposed_at() in the future. + */ + int exposed_at(RasterIndex row, RasterIndex col) const + { + // Future code could remove total exposed and compute that on the fly. + return total_exposed_(row, col); + } + + /** + * @brief Get number of resistant hosts at a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Number of resistant hosts + * + * @note The number is computed from all cohorts. This will be merged with + * exposed_at() in the future. + */ + int computed_exposed_at(RasterIndex row, RasterIndex col) const + { + int sum = 0; + for (const auto& raster : exposed_) + sum += raster(row, col); + return sum; + } + + /** + * @brief Get exposed hosts in each cohort at a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Vector with number of exposed hosts per cohort + */ + std::vector exposed_by_group_at(RasterIndex row, RasterIndex col) const + { + std::vector all; + all.reserve(exposed_.size()); + for (const auto& raster : exposed_) + all.push_back(raster(row, col)); + return all; + } + + /** + * @brief Get infected hosts in each mortality cohort at a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Vector with number of infected hosts per mortality cohort + */ + std::vector mortality_by_group_at(RasterIndex row, RasterIndex col) const + { + std::vector all; + all.reserve(mortality_tracker_vector_.size()); + for (const auto& raster : mortality_tracker_vector_) + all.push_back(raster(row, col)); + return all; + } + + /** + * @brief Get number of resistant hosts at a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Number of resistant hosts + */ + int resistant_at(RasterIndex row, RasterIndex col) const + { + return resistant_(row, col); + } + + /** + * @copydoc pops::HostPoolInterface::total_hosts_at() + * + * @note Computes the total host from susceptible and infected and does not consider + * exposed and resistant. + * + * @note Computes the value on the fly and does not use the raster storage for total + * host. + */ + int total_hosts_at(RasterIndex row, RasterIndex col) const override + { + return susceptible_at(row, col) + infected_at(row, col); + } + + /** + * @brief Return true if cell is outside of the rectangle covered by hosts + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return true if cell is outside, false otherwise + */ + bool is_outside(RasterIndex row, RasterIndex col) + { + return row < 0 || row >= rows_ || col < 0 || col >= cols_; + } + + /** Infect exposed hosts (E to I transition in the SEI model) + * + * Applicable to SEI model, no-operation otherwise, i.e., the state + * is left intact for SI. + * + * Before the first latency period is over, + * the E to I transition won't happen because no item in the exposed + * vector is old enough to become infected. + * + * The position of the items in the exposed vector determines their + * age, i.e., for how long the hosts are exposed. The oldest item + * is at the front and youngest at the end. + * Before the the first latency period is over, items in the front + * are still empty (unused) because no hosts were exposed for the + * given time period. + * After the first latency + * period, this needs to be true before the function is called and + * it is true after the function + * finished with the difference that after the function is called, + * the last item is empty in the sense that it does not contain any + * hosts. + * + * When the E to I transition happens, hosts from the oldest item + * in the exposed vector are moved to the infected (and mortality + * tracker). They are removed from the exposed item and this item + * is moved to the back of the vector. + * + * Like in disperse(), there is no distinction between *infected* + * and *mortality_tracker*, but different usage is expected outside + * of this function. + * + * The raster class used with the simulation class needs to support + * `.fill(value)` method for this function to work. + * + * Step is used to evaluate the latency period. + * + * @param step Step in the simulation (>=0) + */ + void step_forward(unsigned step) + { + if (model_type_ == ModelType::SusceptibleExposedInfected) { + if (step >= latency_period_) { + // Oldest item needs to be in the front + auto& oldest = exposed_.front(); + // Move hosts to infected raster + infected_ += oldest; + mortality_tracker_vector_.back() += oldest; + total_exposed_ += (oldest * (-1)); + // Reset the raster + // (hosts moved from the raster) + oldest.fill(0); + } + // Age the items and the used one to the back + // elements go one position to the left + // new oldest goes to the front + // old oldest goes to the back + rotate_left_by_one(exposed_); + } + else if (model_type_ == ModelType::SusceptibleInfected) { + // no-op + } + else { + throw std::runtime_error( + "Unknown ModelType value in Simulation::infect_exposed()"); + } + } + + /** + * @brief Get suitable cells index + * + * Suitable cells are all cells which need to be modified when all cells with host + * need to be modified. + * + * @return List of cell indices + */ + const std::vector>& suitable_cells() const + { + return suitable_cells_; + } + + std::vector& host_pools() + { + return host_pools_; + } + +private: + std::vector host_pools_ = {this}; // non-owning + /** + * @brief Reset total host value using the individual pools. + * + * This considers susceptible, exposed, infected, and resistant. + * + * This function is needed when the user needs total host as a result and everything + * is provided through individual rasters (otherwise the simplest implementation + * would just compute it on the fly always. + * + * @param row Row index of the cell + * @param col Column index of the cell + */ + void reset_total_host(RasterIndex row, RasterIndex col) + { + total_hosts_(row, col) = susceptible_(row, col) + computed_exposed_at(row, col) + + infected_(row, col) + resistant_(row, col); + } + + IntegerRaster& susceptible_; + IntegerRaster& infected_; + + std::vector& exposed_; + unsigned latency_period_{0}; + IntegerRaster& total_exposed_; + + IntegerRaster& resistant_; + + std::vector& mortality_tracker_vector_; + IntegerRaster& died_; + + IntegerRaster& total_hosts_; + const Environment& environment_; + + ModelType model_type_; + + bool dispersers_stochasticity_{false}; + double reproductive_rate_{0}; + bool establishment_stochasticity_{true}; + double deterministic_establishment_probability_{0}; + + const PestHostUseTable* pest_host_use_table_{nullptr}; + const CompetencyTable* competency_table_{nullptr}; + + RasterIndex rows_{0}; + RasterIndex cols_{0}; + + std::vector>& suitable_cells_; +}; + +} // namespace pops + +#endif // POPS_HOST_POOL_HPP diff --git a/inst/include/host_pool_interface.hpp b/inst/include/host_pool_interface.hpp new file mode 100644 index 00000000..75df0602 --- /dev/null +++ b/inst/include/host_pool_interface.hpp @@ -0,0 +1,50 @@ +/* + * PoPS model - pest or pathogen spread simulation + * + * Copyright (C) 2023 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_HOST_POOL_INTERFACE_HPP +#define POPS_HOST_POOL_INTERFACE_HPP + +namespace pops { + +/** + * Interface declaration for a host pool. + * + * Currently, the interface is providing only total number of hosts because that's the + * only function needed throughout the code. + * + * @note The usefulness of this interface will be evaluated later on, e.g., if we need + * functionally different hosts or if we have dependencies which the interface would + * address. + */ +template +class HostPoolInterface +{ +public: + virtual ~HostPoolInterface() {} + + /** + * @brief Get total number of hosts for a cell + * + * @param row Row index of the cell + * @param col Column index the cell + * + * @return Number of hosts + */ + virtual int total_hosts_at(RasterIndex row, RasterIndex col) const = 0; +}; + +} // namespace pops + +#endif // POPS_HOST_POOL_INTERFACE_HPP diff --git a/inst/include/model.hpp b/inst/include/model.hpp index 161feda6..2137057d 100644 --- a/inst/include/model.hpp +++ b/inst/include/model.hpp @@ -28,7 +28,10 @@ #include "config.hpp" #include "treatments.hpp" #include "spread_rate.hpp" -#include "simulation.hpp" +#include "host_pool.hpp" +#include "pest_pool.hpp" +#include "actions.hpp" +#include "multi_host_pool.hpp" #include "switch_kernel.hpp" #include "kernel.hpp" #include "scheduling.hpp" @@ -57,17 +60,25 @@ class Model UniformDispersalKernel uniform_kernel; DeterministicNeighborDispersalKernel natural_neighbor_kernel; DeterministicNeighborDispersalKernel anthro_neighbor_kernel; - Simulation simulation_; KernelFactory& kernel_factory_; /** - * Surrounding environment (currently used for soils only) + * Surrounding environment */ - Environment environment_; + Environment< + IntegerRaster, + FloatRaster, + RasterIndex, + RandomNumberGeneratorProvider> + environment_; /** * Optionally created soil pool */ - std::shared_ptr> soil_pool_{ - nullptr}; + std::shared_ptr>> + soil_pool_{nullptr}; unsigned last_index{0}; /** @@ -115,6 +126,19 @@ class Model } public: + using StandardSingleHostPool = HostPool< + IntegerRaster, + FloatRaster, + RasterIndex, + RandomNumberGeneratorProvider>; + using StandardMultiHostPool = MultiHostPool< + StandardSingleHostPool, + IntegerRaster, + FloatRaster, + RasterIndex, + RandomNumberGeneratorProvider>; + using StandardPestPool = PestPool; + Model( const Config& config, KernelFactory& kernel_factory = @@ -126,18 +150,8 @@ class Model uniform_kernel(config.rows, config.cols), natural_neighbor_kernel(direction_from_string(config.natural_direction)), anthro_neighbor_kernel(direction_from_string(config.anthro_direction)), - simulation_( - config.rows, - config.cols, - model_type_from_string(config.model_type), - config.latency_period_steps, - config.generate_stochasticity, - config.establishment_stochasticity, - config.movement_stochasticity), kernel_factory_(kernel_factory) - { - simulation_.set_environment(&this->environment()); - } + {} /** * @brief Run one step of the simulation. @@ -149,21 +163,20 @@ class Model * total number of hosts because movement does not support non-host * individuals. * - * No treatment can be applied when movement is active because host movement does - * not support resistant hosts. - * * *dispersers* is for internal use and for tracking dispersers creation. * The input values are ignored and the output is not the current existing * dispersers, but only the number of dispersers generated (and subsequently * used) in this step. There are no dispersers in between simulation steps. * - * @param step Step number in the simulation. + * @param step Step number in the simulation * @param[in,out] infected Infected hosts * @param[in,out] susceptible Susceptible hosts * @param[in,out] total_hosts All host individuals in the area. Is equal to - * infected + exposed + susceptible in the cell. + * infected + exposed + susceptible in the cell * @param[in,out] total_populations All host and non-host individuals in the area * @param[out] dispersers Dispersing individuals (used internally) + * @param[out] established_dispersers Dispersers originating from a given cell which + * established elsewhere * @param[in,out] total_exposed Sum of all exposed hosts (if SEI model is active) * @param[in,out] exposed Exposed hosts (if SEI model is active) * @param[in,out] mortality_tracker Mortality tracker used to generate *died*. @@ -173,14 +186,11 @@ class Model * schedule * @param[in] temperatures Vector of temperatures used to evaluate lethal * temperature - * @param[in] weather_coefficient Weather coefficient (for the current step) - * @param[in,out] treatments Treatments to be applied (also tracks use of - * treatments) + * @param[in] survival_rates Pest survival rates * @param[in,out] resistant Resistant hosts (host temporarily removed from * susceptible hosts) * @param[in,out] outside_dispersers Dispersers escaping the rasters (adds to the * vector) - * @param[in,out] spread_rate Spread rate tracker * @param[in,out] quarantine Quarantine escape tracker * @param[in] quarantine_areas Quarantine areas * @param[in] movements Table of host movements @@ -205,148 +215,200 @@ class Model IntegerRaster& died, const std::vector& temperatures, const std::vector& survival_rates, - Treatments& treatments, IntegerRaster& resistant, std::vector>& outside_dispersers, // out - SpreadRate& spread_rate, // out - QuarantineEscape& quarantine, // out + QuarantineEscapeAction& quarantine, // out const IntegerRaster& quarantine_areas, const std::vector> movements, const Network& network, std::vector>& suitable_cells) + { + StandardSingleHostPool host_pool( + model_type_from_string(config_.model_type), + susceptible, + exposed, + config_.latency_period_steps, + infected, + total_exposed, + resistant, + mortality_tracker, + died, + total_hosts, + environment_, + config_.generate_stochasticity, + config_.reproductive_rate, + config_.establishment_stochasticity, + config_.establishment_probability, + config_.rows, + config_.cols, + suitable_cells); + std::vector host_pools = {&host_pool}; + StandardMultiHostPool multi_host_pool(host_pools, config_); + StandardPestPool pest_pool{ + dispersers, established_dispersers, outside_dispersers}; + SpreadRateAction spread_rate( + multi_host_pool, + config_.rows, + config_.cols, + config_.ew_res, + config_.ns_res, + 0); + Treatments> treatments( + config_.scheduler()); + run_step( + step, + multi_host_pool, + pest_pool, + dispersers, + total_populations, + treatments, + temperatures, + survival_rates, + spread_rate, + quarantine, + quarantine_areas, + movements, + network); + } + + /** + * @brief Run one step of the simulation. + * + * @param step Step number in the simulation + * @param host_pool Host pool + * @param pest_pool Pest pool + * @param[out] dispersers Dispersing individuals (used directly for kernels) + * @param[in,out] total_populations All host and non-host individuals in the area + * @param[in,out] treatments Treatments to be applied (also tracks use of + * treatments) + * @param[in] temperatures Vector of temperatures used to evaluate lethal + * temperature + * @param[in] survival_rates Pest survival rates + * @param spread_rate Spread rate action + * @param[in,out] quarantine Quarantine escape tracker + * @param[in] quarantine_areas Quarantine areas + * @param[in] movements Table of host movements + * @param network Network (initialized or Network::null_network() if unused) + */ + void run_step( + int step, + StandardMultiHostPool& host_pool, + StandardPestPool& pest_pool, + IntegerRaster& dispersers, + IntegerRaster& total_populations, + Treatments& treatments, + const std::vector& temperatures, + const std::vector& survival_rates, + SpreadRateAction& spread_rate, + QuarantineEscapeAction& quarantine, + const IntegerRaster& quarantine_areas, + const std::vector> movements, + const Network& network) { // Soil step is the same as simulation step. if (soil_pool_) soil_pool_->next_step(step); - // removal of dispersers due to lethal temperatures if (config_.use_lethal_temperature && config_.lethal_schedule()[step]) { int lethal_step = simulation_step_to_action_step(config_.lethal_schedule(), step); - simulation_.remove( - infected, - susceptible, - temperatures[lethal_step], - config_.lethal_temperature, - suitable_cells); + this->environment().update_temperature(temperatures[lethal_step]); + RemoveByTemperature< + StandardMultiHostPool, + IntegerRaster, + FloatRaster, + RasterIndex, + RandomNumberGeneratorProvider> + remove(this->environment(), config_.lethal_temperature); + remove.action(host_pool, generator_provider_); } // 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, - generator_provider_); + SurvivalRateAction + survival(survival_rates[survival_step]); + survival.action(host_pool, generator_provider_); } // actual spread if (config_.spread_schedule()[step]) { - simulation_.generate( - dispersers, - established_dispersers, - infected, - config_.weather, - config_.reproductive_rate, - suitable_cells, - generator_provider_); - auto dispersal_kernel = kernel_factory_(config_, dispersers, network); auto overpopulation_kernel = create_overpopulation_movement_kernel(dispersers, network); + SpreadAction< + StandardMultiHostPool, + StandardPestPool, + IntegerRaster, + FloatRaster, + RasterIndex, + decltype(dispersal_kernel), + RandomNumberGeneratorProvider> + spread_action{dispersal_kernel}; - simulation_.disperse_and_infect( - step, - dispersers, - established_dispersers, - susceptible, - exposed, - infected, - mortality_tracker.back(), - total_populations, - total_exposed, - outside_dispersers, - config_.weather, - dispersal_kernel, - suitable_cells, - config_.establishment_probability, - generator_provider_); + environment_.set_total_population(&total_populations); + // Soils are activated by an independent function call for model, but spread + // action is temporary, so it is activated for every step. + if (this->soil_pool_) { + spread_action.activate_soils( + soil_pool_, config_.dispersers_to_soils_percentage); + } + spread_action.action(host_pool, pest_pool, generator_provider_); + host_pool.step_forward(step); if (config_.use_overpopulation_movements) { - simulation_.move_overpopulated_pests( - susceptible, - infected, - total_hosts, - outside_dispersers, - overpopulation_kernel, - suitable_cells, - config_.overpopulation_percentage, - config_.leaving_percentage, - generator_provider_); + MoveOverpopulatedPests< + StandardMultiHostPool, + StandardPestPool, + IntegerRaster, + FloatRaster, + RasterIndex, + decltype(overpopulation_kernel)> + move_pest{ + overpopulation_kernel, + config_.overpopulation_percentage, + config_.leaving_percentage, + config_.rows, + config_.cols}; + move_pest.action(host_pool, pest_pool, generator_provider_); } if (config_.use_movements) { - // to do fix movements to use entire mortality tracker - last_index = simulation_.movement( - infected, - susceptible, - mortality_tracker, - exposed, - resistant, - total_hosts, - total_exposed, - step, - last_index, - movements, - config_.movement_schedule, - suitable_cells, - generator_provider_); + HostMovement< + StandardMultiHostPool, + IntegerRaster, + FloatRaster, + RasterIndex> + host_movement{ + static_cast( + step), // Step is int, but indexing happens with unsigned. + last_index, + movements, + config_.movement_schedule}; + last_index = host_movement.action(host_pool, generator_provider_); } } // treatments if (config_.use_treatments) { - bool managed = treatments.manage( - step, - infected, - exposed, - susceptible, - resistant, - total_hosts, - suitable_cells); - if (managed && config_.use_mortality) { - // treatments apply to all mortality tracker cohorts - for (auto& raster : mortality_tracker) { - treatments.manage_mortality(step, raster, suitable_cells); - } + for (auto& host : host_pool.host_pools()) { + treatments.manage(step, *host); } } if (config_.use_mortality && config_.mortality_schedule()[step]) { // expectation is that mortality tracker is of length (1/mortality_rate // + mortality_time_lag). // TODO: died.zero(); should be done by the caller if needed, document! - simulation_.mortality( - infected, - total_hosts, - config_.mortality_rate, - config_.mortality_time_lag, - died, - mortality_tracker, - suitable_cells); + Mortality mortality; + mortality.action(host_pool); } // compute spread rate if (config_.use_spreadrates && config_.spread_rate_schedule()[step]) { unsigned rates_step = simulation_step_to_action_step(config_.spread_rate_schedule(), step); - spread_rate.compute_step_spread_rate(infected, rates_step, suitable_cells); + spread_rate.action(host_pool, rates_step); } // compute quarantine escape if (config_.use_quarantine && config_.quarantine_schedule()[step]) { unsigned action_step = simulation_step_to_action_step(config_.quarantine_schedule(), step); - quarantine.infection_escape_quarantine( - infected, quarantine_areas, action_step, suitable_cells); + quarantine.action(host_pool, quarantine_areas, action_step); } } @@ -363,7 +425,12 @@ class Model * @brief Get surrounding environment * @return Environment object by reference */ - Environment& environment() + Environment< + IntegerRaster, + FloatRaster, + RasterIndex, + RandomNumberGeneratorProvider>& + environment() { return environment_; } @@ -379,13 +446,15 @@ class Model void activate_soils(std::vector& rasters) { // The soil pool is created again for every new activation. - this->soil_pool_.reset(new SoilPool( + this->soil_pool_.reset(new SoilPool< + IntegerRaster, + FloatRaster, + RasterIndex, + RandomNumberGeneratorProvider>( rasters, this->environment_, config_.generate_stochasticity, config_.establishment_stochasticity)); - this->simulation_.activate_soils( - this->soil_pool_, config_.dispersers_to_soils_percentage); } }; diff --git a/inst/include/model_type.hpp b/inst/include/model_type.hpp new file mode 100644 index 00000000..1ac38182 --- /dev/null +++ b/inst/include/model_type.hpp @@ -0,0 +1,63 @@ +/* + * PoPS model - pest or pathogen spread simulation + * + * Copyright (C) 2023 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_MODEL_TYPE_HPP +#define POPS_MODEL_TYPE_HPP + +namespace pops { + +/** The type of a epidemiological model (SI or SEI) + */ +enum class ModelType +{ + SusceptibleInfected, ///< SI (susceptible - infected) + SusceptibleExposedInfected ///< SEI (susceptible - exposed - infected) +}; + +/*! Get a corresponding enum value for a string which is a model type name. + * + * Throws an std::invalid_argument exception if the value was not + * found or is not supported (which is the same thing). + */ +inline ModelType model_type_from_string(const std::string& text) +{ + if (text == "SI" || text == "SusceptibleInfected" || text == "susceptible-infected" + || text == "susceptible_infected") + return ModelType::SusceptibleInfected; + else if ( + text == "SEI" || text == "SusceptibleExposedInfected" + || text == "susceptible-exposed-infected" + || text == "susceptible_exposed_infected") + return ModelType::SusceptibleExposedInfected; + else + throw std::invalid_argument( + "model_type_from_string: Invalid" + " value '" + + text + "' provided"); +} + +/*! Overload which allows to pass C-style string which is nullptr (NULL) + * + * @see model_type_from_string(const std::string& text) + */ +inline ModelType model_type_from_string(const char* text) +{ + // call the string version + return model_type_from_string(text ? std::string(text) : std::string()); +} + +} // namespace pops + +#endif // POPS_MODEL_TYPE_HPP diff --git a/inst/include/multi_host_pool.hpp b/inst/include/multi_host_pool.hpp new file mode 100644 index 00000000..e7bdf5aa --- /dev/null +++ b/inst/include/multi_host_pool.hpp @@ -0,0 +1,461 @@ +/* + * PoPS model - pest or pathogen spread simulation + * + * Copyright (C) 2023 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_MULTI_HOST_POOL_HPP +#define POPS_MULTI_HOST_POOL_HPP + +#include +#include +#include + +#include "competency_table.hpp" +#include "config.hpp" +#include "pest_host_use_table.hpp" +#include "utils.hpp" + +namespace pops { + +/** + * Host pool for multiple hosts + * + * Keeps the same interface as (single) HostPool given that HostPool has methods to + * behave in a mutli-host way if needed. This allows most external operations such as + * actions to work without distinguishing single- and multi-host pools. + */ +template< + typename HostPool, + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex, + typename GeneratorProvider> +class MultiHostPool +{ +public: + /** + * Standard random number generator to be passed directly to the methods. + */ + using Generator = typename GeneratorProvider::Generator; + + /** + * @brief Create MultiHostPool object with given host pools and configuration + * + * @param host_pools List of host pools to use + * @param config Configuration to use + */ + MultiHostPool(const std::vector& host_pools, const Config& config) + : host_pools_(host_pools), config_(config) + {} + + /** + * @brief Set pest-host-use table for all hosts + * + * The existing object will be used (not copy is performed). + * + * @param table Reference to the table object + */ + void set_pest_host_use_table(const PestHostUseTable& table) + { + for (auto& host_pool : host_pools_) { + host_pool->set_pest_host_use_table(table); + } + } + + /** + * @brief Set competency table for all hosts + * + * The existing object will be used (not copy is performed). + * + * @param table Reference to the table object + */ + void set_competency_table(const CompetencyTable& table) + { + for (auto& host_pool : host_pools_) { + host_pool->set_competency_table(table); + } + } + + /** + * @brief Get suitable cells spatial index + * + * Always uses the index only from the first single-host pool + * assuming that the indices are the same. + * Assumes that at least one single-host pool was added. + * + * @return Const reference to the index + */ + const std::vector>& suitable_cells() const + { + return host_pools_.at(0)->suitable_cells(); + } + + /** + * @brief Check whether the cell is outside of the raster extent + * + * Always uses the only the first single-host pool for the check + * assuming that the extents are the same. + * Assumes that at least one single-host pool was added. + * + * @param row Row index of the cell + * @param col Column index of the cell + * @return true if outside of the raster, false if inside + */ + bool is_outside(RasterIndex row, RasterIndex col) + { + return host_pools_.at(0)->is_outside(row, col); + } + + /** + * @brief Step all hosts forward + * @param step Step in the simulation (>=0) + * + * @see HostPool::step_forward() + */ + void step_forward(unsigned step) + { + for (auto& host_pool : host_pools_) { + host_pool->step_forward(step); + } + } + + /** + * @copydoc HostPool::remove_all_infected_at() + */ + void remove_all_infected_at(RasterIndex row, RasterIndex col, Generator& generator) + { + for (auto& host_pool : host_pools_) { + host_pool->remove_all_infected_at(row, col, generator); + } + } + + /** + * @copydoc HostPool::remove_infection_by_ratio_at() + */ + void remove_infection_by_ratio_at( + RasterIndex row, RasterIndex col, double ratio, Generator& generator) + { + for (auto& host_pool : host_pools_) { + host_pool->remove_infection_by_ratio_at(row, col, ratio, generator); + } + } + + /** + * @copydoc HostPool::apply_mortality_at(RasterIndex, RasterIndex, double, int) + */ + void apply_mortality_at( + RasterIndex row, RasterIndex col, double mortality_rate, int mortality_time_lag) + { + for (auto& host_pool : host_pools_) { + host_pool->apply_mortality_at(row, col, mortality_rate, mortality_time_lag); + } + } + + /** + * @copydoc HostPool::apply_mortality_at(RasterIndex, RasterIndex) + */ + void apply_mortality_at(RasterIndex row, RasterIndex col) + { + for (auto& host_pool : host_pools_) { + host_pool->apply_mortality_at(row, col); + } + } + + /** + * @copydoc HostPool::step_forward_mortality() + */ + void step_forward_mortality() + { + for (auto& host_pool : host_pools_) { + host_pool->step_forward_mortality(); + } + } + + /** + * @brief Total number of all hosts over all host pools + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Number of hosts + */ + int total_hosts_at(RasterIndex row, RasterIndex col) const + { + int sum = 0; + for (auto& host_pool : host_pools_) { + sum += host_pool->total_hosts_at(row, col); + } + return sum; + } + + /** + * @brief Get number of infected hosts over all host pools at a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Number of infected hosts + */ + int infected_at(RasterIndex row, RasterIndex col) const + { + int infected = 0; + for (auto& host_pool : host_pools_) { + infected += host_pool->infected_at(row, col); + } + return infected; + } + + /** + * @copydoc HostPool::dispersers_from() + */ + int dispersers_from(RasterIndex row, RasterIndex col, Generator& generator) const + { + int sum{0}; + for (auto& host_pool : host_pools_) { + sum += host_pool->dispersers_from(row, col, generator); + } + return sum; + } + + /** + * @brief Move pests from a cell + * + * This is a multi-host version of single host pests_from method. + * It randomly distributes the number to un-infect between multiple + * hosts. + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param count Number of pests requested to move from the cell + * @param generator Random number generator + * + * @return Number of pests actually moved from the cell + * + * @note For consitency with the previous implementation, this does not modify + * mortality cohorts nor touches the exposed cohorts. + */ + int pests_from(RasterIndex row, RasterIndex col, int count, Generator& generator) + { + std::vector infected; + int index = 0; + for (auto& host_pool : host_pools_) { + infected.insert(infected.end(), host_pool->infected_at(row, col), index); + index++; + } + + index = 0; + int collect_count = 0; + std::vector draw = draw_n_from_v(infected, count, generator); + for (auto& host_pool : host_pools_) { + count = std::count(draw.begin(), draw.end(), index); + collect_count += host_pool->pests_from(row, col, count, generator); + index++; + } + + return collect_count; + } + /** + * @brief Move pests to a cell + * + * This is a multi-host version of single host pests_to method. + * It randomly distributes the number to infect between multiple + * hosts. + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param count Number of pests requested to move to the cell + * @param generator Random number generator + * + * @return Number of accepted pests + * + * @note For consistency with the previous implementation, this does not make hosts + * exposed in the SEI model. Instead, the hosts are infected right away. This may + * become a feature in the future. + * + * @note For consistency with the previous implementation, this does not modify the + * mortality cohorts. This wil need to be fixed in the future. + * + * @note This may be merged with add_disperser_at() in the future. + */ + int pests_to(RasterIndex row, RasterIndex col, int count, Generator& generator) + { + std::vector susceptible; + int index = 0; + for (auto& host_pool : host_pools_) { + susceptible.insert( + susceptible.end(), host_pool->susceptible_at(row, col), index); + index++; + } + index = 0; + int collect_count = 0; + std::vector draw = draw_n_from_v(susceptible, count, generator); + for (auto& host_pool : host_pools_) { + count = std::count(draw.begin(), draw.end(), index); + collect_count += host_pool->pests_to(row, col, count, generator); + index++; + } + + return collect_count; + } + + /** + * @brief Move a disperser to a cell with establishment test + * + * Config::arrival_behavior() is used to determine if the dispersers should be + * landing in a cell (`"land"`) or infecting a specific host directly (`"infect"`). + * Landing performs the establishment test on a combined probability from all hosts + * in a cell and then uses one host to add the new disperser to. Infecting picks a + * host and lets the host accept or reject the disperser based on its own + * establishment test. + * + * Uses static function HostPool::can_disperser_establish() from HostPool class to + * do the test for multiple hosts for landing, so any host pool class needs to + * implement that besides its disperser methods. It assumes that the configuration + * for multi-host is applicable for the combined establishment probability of + * individual hosts. + * + * For one single host pool only and host pool which implements its `disperser_to()` + * using `can_disperser_establish()` and `add_disperser_at()`, this gives identical + * results for both landing and infecting arrival behavior. + * + * @see HostPool::can_disperser_establish() + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param generator Random number generator + * + * @return 1 if established, 0 otherwise + */ + int disperser_to(RasterIndex row, RasterIndex col, Generator& generator) + { + std::vector probabilities; + double total_s_score = 0; + for (auto& host_pool : host_pools_) { + // The probability accounts for weather and, importantly, number of + // susceptible hosts, so host pool without available hosts in a given cell + // is less likely to be selected over a host pool with available hosts + // (however, it can happen in case zeros are not exactly zeros in the + // discrete distribution used later and code checks can't tell, so we need + // to account for that case later anyway). + double s_for_item = host_pool->establishment_probability_at(row, col); + // The resulting s can be 0-1. While the probabilities are used as weights + // for picking the host, so their absolute range does not matter, the total + // is used as probablity in a stochastic test. The stochastic challenge may + // be against 0.5 + 0.7 = 1.2 which which is fine as long as we are fine + // with probabilities 0.5 and 0.7 from two host translating to 100% + // probability. + probabilities.push_back(s_for_item); + total_s_score += s_for_item; + } + if (total_s_score <= 0) { + // While the score should always be >= 0, it may be == 0 if no hosts are + // present. No hosts present cause all probabilities to be zero which is not + // permissible for the host picking later and it is enough information for + // us to know there won't be any establishment. + return 0; + } + auto host = pick_host_by_probability(host_pools_, probabilities, generator); + if (config_.arrival_behavior() == "land") { + // The operations are ordered so that for single host, this gives an + // identical results to the infect behavior (influenced by usage of random + // numbers and presence of susceptible hosts). + if (host->susceptible_at(row, col) <= 0) + return 0; + bool establish = HostPool::can_disperser_establish( + total_s_score, + config_.establishment_stochasticity, + config_.establishment_probability, + generator); + if (establish) + return host->add_disperser_at(row, col); // simply increases the counts + return 0; + } + else if (config_.arrival_behavior() == "infect") { + return host->disperser_to(row, col, generator); // with establishment test + } + throw std::invalid_argument( + "Unknown value for arrival_behavior: " + config_.arrival_behavior()); + } + /** + * @brief Move hosts from a cell to a cell (multi-host) + * + * Currenty just works for first host. + * + * @param row_from Row index of the source cell + * @param col_from Column index of the source cell + * @param row_to Row index of the target cell + * @param col_to Column index of the target cell + * @param count Number of pests to move + * @param generator Random number generator for distribution of host among pools + * + * @return Number of pests actually moved between cells + * + * @note Mortality is not supported. + */ + int move_hosts_from_to( + RasterIndex row_from, + RasterIndex col_from, + RasterIndex row_to, + RasterIndex col_to, + int count, + Generator& generator) + { + return host_pools_[0]->move_hosts_from_to( + row_from, col_from, row_to, col_to, count, generator); + } + + /** + * @brief Get list of host pools + * @return Reference to host pools + */ + const std::vector& host_pools() + { + return host_pools_; + } + +private: + /** + * @brief Pick a host given a probability for each host + * + * The probabilities are used as weights so they don't need to be 0-1 nor they need + * to add up to 1. However, their sum needs to be >0. + * + * @param hosts List of pointers to host pools + * @param probabilities Probability values for each host pool + * @param generator Random number generator + * + * @return Pointer to selected host pool + */ + static HostPool* pick_host_by_probability( + std::vector& hosts, + const std::vector& probabilities, + Generator& generator) + { + std::discrete_distribution distribution{ + probabilities.begin(), probabilities.end()}; + return hosts.at(distribution(generator)); + } + + /** + * List of non-owning pointers to individual host pools. + */ + std::vector host_pools_; + /** + * Reference to configuration + */ + const Config& config_; +}; + +} // namespace pops + +#endif // POPS_MULTI_HOST_POOL_HPP diff --git a/inst/include/natural_kernel.hpp b/inst/include/natural_kernel.hpp index 4da5f013..880eab6c 100644 --- a/inst/include/natural_kernel.hpp +++ b/inst/include/natural_kernel.hpp @@ -33,7 +33,9 @@ namespace pops { * * Kernel parameters are taken from the configuration. * + * @param config Configuration for the kernel * @param dispersers The disperser raster (reference, for deterministic kernel) + * * @return Created kernel */ template diff --git a/inst/include/network.hpp b/inst/include/network.hpp index 8c2132a8..56caf671 100644 --- a/inst/include/network.hpp +++ b/inst/include/network.hpp @@ -860,7 +860,6 @@ class Network * @brief Read segments from a stream. * * @param stream Input stream with segments - * @param node_ids Container with node IDs to use */ template void load_segments(InputStream& stream) @@ -884,11 +883,11 @@ class Network "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); - // } + 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); + } Segment segment; if (has_probability) { diff --git a/inst/include/pest_host_use_table.hpp b/inst/include/pest_host_use_table.hpp new file mode 100644 index 00000000..28c49e7e --- /dev/null +++ b/inst/include/pest_host_use_table.hpp @@ -0,0 +1,81 @@ +/* + * PoPS model - Pest-host-use table for hosts and pest + * + * Copyright (C) 2023 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_PEST_HOST_USE_TABLE_HPP +#define POPS_PEST_HOST_USE_TABLE_HPP + +#include +#include +#include + +#include "config.hpp" + +namespace pops { + +template +class PestHostUseTable +{ +public: + using Environment = typename HostPool::Environment; + + PestHostUseTable(const Environment& environment) : environment_(environment) {} + PestHostUseTable(const Config& config, const Environment& environment) + : environment_(environment) + { + for (const auto& row : config.pest_host_use_table_data()) { + susceptibilities_.push_back(row.susceptibility); + mortality_rates_.push_back(row.mortality_rate); + mortality_time_lags_.push_back(row.mortality_time_lag); + } + } + + void + add_host_info(double susceptibility, double mortality_rate, int mortality_time_lag) + { + susceptibilities_.push_back(susceptibility); + mortality_rates_.push_back(mortality_rate); + mortality_time_lags_.push_back(mortality_time_lag); + } + + double susceptibility(const HostPool* host) const + { + // This is using index because the environment is part of competency table, + // otherwise a map which would use pointer to host would work here, too. + auto host_index = environment_.host_index(host); + return susceptibilities_.at(host_index); + } + + double mortality_rate(const HostPool* host) const + { + auto host_index = environment_.host_index(host); + return mortality_rates_.at(host_index); + } + + double mortality_time_lag(const HostPool* host) const + { + auto host_index = environment_.host_index(host); + return mortality_time_lags_.at(host_index); + } + +private: + std::vector susceptibilities_; + std::vector mortality_rates_; + std::vector mortality_time_lags_; + const Environment& environment_; +}; + +} // namespace pops + +#endif // POPS_PEST_HOST_USE_TABLE_HPP diff --git a/inst/include/pest_pool.hpp b/inst/include/pest_pool.hpp new file mode 100644 index 00000000..e91feb87 --- /dev/null +++ b/inst/include/pest_pool.hpp @@ -0,0 +1,134 @@ +/* + * PoPS model - pest or pathogen spread simulation + * + * Copyright (C) 2023 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_PEST_POOL_HPP +#define POPS_PEST_POOL_HPP + +#include +#include + +namespace pops { + +/** + * A class to track pests and manage pest which are not on hosts. + * + * This class can be viewed as a pests pool where pests can be + * when they are not on hosts. However, it also tracks other information + * about pests which is not captured by hosts. + */ +template +class PestPool +{ +public: + /** + * @brief Create an object with linked data + * @param dispersers Generated dispersers + * @param established_dispersers Established dispersers from a cell + * @param outside_dispersers Dispersers which left the area + * + * The object provides method-based access to the underlying rasters. + * + * The object does not copy or take ownership of the objects passed in the + * constructor. + */ + PestPool( + IntegerRaster& dispersers, + IntegerRaster& established_dispersers, + std::vector>& outside_dispersers) + : dispersers_(dispersers), + established_dispersers_(established_dispersers), + outside_dispersers_(outside_dispersers) + {} + /** + * @brief Set number of dispersers + * @param row Row number + * @param col Column number + * @param value The new value + */ + void set_dispersers_at(RasterIndex row, RasterIndex col, int value) + { + dispersers_(row, col) = value; + } + /** + * @brief Return number of dispersers + * @param row Row number + * @param col Column number + * @return The current value + */ + int dispersers_at(RasterIndex row, RasterIndex col) const + { + return dispersers_(row, col); + } + /** + * @brief Set number of established dispersers + * + * Established 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 + */ + void remove_established_dispersers_at(RasterIndex row, RasterIndex col, int count) + { + established_dispersers_(row, col) -= count; + } + /** + * @brief Add a disperser which left the study area + * @param row Row number + * @param col Column number + */ + void add_outside_disperser_at(RasterIndex row, RasterIndex col) + { + outside_dispersers_.emplace_back(row, col); + } + /** + * @brief Add a dispersers which left the study area + * @param row Row number + * @param col Column number + * @param count Number of dispersers which left + */ + void add_outside_dispersers_at(RasterIndex row, RasterIndex col, int count) + { + outside_dispersers_.reserve(outside_dispersers_.size() + count); + for (int pest = 0; pest < count; ++pest) + outside_dispersers_.emplace_back(row, col); + } + +private: + /// Generated dispersers + IntegerRaster& dispersers_; + /// Origins of established dispersers + IntegerRaster& established_dispersers_; + /// Destinations of dispersers which left + std::vector>& outside_dispersers_; +}; + +} // namespace pops + +#endif // POPS_PEST_POOL_HPP diff --git a/inst/include/quarantine.hpp b/inst/include/quarantine.hpp index 8db97386..71da68b1 100644 --- a/inst/include/quarantine.hpp +++ b/inst/include/quarantine.hpp @@ -76,7 +76,7 @@ Directions directions_from_string(const std::string& text, char delimiter = ',') * Class storing and computing quarantine escape metrics for one simulation. */ template -class QuarantineEscape +class QuarantineEscapeAction { private: RasterIndex width_; @@ -167,7 +167,7 @@ class QuarantineEscape } public: - QuarantineEscape( + QuarantineEscapeAction( const IntegerRaster& quarantine_areas, double ew_res, double ns_res, @@ -188,26 +188,24 @@ class QuarantineEscape quarantine_boundary(quarantine_areas); } - QuarantineEscape() = delete; + QuarantineEscapeAction() = delete; /** * Computes whether infection in certain step escaped from quarantine areas * and if not, computes and saves minimum distance and direction to quarantine areas * for the specified step. Aggregates over all quarantine areas. */ - void infection_escape_quarantine( - const IntegerRaster& infected, - const IntegerRaster& quarantine_areas, - unsigned step, - const std::vector>& suitable_cells) + template + void + action(const Hosts& hosts, const IntegerRaster& quarantine_areas, unsigned step) { DistDir min_dist_dir = std::make_tuple(std::numeric_limits::max(), Direction::None); - for (auto indices : suitable_cells) { + for (auto indices : hosts.suitable_cells()) { int i = indices[0]; int j = indices[1]; - if (!infected(i, j)) + if (!hosts.infected_at(i, j)) continue; int area = quarantine_areas(i, j); if (area == 0) { @@ -270,7 +268,8 @@ class QuarantineEscape */ template double quarantine_escape_probability( - const std::vector> escape_infos, unsigned step) + const std::vector> escape_infos, + unsigned step) { bool escape; DistDir distdir; @@ -290,7 +289,8 @@ double quarantine_escape_probability( */ template std::vector distance_direction_to_quarantine( - const std::vector> escape_infos, unsigned step) + const std::vector> escape_infos, + unsigned step) { bool escape; DistDir distdir; @@ -312,7 +312,8 @@ std::vector distance_direction_to_quarantine( */ template std::string write_quarantine_escape( - const std::vector> escape_infos, unsigned num_steps) + const std::vector> escape_infos, + unsigned num_steps) { std::stringstream ss; ss << std::setprecision(1) << std::fixed; diff --git a/inst/include/simulation.hpp b/inst/include/simulation.hpp index e9173b14..189b56fa 100644 --- a/inst/include/simulation.hpp +++ b/inst/include/simulation.hpp @@ -1,12 +1,9 @@ /* * PoPS model - pest or pathogen spread simulation * - * Copyright (C) 2015-2022 by the authors. + * Copyright (C) 2015-2023 by the authors. * - * Authors: Zexi Chen (zchen22 ncsu edu) - * Vaclav Petras (wenzeslaus gmail com) - * Anna Petrasova (kratochanna gmail com) - * Chris Jones (cjones1688 gmail com) + * Authors: Vaclav Petras (wenzeslaus gmail com) * * The code contained herein is licensed under the GNU General Public * License. You may obtain a copy of the GNU General Public License @@ -27,51 +24,15 @@ #include #include +#include "actions.hpp" #include "utils.hpp" #include "soils.hpp" +#include "model_type.hpp" +#include "pest_pool.hpp" +#include "host_pool.hpp" namespace pops { -/** The type of a epidemiological model (SI or SEI) - */ -enum class ModelType -{ - SusceptibleInfected, ///< SI (susceptible - infected) - SusceptibleExposedInfected ///< SEI (susceptible - exposed - infected) -}; - -/*! Get a corresponding enum value for a string which is a model type name. - * - * Throws an std::invalid_argument exception if the value was not - * found or is not supported (which is the same thing). - */ -inline ModelType model_type_from_string(const std::string& text) -{ - if (text == "SI" || text == "SusceptibleInfected" || text == "susceptible-infected" - || text == "susceptible_infected") - return ModelType::SusceptibleInfected; - else if ( - text == "SEI" || text == "SusceptibleExposedInfected" - || text == "susceptible-exposed-infected" - || text == "susceptible_exposed_infected") - return ModelType::SusceptibleExposedInfected; - else - throw std::invalid_argument( - "model_type_from_string: Invalid" - " value '" - + text + "' provided"); -} - -/*! Overload which allows to pass C-style string which is nullptr (NULL) - * - * @see model_type_from_string(const std::string& text) - */ -inline ModelType model_type_from_string(const char* text) -{ - // call the string version - return model_type_from_string(text ? std::string(text) : std::string()); -} - /*! The main class to control the spread simulation. * * The Simulation class handles the mechanics of the model, but the @@ -99,7 +60,11 @@ inline ModelType model_type_from_string(const char* text) * type are performed and a signed type might be required in the future. * A default is provided, but it can be changed in the future. */ -template +template< + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex = int, + typename Generator = DefaultSingleGeneratorProvider> class Simulation { private: @@ -111,90 +76,27 @@ class Simulation ModelType model_type_; unsigned latency_period_; /// Non-owning pointer to environment for weather - const Environment* environment_{nullptr}; + // can be const in simulation, except for need to add host now + Environment* environment_{ + nullptr}; /** * Optional soil pool */ - std::shared_ptr> soil_pool_{ - nullptr}; + std::shared_ptr> + soil_pool_{nullptr}; /** * Percentage (0-1 ratio) of disperers to be send to soil */ double to_soil_percentage_{0}; - /** - * @brief Move (add) disperser to a cell in the host pool - * - * Processes event when a disperser lands in a cell potentially establishing on a - * host. The disperser may or may not establish a based on host availability, - * weather, establishment probability, and stochasticity. - * - * Any new dispersers targeting host in the host pool should be added using this - * function. - * - * @param row Row number of the target cell - * @param col Column number of the target cell - * @param susceptible Raster of susceptible hosts - * @param exposed_or_infected Raster of exposed or infected hosts - * @param mortality_tracker Raster tracking hosts for mortality - * @param total_populations Raster of all individuals (hosts and non-hosts) - * @param total_exposed Raster tracking all exposed hosts - * @param weather Whether to use weather - * @param establishment_probability Fixed probability disperser establishment - * - * @return true if disperser has established in the cell, false otherwise - * - * @throw std::runtime_error if model type is unsupported (i.e., not SI or SEI) - */ - template - int disperser_to( - RasterIndex row, - RasterIndex col, - IntegerRaster& susceptible, - IntegerRaster& exposed_or_infected, - IntegerRaster& mortality_tracker, - const IntegerRaster& total_populations, - IntegerRaster& total_exposed, - bool weather, - double establishment_probability, - Generator& generator) - { - std::uniform_real_distribution distribution_uniform(0.0, 1.0); - if (susceptible(row, col) > 0) { - double probability_of_establishment = - (double)(susceptible(row, col)) / total_populations(row, col); - double establishment_tester = 1 - establishment_probability; - if (establishment_stochasticity_) - establishment_tester = distribution_uniform(generator); - - if (weather) - probability_of_establishment *= - environment()->weather_coefficient_at(row, col); - if (establishment_tester < probability_of_establishment) { - exposed_or_infected(row, col) += 1; - susceptible(row, col) -= 1; - if (model_type_ == ModelType::SusceptibleInfected) { - mortality_tracker(row, col) += 1; - } - else if (model_type_ == ModelType::SusceptibleExposedInfected) { - total_exposed(row, col) += 1; - } - else { - throw std::runtime_error( - "Unknown ModelType value in " - "Simulation::disperse()"); - } - return true; - } - } - return false; - } - public: - /** Creates simulation object and seeds the internal random number generator. - * - * The same random number generator is used throughout the simulation - * and is seeded once at the beginning. + // Host pool has the provider from model, but in test, it gets plain engine. + using StandardHostPool = + HostPool; + using StandardPestPool = PestPool; + + /** + * Creates simulation object with the values which are fixed during the simulation. * * The number or rows and columns needs to be the same as the size * of rasters used with the Simulation object @@ -202,7 +104,6 @@ class Simulation * * @param model_type Type of the model (SI or SEI) * @param latency_period Length of the latency period in steps - * @param random_seed Number to seed the random number generator * @param rows Number of rows * @param cols Number of columns * @param dispersers_stochasticity Enable stochasticity in generating of dispersers @@ -234,49 +135,82 @@ class Simulation * * The simulation object does not take ownership of the environment. */ + // parameter and attribute should be const void set_environment( - const Environment* environment) + Environment* environment) { this->environment_ = environment; } /** * @brief Get environment used in the simulation + * + * @param allow_empty if true, empty (non-functional) environment is returned * @return Const pointer to the environment * @throw std::logic_error when environment is not set */ - const Environment* environment() + const Environment* + environment(bool allow_empty = false) { - if (!this->environment_) + static Environment empty; + if (!this->environment_) { + if (allow_empty) + return ∅ throw std::logic_error("Environment used in Simulation, but not provided"); + } return this->environment_; } - /** removes infected based on min or max temperature tolerance + /** Remove infected based on temperature tolerance * * @param infected Currently infected hosts * @param susceptible Currently susceptible hosts - * @param temperature Spatially explicit temperature + * @param exposed Currently exposed hosts + * @param total_exposed Total exposed in all exposed cohorts + * @param mortality_tracker_vector Mortality tracker * @param lethal_temperature temperature at which lethal conditions occur * @param suitable_cells used to run model only where host are known to occur + * @param generator Provider of random number generators */ + template void remove( IntegerRaster& infected, IntegerRaster& susceptible, - const FloatRaster& temperature, + std::vector& exposed, + IntegerRaster& total_exposed, + std::vector& mortality_tracker_vector, double lethal_temperature, - const std::vector>& suitable_cells) + std::vector>& suitable_cells, + GeneratorProvider& generator) { - for (auto indices : suitable_cells) { - int i = indices[0]; - int j = indices[1]; - if (temperature(i, j) < lethal_temperature) { - // move infested/infected host back to susceptible pool - susceptible(i, j) += infected(i, j); - // remove all infestation/infection in the infected class - infected(i, j) = 0; - } - } + IntegerRaster empty; + StandardHostPool hosts( + model_type_, + susceptible, + exposed, + 0, + infected, + total_exposed, + empty, + mortality_tracker_vector, + empty, + empty, + *environment(true), + false, + 0, + false, + 0, + 0, + 0, + suitable_cells); + RemoveByTemperature< + StandardHostPool, + IntegerRaster, + FloatRaster, + RasterIndex, + GeneratorProvider> + remove(*environment(false), lethal_temperature); + remove.action(hosts, generator); } /** Removes percentage of exposed and infected @@ -288,6 +222,7 @@ class Simulation * @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 + * @param generator Provider of random number generators */ template void remove_percentage( @@ -297,57 +232,32 @@ class Simulation std::vector& exposed, IntegerRaster& total_exposed, const FloatRaster& survival_rate, - const std::vector>& suitable_cells, + std::vector>& suitable_cells, GeneratorProvider& generator) { - 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 mortality_draw = draw_n_from_cohorts( - mortality_tracker_vector, - removed_infected, - i, - j, - generator.survival_rate()); - 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 exposed_draw = draw_n_from_cohorts( - exposed, - total_removed_exposed, - i, - j, - generator.survival_rate()); - 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; - } - } + IntegerRaster empty; + StandardHostPool hosts( + model_type_, + susceptible, + exposed, + 0, + infected, + total_exposed, + empty, + mortality_tracker_vector, + empty, + empty, + *environment(true), + false, + 0, + false, + 0, + 0, + 0, + suitable_cells); + SurvivalRateAction survival( + survival_rate); + survival.action(hosts, generator); } /** kills infected hosts based on mortality rate and timing. In the last year @@ -374,39 +284,32 @@ class Simulation int mortality_time_lag, IntegerRaster& died, std::vector& mortality_tracker_vector, - const std::vector>& suitable_cells) + std::vector>& suitable_cells) { - - int max_index = mortality_tracker_vector.size() - mortality_time_lag - 1; - - for (auto indices : suitable_cells) { - int i = indices[0]; - int j = indices[1]; - for (int index = 0; index <= max_index; index++) { - int mortality_in_index = 0; - if (mortality_tracker_vector[index](i, j) > 0) { - // used to ensure that all infected hosts in the last year of - // tracking mortality - if (index == 0) { - mortality_in_index = mortality_tracker_vector[index](i, j); - } - else { - mortality_in_index = - mortality_rate * mortality_tracker_vector[index](i, j); - } - mortality_tracker_vector[index](i, j) -= mortality_in_index; - died(i, j) += mortality_in_index; - if (infected(i, j) > 0) { - infected(i, j) -= mortality_in_index; - } - if (total_hosts(i, j) > 0) { - total_hosts(i, j) -= mortality_in_index; - } - } - } - } - - rotate_left_by_one(mortality_tracker_vector); + IntegerRaster empty; + std::vector empty_vector; + StandardHostPool hosts{ + model_type_, + empty, + empty_vector, + 0, + infected, + empty, + empty, + mortality_tracker_vector, + died, + total_hosts, + *environment(true), + false, + 0, + false, + 0, + 0, + 0, + suitable_cells}; + Mortality mortality( + mortality_rate, mortality_time_lag); + mortality.action(hosts); } /** Moves hosts from one location to another @@ -430,10 +333,10 @@ class Simulation * @param movement_schedule a vector matching movements with the step at which the * movement from movements are applied * @param suitable_cells List of indices of cells with hosts + * @param generator Provider of random number generators * * @note Mortality and non-host individuals are not supported in movements. */ - template unsigned movement( IntegerRaster& infected, IntegerRaster& susceptible, @@ -447,113 +350,45 @@ class Simulation const std::vector>& movements, std::vector movement_schedule, std::vector>& suitable_cells, - GeneratorProvider& generator) + Generator& generator) { - for (unsigned i = last_index; i < movements.size(); i++) { - auto moved = movements[i]; - unsigned move_schedule = movement_schedule[i]; - if (move_schedule != step) { - return i; - } - int infected_moved = 0; - int exposed_moved = 0; - int susceptible_moved = 0; - int resistant_moved = 0; - int total_hosts_moved = 0; - int row_from = moved[0]; - int col_from = moved[1]; - int row_to = moved[2]; - int col_to = moved[3]; - int hosts = moved[4]; - std::vector exposed_categories; - std::vector mortality_categories; - if (hosts > total_hosts(row_from, col_from)) { - total_hosts_moved = total_hosts(row_from, col_from); - } - else { - total_hosts_moved = hosts; - } - auto total_infecteds = infected(row_from, col_from); - auto suscepts = susceptible(row_from, col_from); - auto expose = total_exposed(row_from, col_from); - auto resist = resistant(row_from, col_from); - // set up vector of numeric categories (infected = 1, susceptible = 2, - // exposed = 3) for drawing # moved in each category - std::vector categories(total_infecteds, 1); - categories.insert(categories.end(), suscepts, 2); - categories.insert(categories.end(), expose, 3); - categories.insert(categories.end(), resist, 4); - - std::vector draw = - draw_n_from_v(categories, total_hosts_moved, generator.movement()); - infected_moved = std::count(draw.begin(), draw.end(), 1); - susceptible_moved = std::count(draw.begin(), draw.end(), 2); - exposed_moved = std::count(draw.begin(), draw.end(), 3); - resistant_moved = std::count(draw.begin(), draw.end(), 4); - - if (exposed_moved > 0) { - std::vector exposed_draw = draw_n_from_cohorts( - exposed, exposed_moved, row_from, col_from, generator.movement()); - int index = 0; - for (auto& raster : exposed) { - raster(row_from, col_from) -= exposed_draw[index]; - raster(row_to, col_to) += exposed_draw[index]; - index += 1; - } - } - if (infected_moved > 0) { - std::vector mortality_draw = draw_n_from_cohorts( - mortality_tracker_vector, - infected_moved, - row_from, - col_from, - generator.movement()); - int index = 0; - for (auto& raster : mortality_tracker_vector) { - raster(row_from, col_from) -= mortality_draw[index]; - raster(row_to, col_to) += mortality_draw[index]; - index += 1; - } - } - // check that the location with host movement is in suitable cells. - // Since suitable-cells comes from the total hosts originally. The - // the first check is for total_hosts - if (total_hosts(row_to, col_to) == 0) { - for (auto indices : suitable_cells) { - int i = indices[0]; - int j = indices[1]; - if ((i == row_to) && (j == col_to)) { - std::vector added_index = {row_to, col_to}; - suitable_cells.push_back(added_index); - break; - } - } - } - - infected(row_from, col_from) -= infected_moved; - susceptible(row_from, col_from) -= susceptible_moved; - total_hosts(row_from, col_from) -= total_hosts_moved; - total_exposed(row_from, col_from) -= exposed_moved; - resistant(row_from, col_from) -= resistant_moved; - infected(row_to, col_to) += infected_moved; - susceptible(row_to, col_to) += susceptible_moved; - total_hosts(row_to, col_to) += total_hosts_moved; - total_exposed(row_to, col_to) += exposed_moved; - resistant(row_to, col_to) += resistant_moved; - } - return movements.size(); + HostMovement + host_movement{step, last_index, movements, movement_schedule}; + IntegerRaster empty; + StandardHostPool hosts{ + model_type_, + susceptible, + exposed, + 0, + infected, + total_exposed, + resistant, + mortality_tracker_vector, + empty, + total_hosts, + *environment(true), + false, + 0, + false, + 0, + 0, + 0, + suitable_cells}; + return host_movement.action(hosts, generator); } /** Generates dispersers based on infected * * @param[out] dispersers (existing values are ignored) + * @param[out] established_dispersers Dispersers for a cell established in + * another cell (later modified to final form in disperse()) * @param infected Currently infected hosts * @param weather Whether to use the weather coefficient * @param reproductive_rate reproductive rate (used unmodified when weather * coefficient is not used) * @param[in] suitable_cells List of indices of cells with hosts + * @param generator Provider of random number generators */ - template void generate( IntegerRaster& dispersers, IntegerRaster& established_dispersers, @@ -561,44 +396,44 @@ class Simulation bool weather, double reproductive_rate, const std::vector>& suitable_cells, - GeneratorProvider& generator) + Generator& generator) { - double lambda = reproductive_rate; - for (auto indices : suitable_cells) { - int i = indices[0]; - int j = indices[1]; - if (infected(i, j) > 0) { - if (weather) - lambda = - reproductive_rate * environment()->weather_coefficient_at(i, j); - int dispersers_from_cell = 0; - if (dispersers_stochasticity_) { - std::poisson_distribution distribution(lambda); - for (int k = 0; k < infected(i, j); k++) { - dispersers_from_cell += - distribution(generator.disperser_generation()); - } - } - else { - dispersers_from_cell = lambda * infected(i, j); - } - if (soil_pool_) { - // 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); - soil_pool_->dispersers_to( - dispersers_to_soil, i, j, generator.soil()); - dispersers_from_cell -= dispersers_to_soil; - } - dispersers(i, j) = dispersers_from_cell; - established_dispersers(i, j) = dispersers_from_cell; - } - else { - dispersers(i, j) = 0; - established_dispersers(i, j) = 0; - } - } + IntegerRaster empty; + std::vector empty_vector; + StandardHostPool host_pool{ + model_type_, + empty, + empty_vector, + 0, + const_cast(infected), + empty, + empty, + empty_vector, + empty, + empty, + *environment(!weather), + dispersers_stochasticity_, + reproductive_rate, + false, + 0, + 0, + 0, + const_cast>&>(suitable_cells)}; + std::vector> empty_outside_dispersers; + StandardPestPool pests{ + dispersers, established_dispersers, empty_outside_dispersers}; + std::default_random_engine unused_kernel; + SpreadAction< + StandardHostPool, + StandardPestPool, + IntegerRaster, + FloatRaster, + RasterIndex, + std::default_random_engine, + Generator> + spread_action{unused_kernel}; + spread_action.activate_soils(soil_pool_, to_soil_percentage_); + spread_action.generate(host_pool, pests, generator); } /** Creates dispersal locations for the dispersing individuals @@ -630,10 +465,13 @@ class Simulation * will establish. * * @param[in] dispersers Dispersing individuals ready to be dispersed + * @param[in,out] established_dispersers Dispersers for a cell established in + * another cell * @param[in,out] susceptible Susceptible hosts - * @param[in,out] exposed_or_infected Exposed or infected hosts + * @param[in,out] exposed Exposed hosts + * @param[in,out] infected Infected hosts * @param[in,out] mortality_tracker Newly infected hosts (if applicable) - * @param[in, out] total_exposed Total exposed in all exposed cohorts + * @param[in,out] total_exposed Total exposed in all exposed cohorts * @param[in] total_populations All host and non-host individuals in the area * @param[in,out] outside_dispersers Dispersers escaping the raster * @param weather Whether or not weather coefficients should be used @@ -641,79 +479,141 @@ class Simulation * @param establishment_probability Probability of establishment with no * stochasticity * @param[in] suitable_cells List of indices of cells with hosts + * @param generator Provider of random number generators * * @note If the parameters or their default values don't correspond * with the disperse_and_infect() function, it is a bug. */ - template + template + void disperse( + IntegerRaster& dispersers, + IntegerRaster& established_dispersers, + IntegerRaster& susceptible, + std::vector& exposed, + IntegerRaster& infected, + std::vector& mortality_tracker, + const IntegerRaster& total_populations, + IntegerRaster& total_exposed, + std::vector>& outside_dispersers, + bool weather, + DispersalKernel& dispersal_kernel, + std::vector>& suitable_cells, + double establishment_probability, + Generator& generator) + { + // The interaction does not happen over the member variables yet, use empty + // variables. This requires SI/SEI to be fully resolved in host and not in + // disperse_and_infect. + IntegerRaster empty; + StandardHostPool host_pool{ + model_type_, + susceptible, + exposed, + 0, + infected, + total_exposed, + empty, + mortality_tracker, + empty, + empty, + *environment(!weather), + false, + 0, + establishment_stochasticity_, + establishment_probability, + rows_, + cols_, + suitable_cells}; + // This would be part of the main initialization process. + if (environment_) { + environment_->set_total_population(&total_populations); + } + StandardPestPool pests{dispersers, established_dispersers, outside_dispersers}; + SpreadAction< + StandardHostPool, + StandardPestPool, + IntegerRaster, + FloatRaster, + RasterIndex, + DispersalKernel, + Generator> + spread_action{dispersal_kernel}; + spread_action.activate_soils(soil_pool_, to_soil_percentage_); + spread_action.disperse(host_pool, pests, generator); + } + + // For backwards compatibility for tests (without exposed and mortality) + template void disperse( const IntegerRaster& dispersers, IntegerRaster& established_dispersers, IntegerRaster& susceptible, - IntegerRaster& exposed_or_infected, + IntegerRaster& infected, IntegerRaster& mortality_tracker, const IntegerRaster& total_populations, IntegerRaster& total_exposed, std::vector>& outside_dispersers, bool weather, DispersalKernel& dispersal_kernel, - const std::vector>& suitable_cells, - double establishment_probability, - GeneratorProvider& generator) + std::vector>& suitable_cells, + double establishment_probability = 0.5) { - std::uniform_real_distribution distribution_uniform(0.0, 1.0); - int row; - int col; + std::vector tmp; + tmp.push_back(mortality_tracker); + std::vector empty_vector; + disperse( + dispersers, + established_dispersers, + susceptible, + empty_vector, + infected, + tmp, // mortality + total_populations, + total_exposed, + outside_dispersers, + weather, + dispersal_kernel, + suitable_cells, + establishment_probability); + mortality_tracker = tmp.back(); + } - for (auto indices : suitable_cells) { - int i = indices[0]; - int j = indices[1]; - if (dispersers(i, j) > 0) { - for (int k = 0; k < dispersers(i, j); k++) { - std::tie(row, col) = dispersal_kernel(generator, i, j); - if (row < 0 || row >= rows_ || col < 0 || col >= cols_) { - // export dispersers dispersed outside of modeled area - outside_dispersers.emplace_back(std::make_tuple(row, col)); - established_dispersers(i, j) -= 1; - continue; - } - // Put a disperser to the host pool. - auto dispersed = disperser_to( - row, - col, - susceptible, - exposed_or_infected, - mortality_tracker, - total_populations, - total_exposed, - weather, - establishment_probability, - generator.establishment()); - if (!dispersed) { - established_dispersers(i, j) -= 1; - } - } - } - if (soil_pool_) { - // Get dispersers from the soil if there are any. - auto num_dispersers = - soil_pool_->dispersers_from(i, j, generator.soil()); - // Put each disperser to the host pool. - for (int k = 0; k < num_dispersers; k++) { - disperser_to( - i, - j, - susceptible, - exposed_or_infected, - mortality_tracker, - total_populations, - total_exposed, - weather, - establishment_probability, - generator.soil()); - } - } - } + // For backwards compatibility for tests (without exposed and mortality) + template + void disperse( + IntegerRaster& dispersers, + IntegerRaster& established_dispersers, + IntegerRaster& susceptible, + IntegerRaster& infected, + IntegerRaster& mortality_tracker, + const IntegerRaster& total_populations, + IntegerRaster& total_exposed, + std::vector>& outside_dispersers, + bool weather, + DispersalKernel& dispersal_kernel, + std::vector>& suitable_cells, + double establishment_probability, + Generator& generator) + { + std::vector tmp; + tmp.push_back(mortality_tracker); + std::vector empty_vector; + disperse( + dispersers, + established_dispersers, + susceptible, + empty_vector, + infected, + tmp, // mortality + total_populations, + total_exposed, + outside_dispersers, + weather, + dispersal_kernel, + suitable_cells, + establishment_probability, + generator); + mortality_tracker = tmp.back(); } /** Move overflowing pest population to other hosts. @@ -741,165 +641,62 @@ class Simulation * @param overpopulation_percentage Percentage of occupied hosts when the cell is * considered to be overpopulated * @param leaving_percentage Percentage pests leaving an overpopulated cell + * @param generator Provider of random number generators * * @note Exposed hosts do not count towards total number of pest, * i.e., *total_host* is assumed to be S + E in SEI model. * @note Mortality is not supported by this function, i.e., the mortality rasters * are not modified while the infected are. */ - template + template void move_overpopulated_pests( IntegerRaster& susceptible, IntegerRaster& infected, const IntegerRaster& total_hosts, std::vector>& outside_dispersers, DispersalKernel& dispersal_kernel, - const std::vector>& suitable_cells, + std::vector>& suitable_cells, double overpopulation_percentage, double leaving_percentage, - GeneratorProvider& generator) - { - struct Move - { - RasterIndex row; - RasterIndex col; - int count; - }; - std::vector moves; - - // Identify the moves. Remove from source cells. - for (auto indices : suitable_cells) { - int i = indices[0]; - int j = indices[1]; - int original_count = infected(i, j); - // No move with only one infected host (one unit). - if (original_count <= 1) - continue; - // r = I / (I + S) - // r = I / (I + S + E_1 + E_2 + ...) - double ratio = original_count / double(total_hosts(i, j)); - if (ratio >= overpopulation_percentage) { - int row; - int col; - std::tie(row, col) = dispersal_kernel(generator.overpopulation(), i, j); - // for leaving_percentage == 0.5 - // 2 infected -> 1 leaving - // 3 infected -> 1 leaving - int leaving = original_count * leaving_percentage; - susceptible(i, j) += leaving; - infected(i, j) -= leaving; - if (row < 0 || row >= rows_ || col < 0 || col >= cols_) { - // Collect pests dispersed outside of modeled area. - outside_dispersers.reserve(outside_dispersers.size() + leaving); - for (int pest = 0; pest < leaving; ++pest) - outside_dispersers.emplace_back(row, col); - continue; - } - // Doing the move here would create inconsistent results as some - // target cells would be evaluated after the moved pest arrived, - // possibly triggering another move. - // So, instead, we just collect them and apply later. (The pest is in - // the air when in the moves vector.) - moves.push_back(Move({row, col, leaving})); - } - } - // Perform the moves to target cells. - for (const auto& move : moves) { - RasterIndex row = move.row; - RasterIndex col = move.col; - int leaving = move.count; - // The target cell can accept all. - if (susceptible(row, col) >= leaving) { - susceptible(row, col) -= leaving; - infected(row, col) += leaving; - } - // More pests than the target cell can accept. - // This can happen if there is simply not enough S hosts to accommodate all - // the pests moving from the source or if multiple sources end up in the - // same target cell and there is not enough S hosts to accommodate all of - // them. The pests just disappear in both cases. - else { - leaving = susceptible(row, col); - susceptible(row, col) -= leaving; - infected(row, col) += leaving; - } - } - } - - /** Infect exposed hosts (E to I transition in the SEI model) - * - * Applicable to SEI model, no-operation otherwise, i.e., parameters - * are left intact for other models. - * - * The exposed vector are the hosts exposed in the previous steps. - * The length of the vector is the number of steps of the latency - * period plus one. Before the first latency period is over, - * the E to I transition won't happen because no item in the exposed - * vector is old enough to become infected. - * - * The position of the items in the exposed vector determines their - * age, i.e., for how long the hosts are exposed. The oldest item - * is at the front and youngest at the end. - * Before the the first latency period is over, items in the front - * are still empty (unused) because no hosts were exposed for the - * given time period. - * After the first latency - * period, this needs to be true before the function is called and - * it is true after the function - * finished with the difference that after the function is called, - * the last item is empty in the sense that it does not contain any - * hosts. - * - * When the E to I transition happens, hosts from the oldest item - * in the exposed vector are moved to the infected (and mortality - * tracker). They are removed from the exposed item and this item - * is moved to the back of the vector. - * - * Like in disperse(), there is no distinction between *infected* - * and *mortality_tracker*, but different usage is expected outside - * of this function. - * - * The raster class used with the simulation class needs to support - * `.fill()` method for this function to work. - * - * @param step Step in the simulation (>=0) - * @param exposed Vector of exposed hosts - * @param infected Infected hosts - * @param mortality_tracker Newly infected hosts - * @param total_exposed Total exposed in all exposed cohorts - */ - void infect_exposed( - unsigned step, - std::vector& exposed, - IntegerRaster& infected, - IntegerRaster& mortality_tracker, - IntegerRaster& total_exposed) + Generator& generator) { - if (model_type_ == ModelType::SusceptibleExposedInfected) { - if (step >= latency_period_) { - // Oldest item needs to be in the front - auto& oldest = exposed.front(); - // Move hosts to infected raster - infected += oldest; - mortality_tracker += oldest; - total_exposed += (oldest * (-1)); - // Reset the raster - // (hosts moved from the raster) - oldest.fill(0); - } - // Age the items and the used one to the back - // elements go one position to the left - // new oldest goes to the front - // old oldest goes to the back - rotate_left_by_one(exposed); - } - else if (model_type_ == ModelType::SusceptibleInfected) { - // no-op - } - else { - throw std::runtime_error( - "Unknown ModelType value in Simulation::infect_exposed()"); - } + UNUSED(total_hosts); // Total hosts is computed now. + IntegerRaster empty; + std::vector empty_vector; + StandardHostPool hosts{ + model_type_, + susceptible, + empty_vector, + 0, + infected, + empty, + empty, + empty_vector, + empty, + empty, + *environment(true), + false, + 0, + false, + 0, + 0, + 0, + suitable_cells}; + StandardPestPool pests{empty, empty, outside_dispersers}; + MoveOverpopulatedPests< + StandardHostPool, + StandardPestPool, + IntegerRaster, + FloatRaster, + RasterIndex, + DispersalKernel> + move_pest{ + dispersal_kernel, + overpopulation_percentage, + leaving_percentage, + rows_, + cols_}; + move_pest.action(hosts, pests, generator); } /** Disperse, expose, and infect based on dispersers @@ -928,35 +725,30 @@ class Simulation * has parameter *exposed* which is the same as the one from the * infect_exposed() function. */ - template + template void disperse_and_infect( unsigned step, - const IntegerRaster& dispersers, + IntegerRaster& dispersers, IntegerRaster& established_dispersers, IntegerRaster& susceptible, std::vector& exposed, IntegerRaster& infected, - IntegerRaster& mortality_tracker, + std::vector& mortality_tracker, const IntegerRaster& total_populations, IntegerRaster& total_exposed, std::vector>& outside_dispersers, bool weather, DispersalKernel& dispersal_kernel, - const std::vector>& suitable_cells, + std::vector>& suitable_cells, double establishment_probability, - GeneratorProvider& generator) + Generator& generator) { - auto* infected_or_exposed = &infected; - if (model_type_ == ModelType::SusceptibleExposedInfected) { - // The empty - not yet exposed - raster is in the back - // and will become youngest exposed one. - infected_or_exposed = &exposed.back(); - } this->disperse( dispersers, established_dispersers, susceptible, - *infected_or_exposed, + exposed, + infected, mortality_tracker, total_populations, total_exposed, @@ -966,10 +758,66 @@ class Simulation suitable_cells, establishment_probability, generator); - if (model_type_ == ModelType::SusceptibleExposedInfected) { - this->infect_exposed( - step, exposed, infected, mortality_tracker, total_exposed); - } + IntegerRaster empty; + StandardHostPool host_pool{ + model_type_, + susceptible, + exposed, + latency_period_, + infected, + total_exposed, + empty, + mortality_tracker, + empty, + empty, + *environment(!weather), + false, + 0, + establishment_stochasticity_, + establishment_probability, + rows_, + cols_, + suitable_cells}; + host_pool.step_forward(step); + } + + template + void disperse_and_infect( + unsigned step, + IntegerRaster& dispersers, + IntegerRaster& established_dispersers, + IntegerRaster& susceptible, + std::vector& exposed, + IntegerRaster& infected, + IntegerRaster& mortality_tracker, + const IntegerRaster& total_populations, + IntegerRaster& total_exposed, + std::vector>& outside_dispersers, + bool weather, + DispersalKernel& dispersal_kernel, + std::vector>& suitable_cells, + double establishment_probability, + Generator& generator) + { + std::vector tmp; + tmp.push_back(mortality_tracker); + disperse_and_infect( + step, + dispersers, + established_dispersers, + susceptible, + exposed, + infected, + tmp, + total_populations, + total_exposed, + outside_dispersers, + weather, + dispersal_kernel, + suitable_cells, + establishment_probability, + generator); + mortality_tracker = tmp.back(); } /** @@ -989,7 +837,8 @@ class Simulation * @param dispersers_percentage Percentage of dispersers moving to the soil */ void activate_soils( - std::shared_ptr> soil_pool, + std::shared_ptr> + soil_pool, double dispersers_percentage) { this->soil_pool_ = soil_pool; diff --git a/inst/include/soils.hpp b/inst/include/soils.hpp index 58d62143..8e23435d 100644 --- a/inst/include/soils.hpp +++ b/inst/include/soils.hpp @@ -33,7 +33,11 @@ namespace pops { * * Takes care of adding dispersers to the pool and of taking them out. */ -template +template< + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex, + typename GeneratorProvider> class SoilPool { public: @@ -51,7 +55,8 @@ class SoilPool */ SoilPool( std::vector& rasters, - const Environment& environment, + const Environment& + environment, bool generate_stochasticity = true, bool establishment_stochasticity = true, double fixed_establishment_probability = 0) @@ -63,7 +68,7 @@ class SoilPool { if (rasters.empty()) { throw std::logic_error( - "List of rasters of SpoilPool needs to have at least one item"); + "List of rasters of SoilPool needs to have at least one item"); } } @@ -163,7 +168,8 @@ class SoilPool /** * Surrounding environment */ - const Environment* environment_{nullptr}; + const Environment* + environment_{nullptr}; bool generate_stochasticity_{false}; ///< Outgoing dispersers stochasticity bool establishment_stochasticity_{false}; ///< Incoming dispersers /** diff --git a/inst/include/spread_rate.hpp b/inst/include/spread_rate.hpp index d062b416..b56f85da 100644 --- a/inst/include/spread_rate.hpp +++ b/inst/include/spread_rate.hpp @@ -27,19 +27,87 @@ namespace pops { /** * Class storing and computing step spread rate for one simulation. */ -template -class SpreadRate +template +class SpreadRateAction { +public: + SpreadRateAction( + const Hosts& hosts, + RasterIndex rows, + RasterIndex cols, + double ew_res, + double ns_res, + unsigned num_steps) + : width_(rows), + height_(cols), + west_east_resolution_(ew_res), + north_south_resolution_(ns_res), + num_steps_(num_steps), + boundaries_(num_steps + 1, std::make_tuple(0, 0, 0, 0)), + rates_( + num_steps, + std::make_tuple(std::nan(""), std::nan(""), std::nan(""), std::nan(""))) + { + boundaries_.at(0) = infection_boundary(hosts); + } + + /** + * Returns rate for certain year of simulation + */ + const BBoxFloat& step_rate(unsigned step) const + { + return rates_[step]; + } + + /** + * Computes spread rate in n, s, e, w directions + * for certain simulation year based on provided + * infection raster and bbox of infection in previous year. + * Unit is distance (map units) per year. + * If spread rate is zero and the bbox is touching the edge, + * that means spread is out of bounds and rate is set to NaN. + */ + void action(const Hosts& hosts, unsigned step) + { + BBoxInt bbox = infection_boundary(hosts); + boundaries_.at(step + 1) = bbox; + if (!is_boundary_valid(bbox)) { + rates_.at(step) = + std::make_tuple(std::nan(""), std::nan(""), std::nan(""), std::nan("")); + return; + } + int n1, n2, s1, s2, e1, e2, w1, w2; + std::tie(n1, s1, e1, w1) = boundaries_.at(step); + std::tie(n2, s2, e2, w2) = bbox; + double n_rate = ((n1 - n2) * north_south_resolution_); + double s_rate = ((s2 - s1) * north_south_resolution_); + double e_rate = ((e2 - e1) * west_east_resolution_); + double w_rate = ((w1 - w2) * west_east_resolution_); + + bool bn, bs, be, bw; + std::tie(bn, bs, be, bw) = is_out_of_bounds(bbox); + if (n_rate == 0 && bn) + n_rate = std::nan(""); + if (s_rate == 0 && bs) + s_rate = std::nan(""); + if (e_rate == 0 && be) + e_rate = std::nan(""); + if (w_rate == 0 && bw) + w_rate = std::nan(""); + + rates_.at(step) = std::make_tuple(n_rate, s_rate, e_rate, w_rate); + } + private: - int width; - int height; + int width_; + int height_; // the west-east resolution of the pixel - double west_east_resolution; + double west_east_resolution_; // the north-south resolution of the pixel - double north_south_resolution; - unsigned num_steps; - std::vector boundaries; - std::vector rates; + double north_south_resolution_; + unsigned num_steps_; + std::vector boundaries_; + std::vector rates_; /** * Return tuple of booleans indicating @@ -53,32 +121,32 @@ class SpreadRate std::tie(n, s, e, w) = bbox; if (n == 0) bn = true; - if (s == (height - 1)) + if (s == (height_ - 1)) bs = true; if (w == 0) bw = true; - if (e == (width - 1)) + if (e == (width_ - 1)) be = true; return std::make_tuple(bn, bs, be, bw); } + /** * Finds bbox of infections and returns * north, south, east, west coordinates (as number of rows/cols), * If there is no infection, sets -1 to all directions. */ - BBoxInt infection_boundary( - const Raster& raster, const std::vector>& suitable_cells) + BBoxInt infection_boundary(const Hosts& hosts) { - int n = height - 1; + int n = height_ - 1; int s = 0; int e = 0; - int w = width - 1; + int w = width_ - 1; bool found = false; - for (auto indices : suitable_cells) { + for (auto indices : hosts.suitable_cells()) { int i = indices[0]; int j = indices[1]; - auto value = raster(i, j); + auto value = hosts.infected_at(i, j); if (value > 0) { found = true; if (i < n) @@ -109,78 +177,6 @@ class SpreadRate return false; return true; } - -public: - SpreadRate( - const Raster& raster, - double ew_res, - double ns_res, - unsigned num_steps, - const std::vector>& suitable_cells) - : width(raster.cols()), - height(raster.rows()), - west_east_resolution(ew_res), - north_south_resolution(ns_res), - num_steps(num_steps), - boundaries(num_steps + 1, std::make_tuple(0, 0, 0, 0)), - rates( - num_steps, - std::make_tuple(std::nan(""), std::nan(""), std::nan(""), std::nan(""))) - { - boundaries.at(0) = infection_boundary(raster, suitable_cells); - } - - SpreadRate() = delete; - - /** - * Returns rate for certain year of simulation - */ - const BBoxFloat& step_rate(unsigned step) const - { - return rates[step]; - } - - /** - * Computes spread rate in n, s, e, w directions - * for certain simulation year based on provided - * infection raster and bbox of infection in previous year. - * Unit is distance (map units) per year. - * If spread rate is zero and the bbox is touching the edge, - * that means spread is out of bounds and rate is set to NaN. - */ - void compute_step_spread_rate( - const Raster& raster, - unsigned step, - const std::vector>& suitable_cells) - { - BBoxInt bbox = infection_boundary(raster, suitable_cells); - boundaries.at(step + 1) = bbox; - if (!is_boundary_valid(bbox)) { - rates.at(step) = - std::make_tuple(std::nan(""), std::nan(""), std::nan(""), std::nan("")); - return; - } - int n1, n2, s1, s2, e1, e2, w1, w2; - std::tie(n1, s1, e1, w1) = boundaries.at(step); - std::tie(n2, s2, e2, w2) = bbox; - double n_rate = ((n1 - n2) * north_south_resolution); - double s_rate = ((s2 - s1) * north_south_resolution); - double e_rate = ((e2 - e1) * west_east_resolution); - double w_rate = ((w1 - w2) * west_east_resolution); - - bool bn, bs, be, bw; - std::tie(bn, bs, be, bw) = is_out_of_bounds(bbox); - if (n_rate == 0 && bn) - n_rate = std::nan(""); - if (s_rate == 0 && bs) - s_rate = std::nan(""); - if (e_rate == 0 && be) - e_rate = std::nan(""); - if (w_rate == 0 && bw) - w_rate = std::nan(""); - - rates.at(step) = std::make_tuple(n_rate, s_rate, e_rate, w_rate); - } }; /** @@ -188,9 +184,9 @@ class SpreadRate * from vector of spread rates. * Checks if any rate is nan to not include it in the average. */ -template -BBoxFloat -average_spread_rate(const std::vector>& rates, unsigned step) +template +BBoxFloat average_spread_rate( + const std::vector>& rates, unsigned step) { // loop through stochastic runs double n, s, e, w; diff --git a/inst/include/treatments.hpp b/inst/include/treatments.hpp index 784d2a68..39381a65 100644 --- a/inst/include/treatments.hpp +++ b/inst/include/treatments.hpp @@ -21,6 +21,7 @@ #include "date.hpp" #include "scheduling.hpp" #include "utils.hpp" +#include "host_pool.hpp" #include #include @@ -28,6 +29,9 @@ #include #include +// only temporarily for direct host pool creation +#include + namespace pops { /** @@ -75,7 +79,7 @@ inline TreatmentApplication treatment_app_enum_from_string(const std::string& te * general set of parameters only when * needed for additional classes. */ -template +template class AbstractTreatment { public: @@ -83,20 +87,8 @@ class AbstractTreatment virtual unsigned get_end() = 0; virtual bool should_start(unsigned step) = 0; virtual bool should_end(unsigned step) = 0; - virtual void apply_treatment( - IntegerRaster& infected, - std::vector& exposed, - IntegerRaster& susceptible, - IntegerRaster& resistant, - IntegerRaster& total_hosts, - const std::vector>& spatial_indices) = 0; - virtual void end_treatment( - IntegerRaster& susceptible, - IntegerRaster& resistant, - const std::vector>& spatial_indices) = 0; - virtual void apply_treatment_mortality( - IntegerRaster& infected, - const std::vector>& spatial_indices) = 0; + virtual void apply_treatment(HostPool& host_pool) = 0; + virtual void end_treatment(HostPool& host_pool) = 0; virtual ~AbstractTreatment() {} }; @@ -104,8 +96,8 @@ class AbstractTreatment * Base treatment class. * Holds functions common between all treatment classes. */ -template -class BaseTreatment : public AbstractTreatment +template +class BaseTreatment : public AbstractTreatment { protected: unsigned start_step_; @@ -131,20 +123,23 @@ class BaseTreatment : public AbstractTreatment { return end_step_; } - void apply_treatment_mortality( - IntegerRaster& infected, - const std::vector>& suitable_cells) override + + // returning double allows identical results with the previous version + double get_treated(int i, int j, int count) { - for (auto indices : suitable_cells) { - int i = indices[0]; - int j = indices[1]; - if (application_ == TreatmentApplication::Ratio) { - infected(i, j) = infected(i, j) - (infected(i, j) * map_(i, j)); - } - else if (application_ == TreatmentApplication::AllInfectedInCell) { - infected(i, j) = map_(i, j) ? 0 : infected(i, j); - } + return get_treated(i, j, count, this->application_); + } + + double get_treated(int i, int j, int count, TreatmentApplication application) + { + if (application == TreatmentApplication::Ratio) { + return count * this->map_(i, j); } + else if (application == TreatmentApplication::AllInfectedInCell) { + return this->map_(i, j) ? count : 0; + } + throw std::runtime_error( + "BaseTreatment::get_treated: unknown TreatmentApplication"); } }; @@ -153,15 +148,15 @@ class BaseTreatment : public AbstractTreatment * Removes percentage (given by treatment efficiency) * of infected and susceptible host (e.g. cut down trees). */ -template -class SimpleTreatment : public BaseTreatment +template +class SimpleTreatment : public BaseTreatment { public: SimpleTreatment( const FloatRaster& map, unsigned start, TreatmentApplication treatment_application) - : BaseTreatment(map, start, treatment_application) + : BaseTreatment(map, start, treatment_application) {} bool should_start(unsigned step) override { @@ -173,54 +168,36 @@ class SimpleTreatment : public BaseTreatment { return false; } - void apply_treatment( - IntegerRaster& infected, - std::vector& exposed, - IntegerRaster& susceptible, - IntegerRaster& resistant, - IntegerRaster& total_hosts, - const std::vector>& suitable_cells) override + void apply_treatment(HostPool& host_pool) override { - for (auto indices : suitable_cells) { + for (auto indices : host_pool.suitable_cells()) { int i = indices[0]; int j = indices[1]; - int new_exposed_total = 0; - int new_infected = 0; - int new_susceptible = 0; - int new_exposed_individual = 0; - if (this->application_ == TreatmentApplication::Ratio) { - new_infected = infected(i, j) - (infected(i, j) * this->map_(i, j)); - infected(i, j) = new_infected; - } - else if (this->application_ == TreatmentApplication::AllInfectedInCell) { - new_infected = this->map_(i, j) ? 0 : infected(i, j); - infected(i, j) = new_infected; + 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; + for (int count : host_pool.mortality_by_group_at(i, j)) { + remove_mortality.push_back(this->get_treated(i, j, count)); } - for (auto& raster : exposed) { - if (this->application_ == TreatmentApplication::Ratio) { - new_exposed_individual = - raster(i, j) - (raster(i, j) * this->map_(i, j)); - raster(i, j) = new_exposed_individual; - new_exposed_total += new_exposed_individual; - } - else if ( - this->application_ == TreatmentApplication::AllInfectedInCell) { - new_exposed_individual = - raster(i, j) - (raster(i, j) * this->map_(i, j)); - raster(i, j) = new_exposed_individual; - new_exposed_total += new_exposed_individual; - } + + 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)); } - new_susceptible = - susceptible(i, j) - (susceptible(i, j) * this->map_(i, j)); - susceptible(i, j) = new_susceptible; - total_hosts(i, j) = - new_infected + new_susceptible + new_exposed_total + resistant(i, j); + host_pool.completely_remove_hosts_at( + i, + j, + remove_susceptible, + remove_exposed, + remove_infected, + remove_mortality); } } - void end_treatment( - IntegerRaster&, IntegerRaster&, const std::vector>&) override + void end_treatment(HostPool& host_pool) override { + UNUSED(host_pool); return; } }; @@ -231,8 +208,8 @@ class SimpleTreatment : public BaseTreatment * of infected and susceptible to resistant pool * and after certain number of days back to susceptible. */ -template -class PesticideTreatment : public BaseTreatment +template +class PesticideTreatment : public BaseTreatment { public: PesticideTreatment( @@ -240,7 +217,7 @@ class PesticideTreatment : public BaseTreatment unsigned start, unsigned end, TreatmentApplication treatment_application) - : BaseTreatment(map, start, treatment_application) + : BaseTreatment(map, start, treatment_application) { this->end_step_ = end; } @@ -256,58 +233,41 @@ class PesticideTreatment : public BaseTreatment return true; return false; } - - void apply_treatment( - IntegerRaster& infected, - std::vector& exposed_vector, - IntegerRaster& susceptible, - IntegerRaster& resistant, - IntegerRaster& total_hosts, - const std::vector>& suitable_cells) override + void apply_treatment(HostPool& host_pool) override { - UNUSED(total_hosts); - for (auto indices : suitable_cells) { + for (auto indices : host_pool.suitable_cells()) { int i = indices[0]; int j = indices[1]; - int infected_resistant = 0; - int exposed_resistant_sum = 0; - int susceptible_resistant = susceptible(i, j) * this->map_(i, j); - int current_resistant = resistant(i, j); - if (this->application_ == TreatmentApplication::Ratio) { - infected_resistant = infected(i, j) * this->map_(i, j); - } - else if (this->application_ == TreatmentApplication::AllInfectedInCell) { - infected_resistant = this->map_(i, j) ? infected(i, j) : 0; + // 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); + 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)); } - infected(i, j) -= infected_resistant; - for (auto& exposed : exposed_vector) { - int exposed_resistant = 0; - if (this->application_ == TreatmentApplication::Ratio) { - exposed_resistant = exposed(i, j) * this->map_(i, j); - } - else if ( - this->application_ == TreatmentApplication::AllInfectedInCell) { - exposed_resistant = this->map_(i, j) ? exposed(i, j) : 0; - } - exposed(i, j) -= exposed_resistant; - exposed_resistant_sum += exposed_resistant; + 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(i, j) = infected_resistant + exposed_resistant_sum - + susceptible_resistant + current_resistant; - susceptible(i, j) -= susceptible_resistant; + host_pool.make_resistant_at( + i, + j, + susceptible_resistant, + resistant_exposed_list, + this->get_treated(i, j, host_pool.infected_at(i, j)), + resistant_mortality_list); } } - void end_treatment( - IntegerRaster& susceptible, - IntegerRaster& resistant, - const std::vector>& suitable_cells) override + void end_treatment(HostPool& host_pool) override { - for (auto indices : suitable_cells) { + for (auto indices : host_pool.suitable_cells()) { int i = indices[0]; int j = indices[1]; if (this->map_(i, j) > 0) { - susceptible(i, j) += resistant(i, j); - resistant(i, j) = 0; + host_pool.remove_resistance_at(i, j); } } } @@ -325,11 +285,11 @@ class PesticideTreatment : public BaseTreatment * populations that overlap both spatially and temporally * are returned to the susceptible pool when the first treatment ends. */ -template +template class Treatments { private: - std::vector*> treatments; + std::vector*> treatments; Scheduler scheduler_; public: @@ -362,13 +322,13 @@ class Treatments { unsigned start = scheduler_.schedule_action_date(start_date); if (num_days == 0) - treatments.push_back(new SimpleTreatment( + treatments.push_back(new SimpleTreatment( map, start, treatment_application)); else { Date end_date(start_date); end_date.add_days(num_days); unsigned end = scheduler_.schedule_action_date(end_date); - treatments.push_back(new PesticideTreatment( + treatments.push_back(new PesticideTreatment( map, start, end, treatment_application)); } } @@ -379,64 +339,25 @@ class Treatments * activated/deactivated. * * \param current simulation step - * \param infected raster of infected host - * \param exposed Exposed hosts per cohort - * \param susceptible raster of susceptible host - * \param resistant raster of resistant host - * \param total_hosts All host individuals in the area (I+E+S in the cell) - * \param suitable_cells List of indices of cells with hosts + * \param host_pool Host * * \return true if any management action was necessary */ - bool manage( - unsigned current, - IntegerRaster& infected, - std::vector& exposed, - IntegerRaster& susceptible, - IntegerRaster& resistant, - IntegerRaster& total_hosts, - const std::vector>& suitable_cells) + bool manage(unsigned current, HostPool& host_pool) { bool changed = false; for (unsigned i = 0; i < treatments.size(); i++) { if (treatments[i]->should_start(current)) { - treatments[i]->apply_treatment( - infected, - exposed, - susceptible, - resistant, - total_hosts, - suitable_cells); + treatments[i]->apply_treatment(host_pool); changed = true; } else if (treatments[i]->should_end(current)) { - treatments[i]->end_treatment(susceptible, resistant, suitable_cells); + treatments[i]->end_treatment(host_pool); changed = true; } } return changed; } - /*! - * \brief Separately manage mortality infected cohorts - * \param current simulation step - * \param infected raster of infected host - * \param suitable_cells List of indices of cells with hosts - * - * \return true if any management action was necessary - */ - bool manage_mortality( - unsigned current, - IntegerRaster& infected, - const std::vector>& suitable_cells) - { - bool applied = false; - for (unsigned i = 0; i < treatments.size(); i++) - if (treatments[i]->should_start(current)) { - treatments[i]->apply_treatment_mortality(infected, suitable_cells); - applied = true; - } - return applied; - } /*! * \brief Used to remove treatments after certain step. * Needed for computational steering. diff --git a/inst/include/utils.hpp b/inst/include/utils.hpp index 05f211b7..20653c63 100644 --- a/inst/include/utils.hpp +++ b/inst/include/utils.hpp @@ -49,6 +49,7 @@ #include #include #include +#include /** * Return true if _container_ contains _value_. @@ -56,7 +57,16 @@ template bool container_contains(const Container& container, const Value& value) { - return container.find(value) != container.end(); + return std::find(container.begin(), container.end(), value) != container.end(); +} + +/** + * Return true if _container_ contains _key_. + */ +template +bool container_contains(const std::map& container, const Key& key) +{ + return container.count(key) > 0; } /** @@ -380,4 +390,23 @@ std::vector> find_suitable_cells(const RasterType& rast return cells; } +template +std::vector> +find_suitable_cells(const std::vector& rasters) +{ + std::vector> cells; + // The assumption is that the raster is sparse (otherwise we would not be doing + // this), so we have no number for the reserve method. + for (RasterIndex row = 0; row < rasters[0]->rows(); ++row) { + for (RasterIndex col = 0; col < rasters[0]->cols(); ++col) { + for (const auto& raster : rasters) { + if (raster->operator()(row, col) > 0) { + cells.push_back({row, col}); + break; + } + } + } + } + return cells; +} #endif // POPS_UTILS_HPP diff --git a/src/pops.cpp b/src/pops.cpp index b0768765..f3210198 100644 --- a/src/pops.cpp +++ b/src/pops.cpp @@ -196,6 +196,7 @@ List pops_model_cpp( std::vector> spread_rates_vector; std::tuple spread_rates; + IntegerMatrix dispersers(config.rows, config.cols); IntegerMatrix total_dispersers(config.rows, config.cols); IntegerMatrix established_dispersers(config.rows, config.cols); @@ -218,7 +219,45 @@ List pops_model_cpp( config.create_schedules(); - Treatments treatments(config.scheduler()); + ModelType mt = model_type_from_string(config.model_type); + using PoPSModel = Model; + PoPSModel model(config); + + PestHostUseTable pest_host_use_table( + config, model.environment()); + CompetencyTable competency_table( + config, model.environment()); + PoPSModel::StandardSingleHostPool host_pool( + mt, + susceptible, + exposed, + config.latency_period_steps, + infected, + total_exposed, + resistant, + mortality_tracker, + mortality, + total_hosts, + model.environment(), + config.generate_stochasticity, + config.reproductive_rate, + config.establishment_stochasticity, + config.establishment_probability, + config.rows, + config.cols, + spatial_indices); + + model.environment().add_host(&host_pool); + competency_table.add_host_competencies({1}, 1); + PoPSModel::StandardMultiHostPool multi_host_pool({&host_pool}, config); + multi_host_pool.set_pest_host_use_table(pest_host_use_table); + multi_host_pool.set_competency_table(competency_table); + PoPSModel::StandardPestPool pest_pool( + dispersers, + established_dispersers, + outside_dispersers); + + Treatments treatments(config.scheduler()); for (unsigned t = 0; t < treatment_maps.size(); t++) { treatments.add_treatment( treatment_maps[t], @@ -241,8 +280,13 @@ List pops_model_cpp( else { spread_rate_outputs = 0; } - SpreadRate spreadrate( - infected, config.ew_res, config.ns_res, spread_rate_outputs, spatial_indices); + SpreadRateAction spreadrate( + multi_host_pool, + config.rows, + config.cols, + config.ew_res, + config.ns_res, + spread_rate_outputs); unsigned move_scheduled; if (config.use_movements) { for (unsigned move = 0; move < movements_dates.size(); ++move) { @@ -261,7 +305,7 @@ List pops_model_cpp( quarantine_outputs = 0; } - QuarantineEscape quarantine( + QuarantineEscapeAction quarantine( quarantine_areas, config.ew_res, config.ns_res, @@ -294,13 +338,8 @@ List pops_model_cpp( network->load(network_stream); } - ModelType mt = model_type_from_string(config.model_type); WeatherType weather_typed = weather_type_from_string(config.weather_type); - Simulation simulation( - config.rows, config.cols, mt, config.latency_period_steps); - - Model model(config); if (use_soils) { model.activate_soils(soil_reservoirs); } @@ -308,8 +347,6 @@ List pops_model_cpp( for (unsigned current_index = 0; current_index < config.scheduler().get_num_steps(); ++current_index) { - IntegerMatrix dispersers(config.rows, config.cols); - auto weather_step = config.simulation_step_to_weather_step(current_index); if (weather_typed == WeatherType::Probabilistic) { model.environment().update_weather_from_distribution( @@ -322,27 +359,18 @@ List pops_model_cpp( model.run_step( current_index, - infected, - susceptible, - total_populations, - total_hosts, + multi_host_pool, + pest_pool, dispersers, - established_dispersers, - total_exposed, - exposed, - mortality_tracker, - mortality, + total_populations, + treatments, temperature, survival_rates, - treatments, - resistant, - outside_dispersers, spreadrate, quarantine, quarantine_areas, movements, - network ? *network : Network::null_network(), - spatial_indices); + network ? *network : Network::null_network()); // keeps track of cumulative dispersers or propagules from a site. if (config.spread_schedule()[current_index]) { From 993ac348694de69df4b724a11261b234c9bd8a5c Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Wed, 29 Nov 2023 11:22:24 -0500 Subject: [PATCH 02/12] small changes in API --- inst/cpp/pops-core | 2 +- inst/include/competency_table.hpp | 61 +++++++++++++++-- inst/include/config.hpp | 91 ++++++++++++++++++++++++-- inst/include/environment.hpp | 28 +++++++- inst/include/environment_interface.hpp | 18 +++++ inst/include/host_pool.hpp | 65 ++++++++++++++++-- inst/include/host_pool_interface.hpp | 4 +- inst/include/model.hpp | 11 ++-- inst/include/multi_host_pool.hpp | 2 +- inst/include/pest_host_use_table.hpp | 48 ++++++++++++-- inst/include/pest_pool.hpp | 11 ++++ inst/include/simulation.hpp | 13 +++- inst/include/utils.hpp | 7 +- src/pops.cpp | 6 +- 14 files changed, 331 insertions(+), 36 deletions(-) diff --git a/inst/cpp/pops-core b/inst/cpp/pops-core index bf490a9f..f334c0d8 160000 --- a/inst/cpp/pops-core +++ b/inst/cpp/pops-core @@ -1 +1 @@ -Subproject commit bf490a9f41851ea0e980cd05f386d8a041cff892 +Subproject commit f334c0d8fad399acd93bd7335c92e798f94ef8e8 diff --git a/inst/include/competency_table.hpp b/inst/include/competency_table.hpp index 90aa496b..f98e5034 100644 --- a/inst/include/competency_table.hpp +++ b/inst/include/competency_table.hpp @@ -24,13 +24,30 @@ namespace pops { -template +/** + * Competency table holding combinations of host presences and absences and the + * corresponding competency score. + */ +template class CompetencyTable { public: + using RasterIndex = typename HostPool::RasterIndex; using Environment = typename HostPool::Environment; + /** + * @brief Create an empty competency table + * + * @param environment Reference to the environment + */ CompetencyTable(const Environment& environment) : environment_(environment) {} + + /** + * @brief Create a competency table using values in config + * + * @param config Configuration with competency table data + * @param environment Reference to the environment + */ CompetencyTable(const Config& config, const Environment& environment) : environment_(environment) { @@ -39,12 +56,34 @@ class CompetencyTable } } + /** + * @brief Add competencies for a combination of host presences and absences + * + * Order of presence-absence data needs to be the same as host order in the + * environment. + * + * Order of calls does not matter. + * + * @param presence_absence Presence (true) and absence (false) for each host. + * @param competency Competency for a given presence-absence combination. + */ void add_host_competencies(const std::vector& presence_absence, double competency) { competency_table_.emplace_back(presence_absence, competency); } + /** + * @brief Get competency at a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * @param host Pointer to a specific host pool which is asking about competency + * + * @return Competency score + * + * @see find_competency() + */ double competency_at(RasterIndex row, RasterIndex col, const HostPool* host) const { auto presence_absence = environment_.host_presence_at(row, col); @@ -52,6 +91,18 @@ class CompetencyTable return find_competency(presence_absence, host_index); } +private: + /** + * @brief Find competency score in the table + * + * The *host_index* parameter is used to skip records which don't include the given + * host. + * + * @param presence_absence Actual presence-absence data at a given place (cell) + * @param host_index Index of the specific pool host originating the request + * + * @return Competency score + */ double find_competency(const std::vector& presence_absence, size_t host_index) const { @@ -89,7 +140,9 @@ class CompetencyTable return competency; } -private: + /** + * @brief One row of the comptenecy table + */ struct TableRow { std::vector presence_absence; @@ -100,8 +153,8 @@ class CompetencyTable : presence_absence(presence_absence), competency(competency) {} }; - std::vector competency_table_; - const Environment& environment_; + std::vector competency_table_; ///< Internal table for lookup + const Environment& environment_; ///< Environment (for host index) }; } // namespace pops diff --git a/inst/include/config.hpp b/inst/include/config.hpp index 14c4c770..c51311bf 100644 --- a/inst/include/config.hpp +++ b/inst/include/config.hpp @@ -118,12 +118,20 @@ std::map read_key_value_pairs( class Config { public: + /** + * @brief Row of a table with pest-host-use data + * + * One row is stored for each host. + */ struct PestHostUseTableDataRow { - double susceptibility; - double mortality_rate; - double mortality_time_lag; + double susceptibility; ///< Susceptibility for a given host + double mortality_rate; ///< Mortality rate for a given host + double mortality_time_lag; ///< Time lag of mortality in mortality steps }; + /** + * @brief Row of a table with competency data + */ struct CompetencyTableDataRow { std::vector presence_absence; @@ -182,8 +190,13 @@ class Config bool use_mortality{false}; std::string mortality_frequency; unsigned mortality_frequency_n; + /** Mortality rate if used without pest-host-use table + * + * @see read_pest_host_use_table() + */ double mortality_rate{0}; - int mortality_time_lag{0}; // Time lag of mortality in mortality steps + /** Time lag of mortality in simulation steps if used without pest-host-use table */ + int mortality_time_lag{0}; // Quarantine bool use_quarantine{false}; std::string quarantine_frequency; @@ -450,6 +463,10 @@ class Config season_end_month_ = std::stoi(end); } + /** + * @brief Set disperser arrival behavior for model + * @param value Arrival behavior string "infect" or "land" + */ void set_arrival_behavior(std::string value) { if (value != "infect" && value != "land") { @@ -459,6 +476,10 @@ class Config arrival_behavior_ = value; } + /** + * @brief Get disperser arrival behavior for model + * @return Arrival behavior as string "infect" or "land" + */ const std::string& arrival_behavior() const { return arrival_behavior_; @@ -512,15 +533,36 @@ class Config this->multiple_random_seeds = true; } + /** + * @brief Get data for the pest-host-use table + * @return Reference to the internal table + * + * @see PestHostUseTableDataRow + */ const std::vector& pest_host_use_table_data() const { return pest_host_use_table_data_; } + + /** + * @brief Get data for the competency table + * @return Reference to the internal table + * + * @see CompetencyTableDataRow + */ const std::vector& competency_table_data() const { return competency_table_data_; } + /** + * @brief Read pest-host-use table data from vector of vectors of doubles + * + * The nested vectors need to be of size 3. The order of values is susceptibility, + * mortality rate, and mortality time lag. + * + * @param values Table data + */ void read_pest_host_use_table(const std::vector>& values) { for (const auto& row : values) { @@ -536,6 +578,18 @@ class Config } } + /** + * @brief Use existing config parameters to create pest-host-use table + * + * This will create table with date for the given number of hosts with values for + * all hosts being the same. Susceptibility is set to 1 and mortality is taken from + * existing config attributes *mortality_rate* and *mortality_time_lag*. + * + * @param num_of_hosts Number of hosts + * + * @see #mortality_rate + * @see #mortality_time_lag + */ void create_pest_host_use_table_from_parameters(int num_of_hosts) { for (int i = 0; i < num_of_hosts; ++i) { @@ -547,6 +601,31 @@ class Config } } + /** + * @brief Read competency table from vector of vectors of doubles + * + * The nested vectors are rows of the table which need to have size 2 or higher. + * Each vector contains the combination of hosts and competency score. + * + * First n-1 items are the host presence and absence data which will be converted + * from double to bool, i.e., use 0 for absence, 1 for presence. The number of these + * items should be the number of hosts. Last (the nth) item in each vector is the + * competency score for the given combination and will be used as double. + * + * For 1 host with competency 1, the table is `{{1, 1}, {0, 0}}`. For 2 hosts, the + * table may look like this: + * + * ``` + * { + * {1, 0, 0.1}, + * {0, 1, 0.4}, + * {1, 1, 0.8}, + * {0, 0, 0} + * } + * ``` + * + * @param values Table data + */ void read_competency_table(const std::vector>& values) { for (const auto& row : values) { @@ -576,7 +655,7 @@ class Config Scheduler scheduler_{date_start_, date_end_, step_unit_, step_num_units_}; bool schedules_created_{false}; - std::string arrival_behavior_{"infect"}; + std::string arrival_behavior_{"infect"}; ///< Disperser arrival behavior std::vector spread_schedule_; std::vector output_schedule_; @@ -587,7 +666,9 @@ class Config std::vector quarantine_schedule_; std::vector weather_table_; + /** Storage for the pest-host-use table data */ std::vector pest_host_use_table_data_; + /** Storage for the competency table data */ std::vector competency_table_data_; }; diff --git a/inst/include/environment.hpp b/inst/include/environment.hpp index 1cdedf9b..7cb6da4a 100644 --- a/inst/include/environment.hpp +++ b/inst/include/environment.hpp @@ -200,6 +200,11 @@ class Environment return sum; } + /** + * @copydoc EnvironmentInterface::host_presence_at() + * + * Hosts are in the order of how host pools were registered to the environment. + */ std::vector host_presence_at(RasterIndex row, RasterIndex col) const override { std::vector presence; @@ -219,14 +224,35 @@ class Environment total_population_ = individuals; } + // While copydoc generates correct doc here, it triggers a + // warning with Doxygen 1.8.17, so using see here for now. + /** + * @see EnvironmentInterface::add_host() + * + * @note The function is no-op if host pool already registered. This may throw an + * exception in the future. + */ void add_host(const HostPoolInterface* host) override { - // no-op if already there, may become an error in the future if (container_contains(hosts_, host)) return; hosts_.push_back(host); } + /** + * @brief Remove all hosts from the environment. + * + * This function is useful for reusing an environment object in different contexts, + * especially in tests, but it does not have an epidemiological meaning. + */ + void remove_hosts() + { + hosts_.clear(); + } + + /** + * @copydoc EnvironmentInterface::host_index() + */ size_t host_index(const HostPoolInterface* host) const override { auto it = std::find(hosts_.begin(), hosts_.end(), host); diff --git a/inst/include/environment_interface.hpp b/inst/include/environment_interface.hpp index 62e3e0f2..ea804617 100644 --- a/inst/include/environment_interface.hpp +++ b/inst/include/environment_interface.hpp @@ -38,12 +38,30 @@ class EnvironmentInterface virtual double influence_probability_of_establishment_at( RasterIndex row, RasterIndex col, double value) const = 0; virtual int total_population_at(RasterIndex row, RasterIndex col) const = 0; + /** + * @brief Get presence-absence for host at a given cell + * + * @param row Row index of the cell + * @param col Column index of the cell + * + * @return Vector of presence-absence data + */ virtual std::vector host_presence_at(RasterIndex row, RasterIndex col) const = 0; // in general, we may have other individuals-non const virtual void set_other_individuals(const IntegerRaster* individuals) = 0; virtual void set_total_population(const IntegerRaster* individuals) = 0; + /** + * @brief Register a host pool in the environment + * @param host Host pool to add + */ virtual void add_host(const HostPoolInterface* host) = 0; + /** + * @brief Get index of a host pool in the environment + * @param host Pointer to a specific host pool + * @return Index based on the order within the environment + * @throw std::invalid_argument if host is not present + */ virtual size_t host_index(const HostPoolInterface* host) const = 0; virtual const FloatRaster& weather_coefficient() const = 0; virtual void update_temperature(const FloatRaster& raster) = 0; diff --git a/inst/include/host_pool.hpp b/inst/include/host_pool.hpp index dafb8b20..7c55bcf4 100644 --- a/inst/include/host_pool.hpp +++ b/inst/include/host_pool.hpp @@ -84,6 +84,9 @@ class HostPool : public HostPoolInterface * host infection over time. Expectation is that mortality tracker is of * length (1/mortality_rate + mortality_time_lag). * + * 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 susceptible Raster of susceptible hosts * @param exposed Raster of exposed or infected hosts @@ -114,7 +117,7 @@ class HostPool : public HostPoolInterface std::vector& mortality_tracker_vector, IntegerRaster& died, IntegerRaster& total_hosts, - const Environment& environment, + Environment& environment, bool dispersers_stochasticity, double reproductive_rate, bool establishment_stochasticity, @@ -140,15 +143,39 @@ class HostPool : public HostPoolInterface rows_(rows), cols_(cols), suitable_cells_(suitable_cells) - {} + { + environment.add_host(this); + } + /** + * @brief Set pest-host-use table + * + * If set, apply_mortality_at(RasterIndex, RasterIndex) can be used instead of the + * version which takes mortality parameters directly. Susceptibility is used + * automatically if the table is set. + * + * Pointer to the existing object is stored and used. So, the table can be modified, + * but the table object needs to exists during the lifetime of this object. + * + * @param pest_host_use_table Reference to the table + */ void set_pest_host_use_table(const PestHostUseTable& pest_host_use_table) { this->pest_host_use_table_ = &pest_host_use_table; } - void - set_competency_table(const CompetencyTable& competency_table) + /** + * @brief Set competency table + * + * Competency is used automatically if the table is set. No competency is considered + * if the table is not set. + * + * Pointer to the existing object is stored and used. So, the table can be modified, + * but the table object needs to exists during the lifetime of this object. + * + * @param competency_table Reference to the table + */ + void set_competency_table(const CompetencyTable& competency_table) { this->competency_table_ = &competency_table; } @@ -186,6 +213,19 @@ class HostPool : public HostPoolInterface return 0; } + /** + * @brief Test whether a disperser can establish + * + * This static (object-idependent) function to test disperser establishement allows + * code reuse between a single host and a multi-host case. + * + * @param probability_of_establishment Probability of establishment + * @param establishment_stochasticity true if establishment stochasticity is enabled + * @param deterministic_establishment_probability Establishment probability for + * deterministic establishment + * @param generator Random number generator (used with enabled stochasticity) + * @return true if disperser can establish, false otherwise + */ static bool can_disperser_establish( double probability_of_establishment, bool establishment_stochasticity, @@ -605,9 +645,10 @@ class HostPool : public HostPoolInterface } /** - * @brief Remove percentage of infestation/infection + * @brief Remove percentage of infestation/infection (I->S) * - * remove the same percentage for total exposed and remove randomly from each cohort + * Besides removing percentage of infected, it removes the same percentage for total + * exposed and remove individuals randomly from each cohort. * * @param row Row index of the cell * @param col Column index of the cell @@ -1061,6 +1102,14 @@ class HostPool : public HostPoolInterface return suitable_cells_; } + /** + * @brief Get list with this host pools + * + * This is for compatibility with multi-host pool. In case of this host pool, it + * returns one item which is a pointer to this host pool. + * + * @return Read-write reference to a vector of size 1. + */ std::vector& host_pools() { return host_pools_; @@ -1108,8 +1157,10 @@ class HostPool : public HostPoolInterface bool establishment_stochasticity_{true}; double deterministic_establishment_probability_{0}; + /** Pest-host-use table */ const PestHostUseTable* pest_host_use_table_{nullptr}; - const CompetencyTable* competency_table_{nullptr}; + /** Competency table */ + const CompetencyTable* competency_table_{nullptr}; RasterIndex rows_{0}; RasterIndex cols_{0}; diff --git a/inst/include/host_pool_interface.hpp b/inst/include/host_pool_interface.hpp index 75df0602..877e8d80 100644 --- a/inst/include/host_pool_interface.hpp +++ b/inst/include/host_pool_interface.hpp @@ -28,10 +28,12 @@ namespace pops { * functionally different hosts or if we have dependencies which the interface would * address. */ -template +template class HostPoolInterface { public: + using RasterIndex = RasterIndexType; + virtual ~HostPoolInterface() {} /** diff --git a/inst/include/model.hpp b/inst/include/model.hpp index 2137057d..89b36e1d 100644 --- a/inst/include/model.hpp +++ b/inst/include/model.hpp @@ -126,17 +126,20 @@ class Model } public: + /** Type for single-host pool */ using StandardSingleHostPool = HostPool< IntegerRaster, FloatRaster, RasterIndex, RandomNumberGeneratorProvider>; + /** Type for multi-host pool */ using StandardMultiHostPool = MultiHostPool< StandardSingleHostPool, IntegerRaster, FloatRaster, RasterIndex, RandomNumberGeneratorProvider>; + /** Type for pest pool */ using StandardPestPool = PestPool; Model( @@ -259,7 +262,6 @@ class Model step, multi_host_pool, pest_pool, - dispersers, total_populations, treatments, temperatures, @@ -277,7 +279,6 @@ class Model * @param step Step number in the simulation * @param host_pool Host pool * @param pest_pool Pest pool - * @param[out] dispersers Dispersing individuals (used directly for kernels) * @param[in,out] total_populations All host and non-host individuals in the area * @param[in,out] treatments Treatments to be applied (also tracks use of * treatments) @@ -294,7 +295,6 @@ class Model int step, StandardMultiHostPool& host_pool, StandardPestPool& pest_pool, - IntegerRaster& dispersers, IntegerRaster& total_populations, Treatments& treatments, const std::vector& temperatures, @@ -332,9 +332,10 @@ class Model } // actual spread if (config_.spread_schedule()[step]) { - auto dispersal_kernel = kernel_factory_(config_, dispersers, network); + auto dispersal_kernel = + kernel_factory_(config_, pest_pool.dispersers(), network); auto overpopulation_kernel = - create_overpopulation_movement_kernel(dispersers, network); + create_overpopulation_movement_kernel(pest_pool.dispersers(), network); SpreadAction< StandardMultiHostPool, StandardPestPool, diff --git a/inst/include/multi_host_pool.hpp b/inst/include/multi_host_pool.hpp index e7bdf5aa..280bc70d 100644 --- a/inst/include/multi_host_pool.hpp +++ b/inst/include/multi_host_pool.hpp @@ -79,7 +79,7 @@ class MultiHostPool * * @param table Reference to the table object */ - void set_competency_table(const CompetencyTable& table) + void set_competency_table(const CompetencyTable& table) { for (auto& host_pool : host_pools_) { host_pool->set_competency_table(table); diff --git a/inst/include/pest_host_use_table.hpp b/inst/include/pest_host_use_table.hpp index 28c49e7e..d2879a11 100644 --- a/inst/include/pest_host_use_table.hpp +++ b/inst/include/pest_host_use_table.hpp @@ -24,13 +24,29 @@ namespace pops { +/** + * Pest-host-use table holding susceptibilities, mortality rates, and mortality time + * lags for multiple hosts. + */ template class PestHostUseTable { public: using Environment = typename HostPool::Environment; + /** + * @brief Create an empty pest-host-use table + * + * @param environment Reference to the environment + */ PestHostUseTable(const Environment& environment) : environment_(environment) {} + + /** + * @brief Create a pest-host-use table using values in config + * + * @param config Configuration with pest-host-use table data + * @param environment Reference to the environment + */ PestHostUseTable(const Config& config, const Environment& environment) : environment_(environment) { @@ -41,6 +57,15 @@ class PestHostUseTable } } + /** + * @brief Add info for one host to the table + * + * Order of addition matters and should be the same as additions to the environment. + * + * @param susceptibility Host susceptibility + * @param mortality_rate Host mortality rate + * @param mortality_time_lag Host mortality time lag + */ void add_host_info(double susceptibility, double mortality_rate, int mortality_time_lag) { @@ -49,6 +74,11 @@ class PestHostUseTable mortality_time_lags_.push_back(mortality_time_lag); } + /** + * @brief Get susceptibility for the given host + * @param host Pointer to the host to get the information for + * @return Susceptibility score + */ double susceptibility(const HostPool* host) const { // This is using index because the environment is part of competency table, @@ -57,12 +87,22 @@ class PestHostUseTable return susceptibilities_.at(host_index); } + /** + * @brief Get mortality rate for the given host + * @param host Pointer to the host to get the information for + * @return Mortality rate value + */ double mortality_rate(const HostPool* host) const { auto host_index = environment_.host_index(host); return mortality_rates_.at(host_index); } + /** + * @brief Get mortality time lag for the given host + * @param host Pointer to the host to get the information for + * @return Mortality time lag value + */ double mortality_time_lag(const HostPool* host) const { auto host_index = environment_.host_index(host); @@ -70,10 +110,10 @@ class PestHostUseTable } private: - std::vector susceptibilities_; - std::vector mortality_rates_; - std::vector mortality_time_lags_; - const Environment& environment_; + std::vector susceptibilities_; ///< List of susceptibilities for hosts + std::vector mortality_rates_; ///< List of mortality_rates for hosts + std::vector mortality_time_lags_; ///< List of mortality time lags for hosts + const Environment& environment_; ///< Environment used for host indexing }; } // namespace pops diff --git a/inst/include/pest_pool.hpp b/inst/include/pest_pool.hpp index e91feb87..7eb3e6bb 100644 --- a/inst/include/pest_pool.hpp +++ b/inst/include/pest_pool.hpp @@ -71,6 +71,17 @@ class PestPool { return dispersers_(row, col); } + /** + * @brief Get raster of dispersers + * + * @note Most interaction should be done by dispersers_at(). + * + * @return Read-only reference to dispersers raster + */ + const IntegerRaster& dispersers() + { + return dispersers_; + } /** * @brief Set number of established dispersers * diff --git a/inst/include/simulation.hpp b/inst/include/simulation.hpp index 189b56fa..d01fc301 100644 --- a/inst/include/simulation.hpp +++ b/inst/include/simulation.hpp @@ -149,7 +149,7 @@ class Simulation * @return Const pointer to the environment * @throw std::logic_error when environment is not set */ - const Environment* + Environment* environment(bool allow_empty = false) { static Environment empty; @@ -211,6 +211,7 @@ class Simulation GeneratorProvider> remove(*environment(false), lethal_temperature); remove.action(hosts, generator); + this->environment(true)->remove_hosts(); } /** Removes percentage of exposed and infected @@ -258,6 +259,7 @@ class Simulation SurvivalRateAction survival( survival_rate); survival.action(hosts, generator); + this->environment(true)->remove_hosts(); } /** kills infected hosts based on mortality rate and timing. In the last year @@ -310,6 +312,7 @@ class Simulation Mortality mortality( mortality_rate, mortality_time_lag); mortality.action(hosts); + this->environment(true)->remove_hosts(); } /** Moves hosts from one location to another @@ -374,7 +377,9 @@ class Simulation 0, 0, suitable_cells}; - return host_movement.action(hosts, generator); + auto ret = host_movement.action(hosts, generator); + this->environment(true)->remove_hosts(); + return ret; } /** Generates dispersers based on infected @@ -434,6 +439,7 @@ class Simulation spread_action{unused_kernel}; spread_action.activate_soils(soil_pool_, to_soil_percentage_); spread_action.generate(host_pool, pests, generator); + this->environment(!weather)->remove_hosts(); } /** Creates dispersal locations for the dispersing individuals @@ -540,6 +546,7 @@ class Simulation spread_action{dispersal_kernel}; spread_action.activate_soils(soil_pool_, to_soil_percentage_); spread_action.disperse(host_pool, pests, generator); + this->environment(!weather)->remove_hosts(); } // For backwards compatibility for tests (without exposed and mortality) @@ -697,6 +704,7 @@ class Simulation rows_, cols_}; move_pest.action(hosts, pests, generator); + this->environment(true)->remove_hosts(); } /** Disperse, expose, and infect based on dispersers @@ -779,6 +787,7 @@ class Simulation cols_, suitable_cells}; host_pool.step_forward(step); + this->environment(!weather)->remove_hosts(); } template diff --git a/inst/include/utils.hpp b/inst/include/utils.hpp index 20653c63..0275e1f3 100644 --- a/inst/include/utils.hpp +++ b/inst/include/utils.hpp @@ -362,7 +362,7 @@ std::string quarantine_enum_to_string(Direction type) } /** - * Create a list of suitable cells in from a host raster. + * Create a list of suitable cells from a host raster. * * Suitable cell is defined as cell with value higher than zero, i.e., there is at least * one host. @@ -390,6 +390,11 @@ std::vector> find_suitable_cells(const RasterType& rast return cells; } +/** + * Create a list of suitable cells from a list of host rasters. + * + * @see find_suitable_cells(const RasterType&) + */ template std::vector> find_suitable_cells(const std::vector& rasters) diff --git a/src/pops.cpp b/src/pops.cpp index f3210198..c8c0f40b 100644 --- a/src/pops.cpp +++ b/src/pops.cpp @@ -196,7 +196,7 @@ List pops_model_cpp( std::vector> spread_rates_vector; std::tuple spread_rates; - IntegerMatrix dispersers(config.rows, config.cols); + IntegerMatrix dispersers(config.rows, config.cols, 0); IntegerMatrix total_dispersers(config.rows, config.cols); IntegerMatrix established_dispersers(config.rows, config.cols); @@ -225,7 +225,7 @@ List pops_model_cpp( PestHostUseTable pest_host_use_table( config, model.environment()); - CompetencyTable competency_table( + CompetencyTable competency_table( config, model.environment()); PoPSModel::StandardSingleHostPool host_pool( mt, @@ -247,7 +247,6 @@ List pops_model_cpp( config.cols, spatial_indices); - model.environment().add_host(&host_pool); competency_table.add_host_competencies({1}, 1); PoPSModel::StandardMultiHostPool multi_host_pool({&host_pool}, config); multi_host_pool.set_pest_host_use_table(pest_host_use_table); @@ -361,7 +360,6 @@ List pops_model_cpp( current_index, multi_host_pool, pest_pool, - dispersers, total_populations, treatments, temperature, From 4e97ffa0b1aa4f97a2e60f6bf8f679a1bea8bf66 Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Thu, 30 Nov 2023 19:03:23 -0500 Subject: [PATCH 03/12] fix --- src/pops.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pops.cpp b/src/pops.cpp index c8c0f40b..a00d4f4d 100644 --- a/src/pops.cpp +++ b/src/pops.cpp @@ -196,7 +196,7 @@ List pops_model_cpp( std::vector> spread_rates_vector; std::tuple spread_rates; - IntegerMatrix dispersers(config.rows, config.cols, 0); + IntegerMatrix dispersers(config.rows, config.cols); IntegerMatrix total_dispersers(config.rows, config.cols); IntegerMatrix established_dispersers(config.rows, config.cols); From b6f5d24c69703e072615c41e9252ce402a818d25 Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Thu, 30 Nov 2023 21:44:53 -0500 Subject: [PATCH 04/12] missing pest host table configuration --- src/pops.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/pops.cpp b/src/pops.cpp index a00d4f4d..f8a9b489 100644 --- a/src/pops.cpp +++ b/src/pops.cpp @@ -194,6 +194,8 @@ List pops_model_cpp( bool use_soils = bool_config["use_soils"]; config.dispersers_to_soils_percentage = dispersers_to_soils_percentage; + config.create_pest_host_use_table_from_parameters(1); + std::vector> spread_rates_vector; std::tuple spread_rates; IntegerMatrix dispersers(config.rows, config.cols); From e6591933bb2fd5fd0bae31ce099e40cc11be57f5 Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Fri, 1 Dec 2023 16:10:00 -0500 Subject: [PATCH 05/12] move code changing config before model is created --- src/pops.cpp | 80 +++++++++++++++++++++++++++------------------------- 1 file changed, 41 insertions(+), 39 deletions(-) diff --git a/src/pops.cpp b/src/pops.cpp index f8a9b489..b4fdd30b 100644 --- a/src/pops.cpp +++ b/src/pops.cpp @@ -223,6 +223,47 @@ List pops_model_cpp( ModelType mt = model_type_from_string(config.model_type); using PoPSModel = Model; + + Treatments treatments(config.scheduler()); + for (unsigned t = 0; t < treatment_maps.size(); t++) { + treatments.add_treatment( + treatment_maps[t], + pops::Date(treatment_dates[t]), + pesticide_duration[t], + treatment_application); + config.use_treatments = true; + } + + unsigned move_scheduled; + if (config.use_movements) { + for (unsigned move = 0; move < movements_dates.size(); ++move) { + pops::Date movement_date(movements_dates[move]); + move_scheduled = + unsigned(config.scheduler().schedule_action_date(movement_date)); + config.movement_schedule.push_back(move_scheduled); + } + } + + std::unique_ptr> network{nullptr}; + if (network_config.isNotNull() && network_data_config.isNotNull()) { + // The best place for bbox handling would be with rows, cols, and + // resolution, but since it is required only for network, it is here. + config.bbox.north = bbox["north"]; + config.bbox.south = bbox["south"]; + config.bbox.east = bbox["east"]; + config.bbox.west = bbox["west"]; + List net_config(network_config); + config.network_min_distance = net_config["network_min_distance"]; + config.network_max_distance = net_config["network_max_distance"]; + std::string network_movement = net_config["network_movement"]; + config.network_movement = network_movement; + network.reset(new Network(config.bbox, config.ew_res, config.ns_res)); + List net_data_config(network_data_config); + std::ifstream network_stream{ + Rcpp::as(net_data_config["network_filename"])}; + network->load(network_stream); + } + PoPSModel model(config); PestHostUseTable pest_host_use_table( @@ -258,16 +299,6 @@ List pops_model_cpp( established_dispersers, outside_dispersers); - Treatments treatments(config.scheduler()); - for (unsigned t = 0; t < treatment_maps.size(); t++) { - treatments.add_treatment( - treatment_maps[t], - pops::Date(treatment_dates[t]), - pesticide_duration[t], - treatment_application); - config.use_treatments = true; - } - if (config.use_lethal_temperature) { if (config.num_lethal() > temperature.size()) { Rcerr << "Not enough years of temperature data" << std::endl; @@ -288,15 +319,6 @@ List pops_model_cpp( config.ew_res, config.ns_res, spread_rate_outputs); - unsigned move_scheduled; - if (config.use_movements) { - for (unsigned move = 0; move < movements_dates.size(); ++move) { - pops::Date movement_date(movements_dates[move]); - move_scheduled = - unsigned(config.scheduler().schedule_action_date(movement_date)); - config.movement_schedule.push_back(move_scheduled); - } - } unsigned quarantine_outputs; if (config.use_quarantine) { @@ -319,26 +341,6 @@ List pops_model_cpp( Direction escape_direction; std::vector escape_directions; - std::unique_ptr> network{nullptr}; - if (network_config.isNotNull() && network_data_config.isNotNull()) { - // The best place for bbox handling would be with rows, cols, and - // resolution, but since it is required only for network, it is here. - config.bbox.north = bbox["north"]; - config.bbox.south = bbox["south"]; - config.bbox.east = bbox["east"]; - config.bbox.west = bbox["west"]; - List net_config(network_config); - config.network_min_distance = net_config["network_min_distance"]; - config.network_max_distance = net_config["network_max_distance"]; - std::string network_movement = net_config["network_movement"]; - config.network_movement = network_movement; - network.reset(new Network(config.bbox, config.ew_res, config.ns_res)); - List net_data_config(network_data_config); - std::ifstream network_stream{ - Rcpp::as(net_data_config["network_filename"])}; - network->load(network_stream); - } - WeatherType weather_typed = weather_type_from_string(config.weather_type); if (use_soils) { From 35889ac35a8793df11df7454b353e68212c2e227 Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Fri, 1 Dec 2023 16:33:08 -0500 Subject: [PATCH 06/12] update to latest pops-core --- inst/cpp/pops-core | 2 +- inst/include/competency_table.hpp | 36 +++++++-- inst/include/config.hpp | 79 ++++++++++++++----- inst/include/host_pool.hpp | 34 ++++---- inst/include/multi_host_pool.hpp | 8 +- ...host_use_table.hpp => pest_host_table.hpp} | 24 +++--- src/pops.cpp | 11 ++- 7 files changed, 129 insertions(+), 65 deletions(-) rename inst/include/{pest_host_use_table.hpp => pest_host_table.hpp} (82%) diff --git a/inst/cpp/pops-core b/inst/cpp/pops-core index f334c0d8..237f5770 160000 --- a/inst/cpp/pops-core +++ b/inst/cpp/pops-core @@ -1 +1 @@ -Subproject commit f334c0d8fad399acd93bd7335c92e798f94ef8e8 +Subproject commit 237f577010b369b6c17c28462921672792c9f6a6 diff --git a/inst/include/competency_table.hpp b/inst/include/competency_table.hpp index f98e5034..aef1f0f4 100644 --- a/inst/include/competency_table.hpp +++ b/inst/include/competency_table.hpp @@ -16,6 +16,7 @@ #ifndef POPS_COMPETENCY_TABLE_HPP #define POPS_COMPETENCY_TABLE_HPP +#include #include #include #include @@ -45,14 +46,28 @@ class CompetencyTable /** * @brief Create a competency table using values in config * + * If the table is complete as identified by Config::competency_table_is_complete(), + * it will use lookup using (ordered) std::map instead of the default linear search + * through the provided options. + * * @param config Configuration with competency table data * @param environment Reference to the environment */ CompetencyTable(const Config& config, const Environment& environment) : environment_(environment) { - for (const auto& row : config.competency_table_data()) { - competency_table_.emplace_back(row.presence_absence, row.competency); + if (config.competency_table_is_complete()) { + for (const auto& row : config.competency_table_data()) { + complete_competency_table_[row.presence_absence] = row.competency; + } + complete_ = true; + } + else { + for (const auto& row : config.competency_table_data()) { + partial_competency_table_.emplace_back( + row.presence_absence, row.competency); + } + complete_ = false; } } @@ -70,7 +85,8 @@ class CompetencyTable void add_host_competencies(const std::vector& presence_absence, double competency) { - competency_table_.emplace_back(presence_absence, competency); + partial_competency_table_.emplace_back(presence_absence, competency); + complete_ = false; } /** @@ -87,6 +103,9 @@ class CompetencyTable double competency_at(RasterIndex row, RasterIndex col, const HostPool* host) const { auto presence_absence = environment_.host_presence_at(row, col); + if (complete_) { + return complete_competency_table_.at(presence_absence); + } auto host_index = environment_.host_index(host); return find_competency(presence_absence, host_index); } @@ -109,7 +128,7 @@ class CompetencyTable // Go over all the rows and find the highest competency which fulfilled the // presence criteria. double competency = 0; - for (const auto& table_row : competency_table_) { + for (const auto& table_row : partial_competency_table_) { // probably faster if we just wait for the iteration over the whole thing if (!table_row.presence_absence[host_index]) continue; @@ -141,7 +160,7 @@ class CompetencyTable } /** - * @brief One row of the comptenecy table + * @brief One row of the competency table */ struct TableRow { @@ -153,7 +172,12 @@ class CompetencyTable : presence_absence(presence_absence), competency(competency) {} }; - std::vector competency_table_; ///< Internal table for lookup + /// Internal table for lookup when data is partial + std::vector partial_competency_table_; + /// Internal table for lookup when data is complete + std::map, double> complete_competency_table_; + /// true if the complete table should be used, false for the partial one + bool complete_{false}; const Environment& environment_; ///< Environment (for host index) }; diff --git a/inst/include/config.hpp b/inst/include/config.hpp index c51311bf..79ea58b1 100644 --- a/inst/include/config.hpp +++ b/inst/include/config.hpp @@ -27,6 +27,7 @@ #include "scheduling.hpp" #include "utils.hpp" +#include #include #include #include @@ -119,11 +120,11 @@ class Config { public: /** - * @brief Row of a table with pest-host-use data + * @brief Row of a table with pest-host data * * One row is stored for each host. */ - struct PestHostUseTableDataRow + struct PestHostTableDataRow { double susceptibility; ///< Susceptibility for a given host double mortality_rate; ///< Mortality rate for a given host @@ -190,12 +191,12 @@ class Config bool use_mortality{false}; std::string mortality_frequency; unsigned mortality_frequency_n; - /** Mortality rate if used without pest-host-use table + /** Mortality rate if used without pest-host table * - * @see read_pest_host_use_table() + * @see read_pest_host_table() */ double mortality_rate{0}; - /** Time lag of mortality in simulation steps if used without pest-host-use table */ + /** Time lag of mortality in simulation steps if used without pest-host table */ int mortality_time_lag{0}; // Quarantine bool use_quarantine{false}; @@ -534,14 +535,14 @@ class Config } /** - * @brief Get data for the pest-host-use table + * @brief Get data for the pest-host table * @return Reference to the internal table * - * @see PestHostUseTableDataRow + * @see PestHostTableDataRow */ - const std::vector& pest_host_use_table_data() const + const std::vector& pest_host_table_data() const { - return pest_host_use_table_data_; + return pest_host_table_data_; } /** @@ -556,30 +557,30 @@ class Config } /** - * @brief Read pest-host-use table data from vector of vectors of doubles + * @brief Read pest-host table data from vector of vectors of doubles * * The nested vectors need to be of size 3. The order of values is susceptibility, * mortality rate, and mortality time lag. * * @param values Table data */ - void read_pest_host_use_table(const std::vector>& values) + void read_pest_host_table(const std::vector>& values) { for (const auto& row : values) { if (row.size() < 3) { throw std::invalid_argument( - "3 values are required for each pest-host-use table row"); + "3 values are required for each pest-host table row"); } - PestHostUseTableDataRow resulting_row; + PestHostTableDataRow resulting_row; resulting_row.susceptibility = row[0]; resulting_row.mortality_rate = row[1]; resulting_row.mortality_time_lag = row[2]; - pest_host_use_table_data_.push_back(std::move(resulting_row)); + pest_host_table_data_.push_back(std::move(resulting_row)); } } /** - * @brief Use existing config parameters to create pest-host-use table + * @brief Use existing config parameters to create pest-host table * * This will create table with date for the given number of hosts with values for * all hosts being the same. Susceptibility is set to 1 and mortality is taken from @@ -590,14 +591,14 @@ class Config * @see #mortality_rate * @see #mortality_time_lag */ - void create_pest_host_use_table_from_parameters(int num_of_hosts) + void create_pest_host_table_from_parameters(int num_of_hosts) { for (int i = 0; i < num_of_hosts; ++i) { - PestHostUseTableDataRow resulting_row; + PestHostTableDataRow resulting_row; resulting_row.susceptibility = 1; resulting_row.mortality_rate = this->mortality_rate; resulting_row.mortality_time_lag = this->mortality_time_lag; - pest_host_use_table_data_.push_back(std::move(resulting_row)); + pest_host_table_data_.push_back(std::move(resulting_row)); } } @@ -625,13 +626,29 @@ class Config * ``` * * @param values Table data + * + * @throw std::invalid_argument when rows have different sizes + * @throw std::invalid_argument when row size is less then 2 */ void read_competency_table(const std::vector>& values) { + size_t first_row_size{0}; + bool first_row{true}; for (const auto& row : values) { + if (!first_row && row.size() != first_row_size) { + throw std::invalid_argument( + "All competency table rows must be the same size (" + + std::to_string(row.size()) + + " != " + std::to_string(first_row_size) + ")"); + } + else { + first_row_size = row.size(); + first_row = false; + } if (row.size() < 2) { throw std::invalid_argument( - "At least 2 values are required for each competency table row"); + "At least 2 values are required for each competency table row (not " + + std::to_string(row.size()) + ")"); } CompetencyTableDataRow resulting_row; for (auto it = row.begin(); it < std::prev(row.end()); ++it) { @@ -642,6 +659,26 @@ class Config } } + /** + * @brief Test if the competency table is complete + * + * Complete comptenecy table has 2^N rows where N is number of hosts, i.e., number + * of columns used for the host presence-absence information. + * + * @return true if complete, false otherwise + */ + bool competency_table_is_complete() const + { + size_t num_of_rows{competency_table_data().size()}; + if (num_of_rows == 0) + return false; + // Number of presence-absence records is assumed to be number of hosts. + size_t presence_size{competency_table_data().at(0).presence_absence.size()}; + if (num_of_rows == std::pow(2, presence_size)) + return true; + return false; + } + private: Date date_start_{"0-01-01"}; Date date_end_{"0-01-02"}; @@ -666,8 +703,8 @@ class Config std::vector quarantine_schedule_; std::vector weather_table_; - /** Storage for the pest-host-use table data */ - std::vector pest_host_use_table_data_; + /** Storage for the pest-host table data */ + std::vector pest_host_table_data_; /** Storage for the competency table data */ std::vector competency_table_data_; }; diff --git a/inst/include/host_pool.hpp b/inst/include/host_pool.hpp index 7c55bcf4..649ce2cc 100644 --- a/inst/include/host_pool.hpp +++ b/inst/include/host_pool.hpp @@ -25,7 +25,7 @@ #include "model_type.hpp" #include "environment_interface.hpp" #include "competency_table.hpp" -#include "pest_host_use_table.hpp" +#include "pest_host_table.hpp" namespace pops { @@ -148,7 +148,7 @@ class HostPool : public HostPoolInterface } /** - * @brief Set pest-host-use table + * @brief Set pest-host table * * If set, apply_mortality_at(RasterIndex, RasterIndex) can be used instead of the * version which takes mortality parameters directly. Susceptibility is used @@ -157,11 +157,11 @@ class HostPool : public HostPoolInterface * Pointer to the existing object is stored and used. So, the table can be modified, * but the table object needs to exists during the lifetime of this object. * - * @param pest_host_use_table Reference to the table + * @param table Reference to the table */ - void set_pest_host_use_table(const PestHostUseTable& pest_host_use_table) + void set_pest_host_table(const PestHostTable& table) { - this->pest_host_use_table_ = &pest_host_use_table; + this->pest_host_table_ = &table; } /** @@ -173,11 +173,11 @@ class HostPool : public HostPoolInterface * Pointer to the existing object is stored and used. So, the table can be modified, * but the table object needs to exists during the lifetime of this object. * - * @param competency_table Reference to the table + * @param table Reference to the table */ - void set_competency_table(const CompetencyTable& competency_table) + void set_competency_table(const CompetencyTable& table) { - this->competency_table_ = &competency_table; + this->competency_table_ = &table; } /** @@ -324,8 +324,8 @@ class HostPool : public HostPoolInterface double probability_of_establishment = (double)(susceptible_(row, col)) / environment_.total_population_at(row, col); - if (pest_host_use_table_) { - probability_of_establishment *= pest_host_use_table_->susceptibility(this); + if (pest_host_table_) { + probability_of_establishment *= pest_host_table_->susceptibility(this); } return environment_.influence_probability_of_establishment_at( row, col, probability_of_establishment); @@ -851,7 +851,7 @@ class HostPool : public HostPoolInterface /** * @brief Apply mortality at a given cell * - * Uses pest-host-use table for mortality parameters. + * Uses pest-host table for mortality parameters. * * @param row Row index of the cell * @param col Column index of the cell @@ -860,16 +860,16 @@ class HostPool : public HostPoolInterface */ void apply_mortality_at(RasterIndex row, RasterIndex col) { - if (!pest_host_use_table_) { + if (!pest_host_table_) { throw std::invalid_argument( - "Set pest-host-use table before calling apply_mortality_at " + "Set pest-host table before calling apply_mortality_at " "or provide parameters in the function call"); } this->apply_mortality_at( row, col, - pest_host_use_table_->mortality_rate(this), - pest_host_use_table_->mortality_time_lag(this)); + pest_host_table_->mortality_rate(this), + pest_host_table_->mortality_time_lag(this)); } /** @@ -1157,8 +1157,8 @@ class HostPool : public HostPoolInterface bool establishment_stochasticity_{true}; double deterministic_establishment_probability_{0}; - /** Pest-host-use table */ - const PestHostUseTable* pest_host_use_table_{nullptr}; + /** pest-host table */ + const PestHostTable* pest_host_table_{nullptr}; /** Competency table */ const CompetencyTable* competency_table_{nullptr}; diff --git a/inst/include/multi_host_pool.hpp b/inst/include/multi_host_pool.hpp index 280bc70d..548d9a8d 100644 --- a/inst/include/multi_host_pool.hpp +++ b/inst/include/multi_host_pool.hpp @@ -22,7 +22,7 @@ #include "competency_table.hpp" #include "config.hpp" -#include "pest_host_use_table.hpp" +#include "pest_host_table.hpp" #include "utils.hpp" namespace pops { @@ -59,16 +59,16 @@ class MultiHostPool {} /** - * @brief Set pest-host-use table for all hosts + * @brief Set pest-host table for all hosts * * The existing object will be used (not copy is performed). * * @param table Reference to the table object */ - void set_pest_host_use_table(const PestHostUseTable& table) + void set_pest_host_table(const PestHostTable& table) { for (auto& host_pool : host_pools_) { - host_pool->set_pest_host_use_table(table); + host_pool->set_pest_host_table(table); } } diff --git a/inst/include/pest_host_use_table.hpp b/inst/include/pest_host_table.hpp similarity index 82% rename from inst/include/pest_host_use_table.hpp rename to inst/include/pest_host_table.hpp index d2879a11..0d70ba3e 100644 --- a/inst/include/pest_host_use_table.hpp +++ b/inst/include/pest_host_table.hpp @@ -1,5 +1,5 @@ /* - * PoPS model - Pest-host-use table for hosts and pest + * PoPS model - Pest-host table for hosts and pest * * Copyright (C) 2023 by the authors. * @@ -13,8 +13,8 @@ * http://www.gnu.org/copyleft/gpl.html */ -#ifndef POPS_PEST_HOST_USE_TABLE_HPP -#define POPS_PEST_HOST_USE_TABLE_HPP +#ifndef POPS_PEST_HOST_TABLE_HPP +#define POPS_PEST_HOST_TABLE_HPP #include #include @@ -25,32 +25,32 @@ namespace pops { /** - * Pest-host-use table holding susceptibilities, mortality rates, and mortality time + * Pest-host table holding susceptibilities, mortality rates, and mortality time * lags for multiple hosts. */ template -class PestHostUseTable +class PestHostTable { public: using Environment = typename HostPool::Environment; /** - * @brief Create an empty pest-host-use table + * @brief Create an empty pest-host table * * @param environment Reference to the environment */ - PestHostUseTable(const Environment& environment) : environment_(environment) {} + PestHostTable(const Environment& environment) : environment_(environment) {} /** - * @brief Create a pest-host-use table using values in config + * @brief Create a pest-host table using values in config * - * @param config Configuration with pest-host-use table data + * @param config Configuration with pest-host table data * @param environment Reference to the environment */ - PestHostUseTable(const Config& config, const Environment& environment) + PestHostTable(const Config& config, const Environment& environment) : environment_(environment) { - for (const auto& row : config.pest_host_use_table_data()) { + for (const auto& row : config.pest_host_table_data()) { susceptibilities_.push_back(row.susceptibility); mortality_rates_.push_back(row.mortality_rate); mortality_time_lags_.push_back(row.mortality_time_lag); @@ -118,4 +118,4 @@ class PestHostUseTable } // namespace pops -#endif // POPS_PEST_HOST_USE_TABLE_HPP +#endif // POPS_PEST_HOST_TABLE_HPP diff --git a/src/pops.cpp b/src/pops.cpp index b4fdd30b..28158f9e 100644 --- a/src/pops.cpp +++ b/src/pops.cpp @@ -194,7 +194,7 @@ List pops_model_cpp( bool use_soils = bool_config["use_soils"]; config.dispersers_to_soils_percentage = dispersers_to_soils_percentage; - config.create_pest_host_use_table_from_parameters(1); + config.create_pest_host_table_from_parameters(1); std::vector> spread_rates_vector; std::tuple spread_rates; @@ -263,10 +263,14 @@ List pops_model_cpp( Rcpp::as(net_data_config["network_filename"])}; network->load(network_stream); } + std::vector> competency_table_data; + competency_table_data.push_back({1, 1}); + competency_table_data.push_back({0, 0}); + config.read_competency_table(competency_table_data); PoPSModel model(config); - PestHostUseTable pest_host_use_table( + PestHostTable pest_host_table( config, model.environment()); CompetencyTable competency_table( config, model.environment()); @@ -290,9 +294,8 @@ List pops_model_cpp( config.cols, spatial_indices); - competency_table.add_host_competencies({1}, 1); PoPSModel::StandardMultiHostPool multi_host_pool({&host_pool}, config); - multi_host_pool.set_pest_host_use_table(pest_host_use_table); + multi_host_pool.set_pest_host_table(pest_host_table); multi_host_pool.set_competency_table(competency_table); PoPSModel::StandardPestPool pest_pool( dispersers, From fb9604d1f64971116f351af2341808135d748b7c Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Thu, 14 Dec 2023 09:43:12 -0500 Subject: [PATCH 07/12] update to latest core with random number fix, remove left over file, fix indentation --- inst/cpp/pops-core | 2 +- inst/include/mortality.hpp | 82 -------------------------------- inst/include/multi_host_pool.hpp | 11 +++++ src/pops.cpp | 6 +-- 4 files changed, 15 insertions(+), 86 deletions(-) delete mode 100644 inst/include/mortality.hpp diff --git a/inst/cpp/pops-core b/inst/cpp/pops-core index 237f5770..e791d971 160000 --- a/inst/cpp/pops-core +++ b/inst/cpp/pops-core @@ -1 +1 @@ -Subproject commit 237f577010b369b6c17c28462921672792c9f6a6 +Subproject commit e791d971f65136e11695f1adc5660efd24aa65a8 diff --git a/inst/include/mortality.hpp b/inst/include/mortality.hpp deleted file mode 100644 index 2f646593..00000000 --- a/inst/include/mortality.hpp +++ /dev/null @@ -1,82 +0,0 @@ -/* - * PoPS model - mortality - * - * Copyright (C) 2015-2020 by the authors. - * - * Authors: Chris Jones >cjones1688 gmail com - * Anna Petrasova - * Vaclav Petras - * - * The code contained herein is licensed under the GNU General Public - * License. You may obtain a copy of the GNU General Public License - * Version 2 or later at the following locations: - * - * http://www.opensource.org/licenses/gpl-license.html - * http://www.gnu.org/copyleft/gpl.html -*/ - -#ifndef POPS_MORTALITY_HPP -#define POPS_MORTALITY_HPP - -#include "raster.hpp" -#include "date.hpp" -#include "scheduling.hpp" - -#include -#include -#include -#include - -namespace pops { - - -void mortality( - IntegerRaster& infected, - double mortality_rate, - int current_year, - int first_mortality_year, - IntegerRaster& mortality, - std::vector& mortality_tracker_vector) -{ - if (current_step >= (first_mortality_year)) { - - for (int i = 0; i < rows_; i++) { - for (int j = 0; j < cols_; j++) { - for (int index = 0; index <= max_index; - index++) { - int mortality_in_index = 0; - if (mortality_tracker_vector[index](i, j) > 0) { - mortality_in_index = - mortality_rate - * mortality_tracker_vector[index](i, j); - mortality_tracker_vector[index](i, j) -= - mortality_in_index; - mortality(i, j) += mortality_in_index; - mortality_current_year += mortality_in_index; - if (infected(i, j) > 0) { - infected(i, j) -= mortality_in_index; - } - } - } - } - } - } -} - -if (model_type_ == ModelType::SusceptibleExposedInfected) { - if (step >= latency_period_) { - // Oldest item needs to be in the front - auto& oldest = exposed.front(); - // Move hosts to infected raster - infected += oldest; - mortality_tracker += oldest; - // Reset the raster - // (hosts moved from the raster) - oldest.fill(0); - } - // Age the items and the used one to the back - // elements go one position to the left - // new oldest goes to the front - // old oldest goes to the back - rotate_left_by_one(exposed); -} \ No newline at end of file diff --git a/inst/include/multi_host_pool.hpp b/inst/include/multi_host_pool.hpp index 548d9a8d..e325312d 100644 --- a/inst/include/multi_host_pool.hpp +++ b/inst/include/multi_host_pool.hpp @@ -430,17 +430,28 @@ class MultiHostPool * The probabilities are used as weights so they don't need to be 0-1 nor they need * to add up to 1. However, their sum needs to be >0. * + * If there is only one hosts, it returns that host without using the random number + * generator. + * * @param hosts List of pointers to host pools * @param probabilities Probability values for each host pool * @param generator Random number generator * * @return Pointer to selected host pool + * + * @throw std::invalid_argument if *hosts* is empty. */ static HostPool* pick_host_by_probability( std::vector& hosts, const std::vector& probabilities, Generator& generator) { + if (!hosts.size()) { + std::invalid_argument("List of hosts is empty in multi host pool"); + } + if (hosts.size() == 1) { + return hosts[0]; + } std::discrete_distribution distribution{ probabilities.begin(), probabilities.end()}; return hosts.at(distribution(generator)); diff --git a/src/pops.cpp b/src/pops.cpp index 28158f9e..6764dbec 100644 --- a/src/pops.cpp +++ b/src/pops.cpp @@ -228,9 +228,9 @@ List pops_model_cpp( for (unsigned t = 0; t < treatment_maps.size(); t++) { treatments.add_treatment( treatment_maps[t], - pops::Date(treatment_dates[t]), - pesticide_duration[t], - treatment_application); + pops::Date(treatment_dates[t]), + pesticide_duration[t], + treatment_application); config.use_treatments = true; } From 558146fc97e2970c6913396d1a66f303e3abe8a2 Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Thu, 14 Dec 2023 14:40:57 -0500 Subject: [PATCH 08/12] sync to main pops-core --- inst/cpp/pops-core | 2 +- inst/include/environment.hpp | 2 +- inst/include/environment_interface.hpp | 4 +- inst/include/host_pool.hpp | 32 ++++++------ inst/include/multi_host_pool.hpp | 67 +++++++++++++------------- inst/include/simulation.hpp | 7 ++- 6 files changed, 59 insertions(+), 55 deletions(-) diff --git a/inst/cpp/pops-core b/inst/cpp/pops-core index e791d971..19eab166 160000 --- a/inst/cpp/pops-core +++ b/inst/cpp/pops-core @@ -1 +1 @@ -Subproject commit e791d971f65136e11695f1adc5660efd24aa65a8 +Subproject commit 19eab166224bb5f75a285ec2d96920ec8e2d3f4c diff --git a/inst/include/environment.hpp b/inst/include/environment.hpp index 7cb6da4a..97c66035 100644 --- a/inst/include/environment.hpp +++ b/inst/include/environment.hpp @@ -179,7 +179,7 @@ class Environment return value * weather_coefficient_at(row, col); } - double influence_probability_of_establishment_at( + double influence_suitability_at( RasterIndex row, RasterIndex col, double value) const override { if (!weather_) diff --git a/inst/include/environment_interface.hpp b/inst/include/environment_interface.hpp index ea804617..8990e9e3 100644 --- a/inst/include/environment_interface.hpp +++ b/inst/include/environment_interface.hpp @@ -35,8 +35,8 @@ class EnvironmentInterface virtual double weather_coefficient_at(RasterIndex row, RasterIndex col) const = 0; virtual double influence_reproductive_rate_at( RasterIndex row, RasterIndex col, double value) const = 0; - virtual double influence_probability_of_establishment_at( - RasterIndex row, RasterIndex col, double value) const = 0; + virtual double + influence_suitability_at(RasterIndex row, RasterIndex col, double value) const = 0; virtual int total_population_at(RasterIndex row, RasterIndex col) const = 0; /** * @brief Get presence-absence for host at a given cell diff --git a/inst/include/host_pool.hpp b/inst/include/host_pool.hpp index 649ce2cc..2c600b31 100644 --- a/inst/include/host_pool.hpp +++ b/inst/include/host_pool.hpp @@ -37,7 +37,7 @@ namespace pops { * @tparam RasterIndex Type for indexing the rasters * @tparam GeneratorProvider Provider of random number generators * - * GeneratorProvider needs to provide Generator memeber which is the type of the + * GeneratorProvider needs to provide Generator member which is the type of the * underlying random number generators. */ template< @@ -202,7 +202,7 @@ class HostPool : public HostPoolInterface { if (susceptible_(row, col) <= 0) return 0; - double probability_of_establishment = establishment_probability_at(row, col); + double probability_of_establishment = suitability_at(row, col); bool establish = can_disperser_establish( probability_of_establishment, establishment_stochasticity_, @@ -216,7 +216,7 @@ class HostPool : public HostPoolInterface /** * @brief Test whether a disperser can establish * - * This static (object-idependent) function to test disperser establishement allows + * This static (object-independent) function to test disperser establishment allows * code reuse between a single host and a multi-host case. * * @param probability_of_establishment Probability of establishment @@ -312,23 +312,21 @@ class HostPool : public HostPoolInterface } /** - * @brief Get establishment probability for a cell + * @brief Get suitability score for a cell * * @param row Row index of the cell * @param col Column index of the cell * - * @return Establishment probability + * @return suitability score */ - double establishment_probability_at(RasterIndex row, RasterIndex col) const + double suitability_at(RasterIndex row, RasterIndex col) const { - double probability_of_establishment = - (double)(susceptible_(row, col)) - / environment_.total_population_at(row, col); + double suitability = (double)(susceptible_(row, col)) + / environment_.total_population_at(row, col); if (pest_host_table_) { - probability_of_establishment *= pest_host_table_->susceptibility(this); + suitability *= pest_host_table_->susceptibility(this); } - return environment_.influence_probability_of_establishment_at( - row, col, probability_of_establishment); + return environment_.influence_suitability_at(row, col, suitability); } /** @@ -349,7 +347,7 @@ class HostPool : public HostPoolInterface * * @return Number of pests actually moved from the cell * - * @note For consitency with the previous implementation, this does not modify + * @note For consistency with the previous implementation, this does not modify * mortality cohorts nor touches the exposed cohorts. */ int pests_from(RasterIndex row, RasterIndex col, int count, Generator& generator) @@ -505,7 +503,7 @@ class HostPool : public HostPoolInterface // Returned total hosts actually moved is based only on the total host and no // other checks are performed. This assumes that the counts are correct in the - // object (precodition). + // object (precondition). return total_hosts_moved; } @@ -522,7 +520,7 @@ class HostPool : public HostPoolInterface * @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 threatment code. + * here in the same way as in the original treatment code. * * @note This does not remove resistant just like the original implementation in * treatments. @@ -763,7 +761,7 @@ class HostPool : public HostPoolInterface // fail. if (false && infected != mortality_total) { throw std::invalid_argument( - "Total of mortality values differs from formely infected, now resistant " + "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) @@ -1103,7 +1101,7 @@ class HostPool : public HostPoolInterface } /** - * @brief Get list with this host pools + * @brief Get list which contains this host pool * * This is for compatibility with multi-host pool. In case of this host pool, it * returns one item which is a pointer to this host pool. diff --git a/inst/include/multi_host_pool.hpp b/inst/include/multi_host_pool.hpp index e325312d..9519a702 100644 --- a/inst/include/multi_host_pool.hpp +++ b/inst/include/multi_host_pool.hpp @@ -31,7 +31,7 @@ namespace pops { * Host pool for multiple hosts * * Keeps the same interface as (single) HostPool given that HostPool has methods to - * behave in a mutli-host way if needed. This allows most external operations such as + * behave in a multi-host way if needed. This allows most external operations such as * actions to work without distinguishing single- and multi-host pools. */ template< @@ -313,7 +313,7 @@ class MultiHostPool * * Config::arrival_behavior() is used to determine if the dispersers should be * landing in a cell (`"land"`) or infecting a specific host directly (`"infect"`). - * Landing performs the establishment test on a combined probability from all hosts + * Landing performs the establishment test on a combined suitability from all hosts * in a cell and then uses one host to add the new disperser to. Infecting picks a * host and lets the host accept or reject the disperser based on its own * establishment test. @@ -321,8 +321,7 @@ class MultiHostPool * Uses static function HostPool::can_disperser_establish() from HostPool class to * do the test for multiple hosts for landing, so any host pool class needs to * implement that besides its disperser methods. It assumes that the configuration - * for multi-host is applicable for the combined establishment probability of - * individual hosts. + * for multi-host is applicable for the combined suitability of individual hosts. * * For one single host pool only and host pool which implements its `disperser_to()` * using `can_disperser_establish()` and `add_disperser_at()`, this gives identical @@ -338,33 +337,36 @@ class MultiHostPool */ int disperser_to(RasterIndex row, RasterIndex col, Generator& generator) { - std::vector probabilities; - double total_s_score = 0; + std::vector suitabilities; + double total_suitability_score = 0; for (auto& host_pool : host_pools_) { - // The probability accounts for weather and, importantly, number of + // The suitability accounts for weather and, importantly, number of // susceptible hosts, so host pool without available hosts in a given cell // is less likely to be selected over a host pool with available hosts // (however, it can happen in case zeros are not exactly zeros in the // discrete distribution used later and code checks can't tell, so we need // to account for that case later anyway). - double s_for_item = host_pool->establishment_probability_at(row, col); - // The resulting s can be 0-1. While the probabilities are used as weights - // for picking the host, so their absolute range does not matter, the total - // is used as probablity in a stochastic test. The stochastic challenge may - // be against 0.5 + 0.7 = 1.2 which which is fine as long as we are fine - // with probabilities 0.5 and 0.7 from two host translating to 100% - // probability. - probabilities.push_back(s_for_item); - total_s_score += s_for_item; + double suitability = host_pool->suitability_at(row, col); + // The resulting individual suitability can be 0-1. The individual + // suitabilities are used as weights for picking the host, so their absolute + // range does not matter. The total is used in a stochastic test and should + // be <=1 which should be ensured in the input data. + suitabilities.push_back(suitability); + total_suitability_score += suitability; } - if (total_s_score <= 0) { - // While the score should always be >= 0, it may be == 0 if no hosts are - // present. No hosts present cause all probabilities to be zero which is not - // permissible for the host picking later and it is enough information for - // us to know there won't be any establishment. + if (total_suitability_score <= 0) { + // While the score should always be >= 0, it may be == 0 if no hosts (host + // individuals) are present. No hosts present cause all suitabilities to be + // zero which is not permissible for the host picking later and it is enough + // information for us to know there won't be any establishment. return 0; } - auto host = pick_host_by_probability(host_pools_, probabilities, generator); + if (total_suitability_score > 1) { + throw std::invalid_argument( + "Total suitability score is " + std::to_string(total_suitability_score) + + " but it needs to be <=1"); + } + auto host = pick_host_by_weight(host_pools_, suitabilities, generator); if (config_.arrival_behavior() == "land") { // The operations are ordered so that for single host, this gives an // identical results to the infect behavior (influenced by usage of random @@ -372,7 +374,7 @@ class MultiHostPool if (host->susceptible_at(row, col) <= 0) return 0; bool establish = HostPool::can_disperser_establish( - total_s_score, + total_suitability_score, config_.establishment_stochasticity, config_.establishment_probability, generator); @@ -389,7 +391,7 @@ class MultiHostPool /** * @brief Move hosts from a cell to a cell (multi-host) * - * Currenty just works for first host. + * Currently just works for first host. * * @param row_from Row index of the source cell * @param col_from Column index of the source cell @@ -425,25 +427,25 @@ class MultiHostPool private: /** - * @brief Pick a host given a probability for each host + * @brief Pick a host given a weight for each host * - * The probabilities are used as weights so they don't need to be 0-1 nor they need - * to add up to 1. However, their sum needs to be >0. + * The weights don't need to be 0-1 nor they need to add up to 1. However, their sum + * needs to be >0. * - * If there is only one hosts, it returns that host without using the random number + * If there is only one host, it returns that host without using the random number * generator. * * @param hosts List of pointers to host pools - * @param probabilities Probability values for each host pool + * @param weights Weight values for each host pool * @param generator Random number generator * * @return Pointer to selected host pool * * @throw std::invalid_argument if *hosts* is empty. */ - static HostPool* pick_host_by_probability( + static HostPool* pick_host_by_weight( std::vector& hosts, - const std::vector& probabilities, + const std::vector& weights, Generator& generator) { if (!hosts.size()) { @@ -452,8 +454,7 @@ class MultiHostPool if (hosts.size() == 1) { return hosts[0]; } - std::discrete_distribution distribution{ - probabilities.begin(), probabilities.end()}; + std::discrete_distribution distribution{weights.begin(), weights.end()}; return hosts.at(distribution(generator)); } diff --git a/inst/include/simulation.hpp b/inst/include/simulation.hpp index d01fc301..582185bb 100644 --- a/inst/include/simulation.hpp +++ b/inst/include/simulation.hpp @@ -33,7 +33,12 @@ namespace pops { -/*! The main class to control the spread simulation. +/*! 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. * * The Simulation class handles the mechanics of the model, but the * timing of the events or steps should be handled outside of this From 3c7cf632c29d76d95a6accf712e859799bdeb9c6 Mon Sep 17 00:00:00 2001 From: ChrisJones687 Date: Fri, 15 Dec 2023 15:43:16 -0500 Subject: [PATCH 09/12] Add check that hosts aren't greater than total populations when being drawn from a distribution. --- R/calibrate.R | 3 +++ R/helpers.R | 2 +- R/pops.r | 3 +++ R/pops_multirun.R | 3 +++ R/validate.R | 3 +++ 5 files changed, 13 insertions(+), 1 deletion(-) diff --git a/R/calibrate.R b/R/calibrate.R index cb068681..cefd03bf 100644 --- a/R/calibrate.R +++ b/R/calibrate.R @@ -341,6 +341,9 @@ calibrate <- function(infected_years_file, if (config$use_host_uncertainty) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) + while (all(config$host > config$total_populations)) { + config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) + } } else { config$host <- config$host_mean } diff --git a/R/helpers.R b/R/helpers.R index b5059f1c..b702de46 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -247,7 +247,7 @@ create_random_seeds <- function(n) { matrix_norm_distribution <- function(mean_matrix, sd_matrix) { new_matrix <- round(matrix(mapply(function(x, y) {rnorm(x, y, n = 1)}, x = mean_matrix, y = sd_matrix), - nrow(mean_matrix), ncol(mean_matrix))) + nrow(mean_matrix), ncol(mean_matrix)), digits = 0) new_matrix[is.na(new_matrix)] <- 0 new_matrix[new_matrix < 0] <- 0 return(new_matrix) diff --git a/R/pops.r b/R/pops.r index e6d71aaf..60ecffdb 100644 --- a/R/pops.r +++ b/R/pops.r @@ -384,6 +384,9 @@ pops <- function(infected_file, if (config$use_host_uncertainty) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) + while (all(config$host > config$total_populations)) { + config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) + } } else { config$host <- config$host_mean } diff --git a/R/pops_multirun.R b/R/pops_multirun.R index 96c4cb66..6abf0696 100644 --- a/R/pops_multirun.R +++ b/R/pops_multirun.R @@ -239,6 +239,9 @@ pops_multirun <- function(infected_file, if (config$use_host_uncertainty) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) + while (all(config$host > config$total_populations)) { + config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) + } } else { config$host <- config$host_mean } diff --git a/R/validate.R b/R/validate.R index de8d840a..55d691f7 100644 --- a/R/validate.R +++ b/R/validate.R @@ -258,6 +258,9 @@ validate <- function(infected_years_file, if (config$use_host_uncertainty) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) + while (all(config$host > config$total_populations)) { + config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) + } } else { config$host <- config$host_mean } From 6ba83d0d2adbfeabea6139d0f46f1636ded85123 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 18 Dec 2023 10:23:20 -0500 Subject: [PATCH 10/12] fixed hosts being drawn greater than total_populations --- R/calibrate.R | 8 +++++++- R/pops.r | 8 +++++++- R/pops_multirun.R | 8 +++++++- R/validate.R | 8 +++++++- 4 files changed, 28 insertions(+), 4 deletions(-) diff --git a/R/calibrate.R b/R/calibrate.R index cefd03bf..506fe866 100644 --- a/R/calibrate.R +++ b/R/calibrate.R @@ -327,7 +327,13 @@ calibrate <- function(infected_years_file, random_seeds <- create_random_seeds(1) if (config$use_initial_condition_uncertainty) { config$infected <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + while (any(config$infected < 0)) { + config$infected <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + } exposed2 <- matrix_norm_distribution(config$exposed_mean, config$exposed_sd) + while (any(exposed2 < 0)) { + exposed2 <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + } exposed <- config$exposed exposed[[config$latency_period + 1]] <- exposed2 config$exposed <- exposed @@ -341,7 +347,7 @@ calibrate <- function(infected_years_file, if (config$use_host_uncertainty) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) - while (all(config$host > config$total_populations)) { + while (any(config$host > config$total_populations)) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) } } else { diff --git a/R/pops.r b/R/pops.r index 60ecffdb..d4a61a97 100644 --- a/R/pops.r +++ b/R/pops.r @@ -370,7 +370,13 @@ pops <- function(infected_file, if (config$use_initial_condition_uncertainty) { config$infected <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + while (any(config$infected < 0)) { + config$infected <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + } exposed2 <- matrix_norm_distribution(config$exposed_mean, config$exposed_sd) + while (any(exposed2 < 0)) { + exposed2 <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + } exposed <- config$exposed exposed[[config$latency_period + 1]] <- exposed2 config$exposed <- exposed @@ -384,7 +390,7 @@ pops <- function(infected_file, if (config$use_host_uncertainty) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) - while (all(config$host > config$total_populations)) { + while (any(config$host > config$total_populations)) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) } } else { diff --git a/R/pops_multirun.R b/R/pops_multirun.R index 6abf0696..ff723b58 100644 --- a/R/pops_multirun.R +++ b/R/pops_multirun.R @@ -225,7 +225,13 @@ pops_multirun <- function(infected_file, if (config$use_initial_condition_uncertainty) { config$infected <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + while (any(config$infected < 0)) { + config$infected <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + } exposed2 <- matrix_norm_distribution(config$exposed_mean, config$exposed_sd) + while (any(exposed2 < 0)) { + exposed2 <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + } exposed <- config$exposed exposed[[config$latency_period + 1]] <- exposed2 config$exposed <- exposed @@ -239,7 +245,7 @@ pops_multirun <- function(infected_file, if (config$use_host_uncertainty) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) - while (all(config$host > config$total_populations)) { + while (any(config$host > config$total_populations)) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) } } else { diff --git a/R/validate.R b/R/validate.R index 55d691f7..8fb4e7bb 100644 --- a/R/validate.R +++ b/R/validate.R @@ -244,7 +244,13 @@ validate <- function(infected_years_file, if (config$use_initial_condition_uncertainty) { config$infected <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + while (any(config$infected < 0)) { + config$infected <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + } exposed2 <- matrix_norm_distribution(config$exposed_mean, config$exposed_sd) + while (any(exposed2 < 0)) { + exposed2 <- matrix_norm_distribution(config$infected_mean, config$infected_sd) + } exposed <- config$exposed exposed[[config$latency_period + 1]] <- exposed2 config$exposed <- exposed @@ -258,7 +264,7 @@ validate <- function(infected_years_file, if (config$use_host_uncertainty) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) - while (all(config$host > config$total_populations)) { + while (any(config$host > config$total_populations)) { config$host <- matrix_norm_distribution(config$host_mean, config$host_sd) } } else { From 4602b877fd0fbbedbfe1dcaeb654c50e087e9653 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 18 Dec 2023 10:25:15 -0500 Subject: [PATCH 11/12] update multihost test --- tests/testthat/test-pops-multirun.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-pops-multirun.R b/tests/testthat/test-pops-multirun.R index a874a486..471dcd44 100644 --- a/tests/testthat/test-pops-multirun.R +++ b/tests/testthat/test-pops-multirun.R @@ -507,7 +507,7 @@ test_that("Multirun model outputs work with mask", { end_date <- "2019-12-31" lethal_temperature <- -35 lethal_temperature_month <- 1 - random_seed <- 42 + random_seed <- NA treatments_file <- "" treatment_dates <- c("2019-11-01") treatment_method <- "ratio" @@ -531,7 +531,7 @@ test_that("Multirun model outputs work with mask", { number_of_cores <- 2 model_type <- "SI" latency_period <- 0 - parameter_means <- c(0, 21, 1, 500, 0, 0, 0, 0) + parameter_means <- c(2, 21, 1, 500, 0, 0, 0, 0) parameter_cov_matrix <- matrix(0, nrow = 8, ncol = 8) start_exposed <- FALSE generate_stochasticity <- TRUE From cff1f132617afafba383dcf14d3ff70e688ab202 Mon Sep 17 00:00:00 2001 From: Chris Jones Date: Mon, 18 Dec 2023 11:16:15 -0500 Subject: [PATCH 12/12] don't change the test value --- tests/testthat/test-pops-multirun.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-pops-multirun.R b/tests/testthat/test-pops-multirun.R index 471dcd44..cee9b685 100644 --- a/tests/testthat/test-pops-multirun.R +++ b/tests/testthat/test-pops-multirun.R @@ -531,7 +531,7 @@ test_that("Multirun model outputs work with mask", { number_of_cores <- 2 model_type <- "SI" latency_period <- 0 - parameter_means <- c(2, 21, 1, 500, 0, 0, 0, 0) + parameter_means <- c(0, 21, 1, 500, 0, 0, 0, 0) parameter_cov_matrix <- matrix(0, nrow = 8, ncol = 8) start_exposed <- FALSE generate_stochasticity <- TRUE