metapopulation_mobility_instant.h

Go to the documentation of this file.

/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: 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 METAPOPULATION_MOBILITY_INSTANT_H
#define METAPOPULATION_MOBILITY_INSTANT_H
 
#include "memilio/compartments/simulation.h"
#include "memilio/epidemiology/contact_matrix.h"
#include "memilio/epidemiology/dynamic_npis.h"
#include "memilio/math/eigen.h"
#include "memilio/math/euler.h"
#include "memilio/math/eigen_util.h"
#include "memilio/mobility/graph_simulation.h"
#include "memilio/utils/compiler_diagnostics.h"
#include "memilio/utils/time_series.h"
 
namespace mio
{
 
template <typename FP, class Sim>
class SimulationNode
{
public:
    using Simulation = Sim;
 
    template <class... Args>
        requires std::is_constructible_v<Sim, Args...>
    SimulationNode(Args&&... args)
        : m_simulation(std::forward<Args>(args)...)
        , m_last_state(m_simulation.get_result().get_last_value())
        , m_t0(m_simulation.get_result().get_last_time())
    {
    }
 
    decltype(auto) get_result() const
    {
        return m_simulation.get_result();
    }
    decltype(auto) get_result()
    {
        return m_simulation.get_result();
    }
    Sim& get_simulation()
    {
        return m_simulation;
    }
    const Sim& get_simulation() const
    {
        return m_simulation;
    }
    Eigen::Ref<const Eigen::VectorX<FP>> get_last_state() const
    {
        return m_last_state;
    }
 
    FP get_t0() const
    {
        return m_t0;
    }
 
    void advance(FP t, FP dt)
    {
        m_simulation.advance(t + dt);
        m_last_state = m_simulation.get_result().get_last_value();
    }
 
private:
    Sim m_simulation;
    Eigen::VectorX<FP> m_last_state;
    FP m_t0;
};
 
template <typename FP>
using MobilityCoefficients = DampingMatrixExpression<FP, VectorDampings<FP>>;
 
template <typename FP>
using MobilityCoefficientGroup = DampingMatrixExpressionGroup<FP, MobilityCoefficients<FP>>;
 
template <typename FP>
class MobilityParameters
{
public:
    MobilityParameters(const MobilityCoefficientGroup<FP>& coeffs)
        : m_coefficients(coeffs)
        , m_saved_compartment_indices(0)
    {
    }
 
    MobilityParameters(const Eigen::VectorX<FP>& coeffs)
        : m_coefficients({MobilityCoefficients<FP>(coeffs)})
        , m_saved_compartment_indices(0)
    {
    }
 
    MobilityParameters(const MobilityCoefficientGroup<FP>& coeffs, const std::vector<std::vector<size_t>>& save_indices)
        : m_coefficients(coeffs)
        , m_saved_compartment_indices(save_indices)
    {
    }
 
    MobilityParameters(const Eigen::VectorX<FP>& coeffs, const std::vector<std::vector<size_t>>& save_indices)
        : m_coefficients({MobilityCoefficients<FP>(coeffs)})
        , m_saved_compartment_indices(save_indices)
    {
    }
 
    bool operator==(const MobilityParameters& other) const
    {
        return m_coefficients == other.m_coefficients;
    }
    bool operator!=(const MobilityParameters& other) const
    {
        return m_coefficients != other.m_coefficients;
    }
 
    const MobilityCoefficientGroup<FP>& get_coefficients() const
    {
        return m_coefficients;
    }
    MobilityCoefficientGroup<FP>& get_coefficients()
    {
        return m_coefficients;
    }
    void set_coefficients(const MobilityCoefficientGroup<FP>& coeffs)
    {
        m_coefficients = coeffs;
    }
 
    const auto& get_save_indices() const
    {
        return m_saved_compartment_indices;
    }
 
    const DynamicNPIs<FP>& get_dynamic_npis_infected() const
    {
        return m_dynamic_npis;
    }
    DynamicNPIs<FP>& get_dynamic_npis_infected()
    {
        return m_dynamic_npis;
    }
    void set_dynamic_npis_infected(const DynamicNPIs<FP>& v)
    {
        m_dynamic_npis = v;
    }
    template <class IOContext>
    void serialize(IOContext& io) const
    {
        auto obj = io.create_object("MobilityParameters");
        obj.add_element("Coefficients", m_coefficients);
        obj.add_element("DynamicNPIs", m_dynamic_npis);
    }
 
    template <class IOContext>
    static IOResult<MobilityParameters> deserialize(IOContext& io)
    {
        auto obj = io.expect_object("MobilityParameters");
        auto c   = obj.expect_element("Coefficients", Tag<MobilityCoefficientGroup<FP>>{});
        auto d   = obj.expect_element("DynamicNPIs", Tag<DynamicNPIs<FP>>{});
        return apply(
            io,
            [](auto&& c_, auto&& d_) {
                MobilityParameters params(c_);
                params.set_dynamic_npis_infected(d_);
                return params;
            },
            c, d);
    }
 
private:
    MobilityCoefficientGroup<FP> m_coefficients; //one per group and compartment
    DynamicNPIs<FP> m_dynamic_npis;
    std::vector<std::vector<size_t>> m_saved_compartment_indices; // groups of indices from compartments to save
};
 
template <typename FP>
class MobilityEdge
{
public:
    MobilityEdge(const MobilityParameters<FP>& params)
        : m_parameters(params)
        , m_mobile_population(params.get_coefficients().get_shape().rows())
        , m_return_times(0)
        , m_return_mobile_population(false)
        , m_saved_compartment_indices(params.get_save_indices())
        , m_mobility_results(m_saved_compartment_indices.size() + 1)
    {
    }
 
    MobilityEdge(const Eigen::VectorX<FP>& coeffs)
        : m_parameters(coeffs)
        , m_mobile_population(coeffs.rows())
        , m_return_times(0)
        , m_return_mobile_population(false)
        , m_saved_compartment_indices(0)
        , m_mobility_results(m_saved_compartment_indices.size() + 1)
    {
    }
 
    MobilityEdge(const MobilityParameters<FP>& params, const std::vector<std::vector<size_t>>& save_indices)
        : m_parameters(params)
        , m_mobile_population(params.get_coefficients().get_shape().rows())
        , m_return_times(0)
        , m_return_mobile_population(false)
        , m_saved_compartment_indices(save_indices)
        , m_mobility_results(m_saved_compartment_indices.size() + 1)
    {
    }
 
    MobilityEdge(const Eigen::VectorX<FP>& coeffs, const std::vector<std::vector<size_t>>& save_indices)
        : m_parameters(coeffs)
        , m_mobile_population(coeffs.rows())
        , m_return_times(0)
        , m_return_mobile_population(false)
        , m_saved_compartment_indices(save_indices)
        , m_mobility_results(m_saved_compartment_indices.size() + 1)
    {
    }
 
    const MobilityParameters<FP>& get_parameters() const
    {
        return m_parameters;
    }
 
    TimeSeries<FP>& get_mobility_results()
    {
        return m_mobility_results;
    }
    const TimeSeries<FP>& get_mobility_results() const
    {
        return m_mobility_results;
    }
    template <class Sim>
    void apply_mobility(FP t, FP dt, SimulationNode<FP, Sim>& node_from, SimulationNode<FP, Sim>& node_to);
 
private:
    MobilityParameters<FP> m_parameters;
    TimeSeries<FP> m_mobile_population;
    TimeSeries<FP> m_return_times;
    bool m_return_mobile_population;
    FP m_t_last_dynamic_npi_check                   = -std::numeric_limits<FP>::infinity();
    std::pair<FP, SimulationTime<FP>> m_dynamic_npi = {-std::numeric_limits<FP>::max(), SimulationTime<FP>(0)};
    std::vector<std::vector<size_t>> m_saved_compartment_indices; // groups of indices from compartments to save
    TimeSeries<FP> m_mobility_results; // save results from edges + entry for the total number of commuters
 
    void add_mobility_result_time_point(const FP t);
};
 
template <typename FP>
void MobilityEdge<FP>::add_mobility_result_time_point(const FP t)
{
    const size_t save_indices_size = this->m_saved_compartment_indices.size();
    if (save_indices_size > 0) {
        const auto& last_value = m_mobile_population.get_last_value();
 
        Eigen::VectorX<FP> condensed_values = Eigen::VectorX<FP>::Zero(save_indices_size + 1);
 
        // sum up the values of m_saved_compartment_indices for each group (e.g. Age groups)
        for (size_t i = 0; i < save_indices_size; ++i) {
            FP sum = 0.0;
            for (auto index : this->m_saved_compartment_indices[i]) {
                sum += last_value[index];
            }
            condensed_values[i] = sum;
        }
 
        // the last value is the sum of commuters
        condensed_values[save_indices_size] = m_mobile_population.get_last_value().sum();
 
        // Move the condensed values to the m_mobility_results time series
        m_mobility_results.add_time_point(t, std::move(condensed_values));
    }
}
 
template <typename FP, IsCompartmentalModelSimulation<FP> Sim>
void calculate_mobility_returns(Eigen::Ref<typename TimeSeries<FP>::Vector> mobile_population, const Sim& sim,
                                Eigen::Ref<const typename TimeSeries<FP>::Vector> total, FP t, FP dt)
{
    auto y0 = mobile_population.eval();
    auto y1 = mobile_population;
    EulerIntegratorCore<FP>().step(
        [&](auto&& y, auto&& t_, auto&& dydt) {
            sim.get_model().get_derivatives(total, y, t_, dydt);
        },
        y0, t, dt, y1);
}
 
template <typename FP, class Sim>
FP get_infections_relative(const SimulationNode<FP, Sim>& node, FP t, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
{
    if constexpr (requires { get_infections_relative<FP>(node.get_simulation(), t, y); }) {
        return get_infections_relative<FP>(node.get_simulation(), t, y);
    }
    else {
        mio::unused(node, t, y);
        mio::log_debug("get_infections_relative called without specialization for this simulation type. "
                       "Implement an overload in your model if you intend to use dynamic NPIs.");
        return FP{0};
    }
}
 
template <typename FP, class Sim>
auto get_mobility_factors(const SimulationNode<FP, Sim>& node, FP t, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
{
    if constexpr (requires { get_mobility_factors<FP>(node.get_simulation(), t, y); }) {
        return get_mobility_factors<FP>(node.get_simulation(), t, y);
    }
    else {
        mio::unused(node, t);
        mio::log_debug("get_mobility_factors called without specialization for this simulation type. "
                       "Using default mobility factor (1.0) for all compartments.");
        return Eigen::VectorX<FP>::Ones(y.rows());
    }
}
 
template <typename FP, class Sim>
void test_commuters(SimulationNode<FP, Sim>& node, Eigen::Ref<Eigen::VectorX<FP>> mobile_population, FP time)
{
    if constexpr (requires { test_commuters<FP>(node.get_simulation(), mobile_population, time); }) {
        test_commuters<FP>(node.get_simulation(), mobile_population, time);
    }
    else {
        mio::unused(node, mobile_population, time);
        mio::log_debug("test_commuters called without specialization for this simulation type. "
                       "Implement an overload in your model if you intend to use commuter testing.");
    }
}
 
template <typename FP>
template <class Sim>
void mio::MobilityEdge<FP>::apply_mobility(FP t, FP dt, SimulationNode<FP, Sim>& node_from,
                                           SimulationNode<FP, Sim>& node_to)
{
    //check dynamic npis
    if (m_t_last_dynamic_npi_check == -std::numeric_limits<FP>::infinity()) {
        m_t_last_dynamic_npi_check = node_from.get_t0();
    }
 
    auto& dyn_npis = m_parameters.get_dynamic_npis_infected();
    if (dyn_npis.get_thresholds().size() > 0 &&
        floating_point_greater_equal<FP>(t, m_t_last_dynamic_npi_check + dyn_npis.get_interval().get())) {
        auto inf_rel =
            get_infections_relative<FP>(node_from, t, node_from.get_last_state()) * dyn_npis.get_base_value();
        auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel);
        if (exceeded_threshold != dyn_npis.get_thresholds().end() &&
            (exceeded_threshold->first > m_dynamic_npi.first ||
             t > FP(m_dynamic_npi.second))) { //old NPI was weaker or is expired
            auto t_end    = SimulationTime<FP>(t + dyn_npis.get_duration().get());
            m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end);
            implement_dynamic_npis<FP>(m_parameters.get_coefficients(), exceeded_threshold->second,
                                       SimulationTime<FP>(t), t_end, [this](auto& g) {
                                           return make_mobility_damping_vector<FP>(
                                               m_parameters.get_coefficients().get_shape(), g);
                                       });
        }
        m_t_last_dynamic_npi_check = t;
    }
 
    //returns
    for (Eigen::Index i = m_return_times.get_num_time_points() - 1; i >= 0; --i) {
        if (m_return_times.get_time(i) <= t) {
            auto v0 = find_value_reverse<FP>(node_to.get_result(), m_mobile_population.get_time(i), 1e-10, 1e-10);
            assert(v0 != node_to.get_result().rend() && "unexpected error.");
            calculate_mobility_returns<FP>(m_mobile_population[i], node_to.get_simulation(), *v0,
                                           m_mobile_population.get_time(i), dt);
 
            //the lower-order return calculation may in rare cases produce negative compartments,
            //especially at the beginning of the simulation.
            //fix by subtracting the supernumerous returns from the biggest compartment of the age group.
            Eigen::VectorX<FP> remaining_after_return =
                (node_to.get_result().get_last_value() - m_mobile_population[i]).eval();
            for (Eigen::Index j = 0; j < node_to.get_result().get_last_value().size(); ++j) {
                if (remaining_after_return(j) < 0) {
                    auto num_comparts = (Eigen::Index)Sim::Model::Compartments::Count;
                    auto group        = Eigen::Index(j / num_comparts);
                    auto compart      = j % num_comparts;
                    log(remaining_after_return(j) < -1e-3 ? LogLevel::warn : LogLevel::info,
                        "Underflow during mobility returns at time {}, compartment {}, age group {}: {}", t, compart,
                        group, remaining_after_return(j));
                    Eigen::Index max_index;
                    slice(remaining_after_return, {group * num_comparts, num_comparts}).maxCoeff(&max_index);
                    log_info("Transferring to compartment {}", max_index);
                    max_index += group * num_comparts;
                    m_mobile_population[i](max_index) -= remaining_after_return(j);
                    m_mobile_population[i](j) += remaining_after_return(j);
                }
            }
            node_from.get_result().get_last_value() += m_mobile_population[i];
            node_to.get_result().get_last_value() -= m_mobile_population[i];
            add_mobility_result_time_point(t);
            m_mobile_population.remove_time_point(i);
            m_return_times.remove_time_point(i);
        }
    }
 
    if (!m_return_mobile_population &&
        (m_parameters.get_coefficients().get_matrix_at(SimulationTime<FP>(t)).array() > 0.0).any()) {
        //normal daily mobility
        m_mobile_population.add_time_point(
            t, (node_from.get_last_state().array() *
                m_parameters.get_coefficients().get_matrix_at(SimulationTime<FP>(t)).array() *
                get_mobility_factors<FP>(node_from, t, node_from.get_last_state()).array())
                   .matrix());
        m_return_times.add_time_point(t + dt);
 
        test_commuters<FP>(node_from, m_mobile_population.get_last_value(), t);
 
        node_to.get_result().get_last_value() += m_mobile_population.get_last_value();
        node_from.get_result().get_last_value() -= m_mobile_population.get_last_value();
 
        add_mobility_result_time_point(t);
    }
    m_return_mobile_population = !m_return_mobile_population;
}
 
template <typename FP, class Sim>
void advance_model(FP t, FP dt, SimulationNode<FP, Sim>& node)
{
    node.advance(t, dt);
}
 
template <typename FP, class Sim>
void apply_mobility(FP t, FP dt, MobilityEdge<FP>& mobilityEdge, SimulationNode<FP, Sim>& node_from,
                    SimulationNode<FP, Sim>& node_to)
{
    mobilityEdge.apply_mobility(t, dt, node_from, node_to);
}
 
template <typename FP, class Sim>
GraphSimulation<FP, Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>, FP, FP,
                void (*)(FP, FP, mio::MobilityEdge<FP>&, mio::SimulationNode<FP, Sim>&, mio::SimulationNode<FP, Sim>&),
                void (*)(FP, FP, mio::SimulationNode<FP, Sim>&)>
make_mobility_sim(FP t0, FP dt, const Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>& graph)
{
    return make_graph_sim<FP>(
        t0, dt, graph, static_cast<void (*)(FP, FP, SimulationNode<FP, Sim>&)>(&advance_model<FP, Sim>),
        static_cast<void (*)(FP, FP, MobilityEdge<FP>&, SimulationNode<FP, Sim>&, SimulationNode<FP, Sim>&)>(
            &apply_mobility<FP>));
}
 
template <typename FP, class Sim>
GraphSimulation<FP, Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>, FP, FP,
                void (*)(FP, FP, mio::MobilityEdge<FP>&, mio::SimulationNode<FP, Sim>&, mio::SimulationNode<FP, Sim>&),
                void (*)(FP, FP, mio::SimulationNode<FP, Sim>&)>
make_mobility_sim(FP t0, FP dt, Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>&& graph)
{
    return make_graph_sim<FP>(
        t0, dt, std::move(graph), static_cast<void (*)(FP, FP, SimulationNode<FP, Sim>&)>(&advance_model<FP, Sim>),
        static_cast<void (*)(FP, FP, MobilityEdge<FP>&, SimulationNode<FP, Sim>&, SimulationNode<FP, Sim>&)>(
            &apply_mobility<FP>));
}
 
template <typename FP, class Sim>
auto make_no_mobility_sim(FP t0, Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>& graph)
{
    using GraphSim = GraphSimulation<FP, Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>, FP, FP,
                                     void (*)(FP, FP, mio::MobilityEdge<FP>&, mio::SimulationNode<FP, Sim>&,
                                              mio::SimulationNode<FP, Sim>&),
                                     void (*)(FP, FP, mio::SimulationNode<FP, Sim>&)>;
    return GraphSim(t0, std::numeric_limits<FP>::infinity(), graph, &advance_model<FP, Sim>,
                    [](FP, FP, MobilityEdge<FP>&, SimulationNode<FP, Sim>&, SimulationNode<FP, Sim>&) {});
}
 
template <typename FP, class Sim>
auto make_no_mobility_sim(FP t0, Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>&& graph)
{
    using GraphSim = GraphSimulation<FP, Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>, FP, FP,
                                     void (*)(FP, FP, mio::MobilityEdge<FP>&, mio::SimulationNode<FP, Sim>&,
                                              mio::SimulationNode<FP, Sim>&),
                                     void (*)(FP, FP, mio::SimulationNode<FP, Sim>&)>;
    return GraphSim(t0, std::numeric_limits<FP>::infinity(), std::move(graph), &advance_model<FP, Sim>,
                    [](FP, FP, MobilityEdge<FP>&, SimulationNode<FP, Sim>&, SimulationNode<FP, Sim>&) {});
}
 
} // namespace mio
 
#endif //METAPOPULATION_MOBILITY_INSTANT_H