memilio/compartments/simulation.h

Go to the documentation of this file.

/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Jan Kleinert, Daniel Abele
*
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*     http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef MIO_COMPARTMENTS_SIMULATION_H
#define MIO_COMPARTMENTS_SIMULATION_H
 
#include "memilio/compartments/compartmental_model.h"
#include "memilio/compartments/simulation_base.h"
#include "memilio/math/integrator.h"
#include "memilio/math/stepper_wrapper.h"
#include "memilio/utils/time_series.h"
 
namespace mio
{
 
template <typename FP>
using DefaultIntegratorCore = mio::ControlledStepperWrapper<FP, boost::numeric::odeint::runge_kutta_cash_karp54>;
 
template <typename FP, IsCompartmentalModel<FP> M>
class Simulation : public details::SimulationBase<FP, M, OdeIntegrator<FP>>
{
public:
    using Base  = details::SimulationBase<FP, M, OdeIntegrator<FP>>;
    using Model = M;
 
    Simulation(Model const& model, FP t0 = 0., FP dt = 0.1)
        : Base(model, std::make_unique<DefaultIntegratorCore<FP>>(), t0, dt)
    {
    }
 
    Eigen::Ref<Eigen::VectorX<FP>> advance(FP tmax)
    {
        return Base::advance(
            [this](auto&& y, auto&& t, auto&& dydt) {
                Base::get_model().eval_right_hand_side(y, y, t, dydt);
            },
            tmax, Base::get_result());
    }
};
 
template <class Simulation, typename FP>
concept IsCompartmentalModelSimulation = requires(Simulation simulation, FP t) {
    requires IsCompartmentalModel<typename Simulation::Model, FP>;
    simulation.advance(t);
};
 
template <typename FP, class Model, class Sim = mio::Simulation<FP, Model>>
TimeSeries<FP> simulate(FP t0, FP tmax, FP dt, Model const& model,
                        std::unique_ptr<OdeIntegratorCore<FP>>&& integrator_core = nullptr)
{
    model.check_constraints();
    Sim sim(model, t0, dt);
    if (integrator_core) {
        sim.set_integrator_core(std::move(integrator_core));
    }
    sim.advance(tmax);
    return sim.get_result();
}
 
} // namespace mio
 
#endif // MIO_COMPARTMENTS_SIMULATION_H