Skip to content

Commit b44ce1d

Browse files
annawendlerjubicker
annawendler
and
jubicker
authored
1392 Improve examples and some documentation of IDE-SECIR model (#1402)
- Merge examples ide_secir and ide_secir_ageres to avoid code duplication - Add documentation in IDE-SECIR examples, for implemented formulas in general and in get_mean() of StateAgeFunction Co-authored-by: jubicker <[email protected]>
1 parent 9da2316 commit b44ce1d

File tree

7 files changed

+145
-232
lines changed

7 files changed

+145
-232
lines changed

‎cpp/examples/CMakeLists.txt‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -120,10 +120,6 @@ add_executable(ide_secir_example ide_secir.cpp)
120120
target_link_libraries(ide_secir_example PRIVATE memilio ide_secir)
121121
target_compile_options(ide_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
122122

123-
add_executable(ide_secir_ageres_example ide_secir_ageres.cpp)
124-
target_link_libraries(ide_secir_ageres_example PRIVATE memilio ide_secir)
125-
target_compile_options(ide_secir_ageres_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
126-
127123
add_executable(ide_secir_graph_example ide_secir_graph.cpp)
128124
target_link_libraries(ide_secir_graph_example PRIVATE memilio ide_secir)
129125
target_compile_options(ide_secir_graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/examples/ide_initialization.cpp‎

Lines changed: 27 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -54,41 +54,51 @@ std::string setup(int argc, char** argv)
5454

5555
int main(int argc, char** argv)
5656
{
57-
// This is a simple example to demonstrate how to set initial data for the IDE-SECIR model using real data.
57+
// This is a simple example to demonstrate how to set initial data for the IDE-SECIR model using reported data.
5858
// A default initialization is used if no filename is provided in the command line.
5959
// Have a look at the documentation of the set_initial_flows() function in models/ide_secir/parameters_io.h for a
6060
// description of how to download suitable data.
61-
// A valid filename could be for example "../../data/pydata/Germany/cases_all_germany_ma7.json" if the functionality to download real data is used.
62-
// The default parameters of the IDE-SECIR model are used, so that the simulation results are not realistic and are for demonstration purpose only.
61+
// A valid filename could be for example "../../data/pydata/Germany/cases_all_germany_ma7.json" if the functionality
62+
// to download reported data is used.
63+
// The default parameters of the IDE-SECIR model are used, so note that the simulation results are not realistic
64+
// and are for demonstration purpose only.
65+
66+
// Define simulation parameters.
67+
ScalarType t0 = 0.;
68+
ScalarType tmax = 2.;
69+
ScalarType dt = 0.01;
6370

6471
// Initialize model.
6572
size_t num_agegroups = 1;
66-
mio::CustomIndexArray<ScalarType, mio::AgeGroup> total_population =
73+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> total_population_init =
6774
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 80 * 1e6);
68-
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths = mio::CustomIndexArray<ScalarType, mio::AgeGroup>(
75+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths_init = mio::CustomIndexArray<ScalarType, mio::AgeGroup>(
6976
mio::AgeGroup(num_agegroups),
70-
0.); // The number of deaths will be overwritten if real data is used for initialization.
71-
ScalarType dt = 0.5;
72-
mio::isecir::Model model(mio::TimeSeries<ScalarType>((int)mio::isecir::InfectionTransition::Count),
73-
total_population, deaths, num_agegroups);
77+
0.); // The number of deaths will be overwritten if reported data is used for initialization.
78+
79+
mio::isecir::Model model(mio::TimeSeries<ScalarType>((size_t)mio::isecir::InfectionTransition::Count),
80+
total_population_init, deaths_init, num_agegroups);
7481

7582
// Check provided parameters.
7683
std::string filename = setup(argc, argv);
7784
if (filename.empty()) {
7885
std::cout << "You did not provide a valid filename. A default initialization is used." << std::endl;
7986

87+
// Create TimeSeries with num_transitions * num_agegroups elements where initial transitions needed for simulation
88+
// will be stored. We require values for the transitions for a sufficient number of time points before the start of
89+
// the simulation to initialize our model.
8090
using Vec = Eigen::VectorX<ScalarType>;
81-
mio::TimeSeries<ScalarType> init(num_agegroups * (size_t)mio::isecir::InfectionTransition::Count);
82-
init.add_time_point<Eigen::VectorX<ScalarType>>(
91+
mio::TimeSeries<ScalarType> transitions_init((size_t)mio::isecir::InfectionTransition::Count * num_agegroups);
92+
transitions_init.add_time_point<Eigen::VectorX<ScalarType>>(
8393
-7., Vec::Constant((size_t)mio::isecir::InfectionTransition::Count, 1. * dt));
84-
while (init.get_last_time() < -dt / 2.) {
85-
init.add_time_point(init.get_last_time() + dt,
86-
Vec::Constant((int)mio::isecir::InfectionTransition::Count, 1. * dt));
94+
while (transitions_init.get_last_time() < t0 - dt / 2.) {
95+
transitions_init.add_time_point(transitions_init.get_last_time() + dt,
96+
Vec::Constant((size_t)mio::isecir::InfectionTransition::Count, 1. * dt));
8797
}
88-
model.transitions = init;
98+
model.transitions = transitions_init;
8999
}
90100
else {
91-
// Use the real data for initialization.
101+
// Use reported data for initialization.
92102
// Here we assume that the file contains data without age resolution, hence we use read_confirmed_cases_noage()
93103
// for reading the data and mio::ConfirmedCasesNoAgeEntry as EntryType in set_initial_flows().
94104

@@ -112,7 +122,7 @@ int main(int argc, char** argv)
112122

113123
// Carry out simulation.
114124
mio::isecir::Simulation sim(model, dt);
115-
sim.advance(2.);
125+
sim.advance(tmax);
116126

117127
// Print results.
118128
sim.get_transitions().print_table({"S->E", "E->C", "C->I", "C->R", "I->H", "I->R", "H->U", "H->R", "U->D", "U->R"},

‎cpp/examples/ide_secir.cpp‎

Lines changed: 99 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
/*
22
* Copyright (C) 2020-2025 MEmilio
33
*
4-
* Authors: Anna Wendler, Lena Ploetzke
4+
* Authors: Anna Wendler, Lena Ploetzke, Hannah Tritzschak
55
*
66
* Contact: Martin J. Kuehn <[email protected]>
77
*
@@ -17,7 +17,6 @@
1717
* See the License for the specific language governing permissions and
1818
* limitations under the License.
1919
*/
20-
2120
#include "ide_secir/model.h"
2221
#include "ide_secir/infection_state.h"
2322
#include "ide_secir/simulation.h"
@@ -32,96 +31,138 @@
3231

3332
int main()
3433
{
34+
// This is a simple example to demonstrate how use the IDE-SECIR model.
35+
3536
using Vec = Eigen::VectorX<ScalarType>;
3637

37-
size_t num_agegroups = 1;
38+
// Define simulation parameters.
39+
ScalarType t0 = 0.;
40+
ScalarType tmax = 5.;
41+
ScalarType dt = 0.01; // The step size will stay constant throughout the simulation.
3842

39-
ScalarType tmax = 10;
40-
mio::CustomIndexArray<ScalarType, mio::AgeGroup> N =
41-
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 10000.);
42-
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths =
43-
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 13.10462213);
44-
ScalarType dt = 1.;
43+
// Define number of age groups.
44+
size_t num_agegroups = 2;
4545

46-
int num_transitions = (int)mio::isecir::InfectionTransition::Count;
46+
// Define initial values for the total population and number of deaths per age group.
47+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> total_population_init =
48+
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 1000.);
49+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths_init =
50+
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 6.);
4751

48-
// Create TimeSeries with num_transitions * num_agegroups elements where transitions needed for simulation will be
49-
// stored.
50-
mio::TimeSeries<ScalarType> init(num_transitions * num_agegroups);
52+
// Create TimeSeries with num_transitions * num_agegroups elements where initial transitions needed for simulation
53+
// will be stored. We require values for the transitions for a sufficient number of time points before the start of
54+
// the simulation to initialize our model.
55+
size_t num_transitions = (size_t)mio::isecir::InfectionTransition::Count;
56+
mio::TimeSeries<ScalarType> transitions_init(num_transitions * num_agegroups);
5157

52-
// Add time points for initialization of transitions.
58+
// Define vector of transitions that will be added as values to the time points of the TimeSeries transitions_init.
5359
Vec vec_init(num_transitions * num_agegroups);
54-
vec_init[(int)mio::isecir::InfectionTransition::SusceptibleToExposed] = 25.0;
55-
vec_init[(int)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 15.0;
56-
vec_init[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 8.0;
57-
vec_init[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered] = 4.0;
58-
vec_init[(int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] = 1.0;
59-
vec_init[(int)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered] = 4.0;
60-
vec_init[(int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] = 1.0;
61-
vec_init[(int)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = 1.0;
62-
vec_init[(int)mio::isecir::InfectionTransition::InfectedCriticalToDead] = 1.0;
63-
vec_init[(int)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = 1.0;
64-
60+
for (size_t group = 0; group < num_agegroups; ++group) {
61+
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::SusceptibleToExposed] = 25.0;
62+
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] =
63+
15.0;
64+
vec_init[group * num_transitions +
65+
(size_t)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 8.0;
66+
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered] =
67+
4.0;
68+
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] =
69+
1.0;
70+
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered] = 4.0;
71+
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] =
72+
1.0;
73+
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = 1.0;
74+
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedCriticalToDead] = 1.0;
75+
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = 1.0;
76+
}
77+
// Multiply vec_init with dt so that within a time interval of length 1, always the above number of
78+
// individuals are transitioning from one compartment to another, irrespective of the chosen time step size.
6579
vec_init = vec_init * dt;
66-
// Add initial time point to time series.
67-
init.add_time_point(-10, vec_init);
68-
// Add further time points until time 0.
69-
while (init.get_last_time() < -dt / 2) {
70-
init.add_time_point(init.get_last_time() + dt, vec_init);
80+
81+
// In this example, we will set the TransitionDistributions below. For these distributions, setting the initial time
82+
// point of the TimeSeries transitions_init at time -10 will give us a sufficient number of time points before t0=0.
83+
// For more information on this, we refer to the documentation of TransitionDistributions in
84+
// models/ide_secir/parameters.h.
85+
transitions_init.add_time_point(-10, vec_init);
86+
// Add further time points with distance dt until time t0.
87+
while (transitions_init.get_last_time() < t0 - dt / 2) {
88+
transitions_init.add_time_point(transitions_init.get_last_time() + dt, vec_init);
7189
}
7290

7391
// Initialize model.
74-
mio::isecir::Model model(std::move(init), N, deaths, num_agegroups);
92+
mio::isecir::Model model(std::move(transitions_init), total_population_init, deaths_init, num_agegroups);
93+
94+
// Uncomment one of the code blocks below to use a different method to initialize the model, based on a
95+
// given number of either Susceptibles or Recovered instead of using the TimeSeries transitions_init from above.
96+
97+
// Initialization method with given Susceptibles.
98+
// size_t num_infstates = (size_t)mio::isecir::InfectionState::Count;
99+
// for (size_t group = 0; group < num_agegroups; ++group) {
100+
// model.populations.get_last_value()[group * num_infstates + (size_t)mio::isecir::InfectionState::Susceptible] =
101+
// 900;
102+
// }
75103

76-
// Uncomment one of the two lines to use a different method to initialize the model using the TimeSeries init.
77-
// Initialization method with Susceptibles.
78-
// model.populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 1000;
79-
// Initialization method with Recovered.
80-
// model.populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0;
104+
// Initialization method with given Recovered.
105+
// size_t num_infstates = (size_t)mio::isecir::InfectionState::Count;
106+
// for (size_t group = 0; group < num_agegroups; ++group) {
107+
// model.populations.get_last_value()[group * num_infstates + (size_t)mio::isecir::InfectionState::Recovered] = 10;
108+
// }
81109

82110
// Set working parameters.
83-
mio::SmootherCosine<ScalarType> smoothcos(2.0);
84-
mio::StateAgeFunctionWrapper<ScalarType> delaydistribution(smoothcos);
85-
std::vector<mio::StateAgeFunctionWrapper<ScalarType>> vec_delaydistrib(num_transitions, delaydistribution);
86-
// TransitionDistribution is not used for SusceptibleToExposed. Therefore, the parameter can be set to any value.
87-
vec_delaydistrib[(int)mio::isecir::InfectionTransition::SusceptibleToExposed].set_distribution_parameter(-1.);
88-
vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]
89-
.set_distribution_parameter(4.0);
90111

91-
model.parameters.get<mio::isecir::TransitionDistributions>()[mio::AgeGroup(0)] = vec_delaydistrib;
112+
// TransitionDistributions
113+
// In the following, we explicitly set the TransitionDistributions for the first age group. If the model contains
114+
// more age groups, the default distributions are used for these age groups.
115+
mio::SmootherCosine<ScalarType> smoothcos1(3.0);
116+
mio::StateAgeFunctionWrapper<ScalarType> delaydistribution1(smoothcos1);
117+
std::vector<mio::StateAgeFunctionWrapper<ScalarType>> vec_delaydistrib1(num_transitions, delaydistribution1);
118+
// TransitionDistribution is not used for SusceptibleToExposed. Therefore, the parameter can be set to any value.
119+
vec_delaydistrib1[(size_t)mio::isecir::InfectionTransition::SusceptibleToExposed].set_distribution_parameter(-1.);
120+
model.parameters.get<mio::isecir::TransitionDistributions>()[mio::AgeGroup(0)] = vec_delaydistrib1;
92121

122+
// TransitionProbabilities
93123
std::vector<ScalarType> vec_prob(num_transitions, 0.5);
94124
// The following probabilities must be 1, as there is no other way to go.
95-
vec_prob[Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)] = 1;
96-
vec_prob[Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)] = 1;
97-
model.parameters.get<mio::isecir::TransitionProbabilities>()[mio::AgeGroup(0)] = vec_prob;
125+
vec_prob[(size_t)mio::isecir::InfectionTransition::SusceptibleToExposed] = 1;
126+
vec_prob[(size_t)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 1;
127+
for (mio::AgeGroup group = mio::AgeGroup(0); group < mio::AgeGroup(num_agegroups); ++group) {
128+
model.parameters.get<mio::isecir::TransitionProbabilities>()[group] = vec_prob;
129+
}
98130

131+
// Contact patterns
99132
mio::ContactMatrixGroup<ScalarType> contact_matrix = mio::ContactMatrixGroup<ScalarType>(1, num_agegroups);
100133
contact_matrix[0] =
101134
mio::ContactMatrix<ScalarType>(Eigen::MatrixX<ScalarType>::Constant(num_agegroups, num_agegroups, 10.));
102-
model.parameters.get<mio::isecir::ContactPatterns>() = mio::UncertainContactMatrix<ScalarType>(contact_matrix);
135+
model.parameters.get<mio::isecir::ContactPatterns>() = mio::UncertainContactMatrix(contact_matrix);
103136

137+
// Furhter epidemiological parameters
104138
mio::ExponentialSurvivalFunction<ScalarType> exponential(0.5);
105139
mio::StateAgeFunctionWrapper<ScalarType> prob(exponential);
106-
107-
model.parameters.get<mio::isecir::TransmissionProbabilityOnContact>()[mio::AgeGroup(0)] = prob;
108-
model.parameters.get<mio::isecir::RelativeTransmissionNoSymptoms>()[mio::AgeGroup(0)] = prob;
109-
model.parameters.get<mio::isecir::RiskOfInfectionFromSymptomatic>()[mio::AgeGroup(0)] = prob;
110-
140+
for (mio::AgeGroup group = mio::AgeGroup(0); group < mio::AgeGroup(num_agegroups); ++group) {
141+
model.parameters.get<mio::isecir::TransmissionProbabilityOnContact>()[group] = prob;
142+
model.parameters.get<mio::isecir::RelativeTransmissionNoSymptoms>()[group] = prob;
143+
model.parameters.get<mio::isecir::RiskOfInfectionFromSymptomatic>()[group] = prob;
144+
}
111145
model.parameters.set<mio::isecir::Seasonality>(0.1);
112-
// Start the simulation on the 40th day of a year (i.e. in February).
113-
model.parameters.set<mio::isecir::StartDay>(40);
146+
model.parameters.set<mio::isecir::StartDay>(
147+
40); // Start the simulation on the 40th day of a year (i.e. in February).
114148

149+
// Check if all model constraints regarding initial values and parameters are satisfied before simulating.
115150
model.check_constraints(dt);
116151

117152
// Carry out simulation.
118153
mio::isecir::Simulation sim(model, dt);
119154
sim.advance(tmax);
120155

156+
// Interpolate results to days.
121157
auto interpolated_results = mio::interpolate_simulation_result(sim.get_result(), dt / 2.);
122158

123-
interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
159+
// Print results. Note that the column labels are suitable for a simulation with two age groups and may need to be
160+
// adapted when the number of age groups is changed.
161+
// interpolated_results.print_table(
162+
// {"S1", "E1", "C1", "I1", "H1", "U1", "R1", "D1 ", "S2", "E2", "C2", "I2", "H2", "U2", "R2", "D2 "}, 16, 8);
124163
// Uncomment this line to print the transitions.
125-
// sim.get_transitions().print_table(
126-
// {"S->E 1", "E->C 1", "C->I 1", "C->R 1", "I->H 1", "I->R 1", "H->U 1", "H->R 1", "U->D 1", "U->R 1"}, 16, 8);
164+
// sim.get_transitions().print_table({"S->E 1", "E->C 1", "C->I 1", "C->R 1", "I->H 1", "I->R 1", "H->U 1",
165+
// "H->R 1", "U->D 1", "U->R 1", "S->E 2", "E->C 2", "C->I 2", "C->R 2",
166+
// "I->H 2", "I->R 2", "H->U 2", "H->R 2", "U->D 2", "U->R 2"},
167+
// 16, 8);
127168
}

0 commit comments

Comments
 (0)