damping_sampling.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 MIO_EPI_DAMPING_SAMPLING_H
#define MIO_EPI_DAMPING_SAMPLING_H
 
#include "memilio/epidemiology/damping.h"
#include "memilio/utils/uncertain_value.h"
 
namespace mio
{
 
template <typename FP>
class DampingSampling
{
public:
    template <class V>
    DampingSampling(const UncertainValue<FP>& value, DampingLevel level, DampingType type, SimulationTime<FP> time,
                    const std::vector<size_t> matrices, const Eigen::MatrixBase<V>& groups)
        : m_value(value)
        , m_level(level)
        , m_type(type)
        , m_time(time)
        , m_matrices(matrices)
        , m_groups(groups)
    {
    }
 
    const UncertainValue<FP>& get_value() const
    {
        return m_value;
    }
    UncertainValue<FP>& get_value()
    {
        return m_value;
    }
    void set_value(const UncertainValue<FP>& v)
    {
        m_value = v;
    }
 
    DampingLevel get_level() const
    {
        return m_level;
    }
 
    void set_level(DampingLevel l)
    {
        m_level = l;
    }
 
    DampingType get_type() const
    {
        return m_type;
    }
 
    void set_type(DampingType t)
    {
        m_type = t;
    }
 
    SimulationTime<FP> get_time() const
    {
        return m_time;
    }
 
    void set_time(SimulationTime<FP> t)
    {
        m_time = t;
    }
 
    const std::vector<size_t>& get_matrix_indices() const
    {
        return m_matrices;
    }
 
    void set_matrix_indices(const std::vector<size_t>& v)
    {
        m_matrices = v;
    }
 
    const Eigen::VectorX<FP>& get_group_weights() const
    {
        return m_groups;
    }
    template <class V>
    void set_group_weights(const Eigen::MatrixBase<V>& v)
    {
        m_groups = v;
    }
 
    void draw_sample()
    {
        m_value.draw_sample();
    }
 
    bool operator==(const DampingSampling& other) const
    {
        return m_value == other.m_value && m_level == other.m_level && m_type == other.m_type &&
               m_time == other.m_time && m_matrices == other.m_matrices && m_groups == other.m_groups;
    }
    bool operator!=(const DampingSampling& other) const
    {
        return !(*this == other);
    }
    template <class IOContext>
    void serialize(IOContext& io) const
    {
        auto obj = io.create_object("DampingSampling");
        obj.add_element("Time", get_time());
        obj.add_element("Type", get_type());
        obj.add_element("Level", get_level());
        obj.add_element("Value", get_value());
        obj.add_list("MatrixIndices", get_matrix_indices().begin(), get_matrix_indices().end());
        obj.add_element("GroupWeights", get_group_weights());
    }
 
    template <class IOContext>
    static IOResult<DampingSampling> deserialize(IOContext& io)
    {
        auto obj = io.expect_object("DampingSampling");
        auto ti  = obj.expect_element("Time", Tag<SimulationTime<FP>>{});
        auto ty  = obj.expect_element("Type", Tag<DampingType>{});
        auto l   = obj.expect_element("Level", Tag<DampingLevel>{});
        auto v   = obj.expect_element("Value", Tag<UncertainValue<FP>>{});
        auto m   = obj.expect_list("MatrixIndices", Tag<size_t>{});
        auto g   = obj.expect_element("GroupWeights", Tag<Eigen::VectorX<FP>>{});
        return apply(
            io,
            [](auto&& ti_, auto&& ty_, auto&& l_, auto&& v_, auto&& m_, auto&& g_) {
                return DampingSampling(v_, l_, ty_, ti_, m_, g_);
            },
            ti, ty, l, v, m, g);
    }
 
private:
    UncertainValue<FP> m_value;
    DampingLevel m_level;
    DampingType m_type;
    SimulationTime<FP> m_time;
    std::vector<size_t> m_matrices;
    Eigen::VectorX<FP> m_groups;
};
 
template <class DampingExpression, class DampingSamplings, class F>
void apply_dampings(DampingExpression& damping_expression, const DampingSamplings& dampings, F make_matrix)
{
    damping_expression.set_automatic_cache_update(false);
    for (auto& d : dampings) {
        for (auto& i : d.get_matrix_indices()) {
            auto m = make_matrix(d.get_value().value() * d.get_group_weights());
            damping_expression[i].add_damping(m, d.get_level(), d.get_type(), d.get_time());
        }
    }
    damping_expression.set_automatic_cache_update(true);
}
 
template <class V>
auto make_contact_damping_matrix(V&& groups)
{
    auto ones_v = std::decay_t<V>::PlainObject::Constant(groups.size(), 1.0);
    auto prod   = (ones_v - groups) * (ones_v.transpose() - groups.transpose());
    auto ones_m = decltype(prod)::PlainObject::Constant(groups.size(), groups.size(), 1.0);
    return ones_m - sqrt(prod.array()).matrix();
}
 
template <typename FP, class V>
auto make_mobility_damping_vector(ColumnVectorShape<FP> shape, V&& groups)
{
    return Eigen::VectorX<FP>::NullaryExpr(shape.size(), [shape, groups = std::forward<V>(groups)](Eigen::Index i) {
        auto num_groups       = groups.size();
        auto num_compartments = size_t(shape.size()) / num_groups;
        return groups[size_t(i) / num_compartments];
    });
}
 
} // namespace mio
 
#endif // MIO_EPI_DAMPING_SAMPLING_H