flow_model.h

Go to the documentation of this file.

/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Rene Schmieding, Henrik Zunker
*
* 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_FLOW_MODEL_H
#define MIO_COMPARTMENTS_FLOW_MODEL_H
 
#include "memilio/compartments/compartmental_model.h"
#include "memilio/utils/index_range.h"
#include "memilio/utils/flow.h"
#include "memilio/utils/type_list.h" // IWYU pragma: keep for easier flow definitions
 
namespace mio
{
 
namespace details
{
// The following functions are not defined anywhere. Their use is to provide type conversions via their return type.
 
// Function declaration used to remove OmittedTag from the type list of a tuple.
// First a list of tuples is generated for each Tag in Tags, where the tuple is either of type tuple<Tag>, or if
// Tag == OmittedTag, of type tuple<>. This list is then concatenated, effectively removing OmittedTag.
template <class OmittedTag, class... Tags>
decltype(std::tuple_cat(
    std::declval<std::conditional_t<std::is_same_v<OmittedTag, Tags>, std::tuple<>, std::tuple<Tags>>>()...))
    filter_tuple(std::tuple<Tags...>);
 
// Function declaration used to replace type T by std::tuple.
template <template <class...> class T, class... Args>
std::tuple<Args...> as_tuple(T<Args...>);
 
// Function declaration used to replace std::tuple by type T.
template <template <class...> class T, class... Args>
T<Args...> as_index(std::tuple<Args...>);
 
// Remove all occurrences of OmittedTag from the types in a std::tuple<types...>.
template <class OmittedTag, class Tuple>
using filtered_tuple_t = decltype(filter_tuple<OmittedTag>(std::declval<Tuple>()));
 
// Remove all occurrences of OmittedTag from the types in an Index = IndexTemplate<types...>.
template <class OmittedTag, template <class...> class IndexTemplate, class Index>
using filtered_index_t = decltype(as_index<IndexTemplate>(
    std::declval<filtered_tuple_t<OmittedTag, decltype(as_tuple(std::declval<Index>()))>>()));
 
} //namespace details
 
template <typename FP, class Comp, class Pop, class Params, class Flows>
class FlowModel : public CompartmentalModel<FP, Comp, Pop, Params>
{
    using PopIndex = typename Pop::Index;
    // FlowIndex is the same as PopIndex without the category Comp. It is used as argument type for
    // get_flat_flow_index, since a flow is used to determine only the compartment (i.e. Index<Comp>) for the
    // population of the model.
    // The remaining indices (those contained in FlowIndex) must still be provided to get an entry of population.
    // This approach only works with exactly one category of type Comp in PopIndex (hence the assertion below).
    using FlowIndex = details::filtered_index_t<Comp, Index, PopIndex>;
    // Enforce that Comp is a unique Category of PopIndex, since we use Flows (via their source/target) to provide
    // Comp indices for the population.
    static_assert(FlowIndex::size == PopIndex::size - 1, "Compartments must be used exactly once as population index.");
    static_assert(!Flows::has_duplicates, "Flow duplicates detected. Flows must be unique.");
 
public:
    using Base = CompartmentalModel<FP, Comp, Pop, Params>;
    template <class... Args>
    FlowModel(Args... args)
        : CompartmentalModel<FP, Comp, Pop, Params>(args...)
        , m_flow_values((this->populations.numel() / static_cast<size_t>(Comp::Count)) * Flows::size())
    {
    }
 
    // Note: use get_flat_flow_index when accessing flows
    // Note: by convention, we compute incoming flows, thus entries in flows must be non-negative
    virtual void get_flows(Eigen::Ref<const Eigen::VectorX<FP>> /*pop*/, Eigen::Ref<const Eigen::VectorX<FP>> /*y*/,
                           FP /*t*/, Eigen::Ref<Eigen::VectorX<FP>> /*flows*/) const = 0;
 
    void get_derivatives(Eigen::Ref<const Eigen::VectorX<FP>> flows, Eigen::Ref<Eigen::VectorX<FP>> dydt) const
    {
        // set dydt to 0, then iteratively add all flow contributions
        dydt.setZero();
        if constexpr (std::is_same_v<FlowIndex, Index<>>) {
            // special case where PopIndex only contains Comp, hence FlowIndex has no dimensions to iterate over
            get_rhs_impl(flows, dydt, Index<>{});
        }
        else {
            for (FlowIndex I : reduce_index<FlowIndex>(this->populations.size())) {
                get_rhs_impl(flows, dydt, I);
            }
        }
    }
 
    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 final
    {
        m_flow_values.setZero();
        get_flows(pop, y, t, m_flow_values);
        get_derivatives(m_flow_values, dydt);
    }
 
    Eigen::VectorX<FP> get_initial_flows() const
    {
        return Eigen::VectorX<FP>::Zero((this->populations.numel() / static_cast<size_t>(Comp::Count)) * Flows::size());
    }
 
    template <Comp Source, Comp Target>
    size_t get_flat_flow_index(const FlowIndex& indices) const
    {
        if constexpr (std::is_same_v<FlowIndex, Index<>>) {
            return get_flat_flow_index<Source, Target>();
        }
        else {
            const FlowIndex flow_index_dimensions = reduce_index<FlowIndex>(this->populations.size());
            return flatten_index(indices, flow_index_dimensions) * Flows::size() +
                   index_of_type_v<Flow<Source, Target>, Flows>;
        }
    }
 
    template <Comp Source, Comp Target>
    constexpr size_t get_flat_flow_index() const
    {
        static_assert(std::is_same_v<FlowIndex, Index<>>, "Other indices must be specified");
        return index_of_type_v<Flow<Source, Target>, Flows>;
    }
 
private:
    mutable Eigen::VectorX<FP> m_flow_values; 
 
    template <size_t I = 0>
    inline void get_rhs_impl(Eigen::Ref<const Eigen::VectorX<FP>> flows, Eigen::Ref<Eigen::VectorX<FP>> rhs,
                             const FlowIndex& index) const
    {
        using Flow                 = type_at_index_t<I, Flows>;
        const auto flat_flow_index = get_flat_flow_index<Flow::source, Flow::target>(index);
        const auto flat_source_population =
            this->populations.get_flat_index(extend_index<PopIndex>(index, (size_t)Flow::source));
        const auto flat_target_population =
            this->populations.get_flat_index(extend_index<PopIndex>(index, (size_t)Flow::target));
        rhs[flat_source_population] -= flows[flat_flow_index]; // subtract outflow from source compartment
        rhs[flat_target_population] += flows[flat_flow_index]; // add outflow to target compartment
        // handle next flow (if there is one)
        if constexpr (I + 1 < Flows::size()) {
            get_rhs_impl<I + 1>(flows, rhs, index);
        }
    }
};
 
template <class Model, typename FP>
concept IsFlowModel =
    requires(Model m, Eigen::Ref<const Eigen::VectorX<FP>> const_vref, Eigen::Ref<Eigen::VectorX<FP>> vref, FP t) {
        requires IsCompartmentalModel<Model, FP>;
        { m.get_initial_flows() } -> std::convertible_to<Eigen::VectorX<FP>>;
        m.get_flows(const_vref, const_vref, t, vref);
        m.get_derivatives(const_vref, vref);
    };
 
} // namespace mio
 
#endif // MIO_COMPARTMENTS_FLOW_MODEL_H