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