PredatorPreyPlant.hh

Go to the documentation of this file.

#ifndef UTOPIA_MODELS_PREDATORPREYPLANT_HH
#define UTOPIA_MODELS_PREDATORPREYPLANT_HH
 
#include <algorithm>
#include <random>
 
#include <utopia/core/apply.hh>
#include <utopia/core/model.hh>
#include <utopia/core/cell_manager.hh>
 
#include <utopia/models/PredatorPrey/PredatorPrey.hh>
 
#include "species.hh"
 
 
namespace Utopia {
namespace Models {
namespace PredatorPreyPlant {
 
// ++ Type definitions ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
struct State {
    PredatorPrey::SpeciesState predator;
 
    PredatorPrey::SpeciesState prey;
 
    PlantState plant;
 
    template<class RNGType>
    State(const DataIO::Config& cfg, const std::shared_ptr<RNGType>& rng)
    :
        predator{},
        prey{},
        plant{}
    {
        std::uniform_real_distribution<double> dist(0., 1.);
 
        // Get the threshold probability value
        const auto p_plant = get_as<double>("p_plant", cfg);
        const auto p_prey = get_as<double>("p_prey", cfg);
        const auto p_predator = get_as<double>("p_predator", cfg);
 
        auto cfg_prey = get_as<DataIO::Config>("prey", cfg);
        auto cfg_predator = get_as<DataIO::Config>("predator", cfg);
 
        // Check if the max resource limit is actually higher than the lower
        // limit.
        const auto min_init_resources_prey = get_as<int>("min_init_resources", 
            cfg_prey);
        const auto max_init_resources_prey = get_as<int>("max_init_resources", 
            cfg_prey);
        const auto min_init_resources_predator = 
            get_as<int>("min_init_resources", cfg_predator);
        const auto max_init_resources_predator = 
            get_as<int>("max_init_resources", cfg_predator);
        if(max_init_resources_predator < 
            min_init_resources_predator){
            throw::std::invalid_argument("The upper limit for the initial "
            "predator resources needs to be higher than the lower limit.");
        }
        if(max_init_resources_prey < 
            min_init_resources_prey){
            throw::std::invalid_argument("The upper limit for the initial "
            "prey resources needs to be higher than the lower limit.");
        }
        // Set a plant on this cell with a given probability
        plant.on_cell = (dist(*rng) < p_plant);
 
        // Set a predator on a cell with the given probability and the 
        // resources in the given range.
        if (dist(*rng) < p_predator) {
            predator.on_cell = true;
            std::uniform_int_distribution<> init_res_dist_pred(
                min_init_resources_predator,
                max_init_resources_predator
            );
            predator.resources = init_res_dist_pred(*rng);
        }
 
        // Set a prey on this cell with the desired probability; also set the
        // initial resource amount from an integer distribution.
        if (dist(*rng) < p_prey) {
            prey.on_cell = true;
            std::uniform_int_distribution<> init_res_dist_prey(
                min_init_resources_prey,
                max_init_resources_prey
            );
 
            prey.resources = init_res_dist_prey(*rng);
        }
 
 
    }
};
 
 
 
using CellTraits = Utopia::CellTraits<State, Update::manual>;
 
using ModelTypes = ModelTypes<>;
 
// ++ Model definition ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
class PredatorPreyPlant
    : public Model<PredatorPreyPlant, ModelTypes>
{
public:
    using Base = Model<PredatorPreyPlant, ModelTypes>;
 
    using DataSet = typename Base::DataSet;
 
    using CellManager = Utopia::CellManager<CellTraits, PredatorPreyPlant>;
 
    using Cell = typename CellManager::Cell;
 
    using Rule = typename CellManager::RuleFunc;
 
 
private:
    // Base members: _time, _name, _cfg, _hdfgrp, _rng, _monitor, _log, _space
    // ... but you should definitely check out the documentation ;)
 
    // -- Members -------------------------------------------------------------
    CellManager _cm;
 
    // .. Model parameters ....................................................
    SpeciesParams _params;
 
    std::size_t _num_moves;
 
    // .. Temporary objects ...................................................
    CellContainer<Cell> _prey_cell;
    
    CellContainer<Cell> _empty_cell;
 
    CellContainer<Cell> _repro_cell;
 
    CellContainer<Cell> _resource_cell;
 
 
    std::uniform_real_distribution<double> _prob_distr;
 
    std::uniform_int_distribution<> _cm_dist;
 
 
    // .. Datasets ............................................................
    const std::shared_ptr<DataSet> _dset_prey;
    
    const std::shared_ptr<DataSet> _dset_predator;
    
    const std::shared_ptr<DataSet> _dset_resource_prey;
    
    const std::shared_ptr<DataSet> _dset_resource_predator;
 
    const std::shared_ptr<DataSet> _dset_plant;
 
 
    // .. Rule functions and helper methods ...................................
 
    Rule _cost_of_living = [this](const auto& cell) {
        // Get the state of the cell
        auto& state = cell->state;
 
        // Subtract the cost of living and clamp the resources to the limits:
        // If the resources exceed the maximal resources they are equal to
        // the maximal resources and if they go below 0 they are mapped to 0.
        state.predator.resources =
            std::clamp(  state.predator.resources
                       - _params.predator.cost_of_living,
                       0., _params.predator.resource_max);
        state.prey.resources =
            std::clamp(  state.prey.resources 
                       - _params.prey.cost_of_living, 
                       0., _params.prey.resource_max);
 
        // Remove predators that have no resources.
        if (state.predator.on_cell and state.predator.resources <= 0.)
            state.predator.on_cell = false;
 
        // Remove prey that have no resources.
        if (state.prey.on_cell and state.prey.resources <= 0.) 
            state.prey.on_cell = false;
 
        return state;
    };
 
 
    std::shared_ptr<Cell>
    get_random_neighbor(const CellContainer<Cell>& nbs) const
    {
        return nbs[std::uniform_int_distribution<std::size_t>(0, nbs.size() - 1)(*this->_rng)];
    }
 
 
    void move_predator_to_nb_cell(const std::shared_ptr<Cell>& cell, 
                                  const std::shared_ptr<Cell>& nb_cell) {
        auto& state = cell->state;
        auto& nb_state = nb_cell->state;
 
        nb_state.predator.on_cell = true;
        nb_state.predator.resources = state.predator.resources;
        
        state.predator.on_cell = false;
        state.predator.resources = 0.;
    }
 
 
    void move_prey_to_nb_cell(const std::shared_ptr<Cell>& cell, 
                              const std::shared_ptr<Cell>& nb_cell) {
        auto& state = cell->state;
        auto& nb_state = nb_cell->state;
 
        nb_state.prey.on_cell = true;
        nb_state.prey.resources = state.prey.resources;
 
        state.prey.on_cell = false;
        state.prey.resources = 0.;
    }
 
 
 
    auto move_prey(std::shared_ptr<Cell> cell) {
        _resource_cell.clear();
 
        // Collect empty neighboring cells with available resources, where the
        // prey could move
        for (const auto& nb : this->_cm.neighbors_of(cell)) {
            if (nb->state.plant.on_cell 
                and not nb->state.prey.on_cell 
                and not nb->state.predator.on_cell)
                _resource_cell.push_back(nb);
        }
 
        // Choose a random cell with resources.
        if (_resource_cell.size() > 0) {
            std::uniform_int_distribution<> dist(0, _resource_cell.size() - 1);
            const auto nb_cell = _resource_cell[dist(*this->_rng)];
            move_prey_to_nb_cell(cell, nb_cell);
            return nb_cell;
        }
        // if there isn't one, then try a random walk
        else {
            const std::shared_ptr<Cell> nb_cell = get_random_neighbor(_cm.neighbors_of(cell));
            // The will move even if there is a predator.
            if (not nb_cell->state.prey.on_cell) {
                move_prey_to_nb_cell(cell, nb_cell);
                return nb_cell;
            }
            return cell;
        }
    };
 
 
 
    auto move_predator(std::shared_ptr<Cell> cell) {
        // clear the container for cells that contain prey or empty cells
        // in the neighbourhood
        _prey_cell.clear();
 
        // Collect neighboring cells with preys and without predator
        for (const auto& nb : this->_cm.neighbors_of(cell)) {
            if (    nb->state.prey.on_cell
                and not nb->state.predator.on_cell)
                _prey_cell.push_back(nb);
        }
 
        // now update the cell state and the respective neighbor
        if (_prey_cell.size() > 0) {
            // distribution to choose a random cell for the movement
            std::uniform_int_distribution<> dist_prey(0 ,_prey_cell.size() - 1);
            auto nb_cell = _prey_cell[dist_prey(*this->_rng)];
            // move the predator to the given cell
            move_predator_to_nb_cell(cell, nb_cell);
 
            return nb_cell;
        }
        else {
            // If no prey in the neighborhoods try a random step if possible
            const std::shared_ptr<Cell> nb_cell = get_random_neighbor(_cm.neighbors_of(cell));
            if (not nb_cell->state.predator.on_cell){
                // move the predator to the given cell
                move_predator_to_nb_cell(cell, nb_cell);
                return nb_cell;
            }
        }
        return cell;
    };
 
 
    /*+ Go through all cells. If there are both predator and prey on the cell,
    *   then the prey attempts to flee with a given probability.
    *   If there is only a predator on the cell then it moves until it reaches 
    *   the movement limit or reaches a prey. 
    *   If there is only a prey on the cell and there are no plant resources, 
    *   then the prey moves until it reaches the movement limit or it finds a 
    *   cell with plant resource on it.
    */
    void move_entities(std::shared_ptr<Cell> cell) {
        auto& state = cell->state;
 
        // Marker used to compute how many steps the prey or predator has taken
        unsigned int step = 0;
 
        if (state.predator.on_cell) {
            while (    step++ < _params.predator.move_limit
                   and not cell->state.prey.on_cell)
            {
                cell = move_predator(cell);
            }
        }
        else if (state.prey.on_cell and not cell->state.plant.on_cell){
            while(step++ < _params.prey.move_limit and not cell->state.plant.on_cell)
                cell = move_prey(cell);
        }
    };
 
 
 
 
    Rule _flee_prey = [this](const auto& cell) {
        auto& state = cell->state;
        
        if (state.prey.on_cell and state.predator.on_cell and 
            (_prob_distr(*this->_rng) < _params.prey.p_flee)) {
            // Collect empty neighboring cells to which the prey could flee
            _empty_cell.clear();
            for (const auto& nb : this->_cm.neighbors_of(cell)) {
                if (    not nb->state.prey.on_cell
                    and not nb->state.predator.on_cell)
                    _empty_cell.push_back(nb);
            }
 
            // If there is an empty cell, move there
            if (_empty_cell.size() > 0) {
                // Choose a random cell to move to
                std::uniform_int_distribution<> 
                    dist(0, _empty_cell.size() - 1);
                auto nb_cell = _empty_cell[dist(*this->_rng)];
                move_prey_to_nb_cell(cell, nb_cell);
            }
        }
        return state;
    };
 
 
 
 
    Rule _eat = [this](const auto& cell) {
        auto& state = cell->state;
 
        // Predator eats prey
        if (state.predator.on_cell and state.prey.on_cell) {
            // Increment resources and clamp to [0, resource_max]
            state.predator.resources =
                std::clamp(  state.predator.resources 
                           + _params.predator.resource_intake, 
                           0., _params.predator.resource_max);
 
            // Remove the prey from the cell
            state.prey.on_cell = false;
            state.prey.resources = 0.;
        }
 
        // Prey eats plants
        else if (state.prey.on_cell and state.plant.on_cell) {
            // Increment resources and clamp to [0, resource_max]
            state.prey.resources =
                std::clamp(  state.prey.resources 
                           + _params.prey.resource_intake, 
                           0., _params.prey.resource_max);
            
            if (_params.plant.growth_model != GrowthModel::none) {
                state.plant.on_cell = false;
                state.plant.regeneration_counter = 0;
            }
        }
        return state;
    };
 
 
 
    Rule _reproduce = [this](const auto& cell) {
        auto& state = cell->state;
        
        // Reproduction of predators    
        if (    state.predator.on_cell
            and (this->_prob_distr(*this->_rng) < _params.predator.repro_prob)
            and (   state.predator.resources
                 >= _params.predator.repro_resource_requ))
        {
            const std::shared_ptr<Cell> nb_cell = get_random_neighbor(_cm.neighbors_of(cell));
            if (not nb_cell->state.predator.on_cell) {
                nb_cell->state.predator.on_cell = true;
 
                // transfer resources from parent to offspring
                nb_cell->state.predator.resources = _params.predator.repro_cost;
                state.predator.resources -= _params.predator.repro_cost; 
            }
        }
 
        // Reproduction of preys
        if (    state.prey.on_cell
            and this->_prob_distr(*this->_rng) < _params.prey.repro_prob
            and state.prey.resources >= _params.prey.repro_resource_requ)
        {
            const std::shared_ptr<Cell> nb_cell = get_random_neighbor(_cm.neighbors_of(cell));
            if (not nb_cell->state.prey.on_cell) {
                nb_cell->state.prey.on_cell = true;
 
                // transfer resources from parent to offspring
                nb_cell->state.prey.resources = _params.prey.repro_cost;
                state.prey.resources -= _params.prey.repro_cost;
            }
        }
 
        // If there is no plant, there _may_ be plant growth, depending on the
        // plant growth model
        if (not state.plant.on_cell) {
            // TODO Ideally, the model is decided once in the beginning, and
            //      not for each cell anew.
            if (_params.plant.growth_model == GrowthModel::deterministic) {
                // If the regeneration counter if not high enough, increment it
                if (   state.plant.regeneration_counter
                    >= _params.plant.regen_time)
                {
                    // Grow a plant :)
                    state.plant.on_cell = true;
                }
                else {
                    state.plant.regeneration_counter++;
                }
            }
            else if (_params.plant.growth_model == GrowthModel::stochastic) {
                // Regrow with a certain probability
                if (_prob_distr(*this->_rng) < _params.plant.regen_prob) {
                    state.plant.on_cell = true;
                }
            }
            // else: no regrowth
        }
 
        return state;
    };
 
 
public:
    // -- Model setup ---------------------------------------------------------
 
    template <class ParentModel>
    PredatorPreyPlant (
        const std::string& name,
        ParentModel& parent_model,
        const DataIO::Config& custom_cfg = {}
    )
    : 
        Base(name, parent_model, custom_cfg),
        
        // Initialize the cell manager, binding it to this model
        _cm(*this),
 
        // Extract model parameters
        _params(this->_cfg),
        _num_moves([&](){
            const auto f = get_as<double>("num_moves_fraction", this->_cfg);
            return f * this->_cm.cells().size();
        }()),
 
        // Temporary cell containers
        _prey_cell(),
        _empty_cell(),
        _repro_cell(),
 
        // create random distributions
        _prob_distr(0., 1.),
        _cm_dist(0, _cm.cells().size() - 1),
 
        // create datasets
        _dset_prey(this->create_cm_dset("prey", _cm)),
        _dset_predator(this->create_cm_dset("predator", _cm)),
        _dset_resource_prey(this->create_cm_dset("resource_prey", _cm)),
        _dset_resource_predator(this->create_cm_dset("resource_predator", _cm)),
        _dset_plant(this->create_cm_dset("plant", _cm))
    {
        // Inform about number of movements per iteration step:
        this->_log->info("The movement rule will be applied {} times each "
                         "time step.", _num_moves);
 
        // Load the cell state from a file, overwriting the current state
        if (this->_cfg["cell_states_from_file"]) {
            setup_cell_states_from_file(this->_cfg["cell_states_from_file"]);
        }
 
        // Reserve memory in the size of the neighborhood for the temp. vectors
        const auto nb_size = _cm.nb_size();
        _prey_cell.reserve(nb_size);
        _empty_cell.reserve(nb_size);
        _repro_cell.reserve(nb_size);
 
        // Initialization finished
        this->_log->info("{} model fully set up.", this->_name);
    }
 
private:
    // .. Setup functions .....................................................
 
    void setup_cell_states_from_file(const Config& cs_cfg) {
        const auto hdf5_file = get_as<std::string>("hdf5_file", cs_cfg);
 
        if (get_as<bool>("load_predator", cs_cfg)) {
            this->_log->info("Loading predator positions from file ...");
 
            // Use the CellManager to set the cell state from the data
            // given by the `predator` dataset. Load as int to be able to
            // detect that a user supplied invalid values (better than
            // failing silently, which would happen with booleans).
            _cm.set_cell_states<int>(hdf5_file, "predator",
                [this](auto& cell, const int on_cell){
                    if (on_cell == 0 or on_cell == 1) {
                        // Place predator state on cell.
                        cell->state.predator.on_cell = on_cell;
                        // Take care of species resources.
                        if(on_cell){
                            const auto& predator_cfg = 
                                get_as<Config>("predator", 
                                _cfg["cell_manager"]["cell_params"]);
                            int min_init_resources_predator = 
                                get_as<int>("min_init_resources", 
                                predator_cfg);
                            int max_init_resources_predator = 
                                get_as<int>("max_init_resources", 
                                predator_cfg);
                            std::uniform_int_distribution<> 
                                init_res_dist_predator(
                                min_init_resources_predator,
                                max_init_resources_predator
                            );
                            cell->state.predator.resources = 
                                init_res_dist_predator(*this->_rng);
                        }
                        else{
                            cell->state.predator.resources = 0;
                        }
                        return;
                    }
                    throw std::invalid_argument("While setting predator "
                        "positions, encountered an invalid value: "
                        + std::to_string(on_cell) + ". Allowed: 0 or 1.");
                }
            );
            this->_log->info("Predator positions loaded.");
        }
 
        if (get_as<bool>("load_prey", cs_cfg)) {
            this->_log->info("Loading prey positions from file ...");
 
            _cm.set_cell_states<int>(hdf5_file, "prey",
                [this](auto& cell, const int on_cell){
                    if (on_cell == 0 or on_cell == 1) {
                        // Place prey state on cell.
                        cell->state.prey.on_cell = on_cell;
                        // Take care of species resources.
                        if(on_cell){
                            const auto& prey_cfg = 
                                get_as<Config>("prey", 
                                _cfg["cell_manager"]["cell_params"]);
                            int min_init_resources_prey = 
                                get_as<int>("min_init_resources", prey_cfg);
                            int max_init_resources_prey = 
                                get_as<int>("max_init_resources", prey_cfg);
                            std::uniform_int_distribution<> 
                                init_res_dist_prey(
                                min_init_resources_prey,
                                max_init_resources_prey
                            );
                            cell->state.prey.resources = 
                                init_res_dist_prey(*this->_rng);
                        }
                        else{
                            cell->state.prey.resources = 0;
                        }
                        return;
                    }
                    throw std::invalid_argument("While setting prey "
                        "positions, encountered an invalid value: "
                        + std::to_string(on_cell) + ". Allowed: 0 or 1.");
                }
            );
            this->_log->info("Prey positions loaded.");
        }
 
        if (get_as<bool>("load_plant", cs_cfg)) {
            this->_log->info("Loading plant positions from file ...");
 
            _cm.set_cell_states<int>(hdf5_file, "plant",
                [](auto& cell, const int on_cell){
                    if (on_cell == 0 or on_cell == 1) {
                        cell->state.plant.on_cell = on_cell;
                        return;
                    }
                    throw std::invalid_argument("While setting plant "
                        "positions, encountered an invalid value: "
                        + std::to_string(on_cell) + ". Allowed: 0 or 1.");
                }
            );
            this->_log->info("Plant positions loaded.");
        }
    }
 
public:
    // -- Public Interface ----------------------------------------------------
    // .. Simulation Control ..................................................
 
 
    void perform_step() {
        apply_rule<Update::async, Shuffle::off>(_cost_of_living, _cm.cells());
 
        // Choose _num_moves cells randomly and apply the movement rule to them
        // NOTE Be aware that cells can be selected multiple times.
        for (std::size_t i = 0; i < _num_moves; ++i) {
            move_entities(_cm.cells()[_cm_dist(*this->_rng)]);
        }
        apply_rule<Update::async, Shuffle::on>(_flee_prey, _cm.cells(), 
                                                *this->_rng);
        apply_rule<Update::async, Shuffle::off>(_eat, _cm.cells());
        apply_rule<Update::async, Shuffle::on>(_reproduce, _cm.cells(),
                                               *this->_rng);
    }
 
    void monitor () {
        // Calculate the densities for all species
        auto [pred_density, prey_density, plant_density] = [this](){
            double predator_sum = 0.;
            double prey_sum = 0.;
            double plant_sum = 0.;
            double num_cells = this->_cm.cells().size();
 
            for (const auto& cell : this->_cm.cells()) {
                auto state = cell->state;
 
                if (state.prey.on_cell) prey_sum++;
                if (state.predator.on_cell) predator_sum++;
                if (state.plant.on_cell) plant_sum++;
            }
            return std::tuple{predator_sum / num_cells,
                              prey_sum / num_cells,
                              plant_sum /  num_cells};
        }();
 
        this->_monitor.set_entry("predator_density", pred_density);
        this->_monitor.set_entry("prey_density", prey_density);
        this->_monitor.set_entry("plant_density", plant_density);
    }
 
 
    void write_data() {
        // Predator
        _dset_predator->write(_cm.cells().begin(), _cm.cells().end(), 
            [](const auto& cell) {                
                return static_cast<char>(cell->state.predator.on_cell);
            }
        );
 
        // Prey
        _dset_prey->write(_cm.cells().begin(), _cm.cells().end(), 
            [](const auto& cell) {                
                return static_cast<char>(cell->state.prey.on_cell);
            }
        );
 
        // Plant
        _dset_plant->write(_cm.cells().begin(), _cm.cells().end(), 
            [](const auto& cell) {
                return static_cast<char>(cell->state.plant.on_cell);
            }
        );
 
        // resource of predator
        _dset_resource_predator->write(_cm.cells().begin(), _cm.cells().end(), 
            [](const auto& cell) {
                return cell->state.predator.resources;
            }
        );
 
        // resource of prey
        _dset_resource_prey->write(_cm.cells().begin(), _cm.cells().end(), 
            [](const auto& cell) {
                return cell->state.prey.resources;
            }
        );
    }
};
 
} // namespace PredatorPreyPlant
} // namespace Models
} // namespace Utopia
 
#endif // UTOPIA_MODELS_PREDATORPREYPLANT_HH