simulation.hpp

Go to the documentation of this file.

/*
 * PoPS model - pest or pathogen spread simulation
 *
 * Copyright (C) 2015-2020 by the authors.
 *
 * Authors: Zexi Chen (zchen22 ncsu edu)
 *          Vaclav Petras (wenzeslaus gmail com)
 *          Anna Petrasova (kratochanna gmail com)
 *          Chris Jones (cjones1688 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_SIMULATION_HPP
#define POPS_SIMULATION_HPP
 
#include <cmath>
#include <tuple>
#include <vector>
#include <random>
#include <string>
#include <stdexcept>
 
#include "utils.hpp"
 
namespace pops {
 
template<typename Container>
void rotate_left_by_one(Container& container)
{
    std::rotate(container.begin(), container.begin() + 1, container.end());
}
 
template<typename Generator>
std::vector<int> draw_n_from_v(std::vector<int> v, unsigned n, Generator& generator)
{
    if (n > v.size())
        n = v.size();
 
    std::shuffle(v.begin(), v.end(), generator);
    v.erase(v.begin() + n, v.end());
    return v;
}
 
enum class ModelType
{
    SusceptibleInfected,  
    SusceptibleExposedInfected  
};
 
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");
}
 
inline ModelType model_type_from_string(const char* text)
{
    // call the string version
    return model_type_from_string(text ? std::string(text) : std::string());
}
 
template<typename IntegerRaster, typename FloatRaster, typename RasterIndex = int>
class Simulation
{
private:
    RasterIndex rows_;
    RasterIndex cols_;
    bool dispersers_stochasticity_;
    bool establishment_stochasticity_;
    bool movement_stochasticity_;
    ModelType model_type_;
    unsigned latency_period_;
    std::default_random_engine generator_;
 
public:
    Simulation(
        unsigned random_seed,
        RasterIndex rows,
        RasterIndex cols,
        ModelType model_type = ModelType::SusceptibleInfected,
        unsigned latency_period = 0,
        bool dispersers_stochasticity = true,
        bool establishment_stochasticity = true,
        bool movement_stochasticity = true)
        : rows_(rows),
          cols_(cols),
          dispersers_stochasticity_(dispersers_stochasticity),
          establishment_stochasticity_(establishment_stochasticity),
          movement_stochasticity_(movement_stochasticity),
          model_type_(model_type),
          latency_period_(latency_period)
    {
        generator_.seed(random_seed);
    }
 
    Simulation() = delete;
 
    void remove(
        IntegerRaster& infected,
        IntegerRaster& susceptible,
        const FloatRaster& temperature,
        double lethal_temperature,
        const std::vector<std::vector<int>>& suitable_cells)
    {
        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;
            }
        }
    }
 
    void mortality(
        IntegerRaster& infected,
        IntegerRaster& total_hosts,
        double mortality_rate,
        int mortality_time_lag,
        IntegerRaster& died,
        std::vector<IntegerRaster>& mortality_tracker_vector,
        const std::vector<std::vector<int>>& 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);
    }
 
    unsigned movement(
        IntegerRaster& infected,
        IntegerRaster& susceptible,
        std::vector<IntegerRaster>& mortality_tracker_vector,
        std::vector<IntegerRaster>& exposed,
        IntegerRaster& resistant,
        IntegerRaster& total_hosts,
        IntegerRaster& total_exposed,
        unsigned step,
        unsigned last_index,
        const std::vector<std::vector<int>>& movements,
        std::vector<unsigned> movement_schedule,
        std::vector<std::vector<int>>& suitable_cells)
    {
        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<int> exposed_categories;
            std::vector<int> 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<int> categories(total_infecteds, 1);
            categories.insert(categories.end(), suscepts, 2);
            categories.insert(categories.end(), expose, 3);
            categories.insert(categories.end(), resist, 4);
 
            std::vector<int> draw =
                draw_n_from_v(categories, total_hosts_moved, generator_);
            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) {
                int index = 0;
                for (auto& raster : exposed) {
                    auto exposed_count = raster(row_from, col_from);
                    exposed_categories.insert(
                        exposed_categories.end(), exposed_count, index);
 
                    index += 1;
                }
                std::vector<int> exposed_draw =
                    draw_n_from_v(exposed_categories, exposed_moved, generator_);
                index = 0;
                for (auto& raster : exposed) {
                    auto exposed_moved_in_cohort =
                        std::count(exposed_draw.begin(), exposed_draw.end(), index);
                    raster(row_from, col_from) -= exposed_moved_in_cohort;
                    raster(row_to, col_to) += exposed_moved_in_cohort;
                    index += 1;
                }
            }
 
            if (infected_moved > 0) {
                std::vector<int> mortality_categories;
                int index = 0;
                for (auto& raster : mortality_tracker_vector) {
                    auto mortality_count = raster(row_from, col_from);
                    mortality_categories.insert(
                        mortality_categories.end(), mortality_count, index);
                    index += 1;
                }
                std::vector<int> mortality_draw =
                    draw_n_from_v(mortality_categories, infected_moved, generator_);
                index = 0;
                for (auto& raster : mortality_tracker_vector) {
                    auto mortality_moved_in_cohort =
                        std::count(mortality_draw.begin(), mortality_draw.end(), index);
                    raster(row_from, col_from) -= mortality_moved_in_cohort;
                    raster(row_to, col_to) += mortality_moved_in_cohort;
                    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<int> 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();
    }
 
    void generate(
        IntegerRaster& dispersers,
        const IntegerRaster& infected,
        bool weather,
        const FloatRaster& weather_coefficient,
        double reproductive_rate,
        const std::vector<std::vector<int>>& suitable_cells)
    {
        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 * weather_coefficient(i, j);
                int dispersers_from_cell = 0;
                if (dispersers_stochasticity_) {
                    std::poisson_distribution<int> distribution(lambda);
                    for (int k = 0; k < infected(i, j); k++) {
                        dispersers_from_cell += distribution(generator_);
                    }
                }
                else {
                    dispersers_from_cell = lambda * infected(i, j);
                }
                dispersers(i, j) = dispersers_from_cell;
            }
            else {
                dispersers(i, j) = 0;
            }
        }
    }
 
    template<typename DispersalKernel>
    void disperse(
        const IntegerRaster& dispersers,
        IntegerRaster& susceptible,
        IntegerRaster& exposed_or_infected,
        IntegerRaster& mortality_tracker,
        const IntegerRaster& total_populations,
        IntegerRaster& total_exposed,
        std::vector<std::tuple<int, int>>& outside_dispersers,
        bool weather,
        const FloatRaster& weather_coefficient,
        DispersalKernel& dispersal_kernel,
        const std::vector<std::vector<int>>& suitable_cells,
        double establishment_probability = 0.5)
    {
        std::uniform_real_distribution<double> distribution_uniform(0.0, 1.0);
        int row;
        int col;
 
        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));
                        continue;
                    }
                    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 *= weather_coefficient(i, j);
                        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()");
                            }
                        }
                    }
                }
            }
        }
    }
 
    template<typename DispersalKernel>
    void move_overpopulated_pests(
        IntegerRaster& susceptible,
        IntegerRaster& infected,
        const IntegerRaster& total_hosts,
        std::vector<std::tuple<int, int>>& outside_dispersers,
        DispersalKernel& dispersal_kernel,
        const std::vector<std::vector<int>>& suitable_cells,
        double overpopulation_percentage,
        double leaving_percentage)
    {
        struct Move
        {
            RasterIndex row;
            RasterIndex col;
            int count;
        };
        std::vector<Move> 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_, 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 accomodate 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 accomodate all of
            // them. The pests just disappear in both cases.
            else {
                leaving = susceptible(row, col);
                susceptible(row, col) -= leaving;
                infected(row, col) += leaving;
            }
        }
    }
 
    void infect_exposed(
        unsigned step,
        std::vector<IntegerRaster>& exposed,
        IntegerRaster& infected,
        IntegerRaster& mortality_tracker,
        IntegerRaster& total_exposed)
    {
        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;
                // 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()");
        }
    }
 
    template<typename DispersalKernel>
    void disperse_and_infect(
        unsigned step,
        const IntegerRaster& dispersers,
        IntegerRaster& susceptible,
        std::vector<IntegerRaster>& exposed,
        IntegerRaster& infected,
        IntegerRaster& mortality_tracker,
        const IntegerRaster& total_populations,
        IntegerRaster& total_exposed,
        std::vector<std::tuple<int, int>>& outside_dispersers,
        bool weather,
        const FloatRaster& weather_coefficient,
        DispersalKernel& dispersal_kernel,
        const std::vector<std::vector<int>>& suitable_cells,
        double establishment_probability = 0.5)
    {
        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,
            susceptible,
            *infected_or_exposed,
            mortality_tracker,
            total_populations,
            total_exposed,
            outside_dispersers,
            weather,
            weather_coefficient,
            dispersal_kernel,
            suitable_cells,
            establishment_probability);
        if (model_type_ == ModelType::SusceptibleExposedInfected) {
            this->infect_exposed(
                step, exposed, infected, mortality_tracker, total_exposed);
        }
    }
};
 
}  // namespace pops
 
#endif  // POPS_SIMULATION_HPP