Skip to content

Commit adba576

Browse files
kilianvolmerkilianvolmer
committed
Merge branch 'main' of github.com:SciCompMod/memilio-tutorials
2 parents 034c8b5 + 81c23e2 commit adba576

6 files changed

Lines changed: 868 additions & 0 deletions

File tree

‎cpp-tutorials/CMakeLists.txt‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -55,12 +55,21 @@ target_link_libraries(tutorial2 PRIVATE memilio ode_secir Boost::filesystem)
5555
add_executable(tutorial3 tutorial3.cpp)
5656
target_link_libraries(tutorial3 PRIVATE memilio ode_secir Boost::filesystem)
5757

58+
add_executable(exercise3 exercises/exercise3.cpp)
59+
target_link_libraries(exercise3 PRIVATE memilio ode_secir Boost::filesystem)
60+
5861
add_executable(tutorial5 tutorial5.cpp)
5962
target_link_libraries(tutorial5 PRIVATE memilio ode_secir Boost::filesystem)
6063

64+
add_executable(exercise5 exercises/exercise5.cpp)
65+
target_link_libraries(exercise5 PRIVATE memilio ode_secir Boost::filesystem)
66+
6167
add_executable(tutorial7 tutorial7.cpp)
6268
target_link_libraries(tutorial7 PRIVATE memilio ode_secir Boost::filesystem)
6369

70+
add_executable(exercise7 exercises/exercise7.cpp)
71+
target_link_libraries(exercise7 PRIVATE memilio ode_secir Boost::filesystem)
72+
6473
add_executable(tutorial10 tutorial10.cpp)
6574
target_link_libraries(tutorial10 PRIVATE memilio ode_secir Boost::filesystem)
6675

‎cpp-tutorials/exercises/exercise3.cpp‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
#include "ode_secir/model.h"
2+
#include "memilio/compartments/simulation.h"
3+
#include "memilio/data/analyze_result.h"
4+
5+
int main()
6+
{
7+
// In the previous tutorial, we created, initialized and simulated MEmilio's ODE-based SECIR-type model with one (age) group. In this tutorial, we will show how to incorporate non-pharmaceutical interventions (NPIs) through the use of `Dampings` in the ODE-based SECIR-type model.
8+
9+
// *** Set up model. ***
10+
// First we create and initialize a SECIR-type model with one age group. For a detailed description on taht, see Tutorial 1.
11+
size_t num_agegroups = 1;
12+
ScalarType total_population = 100000;
13+
ScalarType t0 = 0;
14+
ScalarType tmax = 100;
15+
ScalarType dt = 0.1;
16+
17+
// Create model
18+
mio::osecir::Model<ScalarType> model(num_agegroups);
19+
20+
// Set infection state stay times (in days)
21+
model.parameters.get<mio::osecir::TimeExposed<ScalarType>>() = 3.2;
22+
model.parameters.get<mio::osecir::TimeInfectedNoSymptoms<ScalarType>>() = 2.;
23+
model.parameters.get<mio::osecir::TimeInfectedSymptoms<ScalarType>>() = 6.;
24+
model.parameters.get<mio::osecir::TimeInfectedSevere<ScalarType>>() = 12.;
25+
model.parameters.get<mio::osecir::TimeInfectedCritical<ScalarType>>() = 8.;
26+
// Set infection state transition probabilities
27+
model.parameters.get<mio::osecir::RelativeTransmissionNoSymptoms<ScalarType>>() = 0.67;
28+
model.parameters.get<mio::osecir::TransmissionProbabilityOnContact<ScalarType>>() = 0.1;
29+
model.parameters.get<mio::osecir::RecoveredPerInfectedNoSymptoms<ScalarType>>() = 0.2;
30+
model.parameters.get<mio::osecir::RiskOfInfectionFromSymptomatic<ScalarType>>() = 0.25;
31+
model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<ScalarType>>() = 0.2;
32+
model.parameters.get<mio::osecir::CriticalPerSevere<ScalarType>>() = 0.25;
33+
model.parameters.get<mio::osecir::DeathsPerCritical<ScalarType>>() = 0.3;
34+
//Set contact frequency
35+
ScalarType contact_frequency = 10;
36+
mio::ContactMatrixGroup<ScalarType>& contact_matrix =
37+
model.parameters.get<mio::osecir::ContactPatterns<ScalarType>>();
38+
contact_matrix[0] = mio::ContactMatrix<ScalarType>(Eigen::MatrixX<ScalarType>::Constant(1, 1, contact_frequency));
39+
40+
// EXERCISE: Please set the model's contact frequency to 10.
41+
42+
// After the model initialization, we add a contact reduction (`Damping`) that represents an NPI like e.g. mask wearing or social distancing. Dampings are a factor applied to the contact frequency and can be added to the model at fixed simulation time points before simulating. They have a *Level* and a *Type*. A damping with a given level and type replaces the previously active one with the same level and type, while all currently active dampings of one level and different types are summed up. If two dampings have different levels (independent of the type) they are combined multiplicatively. In the following we apply a `Damping` of 0.9 after 10 days and another damping of 0.6 after 20 days which means that the contacts are reduced by 10% and 40%, respectively. To always retain a minimum level of contacts, a minimum contact frequency can be set that is never deceeded. In our example we set this minimum contact rate to 0.
43+
contact_matrix[0].add_damping(0.9, mio::SimulationTime<ScalarType>(10.));
44+
contact_matrix[0].add_damping(0.6, mio::SimulationTime<ScalarType>(20.));
45+
46+
//EXERCISE: Please create a damping that replaces the 10% reduction after 10 days by a 20% reduction. Additionally, add a damping after 40 days that increases the contact rate by 50%.
47+
48+
// Again, we start with 0.5% of the population initially in `Exposed` and 0.5% initially in `InfectedNoSymptoms` while the remaining 99% is `Susceptible`.
49+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = 0.005 * total_population;
50+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = 0.005 * total_population;
51+
model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
52+
total_population);
53+
54+
// *** Simulate model. ***
55+
mio::TimeSeries<ScalarType> result = mio::osecir::simulate<ScalarType>(t0, tmax, dt, model);
56+
// Interpolate time series to full days.
57+
auto interpolated_result = mio::interpolate_simulation_result(result);
58+
59+
// *** Print results. ***
60+
interpolated_result.print_table({"S", "E", "C", "C_confirmed", "I", "I_confirmed", "H", "U", "R", "D"}, 12, 4);
61+
}

‎cpp-tutorials/exercises/exercise5.cpp‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
#include "ode_secir/model.h"
2+
#include "memilio/compartments/simulation.h"
3+
#include "memilio/data/analyze_result.h"
4+
5+
int main()
6+
{
7+
// In Tutorial 1, we created, initialized and simulated MEmilio's ODE-based SECIR-type model without any sociodemographic resolution. All ODE-based models have the possibility to add an arbitrary number of sociodemographic groups which can represent certain certain risk groups, like vaccination or age groups. Adding those groups can have a relevant impact on the simulation outcome. If for example older people have a higher risk of severe and critical infections, that can have an impact on ICU occupancy.
8+
// In the following, we initialize and simulate an ODE-based SECIR-type model with three age groups.
9+
10+
// *** Set up model. ***
11+
// First we create and initialize a SECIR-type model with three age groups. For a detailed description on that, see Tutorial 1.
12+
size_t num_agegroups = 3;
13+
ScalarType total_population = 100000;
14+
ScalarType t0 = 0;
15+
ScalarType tmax = 100;
16+
ScalarType dt = 0.1;
17+
18+
// Create model with three age groups
19+
mio::osecir::Model<ScalarType> model(num_agegroups);
20+
21+
// Now, we have to set the epidemiological model parameters which are dependent on age group. A list of all parameters can be found at https://memilio.readthedocs.io/en/latest/cpp/models/osecir.html.
22+
// We choose an increasing risk of severe and critical infections for age group 2 and 3 compared to age group 1. The other parameters are equal for all age groups.
23+
24+
for (size_t i = 0; i < num_agegroups; i++) {
25+
// Set infection state stay times (in days)
26+
model.parameters.get<mio::osecir::TimeExposed<ScalarType>>()[mio::AgeGroup(i)] = 3.2;
27+
model.parameters.get<mio::osecir::TimeInfectedNoSymptoms<ScalarType>>()[mio::AgeGroup(i)] = 2.;
28+
model.parameters.get<mio::osecir::TimeInfectedSymptoms<ScalarType>>()[mio::AgeGroup(i)] = 6.;
29+
model.parameters.get<mio::osecir::TimeInfectedSevere<ScalarType>>()[mio::AgeGroup(i)] = 12.;
30+
model.parameters.get<mio::osecir::TimeInfectedCritical<ScalarType>>()[mio::AgeGroup(i)] = 9.;
31+
// Set infection state transition probabilities
32+
model.parameters.get<mio::osecir::TransmissionProbabilityOnContact<ScalarType>>()[mio::AgeGroup(i)] = 0.1;
33+
model.parameters.get<mio::osecir::RelativeTransmissionNoSymptoms<ScalarType>>()[mio::AgeGroup(i)] = 0.67;
34+
model.parameters.get<mio::osecir::RecoveredPerInfectedNoSymptoms<ScalarType>>()[mio::AgeGroup(i)] = 0.2;
35+
model.parameters.get<mio::osecir::RiskOfInfectionFromSymptomatic<ScalarType>>()[mio::AgeGroup(i)] = 0.25;
36+
model.parameters.get<mio::osecir::DeathsPerCritical<ScalarType>>()[mio::AgeGroup(i)] = 0.3;
37+
}
38+
39+
// The groups have an increasing risk of severe and critical infections
40+
model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<ScalarType>>()[mio::AgeGroup(0)] = 0.2;
41+
model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<ScalarType>>()[mio::AgeGroup(1)] = 0.2 * 1.5;
42+
model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<ScalarType>>()[mio::AgeGroup(2)] = 0.2 * 2;
43+
model.parameters.get<mio::osecir::CriticalPerSevere<ScalarType>>()[mio::AgeGroup(0)] = 0.25;
44+
model.parameters.get<mio::osecir::CriticalPerSevere<ScalarType>>()[mio::AgeGroup(1)] = 0.25 * 1.5;
45+
model.parameters.get<mio::osecir::CriticalPerSevere<ScalarType>>()[mio::AgeGroup(2)] = 0.25 * 2;
46+
47+
//EXERCISE: Please increase the percentage of severe cases per symptomatic cases in age group 2 by 10%.
48+
49+
//Set contact frequency
50+
ScalarType contact_frequency = 10;
51+
mio::ContactMatrixGroup<ScalarType>& contact_matrix =
52+
model.parameters.get<mio::osecir::ContactPatterns<ScalarType>>();
53+
contact_matrix[0] = mio::ContactMatrix<ScalarType>(
54+
Eigen::MatrixX<ScalarType>::Constant(num_agegroups, num_agegroups, contact_frequency));
55+
56+
// In this example, 0.5% of the population is initially in `Exposed` and 0.5% is initially in `InfectedNoSymptoms` while the remaining 99% of the population is `Susceptible`. The fractions are equally distributed to all age groups by:
57+
for (size_t i = 0; i < num_agegroups; i++) {
58+
model.populations[{mio::AgeGroup(i), mio::osecir::InfectionState::Exposed}] =
59+
0.005 * total_population / num_agegroups;
60+
model.populations[{mio::AgeGroup(i), mio::osecir::InfectionState::InfectedNoSymptoms}] =
61+
0.005 * total_population / num_agegroups;
62+
model.populations.set_difference_from_group_total<mio::AgeGroup>(
63+
{mio::AgeGroup(i), mio::osecir::InfectionState::Susceptible}, total_population / num_agegroups);
64+
}
65+
66+
// *** Simulate model. ***
67+
// After having initialized the model, dynamics can be simulated. The simulation output is a time series containing the evolution of all compartments per age group over time. In the following we simulate the model from `t0` to `tmax` with initial step size `dt` and subsequently print the time series result:
68+
mio::TimeSeries<ScalarType> result = mio::osecir::simulate<ScalarType>(t0, tmax, dt, model);
69+
// Interpolate time series to full days.
70+
auto interpolated_result = mio::interpolate_simulation_result(result);
71+
72+
// *** Print results. ***
73+
interpolated_result.print_table();
74+
}

‎cpp-tutorials/exercises/exercise7.cpp‎

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
#include "memilio/config.h"
2+
#include "ode_secir/model.h"
3+
#include "ode_secir/infection_state.h"
4+
#include "ode_secir/parameters.h"
5+
#include "memilio/mobility/metapopulation_mobility_instant.h"
6+
#include "memilio/compartments/simulation.h"
7+
#include "memilio/data/analyze_result.h"
8+
9+
int main()
10+
{
11+
// In the previous tutorials, we saw how to set up and run an age-resolved ODE-based SECIR-type model. However, one limiting assumption of simple ODE-based models is the assumption of homogenous mixing within the population. To overcome this limitation and incorporate spatial heterogeneity, in this example we show how to use MEmilio's graph-based metapopulation model. This model realizes mobility between regions via graph edges, while every region is represented by a graph node containing it's own ODE-based model.
12+
13+
// *** Set up model. ***
14+
// We set the simulation start time `t0`, the end time `tmax` and the initial step size `dt` as:
15+
ScalarType t0 = 0;
16+
ScalarType tmax = 100;
17+
ScalarType dt = 0.1;
18+
19+
// Next, we need to specify the parameters. We will initialize a metapopulation model with two regions. The total population as well as the epidemiological parameters will be the same for both regions.
20+
ScalarType total_population_per_region = 100000;
21+
22+
// We use a model with three age groups for both regions:
23+
size_t num_agegroups = 3;
24+
// Create model with three age groups
25+
mio::osecir::Model<ScalarType> model(num_agegroups);
26+
27+
// Now, we have to set the epidemiological model parameters which are dependent on age group. A list of all parameters can be found at https://memilio.readthedocs.io/en/latest/cpp/models/osecir.html.
28+
// We choose an increasing risk of severe and critical infections for age group 2 and 3 compared to age group 1. The other parameters are equal for all age groups.
29+
30+
for (size_t i = 0; i < num_agegroups; i++) {
31+
// Set infection state stay times (in days)
32+
model.parameters.get<mio::osecir::TimeExposed<ScalarType>>()[mio::AgeGroup(i)] = 3.2;
33+
model.parameters.get<mio::osecir::TimeInfectedNoSymptoms<ScalarType>>()[mio::AgeGroup(i)] = 2.;
34+
model.parameters.get<mio::osecir::TimeInfectedSymptoms<ScalarType>>()[mio::AgeGroup(i)] = 6.;
35+
model.parameters.get<mio::osecir::TimeInfectedSevere<ScalarType>>()[mio::AgeGroup(i)] = 12.;
36+
model.parameters.get<mio::osecir::TimeInfectedCritical<ScalarType>>()[mio::AgeGroup(i)] = 9.;
37+
// Set infection state transition probabilities
38+
model.parameters.get<mio::osecir::TransmissionProbabilityOnContact<ScalarType>>()[mio::AgeGroup(i)] = 0.1;
39+
model.parameters.get<mio::osecir::RelativeTransmissionNoSymptoms<ScalarType>>()[mio::AgeGroup(i)] = 0.67;
40+
model.parameters.get<mio::osecir::RecoveredPerInfectedNoSymptoms<ScalarType>>()[mio::AgeGroup(i)] = 0.2;
41+
model.parameters.get<mio::osecir::RiskOfInfectionFromSymptomatic<ScalarType>>()[mio::AgeGroup(i)] = 0.25;
42+
model.parameters.get<mio::osecir::DeathsPerCritical<ScalarType>>()[mio::AgeGroup(i)] = 0.3;
43+
}
44+
45+
// The groups have an increasing risk of severe and critical infections
46+
model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<ScalarType>>()[mio::AgeGroup(0)] = 0.2;
47+
model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<ScalarType>>()[mio::AgeGroup(1)] = 0.2 * 1.5;
48+
model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<ScalarType>>()[mio::AgeGroup(2)] = 0.2 * 2;
49+
model.parameters.get<mio::osecir::CriticalPerSevere<ScalarType>>()[mio::AgeGroup(0)] = 0.25;
50+
model.parameters.get<mio::osecir::CriticalPerSevere<ScalarType>>()[mio::AgeGroup(1)] = 0.25 * 1.5;
51+
model.parameters.get<mio::osecir::CriticalPerSevere<ScalarType>>()[mio::AgeGroup(2)] = 0.25 * 2;
52+
53+
// Set contact frequency
54+
ScalarType contact_frequency = 10;
55+
mio::ContactMatrixGroup<ScalarType>& contact_matrix =
56+
model.parameters.get<mio::osecir::ContactPatterns<ScalarType>>();
57+
contact_matrix[0] = mio::ContactMatrix<ScalarType>(
58+
Eigen::MatrixX<ScalarType>::Constant(num_agegroups, num_agegroups, contact_frequency));
59+
60+
// Next, we create the graph via:
61+
mio::Graph<mio::SimulationNode<ScalarType, mio::osecir::Simulation<ScalarType>>, mio::MobilityEdge<ScalarType>>
62+
graph;
63+
64+
// we want to add two regions (nodes) to the graph, therefore we need two copies of the model
65+
auto model_region1 = model;
66+
auto model_region2 = model;
67+
68+
// In the graph-based metapopulation model, every graph node gets it's own ODE-based model which is copied when adding a graph node and handing the model to it as parameter. Therefore we can choose different initial conditions (as well as differing parameters) for different graph nodes. In our example, we simulate two regions with only one region having initially infected individuals. We choose 1% initially infected for that region while the other region starts with a totally susceptible population.
69+
// The model compartments for the first node are initialized via:
70+
for (size_t i = 0; i < num_agegroups; i++) {
71+
model_region1.populations[{mio::AgeGroup(i), mio::osecir::InfectionState::Exposed}] =
72+
0.005 * total_population_per_region / num_agegroups;
73+
model_region1.populations[{mio::AgeGroup(i), mio::osecir::InfectionState::InfectedNoSymptoms}] =
74+
0.005 * total_population_per_region / num_agegroups;
75+
model_region1.populations.set_difference_from_group_total<mio::AgeGroup>(
76+
{mio::AgeGroup(i), mio::osecir::InfectionState::Susceptible}, total_population_per_region / num_agegroups);
77+
}
78+
79+
// We set the second model's initial populations to totally susceptible
80+
for (size_t i = 0; i < num_agegroups; i++) {
81+
model_region2.populations.set_difference_from_group_total<mio::AgeGroup>(
82+
{mio::AgeGroup(i), mio::osecir::InfectionState::Susceptible}, total_population_per_region / num_agegroups);
83+
}
84+
85+
// After having initialized the models, we add two nodes (regions) to the graph
86+
graph.add_node(0, model_region1, t0, dt);
87+
88+
//EXERCISE: Please add the second node to the graph.
89+
90+
// If we would simulate the graph-based metapopulation model now, we would just have two independent ODE-based SECIR-type models running with different initial conditions. In reality, there is usually exchange between regions through individuals travelling or commuting from one region to another. This can be realized via graph edges.
91+
// We here use a symmetric mobility i.e. we have the same number of individuals that travel from node 0 to node 1 as vice versa. We let 10% of the population commute via the edges twice a day.
92+
graph.add_edge(
93+
0, 1, Eigen::VectorX<ScalarType>::Constant((size_t)mio::osecir::InfectionState::Count * num_agegroups, 0.1));
94+
95+
// EXERCISE: Please add an edge from node 1 to node 0.
96+
97+
// Exchange commuters twice a day
98+
double dt_exchange = 0.5;
99+
100+
// *** Simulate model. ***
101+
// We now have finished initializing the metapopulation model. The graph-based simulation is created and advanced until `tmax`
102+
auto sim = mio::make_mobility_sim<ScalarType>(t0, dt_exchange, std::move(graph));
103+
104+
// EXERCISE: Please advance the simulation until `tmax`.
105+
106+
// As every graph node has its own model, we get one result time series per node. Those can be accessed as follows
107+
auto result_region0 = sim.get_graph().nodes()[0].property.get_result();
108+
auto result_region1 = sim.get_graph().nodes()[1].property.get_result();
109+
// Interpolate time series to full days.
110+
auto interpolated_result_r0 = mio::interpolate_simulation_result(result_region0);
111+
112+
// *** Print results. ***
113+
interpolated_result_r0.print_table();
114+
}

0 commit comments

Comments
 (0)