feedback_simulation.h

Go to the documentation of this file.

/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: 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_FEEDBACK_SIMULATION_H
#define MIO_COMPARTMENTS_FEEDBACK_SIMULATION_H
 
#include "memilio/epidemiology/age_group.h"
#include "memilio/epidemiology/damping_sampling.h"
#include "memilio/utils/parameter_set.h"
#include "memilio/utils/time_series.h"
#include "memilio/utils/uncertain_value.h"
 
namespace mio
{
template <typename FP>
struct ICUOccupancyHistory {
    using Type = mio::TimeSeries<FP>;
    static Type get_default(AgeGroup size)
    {
        return Type(size.get());
    }
    static std::string name()
    {
        return "ICUOccupancyHistory";
    }
};
 
template <typename FP>
struct GammaShapeParameter {
    using Type = FP;
    static Type get_default(AgeGroup)
    {
        return FP(6);
    }
    static std::string name()
    {
        return "GammaShapeParameter";
    }
};
 
template <typename FP>
struct GammaScaleParameter {
    using Type = FP;
    static Type get_default(AgeGroup)
    {
        return FP(0.4);
    }
    static std::string name()
    {
        return "GammaScaleParameter";
    }
};
 
struct GammaCutOff {
    using Type = size_t;
    static Type get_default(AgeGroup)
    {
        return Type(45);
    }
    static std::string name()
    {
        return "GammaCutOff";
    }
};
 
template <typename FP>
struct ContactReductionMax {
    using Type = std::vector<FP>;
    static Type get_default(AgeGroup)
    {
        return Type(1, FP(1.0));
    }
    static std::string name()
    {
        return "ContactReductionMax";
    }
};
 
template <typename FP>
struct ContactReductionMin {
    using Type = std::vector<FP>;
    static Type get_default(AgeGroup)
    {
        return Type(1, FP(0.0));
    }
    static std::string name()
    {
        return "ContactReductionMin";
    }
};
 
template <typename FP>
struct SoftPlusCurvatureParameter {
    using Type = FP;
    static Type get_default(AgeGroup)
    {
        return Type(0.1);
    }
    static std::string name()
    {
        return "SoftPlusCurvatureParameter";
    }
};
 
template <typename FP>
struct NominalICUCapacity {
    using Type = FP;
    static Type get_default(AgeGroup)
    {
        return Type(9.0);
    }
    static std::string name()
    {
        return "NominalICUCapacity";
    }
};
 
template <typename FP>
struct BlendingFactorLocal {
    using Type = FP;
    static Type get_default(AgeGroup)
    {
        return FP(1.0);
    }
    static std::string name()
    {
        return "BlendingFactorLocal";
    }
};
 
template <typename FP>
struct BlendingFactorRegional {
    using Type = FP;
    static Type get_default(AgeGroup)
    {
        return FP(0.0);
    }
    static std::string name()
    {
        return "BlendingFactorRegional";
    }
};
 
template <typename FP>
using FeedbackSimulationParameters =
    ParameterSet<ICUOccupancyHistory<FP>, GammaShapeParameter<FP>, GammaScaleParameter<FP>, GammaCutOff,
                 ContactReductionMax<FP>, ContactReductionMin<FP>, SoftPlusCurvatureParameter<FP>,
                 NominalICUCapacity<FP>, BlendingFactorLocal<FP>, BlendingFactorRegional<FP>>;
 
template <typename FP, typename Sim, typename ContactPatterns>
class FeedbackSimulation
{
public:
    using Model = typename Sim::Model;
 
    explicit FeedbackSimulation(Sim&& sim, const std::vector<size_t>& icu_indices)
        : m_simulation(std::move(sim))
        , m_icu_indices(icu_indices)
        , m_feedback_parameters(m_simulation.get_model().parameters.get_num_groups())
        , m_perceived_risk(static_cast<Eigen::Index>(m_simulation.get_model().parameters.get_num_groups().get()))
        , m_regional_icu_occupancy(
              static_cast<Eigen::Index>(m_simulation.get_model().parameters.get_num_groups().get()))
        , m_global_icu_occupancy(static_cast<Eigen::Index>(m_simulation.get_model().parameters.get_num_groups().get()))
    {
    }
 
    auto advance(const FP tmax, const FP dt_feedback = 1.0)
    {
        using std::min;
        FP t = m_simulation.get_result().get_last_time();
        while (t < tmax) {
            FP dt_eff = min<FP>(dt_feedback, tmax - t);
            apply_feedback(t + dt_eff);
            m_simulation.advance(t + dt_eff);
            add_icu_occupancy(t + dt_eff);
            t += dt_eff;
        }
        return m_simulation.get_result().get_last_value();
    }
 
    auto& get_result()
    {
        return m_simulation.get_result();
    }
 
    auto& get_model()
    {
        return m_simulation.get_model();
    }
 
    const auto& get_model() const
    {
        return m_simulation.get_model();
    }
 
    auto& get_perceived_risk() const
    {
        return m_perceived_risk;
    }
    auto& get_perceived_risk()
    {
        return m_perceived_risk;
    }
 
    auto& get_parameters()
    {
        return m_feedback_parameters;
    }
 
    FP calc_risk_perceived()
    {
        const auto& local_icu_occ = m_feedback_parameters.template get<ICUOccupancyHistory<FP>>();
        size_t num_time_points    = local_icu_occ.get_num_time_points();
        size_t n = std::min(static_cast<size_t>(num_time_points), m_feedback_parameters.template get<GammaCutOff>());
        FP perceived_risk = 0.0;
        const auto& a     = m_feedback_parameters.template get<GammaShapeParameter<FP>>();
        const auto& b     = m_feedback_parameters.template get<GammaScaleParameter<FP>>();
 
        const auto& blending_local    = m_feedback_parameters.template get<BlendingFactorLocal<FP>>();
        const auto& blending_regional = m_feedback_parameters.template get<BlendingFactorRegional<FP>>();
        const auto& blending_global   = 1 - blending_local - blending_regional;
 
        // assert that that each blending factor is positive
        assert(blending_local >= 0 && blending_regional >= 0 && blending_global >= 0 &&
               "Blending factors must be non-negative.");
 
        for (size_t i = num_time_points - n; i < num_time_points; ++i) {
            size_t day   = i - (num_time_points - n);
            FP gamma     = std::pow(b, a) * std::pow(day, a - 1) * std::exp(-b * day) / std::tgamma(a);
            FP perc_risk = 0.0;
 
            if (blending_local > 0) {
                perc_risk += blending_local * local_icu_occ.get_value(i).sum() /
                             m_feedback_parameters.template get<NominalICUCapacity<FP>>();
            }
            if (blending_regional > 0 && m_regional_icu_occupancy.get_num_time_points() > 0) {
                perc_risk += blending_regional * m_regional_icu_occupancy.get_value(i).sum() /
                             m_feedback_parameters.template get<NominalICUCapacity<FP>>();
            }
            if (blending_global > 0 && m_global_icu_occupancy.get_num_time_points() > 0) {
                perc_risk += blending_global * m_global_icu_occupancy.get_value(i).sum() /
                             m_feedback_parameters.template get<NominalICUCapacity<FP>>();
            }
 
            perc_risk = std::min(perc_risk, 1.0);
            perceived_risk += perc_risk * gamma;
        }
        return perceived_risk;
    }
 
    void add_icu_occupancy(FP t)
    {
        auto& model             = m_simulation.get_model();
        size_t num_groups       = static_cast<size_t>(model.parameters.get_num_groups());
        const auto& last_values = m_simulation.get_result().get_last_value();
        Eigen::VectorX<FP> icu_occ(num_groups);
        for (size_t i = 0; i < m_icu_indices.size(); ++i) {
            icu_occ[i] = last_values[m_icu_indices[i]];
        }
        // normalize by total population and scale to 100,000 inhabitants
        icu_occ = icu_occ / model.populations.get_total() * 100000;
        m_feedback_parameters.template get<ICUOccupancyHistory<FP>>().add_time_point(t, icu_occ);
    }
 
    void set_regional_icu_occupancy(const mio::TimeSeries<FP>& icu_regional)
    {
        m_regional_icu_occupancy = icu_regional;
    }
 
    void set_global_icu_occupancy(const mio::TimeSeries<FP>& icu_global)
    {
        m_global_icu_occupancy = icu_global;
    }
 
    void apply_feedback(FP t)
    {
        using std::exp;
        using std::log;
        using std::min;
 
        // get model parameters
        auto& params          = m_simulation.get_model().parameters;
        const auto num_groups = (size_t)params.get_num_groups();
 
        // get feedback parameters
        const auto& contactReductionMax = m_feedback_parameters.template get<ContactReductionMax<FP>>();
        const auto& contactReductionMin = m_feedback_parameters.template get<ContactReductionMin<FP>>();
        const auto epsilon              = m_feedback_parameters.template get<SoftPlusCurvatureParameter<FP>>();
 
        const size_t num_locations = contactReductionMax.size();
 
        // calculate local perceived risk
        FP perceived_risk = calc_risk_perceived();
 
        // add the perceived risk to the time series
        m_perceived_risk.add_time_point(t, Eigen::VectorX<FP>::Constant(num_groups, perceived_risk));
 
        auto group_weights_all = Eigen::VectorX<FP>::Constant(num_groups, 1.0);
 
        // define a lambda function to create a damping sampling object given a reduction factor and location
        auto set_contact_reduction = [=](FP val, size_t location) {
            auto v = mio::UncertainValue<FP>(val);
            return mio::DampingSampling(v, mio::DampingLevel(0), mio::DampingType(0), mio::SimulationTime<FP>(t),
                                        std::vector<size_t>{location}, group_weights_all);
        };
 
        auto& dampings = params.template get<ContactPatterns>().get_dampings();
 
        // for each location, compute the effective contact reduction factor and add the corresponding damping sample
        for (size_t loc = 0; loc < num_locations; ++loc) {
 
            // compute the effective reduction factor using a softplus function.
            FP reduc_fac = (contactReductionMax[loc] - contactReductionMin[loc]) * epsilon *
                               log(exp(perceived_risk / epsilon) + 1.0) +
                           contactReductionMin[loc];
 
            // clamp the reduction factor to the maximum allowed value.
            reduc_fac = min<FP>(reduc_fac, contactReductionMax[loc]);
 
            // generate the damping for the current location.
            auto damping = set_contact_reduction(reduc_fac, loc);
 
            // add the damping to the contact pattern.
            dampings.push_back(damping);
        }
        // update the contact matrices after all dampings have been added.
        params.template get<ContactPatterns>().make_matrix();
    }
 
private:
    Sim m_simulation; 
    std::vector<size_t> m_icu_indices; 
    FeedbackSimulationParameters<FP> m_feedback_parameters; 
    mio::TimeSeries<FP> m_perceived_risk; 
    mio::TimeSeries<FP> m_regional_icu_occupancy; 
    mio::TimeSeries<FP> m_global_icu_occupancy; 
};
 
} // namespace mio
 
#endif // MIO_COMPARTMENTS_FEEDBACK_SIMULATION_H