compartmental_model.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_COMPARTMENTAL_MODEL_H
#define MIO_COMPARTMENTS_COMPARTMENTAL_MODEL_H
 
#include "memilio/config.h" // IWYU pragma: keep
#include "memilio/math/eigen.h" // IWYU pragma: keep
 
#include <concepts>
 
namespace mio
{
 
template <class T>
concept HasCheckConstraints = requires(T t) {
    { t.check_constraints() } -> std::same_as<bool>;
};
 
template <class T>
concept HasApplyConstraints = requires(T t) {
    { t.apply_constraints() } -> std::same_as<bool>;
};
 
template <typename FP, class Comp, class Pop, class Params>
struct CompartmentalModel {
public:
    using Compartments = Comp;
    using Populations  = Pop;
    using ParameterSet = Params;
    CompartmentalModel(Populations const& po, ParameterSet const& pa)
        : populations{std::move(po)}
        , parameters{pa}
    {
    }
 
    CompartmentalModel(const CompartmentalModel&)            = default;
    CompartmentalModel(CompartmentalModel&&)                 = default;
    CompartmentalModel& operator=(const CompartmentalModel&) = default;
    CompartmentalModel& operator=(CompartmentalModel&&)      = default;
    virtual ~CompartmentalModel()                            = default;
 
    // REMARK: Not pure virtual for easier java/python bindings.
    virtual void get_derivatives(Eigen::Ref<const Eigen::VectorX<FP>> /*pop*/,
                                 Eigen::Ref<const Eigen::VectorX<FP>> /*y*/, FP /*t*/,
                                 Eigen::Ref<Eigen::VectorX<FP>> /*dydt*/) const {};
 
    void eval_right_hand_side(Eigen::Ref<const Eigen::VectorX<FP>> pop, Eigen::Ref<const Eigen::VectorX<FP>> y, FP t,
                              Eigen::Ref<Eigen::VectorX<FP>> dydt) const
    {
        dydt.setZero();
        this->get_derivatives(pop, y, t, dydt);
    }
 
    Eigen::VectorX<FP> get_initial_values() const
    {
        return populations.get_compartments();
    }
 
    bool apply_constraints()
    {
        if constexpr (HasApplyConstraints<ParameterSet>) {
            // use bitwise instead of logical or, so that both apply_constraint functions are called
            return ((int)parameters.apply_constraints() | (int)populations.apply_constraints());
        }
        else {
            return populations.check_constraints();
        }
    }
 
    bool check_constraints() const
    {
        if constexpr (HasCheckConstraints<ParameterSet>) {
            return (parameters.check_constraints() || populations.check_constraints());
        }
        else {
            return populations.check_constraints();
        }
    }
 
    Populations populations{};
    ParameterSet parameters{};
};
 
template <class Model, typename FP>
concept IsCompartmentalModel =
    requires(Model m, Eigen::Ref<const Eigen::VectorX<FP>> pop_y, Eigen::Ref<Eigen::VectorX<FP>> dydt, FP t) {
        { m.get_initial_values() } -> std::convertible_to<Eigen::VectorX<FP>>;
        m.eval_right_hand_side(pop_y, pop_y, t, dydt);
    };
 
} // namespace mio
 
#endif // MIO_COMPARTMENTS_COMPARTMENTAL_MODEL_H