Skip to content

Commit 01af90b

Browse files
HenrZujubickerreneSchm
authored
978 Local Feedback Mechanism (#1211)
- Add a local Feedback Mechanism based on Dönges et al. (doi.org/10.3389/fphy.2022.842180) Co-authored-by: jubicker <[email protected]> Co-authored-by: reneSchm <[email protected]>
1 parent 4f6c0d0 commit 01af90b

File tree

6 files changed

+666
-0
lines changed

6 files changed

+666
-0
lines changed

‎cpp/examples/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,10 @@ add_executable(ode_secir_contact_changes ode_secir_contact_changes.cpp)
7272
target_link_libraries(ode_secir_contact_changes PRIVATE memilio ode_secir)
7373
target_compile_options(ode_secir_contact_changes PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
7474

75+
add_executable(ode_secir_feedback ode_secir_feedback.cpp)
76+
target_link_libraries(ode_secir_feedback PRIVATE memilio ode_secir)
77+
target_compile_options(ode_secir_feedback PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
78+
7579
add_executable(ode_secirvvs_example ode_secirvvs.cpp)
7680
target_link_libraries(ode_secirvvs_example PRIVATE memilio ode_secirvvs)
7781
target_compile_options(ode_secirvvs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/examples/ode_secir_feedback.cpp‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
/*
2+
* Copyright (C) 2020-2025 MEmilio
3+
*
4+
* Authors: Henrik Zunker
5+
*
6+
* Contact: Martin J. Kuehn <[email protected]>
7+
*
8+
* Licensed under the Apache License, Version 2.0 (the "License");
9+
* you may not use this file except in compliance with the License.
10+
* You may obtain a copy of the License at
11+
*
12+
* http://www.apache.org/licenses/LICENSE-2.0
13+
*
14+
* Unless required by applicable law or agreed to in writing, software
15+
* distributed under the License is distributed on an "AS IS" BASIS,
16+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17+
* See the License for the specific language governing permissions and
18+
* limitations under the License.
19+
*/
20+
#include "ode_secir/model.h"
21+
#include "memilio/compartments/feedback_simulation.h"
22+
#include "memilio/utils/logging.h"
23+
24+
void initialize_model(mio::osecir::Model<double>& model, int total_population, double cont_freq)
25+
{
26+
model.parameters.set<mio::osecir::StartDay>(60);
27+
model.parameters.set<mio::osecir::Seasonality<double>>(0.2);
28+
29+
// time-related parameters
30+
model.parameters.get<mio::osecir::TimeExposed<double>>() = 3.2;
31+
model.parameters.get<mio::osecir::TimeInfectedNoSymptoms<double>>() = 2.0;
32+
model.parameters.get<mio::osecir::TimeInfectedSymptoms<double>>() = 5.8;
33+
model.parameters.get<mio::osecir::TimeInfectedSevere<double>>() = 9.5;
34+
model.parameters.get<mio::osecir::TimeInfectedCritical<double>>() = 7.1;
35+
36+
// Set transmission and isolation parameters
37+
model.parameters.get<mio::osecir::TransmissionProbabilityOnContact<double>>() = 0.05;
38+
model.parameters.get<mio::osecir::RelativeTransmissionNoSymptoms<double>>() = 0.7;
39+
model.parameters.get<mio::osecir::RecoveredPerInfectedNoSymptoms<double>>() = 0.09;
40+
model.parameters.get<mio::osecir::RiskOfInfectionFromSymptomatic<double>>() = 0.25;
41+
model.parameters.get<mio::osecir::MaxRiskOfInfectionFromSymptomatic<double>>() = 0.45;
42+
model.parameters.get<mio::osecir::TestAndTraceCapacity<double>>() = 35;
43+
model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<double>>() = 0.2;
44+
model.parameters.get<mio::osecir::CriticalPerSevere<double>>() = 0.25;
45+
model.parameters.get<mio::osecir::DeathsPerCritical<double>>() = 0.3;
46+
47+
// contact matrix
48+
mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::osecir::ContactPatterns<double>>();
49+
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
50+
51+
// initial population
52+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = 40;
53+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = 30;
54+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptomsConfirmed}] = 0;
55+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = 20;
56+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 0;
57+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSevere}] = 10;
58+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}] = 5;
59+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Recovered}] = 20;
60+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Dead}] = 0;
61+
model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
62+
total_population);
63+
64+
model.apply_constraints();
65+
}
66+
67+
void initialize_feedback(mio::FeedbackSimulation<double, mio::Simulation<double, mio::osecir::Model<double>>,
68+
mio::osecir::ContactPatterns<double>>& feedback_simulation)
69+
{
70+
// nominal ICU capacity
71+
feedback_simulation.get_parameters().template get<mio::NominalICUCapacity<double>>() = 10;
72+
73+
// ICU occupancy in the past for memory kernel
74+
auto& icu_occupancy = feedback_simulation.get_parameters().template get<mio::ICUOccupancyHistory<double>>();
75+
Eigen::VectorXd icu_day = Eigen::VectorXd::Constant(1, 1);
76+
const auto cutoff = static_cast<int>(feedback_simulation.get_parameters().template get<mio::GammaCutOff>());
77+
for (int t = -cutoff; t <= 0; ++t) {
78+
icu_occupancy.add_time_point(t, icu_day);
79+
}
80+
81+
// bounds for contact reduction measures
82+
feedback_simulation.get_parameters().template get<mio::ContactReductionMin<double>>() = {0.1};
83+
feedback_simulation.get_parameters().template get<mio::ContactReductionMax<double>>() = {0.8};
84+
}
85+
86+
int main()
87+
{
88+
// This example demonstrates the implementation of a feedback mechanism for a ODE SECIR model.
89+
// It shows how the perceived risk dynamically impacts contact reduction measures.
90+
// The feedback mechanism adjusts contact rates during simulation based on the perceived
91+
// risk which is calculated from the ICU occupancy using a memory kernel.
92+
mio::set_log_level(mio::LogLevel::warn);
93+
94+
const double tmax = 35;
95+
const int total_population = 1000;
96+
const double cont_freq = 10;
97+
98+
// create and initialize ODE model for a single age group
99+
mio::osecir::Model model(1);
100+
initialize_model(model, total_population, cont_freq);
101+
102+
// determine the index for the ICU state (InfectedCritical) for feedback mechanism
103+
auto icu_index = std::vector<size_t>{
104+
model.populations.get_flat_index({mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical})};
105+
106+
// create simulation objects: first a secir simulation, then a feedback simulation
107+
auto simulation = mio::osecir::Simulation<double, mio::Simulation<double, mio::osecir::Model<double>>>(model);
108+
auto feedback_simulation =
109+
mio::FeedbackSimulation<double, mio::Simulation<double, mio::osecir::Model<double>>,
110+
mio::osecir::ContactPatterns<double>>(std::move(simulation), icu_index);
111+
112+
// set up the parameters for the feedback simulation
113+
initialize_feedback(feedback_simulation);
114+
115+
// run the simulation with feedback mechanism
116+
feedback_simulation.advance(tmax);
117+
118+
// print the perceived risk and the final total population
119+
auto& perceived_risk = feedback_simulation.get_perceived_risk();
120+
perceived_risk.print_table({"Perceived Risk"});
121+
return 0;
122+
}

‎cpp/memilio/compartments/README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ Classes:
55
- 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.
66
- 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.
77
- 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.
8+
- 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)).
89

910
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.
1011

0 commit comments

Comments
 (0)