diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index e7c5cb99f4..69d45d2be8 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -72,6 +72,10 @@ add_executable(ode_secir_contact_changes ode_secir_contact_changes.cpp) target_link_libraries(ode_secir_contact_changes PRIVATE memilio ode_secir) target_compile_options(ode_secir_contact_changes PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ode_secir_feedback ode_secir_feedback.cpp) +target_link_libraries(ode_secir_feedback PRIVATE memilio ode_secir) +target_compile_options(ode_secir_feedback PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(ode_secirvvs_example ode_secirvvs.cpp) target_link_libraries(ode_secirvvs_example PRIVATE memilio ode_secirvvs) target_compile_options(ode_secirvvs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/examples/ode_secir_feedback.cpp b/cpp/examples/ode_secir_feedback.cpp new file mode 100644 index 0000000000..eea8885993 --- /dev/null +++ b/cpp/examples/ode_secir_feedback.cpp @@ -0,0 +1,122 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* 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. +*/ +#include "ode_secir/model.h" +#include "memilio/compartments/feedback_simulation.h" +#include "memilio/utils/logging.h" + +void initialize_model(mio::osecir::Model& model, int total_population, double cont_freq) +{ + model.parameters.set(60); + model.parameters.set>(0.2); + + // time-related parameters + model.parameters.get>() = 3.2; + model.parameters.get>() = 2.0; + model.parameters.get>() = 5.8; + model.parameters.get>() = 9.5; + model.parameters.get>() = 7.1; + + // Set transmission and isolation parameters + model.parameters.get>() = 0.05; + model.parameters.get>() = 0.7; + model.parameters.get>() = 0.09; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.45; + model.parameters.get>() = 35; + model.parameters.get>() = 0.2; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.3; + + // contact matrix + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); + + // initial population + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = 40; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = 30; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = 20; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 0; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSevere}] = 10; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}] = 5; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Recovered}] = 20; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Dead}] = 0; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, + total_population); + + model.apply_constraints(); +} + +void initialize_feedback(mio::FeedbackSimulation>, + mio::osecir::ContactPatterns>& feedback_simulation) +{ + // nominal ICU capacity + feedback_simulation.get_parameters().template get>() = 10; + + // ICU occupancy in the past for memory kernel + auto& icu_occupancy = feedback_simulation.get_parameters().template get>(); + Eigen::VectorXd icu_day = Eigen::VectorXd::Constant(1, 1); + const auto cutoff = static_cast(feedback_simulation.get_parameters().template get()); + for (int t = -cutoff; t <= 0; ++t) { + icu_occupancy.add_time_point(t, icu_day); + } + + // bounds for contact reduction measures + feedback_simulation.get_parameters().template get>() = {0.1}; + feedback_simulation.get_parameters().template get>() = {0.8}; +} + +int main() +{ + // This example demonstrates the implementation of a feedback mechanism for a ODE SECIR model. + // It shows how the perceived risk dynamically impacts contact reduction measures. + // The feedback mechanism adjusts contact rates during simulation based on the perceived + // risk which is calculated from the ICU occupancy using a memory kernel. + mio::set_log_level(mio::LogLevel::warn); + + const double tmax = 35; + const int total_population = 1000; + const double cont_freq = 10; + + // create and initialize ODE model for a single age group + mio::osecir::Model model(1); + initialize_model(model, total_population, cont_freq); + + // determine the index for the ICU state (InfectedCritical) for feedback mechanism + auto icu_index = std::vector{ + model.populations.get_flat_index({mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical})}; + + // create simulation objects: first a secir simulation, then a feedback simulation + auto simulation = mio::osecir::Simulation>>(model); + auto feedback_simulation = + mio::FeedbackSimulation>, + mio::osecir::ContactPatterns>(std::move(simulation), icu_index); + + // set up the parameters for the feedback simulation + initialize_feedback(feedback_simulation); + + // run the simulation with feedback mechanism + feedback_simulation.advance(tmax); + + // print the perceived risk and the final total population + auto& perceived_risk = feedback_simulation.get_perceived_risk(); + perceived_risk.print_table({"Perceived Risk"}); + return 0; +} diff --git a/cpp/memilio/compartments/README.md b/cpp/memilio/compartments/README.md index 069e1fafef..b603667eae 100644 --- a/cpp/memilio/compartments/README.md +++ b/cpp/memilio/compartments/README.md @@ -5,6 +5,7 @@ Classes: - FlowModel: A CompartmentalModel defined by the flows between populations. Requires the additional template parameter Flows, which is a [TypeList](../utils/type_list.h) of [Flow](../utils/flow.h). Instead of defining `get_derivatives`, the function `get_flows` has to be implemented for a FlowModel. - Simulation: Template class that runs the simulation using a specified compartment model. Can be derived from to implement behavior that cannot be modeled inside the usual compartment flow structure. - FlowSimulation: Template class derived from Simulation, that additionally computes the flows between populations. It requires a FlowModel. While the compartment and flow results can be changed outside of the simulation, it is required that both still have the same number of time points to further advance this simulation. +- FeedbackSimulation: Template class derived from Simulation (or FlowSimulation) that extends the standard simulation functionality with a feedback mechanism. FeedbackSimulation holds its own parameters (using a specialized ParameterSet) and a perceived risk time series. It overrides the simulation advance routine to compute and apply contact adjustments at regular intervals. The contact adjustments are computed by first deriving the perceived risk from the historical ICU occupancy—normalized by the nominal ICU capacity—and then transforming that risk via a softplus function into an effective damping factor to reduce contacts as described in Dönges et al. ([doi.org/10.3389/fphy.2022.842180](https://doi.org/10.3389/fphy.2022.842180)). Note that currently, only the last value of the simulation result (i.e. the TimeSeries obtained from `get_result()`) is used to advance the simulation. Thus, any changes to the population of the simulation's model will have no effect on the simulation result. diff --git a/cpp/memilio/compartments/feedback_simulation.h b/cpp/memilio/compartments/feedback_simulation.h new file mode 100644 index 0000000000..4f0d035175 --- /dev/null +++ b/cpp/memilio/compartments/feedback_simulation.h @@ -0,0 +1,378 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* 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 FEEDBACK_SIMULATION_H +#define FEEDBACK_SIMULATION_H + +#include "memilio/compartments/simulation.h" +#include "memilio/utils/time_series.h" +#include "memilio/utils/parameter_set.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/uncertain_value.h" +#include "memilio/epidemiology/damping_sampling.h" + +namespace mio +{ +/** + * @brief Daily local ICU occupancy per age group. + * + * This parameter stores all historic (local) ICU occupancy values normalized per 100,000 inhabitants. + * The values are extracted from the simulation results based on the provided ICU compartment indices + * and are updated at each feedback step. + */ +template +struct ICUOccupancyHistory { + using Type = mio::TimeSeries; + static Type get_default(AgeGroup size) + { + return Type(size.get()); + } + static std::string name() + { + return "ICUOccupancyHistory"; + } +}; + +/** + * @brief Shape parameter of the gamma distribution. + */ +template +struct GammaShapeParameter { + using Type = FP; + static Type get_default(AgeGroup) + { + return FP(6); + } + static std::string name() + { + return "GammaShapeParameter"; + } +}; + +/** + * @brief Scale parameter of the gamma distribution. + */ +template +struct GammaScaleParameter { + using Type = FP; + static Type get_default(AgeGroup) + { + return FP(0.4); + } + static std::string name() + { + return "GammaScaleParameter"; + } +}; + +/** + * @brief Number of days in the past considered for the gamma distribution. + */ +struct GammaCutOff { + using Type = size_t; + static Type get_default(AgeGroup) + { + return 45; + } + static std::string name() + { + return "GammaCutOff"; + } +}; + +/** + * @brief Maximum allowed contact reduction factors per location. + */ +template +struct ContactReductionMax { + using Type = std::vector; + static Type get_default(AgeGroup) + { + return Type(1, FP(1.0)); + } + static std::string name() + { + return "ContactReductionMax"; + } +}; + +/** + * @brief Minimum allowed contact reduction factors per location. + */ +template +struct ContactReductionMin { + using Type = std::vector; + static Type get_default(AgeGroup) + { + return Type(1, FP(0.0)); + } + static std::string name() + { + return "ContactReductionMin"; + } +}; + +/** + * @brief Soft-plus curvature parameter for contact adjustment. + */ +template +struct SoftPlusCurvatureParameter { + using Type = FP; + static Type get_default(AgeGroup) + { + return FP(0.1); + } + static std::string name() + { + return "SoftPlusCurvatureParameter"; + } +}; + +/** + * @brief Nominal ICU capacity. + */ +template +struct NominalICUCapacity { + using Type = FP; + static Type get_default(AgeGroup) + { + return FP(9.0); + } + static std::string name() + { + return "NominalICUCapacity"; + } +}; + +template +using FeedbackSimulationParameters = + ParameterSet, GammaShapeParameter, GammaScaleParameter, GammaCutOff, + ContactReductionMax, ContactReductionMin, SoftPlusCurvatureParameter, + NominalICUCapacity>; + +/** + * @brief A generic feedback simulation extending existing simulations with a feedback mechanism. + * + * This class wraps any simulation (e.g. Simulation or FlowSimulation) and applies additional + * feedback logic in a model-independent way. Model-specific details—such as the ICU compartments— + * are provided via arguments. + * + * @tparam FP The floating point type. + * @tparam Sim The simulation type. + * @tparam ContactPatterns The model-specific contact patterns type. + */ +template +class FeedbackSimulation +{ +public: + /** + * @brief Constructs the FeedbackSimulation by taking ownership of an existing simulation instance. + * + * @param sim The simulation instance to be extended with feedback mechanism. + * @param icu_indices A vector of indices indicating ICU compartments for specific model. + */ + explicit FeedbackSimulation(Sim&& sim, const std::vector& 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(m_simulation.get_model().parameters.get_num_groups())) + { + } + + /** + * @brief Advances the simulation until tmax while applying feedback at fixed intervals. + * + * The simulation is advanced in steps of dt_feedback. At each step, feedback + * is applied, then the simulation is advanced, and afterwards the current ICU occupancy is stored. + * + * Note that the simulation may make additional substeps depending on its own + * timestep dt. When using fixed-step integrators, dt_feedback should be an integer multiple of + * the simulation timestep dt. + * + * @param tmax The maximum simulation time. + * @param dt_feedback The feedback time step (default 1.0). + * @return The result in the last time step of the simulation. + */ + auto advance(const FP tmax, const FP dt_feedback = 1.0) + { + FP t = m_simulation.get_result().get_last_time(); + while (t < tmax) { + FP dt_eff = std::min(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(); + } + + /** + * @brief Returns the simulation result. + */ + auto& get_result() + { + return m_simulation.get_result(); + } + + /** + * @brief Returns the model used in the simulation. + */ + auto& get_model() + { + return m_simulation.get_model(); + } + + /** + * @brief Returns the perceived risk time series. + */ + auto& get_perceived_risk() const + { + return m_perceived_risk; + } + auto& get_perceived_risk() + { + return m_perceived_risk; + } + + /** + * @brief Returns the local feedback parameters. + */ + auto& get_parameters() + { + return m_feedback_parameters; + } + + /** + * @brief Calculates the perceived risk based on ICU occupancy data. + * + * The risk is computed as a weighted sum of ICU occupancy over recent time points using a gamma-distribution + * acting as memory kernel. + * + * @return The computed (local) perceived risk. + */ + FP calc_risk_perceived() + { + const auto& icu_occ = m_feedback_parameters.template get>(); + size_t num_time_points = icu_occ.get_num_time_points(); + size_t n = std::min(static_cast(num_time_points), m_feedback_parameters.template get()); + FP perceived_risk = 0.0; + const auto& a = m_feedback_parameters.template get>(); + const auto& b = m_feedback_parameters.template get>(); + 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 = icu_occ.get_value(i).sum() / m_feedback_parameters.template get>(); + perc_risk = std::min(perc_risk, 1.0); + perceived_risk += perc_risk * gamma; + } + return perceived_risk; + } + + /** + * @brief Adds the current (local) ICU occupancy into the parameter containing all historical ICU occupancy values. + * + * This function use the latest simulation results and extracts the ICU occupancy + * based on the indices given. The occupancy values are then normalized to 100,000 inhabitants, + * and then stored (as a new time point) in the ICUOccupancyHistory parameter. + * + * @param t The current simulation time at which the ICU occupancy is recorded. + */ + void add_icu_occupancy(FP t) + { + auto& model = m_simulation.get_model(); + size_t num_groups = static_cast(model.parameters.get_num_groups()); + const auto& last_values = m_simulation.get_result().get_last_value(); + Eigen::VectorXd 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>().add_time_point(t, icu_occ); + } + + /** + * @brief Transforms the perceived risk into a contact reduction factor and applies it to the contact patterns. + * + * This function computes a contact reduction factor for each location based on the perceived risk. For smooth transitions, + * we use a softplus function as described in Dönges et al. (doi.org/10.3389/fphy.2022.842180). + * + * @param t The simulation time at which the feedback adjustments are applied. + */ + void apply_feedback(FP t) + { + // 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>(); + const auto& contactReductionMin = m_feedback_parameters.template get>(); + const auto epsilon = m_feedback_parameters.template get>(); + + 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::VectorXd::Constant(num_groups, perceived_risk)); + + auto group_weights_all = Eigen::VectorXd::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(val); + return mio::DampingSampling(v, mio::DampingLevel(0), mio::DampingType(0), mio::SimulationTime(t), + std::vector{location}, group_weights_all); + }; + + auto& dampings = params.template get().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. + ScalarType reduc_fac = (contactReductionMax[loc] - contactReductionMin[loc]) * epsilon * + std::log(std::exp(perceived_risk / epsilon) + 1.0) + + contactReductionMin[loc]; + + // clamp the reduction factor to the maximum allowed value. + reduc_fac = std::min(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().make_matrix(); + } + +private: + Sim m_simulation; ///< The simulation instance. + std::vector m_icu_indices; ///< The ICU compartment indices from the model. + FeedbackSimulationParameters m_feedback_parameters; ///< The feedback parameters. + mio::TimeSeries m_perceived_risk; ///< The perceived risk time series. +}; + +} // namespace mio + +#endif // FEEDBACK_SIMULATION_H diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 44e913f442..74e59412fa 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -49,6 +49,7 @@ set(TESTSOURCES test_custom_index_array.cpp test_d_abm_model.cpp test_flows.cpp + test_feedback.cpp test_parameter_set.cpp test_matrix_shape.cpp test_damping_sampling.cpp diff --git a/cpp/tests/test_feedback.cpp b/cpp/tests/test_feedback.cpp new file mode 100644 index 0000000000..93a78515b7 --- /dev/null +++ b/cpp/tests/test_feedback.cpp @@ -0,0 +1,160 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* 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. +*/ +#include "load_test_data.h" +#include "memilio/compartments/compartmentalmodel.h" +#include "memilio/compartments/flow_simulation.h" +#include "memilio/compartments/simulation.h" +#include "memilio/compartments/feedback_simulation.h" +#include "ode_seir/infection_state.h" +#include "ode_seir/model.h" +#include "ode_seir/parameters.h" + +#include "gtest/gtest.h" + +class TestFeedbackSimulation : public ::testing::Test +{ +protected: + void SetUp() override + { + // Create an ODE SEIR model with one age group. + auto model = mio::oseir::Model(1); + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 1000; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 1000; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 1000; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}, 10000); + model.parameters.set>(1.0); + model.parameters.set>(5.2); + model.parameters.set>(2); + + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); + + // ICU compartment index + icu_indices = {2}; + + // Create a simulation based on the model. + auto sim = mio::Simulation>(model); + + // Create the FeedbackSimulation by moving the simulation and providing the ICU indices. + feedback_sim = + std::make_unique>, + mio::oseir::ContactPatterns>>(std::move(sim), icu_indices); + } + + std::vector icu_indices; + std::unique_ptr>, + mio::oseir::ContactPatterns>> + feedback_sim; +}; + +// Test that advancing the simulation adds ICU occupancy and perceived risk entries. +TEST_F(TestFeedbackSimulation, AdvanceAddICUAndRisk) +{ + // initially, the perceived risk time series should be empty + EXPECT_EQ(feedback_sim->get_perceived_risk().get_num_time_points(), 0); + + // advance the simulation to t = 2.0 in two steps (dt_feedback = 1.0). + feedback_sim->advance(2.0, 1.0); + + // ICU occupancy time series should now have 2 time points. + const auto& icu_occ = feedback_sim->get_parameters().template get>(); + EXPECT_EQ(icu_occ.get_num_time_points(), 2); + + // similarly, the perceived risk time series should also have 2 entries. + EXPECT_EQ(feedback_sim->get_perceived_risk().get_num_time_points(), 2); +} + +TEST_F(TestFeedbackSimulation, CalcPerceivedRisk) +{ + // set GammaShapeParameter to 1 and GammaScaleParameter to 1 so that gamma becomes exp(-day) + auto& fb_params = feedback_sim->get_parameters(); + fb_params.template get>() = 1; + fb_params.template get>() = 1; + fb_params.template get>() = 10; + + // add a single ICU occupancy time points with value 2. + Eigen::VectorXd icu_value(1); + icu_value << 2; + auto& icu_occ = fb_params.template get>(); + icu_occ.add_time_point(0.0, icu_value); + double risk = feedback_sim->calc_risk_perceived(); + // For day 0, we have gamma = exp(0) = 1 and therefore perc_risk = 2 / 10 = 0.2. + EXPECT_NEAR(risk, 0.2, 1e-10); + + // Add another time point with value 2. + icu_occ.add_time_point(1.0, icu_value); + + // Here, we have gamma = exp(-1) = 0.367879441 and therefore expect a risk of 0.2 + 0.07357588823428847. + risk = feedback_sim->calc_risk_perceived(); + EXPECT_NEAR(risk, 0.27357588823428847, 1e-10); +} + +TEST_F(TestFeedbackSimulation, ApplyFeedback) +{ + // bounds for contact reduction measures + auto& feedback_params = feedback_sim->get_parameters(); + feedback_params.template get>() = std::vector{0.8}; + feedback_params.template get>() = std::vector{0.2}; + + // get initial number of damping samples. + auto& contact_patterns = feedback_sim->get_model().parameters.get>(); + size_t initial_dampings = contact_patterns.get_dampings().size(); + + // set all historical ICU occupancy values to 100.0, to have maximum perceived risk. + auto& icu_occ = feedback_params.template get>(); + Eigen::VectorXd icu_value(1); + icu_value << 100.0; + for (int t = -45; t <= 0; ++t) { + icu_occ.add_time_point(t, icu_value); + } + + // apply_feedback at time t = 0. + feedback_sim->apply_feedback(0.0); + + // number of new damping samples added should equal the number of locations, + // as given by the size of the ContactReductionMax vector. + size_t num_locations = feedback_sim->get_parameters().template get>().size(); + EXPECT_EQ(contact_patterns.get_dampings().size(), initial_dampings + num_locations); + + // contact reduction should be near the maximum allowed contact reduction. + EXPECT_NEAR(contact_patterns.get_dampings().back().get_value().value(), + feedback_params.template get>()[0], 1e-3); +} + +TEST_F(TestFeedbackSimulation, AddICUOccupancy) +{ + // check that the ICU occupancy time series is empty + auto& feedback_params = feedback_sim->get_parameters(); + EXPECT_EQ(feedback_params.template get>().get_num_time_points(), 0); + + // add ICU occupancy for t = 1.0. + feedback_sim->add_icu_occupancy(1.0); + + // check that the ICU occupancy time series now has one time point. + EXPECT_EQ(feedback_params.template get>().get_num_time_points(), 1); + // check that stored value is equal to the inital value relative to 100,000 + auto stored_val = feedback_params.template get>().get_value(0); + auto expected_val = + feedback_sim->get_model().populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] / + feedback_sim->get_model().populations.get_total() * 100000; + EXPECT_EQ(stored_val[0], expected_val); +}