dynamic_npis.h

Go to the documentation of this file.

/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Martin J. Kuehn, 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_EPI_DYNAMIC_LOCKDOWN_H
#define MIO_EPI_DYNAMIC_LOCKDOWN_H
 
#include "memilio/epidemiology/damping_sampling.h"
#include "memilio/utils/stl_util.h"
 
#include <cassert>
 
namespace mio
{
 
template <typename FP>
class DynamicNPIs
{
public:
    void set_threshold(FP threshold, const std::vector<DampingSampling<FP>>& dampings)
    {
        insert_sorted_replace(m_thresholds, std::make_pair(threshold, dampings), [](auto& t1, auto& t2) {
            return t1.first > t2.first;
        });
    }
 
    auto get_max_exceeded_threshold(FP value)
    {
        // thresholds are sorted by value descending, so upper_bound returns the first threshold that is smaller using binary search
        auto iter_max_exceeded_threshold =
            std::upper_bound(m_thresholds.begin(), m_thresholds.end(), value, [](auto& val, auto& t2) {
                return val > t2.first;
            });
        return iter_max_exceeded_threshold;
    }
 
    auto get_thresholds() const
    {
        return Range(m_thresholds);
    }
    auto get_thresholds()
    {
        return Range(m_thresholds);
    }
    SimulationTime<FP> get_duration() const
    {
        return m_duration;
    }
 
    void set_duration(SimulationTime<FP> v)
    {
        m_duration = v;
    }
 
    SimulationTime<FP> get_implementation_delay() const
    {
        return m_delay;
    }
    void set_implementation_delay(SimulationTime<FP> delay)
    {
        assert(delay >= SimulationTime<FP>(0.0) && "Implementation delay must be non-negative.");
        m_delay = delay;
    }
    FP get_base_value() const
    {
        return m_base;
    }
    void set_base_value(FP v)
    {
        m_base = v;
    }
    SimulationTime<FP> get_directive_begin() const
    {
        return m_directive_begin;
    }
    void set_directive_begin(SimulationTime<FP> begin)
    {
        assert(begin < m_directive_end && "Directive begin must be before directive end.");
        m_directive_begin = begin;
    }
    SimulationTime<FP> get_directive_end() const
    {
        return m_directive_end;
    }
    void set_directive_end(SimulationTime<FP> end)
    {
        assert(m_directive_begin < end && "Directive end must be strictly after directive begin.");
        m_directive_end = end;
    }
    void draw_sample()
    {
        for (auto&& t : m_thresholds) {
            for (auto&& d : t.second) {
                d.draw_sample();
            }
        }
    }
 
    template <class IOContext>
    void serialize(IOContext& io) const
    {
        auto obj = io.create_object("DynamicNPIs");
        obj.add_list("Thresholds", get_thresholds().begin(), get_thresholds().end());
        obj.add_element("Duration", get_duration());
        obj.add_element("Delay", get_implementation_delay());
        obj.add_element("BaseValue", get_base_value());
        obj.add_element("DirectiveBegin", get_directive_begin());
        obj.add_element("DirectiveEnd", get_directive_end());
    }
 
    template <class IOContext>
    static IOResult<DynamicNPIs> deserialize(IOContext& io)
    {
        auto obj = io.expect_object("DynamicNPIs");
        auto t   = obj.expect_list("Thresholds", Tag<std::pair<FP, std::vector<DampingSampling<FP>>>>{});
        auto d   = obj.expect_element("Duration", Tag<SimulationTime<FP>>{});
        auto i   = obj.expect_element("Delay", Tag<SimulationTime<FP>>{});
        auto b   = obj.expect_element("BaseValue", Tag<FP>{});
        auto f   = obj.expect_element("DirectiveBegin", Tag<SimulationTime<FP>>{});
        auto l   = obj.expect_element("DirectiveEnd", Tag<SimulationTime<FP>>{});
        return apply(
            io,
            [](auto&& t_, auto&& d_, auto&& i_, auto&& b_, auto&& f_, auto&& l_) {
                auto npis = DynamicNPIs();
                npis.set_duration(d_);
                npis.set_implementation_delay(i_);
                npis.set_base_value(b_);
                for (auto&& e : t_) {
                    npis.set_threshold(e.first, e.second);
                }
                npis.set_directive_begin(f_);
                npis.set_directive_end(l_);
                return npis;
            },
            t, d, i, b, f, l);
    }
 
private:
    std::vector<std::pair<FP, std::vector<DampingSampling<FP>>>> m_thresholds;
    SimulationTime<FP> m_duration{14.0};
    SimulationTime<FP> m_delay{0.0};
    SimulationTime<FP> m_directive_begin{SimulationTime<FP>(std::numeric_limits<FP>::lowest())};
    SimulationTime<FP> m_directive_end{SimulationTime<FP>(std::numeric_limits<FP>::max())};
    FP m_base{1.0};
};
 
template <typename FP, class DampingExpr>
std::vector<size_t> get_damping_indices(const DampingExpr& damping_expr, DampingLevel lvl, DampingType type,
                                        SimulationTime<FP> begin, SimulationTime<FP> end)
{
    std::vector<size_t> indices;
    for (size_t i = 0; i < damping_expr.get_dampings().size(); ++i) {
        const auto d = damping_expr.get_dampings()[i];
        if (d.get_level() == lvl && d.get_type() == type && d.get_time() > begin && d.get_time() < end) {
            indices.push_back(i);
        }
    }
    return indices;
}
 
template <typename FP, class DampingExpr>
Eigen::Ref<const typename DampingExpr::Matrix> get_active_damping(const DampingExpr& damping_expr, DampingLevel lvl,
                                                                  DampingType type, SimulationTime<FP> t)
{
    auto ub =
        std::find_if(damping_expr.get_dampings().rbegin(), damping_expr.get_dampings().rend(), [lvl, type, t](auto& d) {
            return d.get_level() == lvl && d.get_type() == type && d.get_time() <= t;
        });
    if (ub != damping_expr.get_dampings().rend()) {
        return ub->get_coeffs();
    }
    return DampingExpr::Matrix::Zero(damping_expr.get_shape().rows(), damping_expr.get_shape().cols());
}
 
template <typename FP, class DampingExprGroup, class MakeMatrix>
void implement_dynamic_npis(DampingExprGroup& damping_expr_group, const std::vector<DampingSampling<FP>>& npis,
                            SimulationTime<FP> begin, SimulationTime<FP> end, MakeMatrix&& make_matrix)
{
    for (auto& npi : npis) {
        for (auto& mat_idx : npi.get_matrix_indices()) {
            auto type          = npi.get_type();
            auto level         = npi.get_level();
            auto& damping_expr = damping_expr_group[mat_idx];
 
            auto active     = get_active_damping<FP>(damping_expr, level, type, begin);
            auto active_end = get_active_damping<FP>(damping_expr, level, type, end)
                                  .eval(); // copy because it may be removed or changed
            auto value = make_matrix(npi.get_value().value() * npi.get_group_weights());
 
            auto npi_implemented = false;
 
            // add begin of npi if not already bigger
            if ((active.array() < value.array()).any()) {
                damping_expr.add_damping(max(value, active), level, type, begin);
                npi_implemented = true;
            }
 
            // replace dampings during the new npi
            auto damping_indices = get_damping_indices<FP>(damping_expr, level, type, begin, end);
            for (auto& i : damping_indices) {
                auto& d = damping_expr.get_dampings()[i];
                damping_expr.add_damping(max(d.get_coeffs(), value), level, type, d.get_time());
                npi_implemented = true;
            }
 
            // add end of npi to restore active dampings if any change was made
            if (npi_implemented) {
                damping_expr.add_damping(active_end, level, type, end);
            }
        }
    }
 
    // remove duplicates that accumulated because of dampings that become active during the time span
    // a damping is obsolete if the most recent damping of the same type and level has the same value
    for (auto& damping_expr : damping_expr_group) {
        // go from the back so indices aren't invalidated when dampings are removed
        // use indices to loop instead of reverse iterators because removing invalidates the current iterator
        for (auto i = int(0); i < int(damping_expr.get_dampings().size()) - 1; ++i) {
            auto it = damping_expr.get_dampings().rbegin() + i;
 
            // look for previous damping of the same type/level
            auto it_prev = std::find_if(it + 1, damping_expr.get_dampings().rend(), [&di = *it](auto& dj) {
                return di.get_level() == dj.get_level() && di.get_type() == dj.get_type();
            });
 
            // remove if match is found and has same value
            if (it_prev != damping_expr.get_dampings().rend() && it->get_coeffs() == it_prev->get_coeffs()) {
                damping_expr.remove_damping(damping_expr.get_dampings().size() - 1 - i);
            }
        }
    }
}
 
} // namespace mio
 
#endif // MIO_EPI_DYNAMIC_LOCKDOWN_H