SimpleFlocking.hh

Go to the documentation of this file.

#ifndef UTOPIA_MODELS_SIMPLEFLOCKING_HH
#define UTOPIA_MODELS_SIMPLEFLOCKING_HH
 
#include <random>
#include <functional>
 
#include <utopia/core/model.hh>
#include <utopia/core/types.hh>
#include <utopia/core/agent_manager.hh>
#include <utopia/core/apply.hh>
 
#include "state.hh"
 
 
namespace Utopia::Models::SimpleFlocking {
 
// ++ Type definitions ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
using ModelTypes = Utopia::ModelTypes<>;
 
 
using AgentTraits = Utopia::AgentTraits<AgentState, Update::sync>;
 
 
 
 
// ++ Model definition ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
class SimpleFlocking:
    public Model<SimpleFlocking, ModelTypes>
{
public:
    using Base = Model<SimpleFlocking, ModelTypes>;
 
    using AgentManager = Utopia::AgentManager<AgentTraits, SimpleFlocking>;
 
    using Agent = AgentManager::Agent;
 
    using AgentPtr = std::shared_ptr<Agent>;
 
    using Rule = typename AgentManager::RuleFunc;
 
    using VoidRule = typename AgentManager::VoidRuleFunc;
 
    using SpaceVec = typename AgentManager::SpaceVec;
 
    using DataSet = typename Base::DataSet;
 
 
private:
    // -- Members -------------------------------------------------------------
 
    AgentManager _am;
 
    // .. Global parameters ...................................................
 
    double _speed;
 
    double _interaction_radius;
 
    double _noise_level;
 
 
    // .. Temporary and helper objects ........................................
 
    std::uniform_real_distribution<double> _noise_distr;
 
 
    // .. Output-related ......................................................
    bool _store_agent_data;
 
    // Agent-specific datasets
    std::shared_ptr<DataSet> _dset_agent_x;
    std::shared_ptr<DataSet> _dset_agent_y;
    std::shared_ptr<DataSet> _dset_agent_orientation;
 
    // Global observables
    std::shared_ptr<DataSet> _dset_orientation_circmean;
    std::shared_ptr<DataSet> _dset_orientation_circstd;
    std::shared_ptr<DataSet> _dset_norm_group_velocity;
 
 
public:
    // -- Model Setup ---------------------------------------------------------
 
 
    template<class ParentModel>
    SimpleFlocking (
        const std::string& name,
        ParentModel& parent_model,
        const DataIO::Config& custom_cfg = {}
    )
    :
        Base(name, parent_model, custom_cfg)
 
    ,   _am(*this)
 
    ,   _speed(get_as<double>("speed", this->_cfg))
    ,   _interaction_radius(get_as<double>("interaction_radius", this->_cfg))
    ,   _noise_level(get_as<double>("noise_level", this->_cfg))
 
    ,   _noise_distr(-_noise_level/2., +_noise_level/2.)
 
    // .. Output-related ......................................................
    ,   _store_agent_data(get_as<bool>("store_agent_data", this->_cfg))
 
    ,   _dset_agent_x(this->create_am_dset("agent/x", _am))
    ,   _dset_agent_y(this->create_am_dset("agent/y", _am))
    ,   _dset_agent_orientation(this->create_am_dset("agent/orientation", _am))
 
    ,   _dset_orientation_circmean(
            this->create_dset("orientation_circmean", {})
        )
    ,   _dset_orientation_circstd(
            this->create_dset("orientation_circstd", {})
        )
    ,   _dset_norm_group_velocity(
            this->create_dset("norm_group_velocity", {})
        )
    {
        set_agent_speed(_speed);
 
        this->_log->info("{} all set up.", this->_name);
        this->_log->info("  Store agent data?  {}", _store_agent_data);
    }
 
 
private:
    // .. Setup functions .....................................................
 
    // .. Helper functions ....................................................
 
public:
    // -- Public Interface ----------------------------------------------------
    // .. Simulation Control ..................................................
 
 
    void perform_step () {
        apply_rule(_adjust_orientation, _am.agents());
        apply_rule(_move, _am.agents());
    }
 
 
 
    void monitor () {
        const auto orientations = get_from_agents([](const auto& agent){
            return agent->state().get_orientation();
        });
 
        const auto [circ_mean, circ_std] = circular_mean_and_std(orientations);
        this->_monitor.set_entry("orientation_mean", circ_mean);
        this->_monitor.set_entry("orientation_std", circ_std);
        this->_monitor.set_entry("norm_group_velocity", norm_group_velocity());
    }
 
 
    void write_data () {
        using WriteT = float;
        const auto& agents = _am.agents();
 
        // -- Global observables
        const auto orientations = get_from_agents([](const auto& agent){
            return agent->state().get_orientation();
        });
 
        const auto [circ_mean, circ_std] = circular_mean_and_std(orientations);
        _dset_orientation_circmean->write(circ_mean);
        _dset_orientation_circstd->write(circ_std);
 
        _dset_norm_group_velocity->write(norm_group_velocity());
 
        // -- Agent-specific data
        // ... only stored optionally
        if (not _store_agent_data) return;
 
        _dset_agent_x->write(
            agents.begin(), agents.end(),
            [](const auto& agent) {
                return static_cast<WriteT>(agent->position()[0]);
        });
 
        _dset_agent_y->write(
            agents.begin(), agents.end(),
            [](const auto& agent) {
                return static_cast<WriteT>(agent->position()[1]);
        });
 
        _dset_agent_orientation->write(
            orientations.begin(), orientations.end(),
            [](const auto& orientation) {
                return static_cast<WriteT>(orientation);
        });
    }
 
 
    // Getters and setters ....................................................
 
    std::size_t num_agents () const {
        return _am.agents().size();
    }
 
    void set_agent_speed (const double new_speed) {
        this->_log->info("Setting all agent's speed to {} ...", new_speed);
        apply_rule(
            [speed=new_speed](const auto& agent){
                auto state = agent->state();
                state.set_speed(speed);
                return state;
            },
            this->_am.agents()
        );
    }
 
 
    double norm_group_velocity () const {
        const auto velocities = get_from_agents([](const auto& agent){
            return agent->state().get_displacement();
        });
        return absolute_group_velocity(velocities) / std::fabs(_speed);
    }
 
    template<
        class Adapter,
        class ValueType = std::invoke_result_t<Adapter, const AgentPtr&>
    >
    std::vector<ValueType> get_from_agents (const Adapter& adapter) const {
        std::vector<ValueType> cont;
        cont.reserve(num_agents());
 
        std::transform(
            this->_am.agents().begin(), this->_am.agents().end(),
            std::back_inserter(cont), adapter
        );
        return cont;
    }
 
 
    // Rules ..................................................................
 
 
    const Rule _adjust_orientation = [this](const auto& agent){
        auto state = agent->state();
 
        // Find all agents within the interaction radius and compute their
        // average orientation, including (!) the current agent.
        // To find an average orientation, we need to separate the x- and y-
        // components and later combine them back into an angle.
        auto agg_sin = std::sin(agent->state().get_orientation());
        auto agg_cos = std::cos(agent->state().get_orientation());
 
        // NOTE The AgentManager::neighbors_of method finds neighbors with
        //      linear complexity in agent number, leading to an overall
        //      quadratic complexity in agent number for this search.
        //      This can be mitigated (on the level of the AgentManager!) by
        //      using a spatially restricted search or a lookup grid.
        for (const auto& nb :
             this->_am.neighbors_of(agent, this->_interaction_radius))
        {
            agg_sin += std::sin(nb->state().get_orientation());
            agg_cos += std::cos(nb->state().get_orientation());
        }
 
        // Can now set the orientation, including noise
        // NOTE Could divide by number of involved agents here, but that is
        //      unnecessary because it cancels out through division anyway.
        state.set_orientation(
            std::atan2(agg_sin, agg_cos) +
            (this->_noise_level > 0. ? _noise_distr(*this->_rng) : 0.)
        );
 
        return state;
    };
 
    const VoidRule _move = [this](const auto& agent){
        this->_am.move_by(agent, agent->state().get_displacement());
    };
 
};
 
 
} // namespace Utopia::Models::SimpleFlocking
 
#endif // UTOPIA_MODELS_SIMPLEFLOCKING_HH