contact_matrix.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_CONTACT_MATRIX_H
#define MIO_EPI_CONTACT_MATRIX_H
 
#include "memilio/epidemiology/damping.h"
#include "memilio/math/math_utils.h"
#include "memilio/utils/stl_util.h"
 
#include <numeric>
#include <ostream>
#include <vector>
 
namespace mio
{
 
template <typename FP, class D>
class DampingMatrixExpression
{
public:
    using DampingsType = D;
    using Shape        = typename DampingsType::Shape;
    using Matrix       = typename Shape::Matrix;
 
    template <class M, class K>
    DampingMatrixExpression(const Eigen::MatrixBase<M>& baseline, const Eigen::MatrixBase<K>& minimum)
        : m_baseline(baseline)
        , m_minimum(minimum)
        , m_dampings(Shape::get_shape_of(m_baseline))
    {
        assert(Shape::get_shape_of(m_minimum) == Shape::get_shape_of(m_baseline));
    }
 
    template <class M>
    explicit DampingMatrixExpression(const Eigen::MatrixBase<M>& baseline)
        : DampingMatrixExpression(
              baseline, Matrix::Zero(Shape::get_shape_of(baseline).rows(), Shape::get_shape_of(baseline).cols()))
    {
    }
 
    template <class... T>
        requires std::is_constructible_v<Shape, T...>
    explicit DampingMatrixExpression(T... shape_args)
        : DampingMatrixExpression(Matrix::Zero(Shape(shape_args...).rows(), Shape(shape_args...).cols()))
    {
    }
 
    template <class... T>
    void add_damping(T&&... t)
    {
        m_dampings.add(std::forward<T>(t)...);
    }
 
    void remove_damping(size_t i)
    {
        m_dampings.remove(i);
    }
 
    void clear_dampings()
    {
        m_dampings.clear();
    }
 
    auto get_dampings() const
    {
        return Range(m_dampings.begin(), m_dampings.end());
    }
    auto get_dampings()
    {
        return Range(m_dampings.begin(), m_dampings.end());
    }
 
    const Matrix& get_baseline() const
    {
        return m_baseline;
    }
    Matrix& get_baseline()
    {
        return m_baseline;
    }
 
    const Matrix& get_minimum() const
    {
        return m_minimum;
    }
    Matrix& get_minimum()
    {
        return m_minimum;
    }
 
    Shape get_shape() const
    {
        return Shape::get_shape_of(m_baseline);
    }
 
    bool operator==(const DampingMatrixExpression& other) const
    {
        return m_baseline == other.m_baseline && m_minimum == other.m_minimum && m_dampings == other.m_dampings;
    }
    bool operator!=(const DampingMatrixExpression& other) const
    {
        return !(*this == other);
    }
 
    void set_automatic_cache_update(bool b)
    {
        m_dampings.set_automatic_cache_update(b);
    }
 
    auto get_matrix_at(SimulationTime<FP> t) const
    {
        assert(Shape::get_shape_of(m_dampings.get_matrix_at(t)) == Shape::get_shape_of(m_baseline));
        return m_baseline - (m_dampings.get_matrix_at(t).array() * (m_baseline - m_minimum).array()).matrix();
    }
 
    friend void PrintTo(const DampingMatrixExpression& self, std::ostream* os)
    {
        *os << '\n' << self.m_baseline;
        *os << '\n' << self.m_minimum;
        PrintTo(self.m_dampings, os);
    }
 
    template <class IOContext>
    void serialize(IOContext& io) const
    {
        auto obj = io.create_object("DampingMatrixExpression");
        obj.add_element("Baseline", get_baseline());
        obj.add_element("Minimum", get_minimum());
        obj.add_list("Dampings", get_dampings().begin(), get_dampings().end());
    }
 
protected:
    template <class IOContext, class Derived>
    static IOResult<Derived> deserialize(IOContext& io, Tag<Derived>)
    {
        auto obj = io.expect_object("DampingMatrixExpression");
        auto b   = obj.expect_element("Baseline", Tag<Matrix>{});
        auto m   = obj.expect_element("Minimum", Tag<Matrix>{});
        auto d   = obj.expect_list("Dampings", Tag<typename DampingsType::value_type>{});
        return apply(
            io,
            [](auto&& b_, auto&& m_, auto&& d_) -> IOResult<Derived> {
                if (Shape::get_shape_of(b_) != Shape::get_shape_of(m_)) {
                    return failure(StatusCode::InvalidValue, "Baseline and Minimum must have the same shape.");
                }
                auto r = Derived(b_, m_);
                for (auto&& e : d_) {
                    if (e.get_shape() != Shape::get_shape_of(b_)) {
                        return failure(StatusCode::InvalidValue, "Dampings must have the same shape as the Baseline.");
                    }
                    r.add_damping(e);
                }
                return success(r);
            },
            b, m, d);
    }
 
public:
    template <class IOContext>
    static IOResult<DampingMatrixExpression> deserialize(IOContext& io)
    {
        return deserialize(io, Tag<DampingMatrixExpression>{});
    }
 
private:
    Matrix m_baseline;
    Matrix m_minimum;
    DampingsType m_dampings;
};
 
template <typename FP, class E>
class DampingMatrixExpressionGroup
{
public:
    using value_type      = E;
    using reference       = value_type&;
    using const_reference = const value_type&;
    using iterator        = typename std::vector<value_type>::iterator;
    using const_iterator  = typename std::vector<value_type>::const_iterator;
 
    using Shape        = typename value_type::Shape;
    using Matrix       = typename value_type::Matrix;
    using DampingsType = typename value_type::DampingsType;
 
    template <class... T>
        requires(std::is_constructible_v<Shape, T...>)
    explicit DampingMatrixExpressionGroup(size_t num_matrices, T... shape_args)
        : m_matrices(num_matrices, value_type{shape_args...})
    {
        assert(num_matrices > 0);
    }
 
    DampingMatrixExpressionGroup(std::initializer_list<value_type> il)
        : m_matrices(il)
    {
        assert(il.size() > 0);
    }
 
    DampingMatrixExpressionGroup(const std::vector<value_type>& v)
        : m_matrices(v)
    {
        assert(v.size() > 0);
    }
 
    reference operator[](size_t i)
    {
        return m_matrices[i];
    }
    const_reference operator[](size_t i) const
    {
        return m_matrices[i];
    }
 
    size_t get_num_matrices() const
    {
        return m_matrices.size();
    }
 
    Shape get_shape() const
    {
        return m_matrices[0].get_shape();
    }
 
    bool operator==(const DampingMatrixExpressionGroup& other) const
    {
        return m_matrices == other.m_matrices;
    }
    bool operator!=(const DampingMatrixExpressionGroup& other) const
    {
        return !(*this == other);
    }
 
    template <class... T>
    void add_damping(T&&... t)
    {
        for (auto& m : *this) {
            m.add_damping(std::forward<T>(t)...);
        }
    }
 
    void clear_dampings()
    {
        for (auto& m : *this) {
            m.clear_dampings();
        }
    }
 
    void set_automatic_cache_update(bool b)
    {
        for (auto& m : *this) {
            m.set_automatic_cache_update(b);
        }
    }
 
    auto get_matrix_at(SimulationTime<FP> t) const
    {
        return Eigen::Matrix<FP, Eigen::Dynamic, Eigen::Dynamic>::NullaryExpr(
            get_shape().rows(), get_shape().cols(), [t, this](Eigen::Index i, Eigen::Index j) {
                return std::accumulate(m_matrices.begin(), m_matrices.end(), FP(0.0), [t, i, j](FP s, auto& m) {
                    return evaluate_intermediate<FP>(s + m.get_matrix_at(t)(i, j));
                });
            });
    }
 
    iterator begin()
    {
        return m_matrices.begin();
    }
    iterator end()
    {
        return m_matrices.end();
    }
    const_iterator begin() const
    {
        return m_matrices.begin();
    }
    const_iterator end() const
    {
        return m_matrices.end();
    }
 
    friend void PrintTo(const DampingMatrixExpressionGroup& self, std::ostream* os)
    {
        for (auto& m : self.m_matrices) {
            PrintTo(m, os);
            *os << '\n';
        }
    }
 
    template <class IOContext>
    void serialize(IOContext& io) const
    {
        auto obj = io.create_object("DampingMatrixExpressionGroup");
        obj.add_list("Matrices", begin(), end());
    }
 
protected:
    template <class IOContext, class Derived>
    static IOResult<Derived> deserialize(IOContext& io, Tag<Derived>)
    {
        auto obj = io.expect_object("DampingMatrixExpressionGroup");
        auto m   = obj.expect_list("Matrices", Tag<value_type>{});
        return apply(
            io,
            [](auto&& m_) -> IOResult<Derived> {
                //validation
                if (m_.empty()) {
                    return failure(StatusCode::InvalidValue,
                                   "DampingMatrixExpressionGroup must have at least one matrix.");
                }
                auto shape = m_[0].get_shape();
                for (size_t i = 1; i < m_.size(); ++i) {
                    if (m_[i].get_shape() != shape) {
                        return failure(StatusCode::InvalidValue,
                                       "Elements of DampingMatrixExpressionGroup must all have the same shape.");
                    }
                }
 
                return success(Derived(m_));
            },
            m);
    }
 
public:
    template <class IOContext>
    static IOResult<DampingMatrixExpressionGroup> deserialize(IOContext& io)
    {
        return deserialize(io, Tag<DampingMatrixExpressionGroup>{});
    }
 
private:
    std::vector<value_type> m_matrices;
};
 
template <typename FP>
class ContactMatrix : public DampingMatrixExpression<FP, SquareDampings<FP>>
{
public:
    using Base = DampingMatrixExpression<FP, SquareDampings<FP>>;
    using Base::Base;
 
    Eigen::Index get_num_groups() const
    {
        return Base::get_shape().rows();
    }
 
    template <class IOContext>
    static IOResult<ContactMatrix> deserialize(IOContext& io)
    {
        return Base::deserialize(io, Tag<ContactMatrix>{});
    }
};
 
template <typename FP>
class ContactMatrixGroup : public DampingMatrixExpressionGroup<FP, ContactMatrix<FP>>
{
public:
    using Base = DampingMatrixExpressionGroup<FP, ContactMatrix<FP>>;
    using Base::Base;
 
    Eigen::Index get_num_groups() const
    {
        return (*this)[0].get_num_groups();
    }
 
    template <class IOContext>
    static IOResult<ContactMatrixGroup> deserialize(IOContext& io)
    {
        return Base::deserialize(io, Tag<ContactMatrixGroup>{});
    }
};
 
} // namespace mio
 
#endif // MIO_EPI_CONTACT_MATRIX_H