Skip to content

Commit 27c9ce8

Browse files
jubickerHenrZu
jubicker
and
HenrZu
authored
1239 implement temporal hybrid model (#1262)
- Implemented a temporal-hybrid model that switches between two models during the simulation according to a predefined condition - Added a singlewell potential for the diffusive ABM - Implemented conversion functions for a temporal dABM-SMM hybrid model - Implemented conversion functions for a temporal dABM-oSECIR hybrid model - Added tests for new functionality - Added example for temporal dABM-oSECIR hybrid model Co-authored-by: Henrik Zunker <69154294+HenrZu@users.noreply.github.com>
1 parent 6c83abd commit 27c9ce8

27 files changed

+1862
-36
lines changed

‎cpp/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,8 @@ set(CMAKE_PDB_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin")
7575
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin")
7676
set(CMAKE_INSTALL_RPATH "${CMAKE_BINARY_DIR}/lib" "${CMAKE_BINARY_DIR}/bin")
7777

78+
file(TO_CMAKE_PATH "${PROJECT_SOURCE_DIR}/.." MEMILIO_BASE_DIR)
79+
7880
# code coverage analysis
7981
# Note: this only works under linux and with make
8082
# Ninja creates different directory names which do not work together with this scrupt
@@ -188,6 +190,7 @@ if(MEMILIO_BUILD_MODELS)
188190
add_subdirectory(models/sde_seirvv)
189191
add_subdirectory(models/graph_abm)
190192
add_subdirectory(models/smm)
193+
add_subdirectory(models/hybrid)
191194
endif()
192195

193196
if(MEMILIO_BUILD_EXAMPLES)

‎cpp/examples/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,10 @@ add_executable(smm_example smm.cpp)
148148
target_link_libraries(smm_example PRIVATE memilio smm)
149149
target_compile_options(smm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
150150

151+
add_executable(temporal_hybrid_example temporal_hybrid_dabm_osecir.cpp)
152+
target_link_libraries(temporal_hybrid_example PRIVATE memilio hybrid ode_secir d_abm)
153+
target_compile_options(temporal_hybrid_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
154+
151155
if(MEMILIO_HAS_JSONCPP)
152156
add_executable(ode_secir_read_graph_example ode_secir_read_graph.cpp)
153157
target_link_libraries(ode_secir_read_graph_example PRIVATE memilio ode_secir)

‎cpp/examples/d_abm.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ enum class InfectionState
4141
int main()
4242
{
4343
//Example how to run a simulation of the diffusive ABM using the quadwell potential
44-
using Model = mio::dabm::Model<QuadWellModel<InfectionState>>;
44+
using Model = mio::dabm::Model<QuadWell<InfectionState>>;
4545
std::vector<Model::Agent> agents(1000);
4646
//Random variables for initialization of agents' position and infection state
4747
auto& pos_rng = mio::UniformDistribution<double>::get_instance();

‎cpp/examples/temporal_hybrid_dabm_osecir.cpp‎

Lines changed: 218 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,218 @@
1+
/*
2+
* Copyright (C) 2020-2025 MEmilio
3+
*
4+
* Authors: Julia Bicker
5+
*
6+
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
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+
21+
#include "d_abm/model.h"
22+
#include "memilio/data/analyze_result.h"
23+
#include "memilio/epidemiology/age_group.h"
24+
#include "memilio/utils/logging.h"
25+
#include "memilio/utils/time_series.h"
26+
#include "models/hybrid/temporal_hybrid_model.h"
27+
#include "models/d_abm/simulation.h"
28+
#include "models/d_abm/single_well.h"
29+
#include "memilio/compartments/simulation.h"
30+
#include "models/ode_secir/model.h"
31+
#include "memilio/utils/random_number_generator.h"
32+
#include "memilio/epidemiology/adoption_rate.h"
33+
#include "memilio/geography/regions.h"
34+
#include "ode_secir/infection_state.h"
35+
#include "models/hybrid/conversion_functions.cpp"
36+
#include <cstddef>
37+
#include <vector>
38+
39+
int main()
40+
{
41+
mio::set_log_level(mio::LogLevel::warn);
42+
// Simple example to demonstrate how to run a simulation using temporal-hybrid model combining the diffusive ABM and the ODE-SECIR-model.
43+
// As condition to switch between models we use a threshold of 20 infected individuals (For <20 Infected the ABM is used and for >=20 Infected the ODE-Model is used).
44+
45+
using ABM = mio::dabm::Model<SingleWell<mio::osecir::InfectionState>>;
46+
using ODE = mio::osecir::Model<double>;
47+
48+
//Initialize ABM population
49+
std::vector<ABM::Agent> agents(1000);
50+
//Random variables used to initialize agents' position and infection state
51+
auto& pos_sampler = mio::UniformDistribution<double>::get_instance();
52+
auto& stat_sampler = mio::DiscreteDistribution<size_t>::get_instance();
53+
//Infection state distribution
54+
std::vector<double> infection_state_dist{0.99, 0.01, 0., 0., 0., 0., 0., 0.};
55+
//Sample agents' position and infection state
56+
for (auto& a : agents) {
57+
//Agents' positions are equally distributed in [-2, 2] x [-2, 2]
58+
a.position = Eigen::Vector2d{pos_sampler(mio::thread_local_rng(), -2., 2.),
59+
pos_sampler(mio::thread_local_rng(), -2., 2.)};
60+
//Agents' infection states are sampled from infection_state_dist
61+
a.status =
62+
static_cast<mio::osecir::InfectionState>(stat_sampler(mio::thread_local_rng(), infection_state_dist));
63+
}
64+
//Transmission parameters used for both models
65+
const double contact_frequency = 10, trans_prob_on_contact = 0.06, time_E = 3., time_Ins = 2.5, time_Isy = 5.2,
66+
time_Isev = 9., time_Icri = 7.2, mu_Ins_R = 0.2, mu_Isy_Isev = 0.1, mu_Isev_Icri = 0.1,
67+
mu_Icri_D = 0.2;
68+
//Initialize ABM adoption rates
69+
std::vector<mio::AdoptionRate<mio::osecir::InfectionState>> adoption_rates;
70+
//Second-order adoption rate (S->E)
71+
adoption_rates.push_back(
72+
{mio::osecir::InfectionState::Susceptible,
73+
mio::osecir::InfectionState::Exposed,
74+
mio::regions::Region(0),
75+
contact_frequency * trans_prob_on_contact,
76+
{{mio::osecir::InfectionState::InfectedNoSymptoms, 1}, {mio::osecir::InfectionState::InfectedSymptoms, 1}}});
77+
//First-order adoption rates
78+
//E->Ins
79+
adoption_rates.push_back({mio::osecir::InfectionState::Exposed,
80+
mio::osecir::InfectionState::InfectedNoSymptoms,
81+
mio::regions::Region(0),
82+
1. / time_E,
83+
{}});
84+
//Ins->Isy
85+
adoption_rates.push_back({mio::osecir::InfectionState::InfectedNoSymptoms,
86+
mio::osecir::InfectionState::InfectedSymptoms,
87+
mio::regions::Region(0),
88+
(1 - mu_Ins_R) / time_Ins,
89+
{}});
90+
//Ins->R
91+
adoption_rates.push_back({mio::osecir::InfectionState::InfectedNoSymptoms,
92+
mio::osecir::InfectionState::Recovered,
93+
mio::regions::Region(0),
94+
mu_Ins_R / time_Ins,
95+
{}});
96+
//Isy->Isev
97+
adoption_rates.push_back({mio::osecir::InfectionState::InfectedSymptoms,
98+
mio::osecir::InfectionState::InfectedSevere,
99+
mio::regions::Region(0),
100+
mu_Isy_Isev / time_Isy,
101+
{}});
102+
//Isy->R
103+
adoption_rates.push_back({mio::osecir::InfectionState::InfectedSymptoms,
104+
mio::osecir::InfectionState::Recovered,
105+
mio::regions::Region(0),
106+
(1 - mu_Isy_Isev) / time_Isy,
107+
{}});
108+
//Isev->Icri
109+
adoption_rates.push_back({mio::osecir::InfectionState::InfectedSevere,
110+
mio::osecir::InfectionState::InfectedCritical,
111+
mio::regions::Region(0),
112+
mu_Isev_Icri / time_Isev,
113+
{}});
114+
//Isev->R
115+
adoption_rates.push_back({mio::osecir::InfectionState::InfectedSevere,
116+
mio::osecir::InfectionState::Recovered,
117+
mio::regions::Region(0),
118+
(1 - mu_Isev_Icri) / time_Isev,
119+
{}});
120+
//Icri->R
121+
adoption_rates.push_back({mio::osecir::InfectionState::InfectedCritical,
122+
mio::osecir::InfectionState::Recovered,
123+
mio::regions::Region(0),
124+
(1 - mu_Icri_D) / time_Icri,
125+
{}});
126+
//Icri->D
127+
adoption_rates.push_back({mio::osecir::InfectionState::InfectedCritical,
128+
mio::osecir::InfectionState::Dead,
129+
mio::regions::Region(0),
130+
mu_Icri_D / time_Icri,
131+
{}});
132+
//Interaction radius and noise
133+
double interaction_radius = 0.4, noise = 0.5;
134+
ABM abm(agents, adoption_rates, interaction_radius, noise,
135+
{mio::osecir::InfectionState::InfectedSevere, mio::osecir::InfectionState::InfectedCritical,
136+
mio::osecir::InfectionState::Dead});
137+
138+
//As we start modeling with the ABM, we don't need to initialize the population for the ODE-model
139+
//Initialize ODE model parameters
140+
ODE ode(1);
141+
ode.parameters.get<mio::osecir::TimeExposed<double>>()[mio::AgeGroup(0)] = time_E;
142+
ode.parameters.get<mio::osecir::TimeInfectedNoSymptoms<double>>()[mio::AgeGroup(0)] = time_Ins;
143+
ode.parameters.get<mio::osecir::TimeInfectedSymptoms<double>>()[mio::AgeGroup(0)] = time_Isy;
144+
ode.parameters.get<mio::osecir::TimeInfectedSevere<double>>()[mio::AgeGroup(0)] = time_Isev;
145+
ode.parameters.get<mio::osecir::TimeInfectedCritical<double>>()[mio::AgeGroup(0)] = time_Icri;
146+
ode.parameters.get<mio::osecir::TransmissionProbabilityOnContact<double>>()[mio::AgeGroup(0)] =
147+
trans_prob_on_contact;
148+
ode.parameters.get<mio::osecir::RecoveredPerInfectedNoSymptoms<double>>()[mio::AgeGroup(0)] = mu_Ins_R;
149+
ode.parameters.get<mio::osecir::SeverePerInfectedSymptoms<double>>()[mio::AgeGroup(0)] = mu_Isy_Isev;
150+
ode.parameters.get<mio::osecir::CriticalPerSevere<double>>()[mio::AgeGroup(0)] = mu_Isev_Icri;
151+
ode.parameters.get<mio::osecir::DeathsPerCritical<double>>()[mio::AgeGroup(0)] = mu_Icri_D;
152+
ode.apply_constraints();
153+
mio::ContactMatrixGroup& contact_matrix = ode.parameters.get<mio::osecir::ContactPatterns<double>>();
154+
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, contact_frequency));
155+
156+
//Set t0 and internal dt for each model
157+
double t0 = 0;
158+
double dt = 0.1;
159+
160+
//Create simulations
161+
auto sim_abm = mio::dabm::Simulation(abm, t0, dt);
162+
auto sim_ode = mio::Simulation(ode, t0, dt);
163+
164+
const auto result_fct_abm = [](const mio::dabm::Simulation<SingleWell<mio::osecir::InfectionState>>& sim,
165+
double /*t*/) {
166+
return sim.get_result();
167+
};
168+
169+
const auto result_fct_ode = [](const mio::Simulation<double, ODE>& sim, double /*t*/) {
170+
return sim.get_result();
171+
};
172+
173+
//Create hybrid simulation
174+
double dt_switch = 0.2;
175+
mio::hybrid::TemporalHybridSimulation<decltype(sim_abm), decltype(sim_ode), mio::TimeSeries<double>,
176+
mio::TimeSeries<double>>
177+
hybrid_sim(std::move(sim_abm), std::move(sim_ode), result_fct_abm, result_fct_ode, true, t0, dt_switch);
178+
179+
//Define switching condition
180+
const auto condition = [](const mio::TimeSeries<double>& result_abm, const mio::TimeSeries<double>& result_ode,
181+
bool abm_used) {
182+
if (abm_used) {
183+
auto& last_value = result_abm.get_last_value().eval();
184+
if ((last_value[(int)mio::osecir::InfectionState::Exposed] +
185+
last_value[(int)mio::osecir::InfectionState::InfectedNoSymptoms] +
186+
last_value[(int)mio::osecir::InfectionState::InfectedNoSymptomsConfirmed] +
187+
last_value[(int)mio::osecir::InfectionState::InfectedSymptoms] +
188+
last_value[(int)mio::osecir::InfectionState::InfectedSymptomsConfirmed] +
189+
last_value[(int)mio::osecir::InfectionState::InfectedSevere] +
190+
last_value[(int)mio::osecir::InfectionState::InfectedCritical]) > 20) {
191+
return true;
192+
}
193+
}
194+
else {
195+
auto& last_value = result_ode.get_last_value().eval();
196+
if ((last_value[(int)mio::osecir::InfectionState::Exposed] +
197+
last_value[(int)mio::osecir::InfectionState::InfectedNoSymptoms] +
198+
last_value[(int)mio::osecir::InfectionState::InfectedNoSymptomsConfirmed] +
199+
last_value[(int)mio::osecir::InfectionState::InfectedSymptoms] +
200+
last_value[(int)mio::osecir::InfectionState::InfectedSymptomsConfirmed] +
201+
last_value[(int)mio::osecir::InfectionState::InfectedSevere] +
202+
last_value[(int)mio::osecir::InfectionState::InfectedCritical]) <= 20) {
203+
return true;
204+
}
205+
}
206+
return false;
207+
};
208+
209+
//Simulate for 10 days
210+
hybrid_sim.advance(10., condition);
211+
212+
auto ts_abm = hybrid_sim.get_result_model1();
213+
auto ts_ode = hybrid_sim.get_result_model2();
214+
215+
//Print result time series
216+
auto ts = mio::interpolate_simulation_result(mio::merge_time_series(ts_abm, ts_ode).value());
217+
ts.print_table({"S", "E", "Ins", "Ins_confirmed", "Isy", "Isy_confirmed", "Isev", "Icri", "R", "D"});
218+
}

‎cpp/memilio/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -109,6 +109,7 @@ add_library(memilio
109109
utils/mioomp.cpp
110110
utils/string_literal.h
111111
utils/type_list.h
112+
utils/base_dir.h
112113
)
113114

114115
target_include_directories(memilio PUBLIC

‎cpp/memilio/config_internal.h.in‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,4 +31,6 @@
3131
#cmakedefine MEMILIO_ENABLE_OPENMP
3232
#cmakedefine MEMILIO_ENABLE_PROFILING
3333

34+
const char* const MEMILIO_BASE_DIR = "${MEMILIO_BASE_DIR}/";
35+
3436
#endif

‎cpp/memilio/data/analyze_result.h‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222

2323
#include "memilio/utils/time_series.h"
2424
#include "memilio/mobility/metapopulation_mobility_instant.h"
25+
#include "memilio/io/io.h"
2526

2627
#include <functional>
2728
#include <vector>
@@ -169,6 +170,97 @@ double result_distance_2norm(const std::vector<mio::TimeSeries<double>>& result1
169170
return std::sqrt(norm_sqr);
170171
}
171172

173+
/**
174+
* @brief This function merges two TimeSeries by copying their time points and values to a new TimeSeries in the correct order.
175+
* If both TimeSeries have values for the same time point, their values are either added or only one value is taken.
176+
* @param[in] ts1 First TimeSeries.
177+
* @param[in] ts2 Second TimeSeries.
178+
* @param[in] add_values Boolean specifying whether the values should be added if both TimeSeries contain the same time point. If false, the value of just the first TimeSeries is taken.
179+
* @tparam FP A floating point type.
180+
* @return A TimeSeries containing all time points and values from both input TimeSeries.
181+
*/
182+
template <class FP>
183+
IOResult<TimeSeries<FP>> merge_time_series(const TimeSeries<FP>& ts1, const TimeSeries<FP>& ts2,
184+
bool add_values = false)
185+
{
186+
TimeSeries<FP> merged_ts(ts1.get_num_elements());
187+
if (ts1.get_num_elements() != ts2.get_num_elements()) {
188+
log_error("TimeSeries have a different number of elements.");
189+
return failure(mio::StatusCode::InvalidValue);
190+
}
191+
else {
192+
Eigen::Index t1_iterator = 0;
193+
Eigen::Index t2_iterator = 0;
194+
bool t1_finished = false;
195+
bool t2_finished = false;
196+
while (!t1_finished || !t2_finished) {
197+
if (!t1_finished) {
198+
if (ts1.get_time(t1_iterator) < ts2.get_time(t2_iterator) ||
199+
t2_finished) { // Current time point of first TimeSeries is smaller than current time point of second TimeSeries or second TimeSeries has already been copied entirely
200+
merged_ts.add_time_point(ts1.get_time(t1_iterator), ts1.get_value(t1_iterator));
201+
t1_iterator += 1;
202+
}
203+
else if (!t2_finished && ts1.get_time(t1_iterator) ==
204+
ts2.get_time(t2_iterator)) { // Both TimeSeries have the current time point
205+
if (add_values) {
206+
merged_ts.add_time_point(ts1.get_time(t1_iterator),
207+
ts1.get_value(t1_iterator) + ts2.get_value(t2_iterator));
208+
}
209+
else {
210+
merged_ts.add_time_point(ts1.get_time(t1_iterator), ts1.get_value(t1_iterator));
211+
log_warning("Both TimeSeries have values for t={}. The value of the first TimeSeries is used",
212+
ts1.get_time(t1_iterator));
213+
}
214+
t1_iterator += 1;
215+
t2_iterator += 1;
216+
if (t2_iterator >=
217+
ts2.get_num_time_points()) { // Check if all values of second TimeSeries have been copied
218+
t2_finished = true;
219+
t2_iterator = ts2.get_num_time_points() - 1;
220+
}
221+
}
222+
if (t1_iterator >=
223+
ts1.get_num_time_points()) { // Check if all values of first TimeSeries have been copied
224+
t1_finished = true;
225+
t1_iterator = ts1.get_num_time_points() - 1;
226+
}
227+
}
228+
if (!t2_finished) {
229+
if (ts2.get_time(t2_iterator) < ts1.get_time(t1_iterator) ||
230+
t1_finished) { // Current time point of second TimeSeries is smaller than current time point of first TimeSeries or first TimeSeries has already been copied entirely
231+
merged_ts.add_time_point(ts2.get_time(t2_iterator), ts2.get_value(t2_iterator));
232+
t2_iterator += 1;
233+
}
234+
else if (!t1_finished && ts2.get_time(t2_iterator) ==
235+
ts1.get_time(t1_iterator)) { // Both TimeSeries have the current time point
236+
if (add_values) {
237+
merged_ts.add_time_point(ts1.get_time(t1_iterator),
238+
ts1.get_value(t1_iterator) + ts2.get_value(t2_iterator));
239+
}
240+
else {
241+
merged_ts.add_time_point(ts1.get_time(t1_iterator), ts1.get_value(t1_iterator));
242+
log_warning("Both TimeSeries have values for t={}. The value of the first TimeSeries is used",
243+
ts1.get_time(t1_iterator));
244+
}
245+
t1_iterator += 1;
246+
t2_iterator += 1;
247+
if (t1_iterator >=
248+
ts1.get_num_time_points()) { // Check if all values of first TimeSeries have been copied
249+
t1_finished = true;
250+
t1_iterator = ts1.get_num_time_points() - 1;
251+
}
252+
}
253+
if (t2_iterator >=
254+
ts2.get_num_time_points()) { // Check if all values of second TimeSeries have been copied
255+
t2_finished = true;
256+
t2_iterator = ts2.get_num_time_points() - 1;
257+
}
258+
}
259+
}
260+
}
261+
return success(merged_ts);
262+
}
263+
172264
} // namespace mio
173265

174266
#endif //MEMILIO_DATA_ANALYZE_RESULT_H

0 commit comments

Comments
 (0)