Vegetation.hh

Go to the documentation of this file.

#ifndef VEGETATION_HH
#define VEGETATION_HH
 
#include <random>
#include <memory>
#include <algorithm>
#include <string>
 
#include <utopia/core/model.hh>
#include <utopia/core/apply.hh>
#include <utopia/core/types.hh>
#include <utopia/core/cell_manager.hh>
 
#include <utopia/data_io/hdfgroup.hh>
#include <utopia/data_io/hdfdataset.hh>
#include <utopia/data_io/cfg_utils.hh>
 
 
namespace Utopia {
namespace Models {
namespace Vegetation {
 
struct CellState {
    double plant_mass;
};
 
using VegetationTypes = Utopia::ModelTypes<>;
 
using CellTraits = Utopia::CellTraits<CellState, Update::sync, true>;
 
 
class Vegetation:
    public Model<Vegetation, VegetationTypes>
{
public:
    using Base = Model<Vegetation, VegetationTypes>;
 
    using Config = typename Base::Config;
 
    using DataSet = typename Base::DataSet;
 
    using CellManager = Utopia::CellManager<CellTraits, Vegetation>;
 
    using RuleFunc = typename CellManager::RuleFunc;
 
private:
    CellManager _cm;
 
    // -- The parameters of the model -- //
    // TODO Consider making these public or implementing setters & getters
 
    std::normal_distribution<double> _rain_dist;
 
    double _growth_rate;
 
    double _seeding_rate;
 
 
    // Datasets -- //
    std::shared_ptr<DataSet> _dset_plant_mass;
 
 
    // -- Rule functions -- //
 
    RuleFunc _growth_seeding = [this](const auto& cell){
        auto state = cell->state();
        auto rain = _rain_dist(*(this->_rng));
        auto mass = state.plant_mass;
 
        if (rain < 1e-16) {
            rain = 0.;
        }
 
        // Distinguish by mass whether to grow or to seed
        // If negative or smaller than machine epsilon at 1.0 for doubles,
        // consider invalid and reseed
        if (not (mass < 1e-16)) {
            // Logistic Growth
            /* Use Beverton-Holt to approximate discretized logistic growth.
             * Be careful that when using the Wikipedia version
             *   [https://en.wikipedia.org/wiki/Beverton–Holt_model], the R0
             * parameter therein is >= 1 -> proliferation rate!
             * This means that, as given therein, we have:
             *   n_{t+1} = (r*n_t)/(1 + (n_t*(r-1)/K))
             *  with r >= 1, which has to be turned into
             *   n_{t+1} =((r+1.)*n_t)/(1 +(n_t*r)/K)
             * when r is given as a growth rate proper.
             */
            state.plant_mass = ((_growth_rate + 1.) * mass)/(1. + (mass * (_growth_rate))/rain);
        }
        else {
            // Seeding
            state.plant_mass = _seeding_rate * rain;
        }
 
        return state;
    };
 
    double calc_mean_mass () const {
        return std::accumulate(_cm.cells().begin(),
                               _cm.cells().end(),
                               0.,
                               [](double sum, const auto& cell) {
                                   return sum + cell->state().plant_mass;
                               }) / _cm.cells().size();
    }
 
public:
 
 
    template<class ParentModel>
    Vegetation (
        const std::string& name,
        const ParentModel& parent_model,
        const DataIO::Config& custom_cfg = {}
    )
    :
        // Construct the base class
        Base(name, parent_model, custom_cfg),
 
        // Initialize the manager, setting the initial state
        _cm(*this, CellState({0.0})),
 
        // Initialize the rain distribution
        _rain_dist{get_as<double>("rain_mean", this->_cfg),
                   get_as<double>("rain_std", this->_cfg)},
 
        // Initialize model parameters from config file
        _growth_rate(get_as<double>("growth_rate", this->_cfg)),
        _seeding_rate(get_as<double>("seeding_rate", this->_cfg)),
 
        // Open dataset for output of cell states
        _dset_plant_mass(this->create_cm_dset("plant_mass", _cm))
    {
        this->_log->info("'{}' model fully set up.", name);
    }
 
    void perform_step ()
    {
        apply_rule(_growth_seeding, _cm.cells());
    }
 
    void write_data ()
    {
        _dset_plant_mass->write(_cm.cells().begin(),
                                _cm.cells().end(),
                                [](auto& cell) {
                                    return cell->state().plant_mass; }
                                );
    }
 
    void monitor () {
        this->_monitor.set_entry("mean_mass", calc_mean_mass());
    }
 
};
 
 
} // namespace Vegetation
} // namespace Models
} // namespace Utopia
 
#endif // VEGETATION_HH