Skip to content

Commit f0191b6

Browse files
johapauHenrZujubicker
authored
905 add age group resolution to seir and sir model (#911)
- Integration of age groups into the seir and sir models - New examples using the age groups Co-authored-by: HenrZu <[email protected]> Co-authored-by: jubicker <[email protected]>
1 parent 3f602b1 commit f0191b6

24 files changed

Lines changed: 1353 additions & 406 deletions

File tree

‎cpp/examples/CMakeLists.txt‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,10 @@ add_executable(ode_seir_example ode_seir.cpp)
1818
target_link_libraries(ode_seir_example PRIVATE memilio ode_seir)
1919
target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
2020

21+
add_executable(ode_seir_ageres_example ode_seir_ageres.cpp)
22+
target_link_libraries(ode_seir_ageres_example PRIVATE memilio ode_seir)
23+
target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
24+
2125
add_executable(ode_sir_example ode_sir.cpp)
2226
target_link_libraries(ode_sir_example PRIVATE memilio ode_sir)
2327
target_compile_options(ode_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
@@ -30,6 +34,10 @@ add_executable(sde_sirs_example sde_sirs.cpp)
3034
target_link_libraries(sde_sirs_example PRIVATE memilio sde_sirs)
3135
target_compile_options(sde_sirs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
3236

37+
add_executable(ode_sir_ageres_example ode_sir_ageres.cpp)
38+
target_link_libraries(ode_sir_ageres_example PRIVATE memilio ode_sir)
39+
target_compile_options(ode_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
40+
3341
add_executable(seir_flows_example ode_seir_flows.cpp)
3442
target_link_libraries(seir_flows_example PRIVATE memilio ode_seir)
3543
target_compile_options(seir_flows_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/examples/graph.cpp‎

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -31,20 +31,30 @@ int main()
3131
const auto tmax = 10.;
3232
const auto dt = 0.5; //time step of migration, daily migration every second step
3333

34-
mio::oseir::Model model;
35-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Susceptible)}] = 10000;
34+
mio::oseir::Model model(1);
35+
36+
// set population
37+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000;
38+
39+
// set transition times
3640
model.parameters.set<mio::oseir::TimeExposed>(1);
37-
model.parameters.get<mio::oseir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
3841
model.parameters.set<mio::oseir::TimeInfected>(1);
3942

43+
// set contact matrix
44+
mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::oseir::ContactPatterns>().get_cont_freq_mat();
45+
contact_matrix[0].get_baseline().setConstant(2.7);
46+
4047
//two mostly identical groups
4148
auto model_group1 = model;
4249
auto model_group2 = model;
50+
4351
//some contact restrictions in group 1
44-
model_group1.parameters.get<mio::oseir::ContactPatterns>().add_damping(0.5, mio::SimulationTime(5));
52+
mio::ContactMatrixGroup& contact_matrix1 = model_group1.parameters.get<mio::oseir::ContactPatterns>().get_cont_freq_mat();
53+
contact_matrix1[0].add_damping(0.5, mio::SimulationTime(5));
54+
4555
//infection starts in group 1
46-
model_group1.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Susceptible)}] = 9990;
47-
model_group1.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}] = 10;
56+
model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9990;
57+
model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 10;
4858

4959
mio::Graph<mio::SimulationNode<mio::Simulation<mio::oseir::Model>>, mio::MigrationEdge> g;
5060
g.add_node(1001, model_group1, t0);

‎cpp/examples/ode_seir.cpp‎

Lines changed: 17 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,8 @@
2323
#include "memilio/compartments/simulation.h"
2424
#include "memilio/utils/logging.h"
2525

26+
#include "memilio/utils/time_series.h"
27+
2628
int main()
2729
{
2830
mio::set_log_level(mio::LogLevel::debug);
@@ -33,26 +35,25 @@ int main()
3335

3436
mio::log_info("Simulating ODE SEIR; t={} ... {} with dt = {}.", t0, tmax, dt);
3537

36-
mio::oseir::Model model;
38+
mio::oseir::Model model(1);
3739

38-
double total_population = 1061000;
39-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}] = 10000;
40-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}] = 1000;
41-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}] = 1000;
42-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Susceptible)}] =
43-
total_population -
44-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}] -
45-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}] -
46-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}];
40+
double total_population = 10000;
41+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100;
42+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100;
43+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100;
44+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] =
45+
total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] -
46+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] -
47+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}];
4748

48-
model.parameters.set<mio::oseir::TransmissionProbabilityOnContact>(0.04);
4949
model.parameters.set<mio::oseir::TimeExposed>(5.2);
50-
model.parameters.set<mio::oseir::TimeInfected>(2);
51-
52-
model.parameters.get<mio::oseir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
50+
model.parameters.set<mio::oseir::TimeInfected>(6);
51+
model.parameters.set<mio::oseir::TransmissionProbabilityOnContact>(0.1);
5352

54-
// contacts increase by 100% after 12.5 days
55-
model.parameters.get<mio::oseir::ContactPatterns>().add_damping(-1., mio::SimulationTime(12.5));
53+
mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::oseir::ContactPatterns>();
54+
contact_matrix[0].get_baseline().setConstant(2.7);
55+
contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.));
56+
5657

5758
model.check_constraints();
5859

‎cpp/examples/ode_seir_ageres.cpp‎

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
#include "ode_seir/model.h"
2+
#include "ode_seir/infection_state.h"
3+
#include "ode_seir/parameters.h"
4+
#include "memilio/compartments/simulation.h"
5+
#include "memilio/utils/logging.h"
6+
7+
int main()
8+
{
9+
mio::set_log_level(mio::LogLevel::debug);
10+
11+
double t0 = 0;
12+
double tmax = 50;
13+
double dt = 0.001;
14+
15+
mio::log_info("Simulating SEIR; t={} ... {} with dt = {}.", t0, tmax, dt);
16+
17+
double cont_freq = 10;
18+
19+
double nb_total_t0 = 10000, nb_exp_t0 = 100, nb_inf_t0 = 50, nb_rec_t0 = 10;
20+
const size_t num_groups = 3;
21+
22+
mio::oseir::Model model(num_groups);
23+
double fact = 1.0 / num_groups;
24+
25+
auto& params = model.parameters;
26+
27+
for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); i++) {
28+
model.populations[{i, mio::oseir::InfectionState::Exposed}] = fact * nb_exp_t0;
29+
model.populations[{i, mio::oseir::InfectionState::Infected}] = fact * nb_inf_t0;
30+
model.populations[{i, mio::oseir::InfectionState::Recovered}] = fact * nb_rec_t0;
31+
model.populations.set_difference_from_group_total<mio::AgeGroup>({i, mio::oseir::InfectionState::Susceptible},
32+
fact * nb_total_t0);
33+
34+
model.parameters.get<mio::oseir::TimeExposed>()[i] = 5.2;
35+
model.parameters.get<mio::oseir::TimeInfected>()[i] = 6;
36+
model.parameters.get<mio::oseir::TransmissionProbabilityOnContact>()[i] = 0.04;
37+
}
38+
39+
mio::ContactMatrixGroup& contact_matrix = params.get<mio::oseir::ContactPatterns>();
40+
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_groups, num_groups, fact * cont_freq));
41+
contact_matrix.add_damping(Eigen::MatrixXd::Constant(num_groups, num_groups, 0.7), mio::SimulationTime(30.));
42+
43+
model.apply_constraints();
44+
45+
mio::TimeSeries<double> seir = simulate(t0, tmax, dt, model);
46+
47+
std::vector<std::string> vars = {"S", "E", "I", "R"};
48+
printf("Number of time points :%d\n", static_cast<int>(seir.get_num_time_points()));
49+
printf("People in\n");
50+
51+
for (size_t k = 0; k < (size_t)mio::oseir::InfectionState::Count; k++) {
52+
double dummy = 0;
53+
54+
for (size_t i = 0; i < (size_t)params.get_num_groups(); i++) {
55+
printf("\t %s[%d]: %.0f", vars[k].c_str(), (int)i,
56+
seir.get_last_value()[k + (size_t)mio::oseir::InfectionState::Count * (int)i]);
57+
dummy += seir.get_last_value()[k + (size_t)mio::oseir::InfectionState::Count * (int)i];
58+
}
59+
60+
printf("\t %s_otal: %.0f\n", vars[k].c_str(), dummy);
61+
}
62+
}

‎cpp/examples/ode_seir_flows.cpp‎

Lines changed: 12 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -38,25 +38,24 @@ int main()
3838

3939
mio::log_info("Simulating SEIR; t={} ... {} with dt = {}.", t0, tmax, dt);
4040

41-
mio::oseir::Model model;
41+
mio::oseir::Model model(1);
42+
43+
constexpr double total_population = 10000;
44+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100;
45+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100;
46+
model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100;
47+
model.populations.set_difference_from_group_total<mio::AgeGroup>(
48+
{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}, total_population);
4249

43-
double total_population = 10000;
44-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}] = 100;
45-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}] = 100;
46-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}] = 100;
47-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Susceptible)}] =
48-
total_population -
49-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}] -
50-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}] -
51-
model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}];
52-
// suscetible now set with every other update
53-
// params.nb_sus_t0 = params.nb_total_t0 - params.nb_exp_t0 - params.nb_inf_t0 - params.nb_rec_t0;
5450
model.parameters.set<mio::oseir::TimeExposed>(5.2);
5551
model.parameters.set<mio::oseir::TimeInfected>(6);
5652
model.parameters.set<mio::oseir::TransmissionProbabilityOnContact>(0.04);
57-
model.parameters.get<mio::oseir::ContactPatterns>().get_baseline()(0, 0) = 10;
53+
54+
mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::oseir::ContactPatterns>().get_cont_freq_mat();
55+
contact_matrix[0].get_baseline().setConstant(10);
5856

5957
model.check_constraints();
58+
6059
auto seir = simulate_flows(t0, tmax, dt, model);
6160

6261
printf("Compartments: \n");

‎cpp/examples/ode_sir.cpp‎

Lines changed: 14 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,8 @@
2424
#include "ode_sir/infection_state.h"
2525
#include "ode_sir/model.h"
2626
#include "ode_sir/parameters.h"
27+
#include <fstream>
28+
#include <stdio.h>
2729

2830
int main()
2931
{
@@ -37,24 +39,23 @@ int main()
3739

3840
mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt);
3941

40-
mio::osir::Model model;
42+
mio::osir::Model model(1);
4143

42-
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Infected)}] = 1000;
43-
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Recovered)}] = 1000;
44-
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Susceptible)}] =
45-
total_population -
46-
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Infected)}] -
47-
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Recovered)}];
44+
model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}] = 1000;
45+
model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}] = 1000;
46+
model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Susceptible}] =
47+
total_population - model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}] -
48+
model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}];
4849
model.parameters.set<mio::osir::TimeInfected>(2);
49-
model.parameters.set<mio::osir::TransmissionProbabilityOnContact>(0.04);
50-
model.parameters.get<mio::osir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
51-
model.parameters.get<mio::osir::ContactPatterns>().add_damping(0.6, mio::SimulationTime(12.5));
52-
53-
auto integrator = std::make_shared<mio::EulerIntegratorCore>();
50+
model.parameters.set<mio::osir::TransmissionProbabilityOnContact>(0.5);
5451

52+
mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::osir::ContactPatterns>().get_cont_freq_mat();
53+
contact_matrix[0].get_baseline().setConstant(2.7);
54+
contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5));
5555
model.check_constraints();
5656

57-
auto sir = simulate(t0, tmax, dt, model, integrator);
57+
auto integrator = std::make_shared<mio::EulerIntegratorCore>();
58+
auto sir = simulate(t0, tmax, dt, model, integrator);
5859

5960
bool print_to_terminal = true;
6061

‎cpp/examples/ode_sir_ageres.cpp‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
#include "memilio/epidemiology/age_group.h"
2+
#include "memilio/utils/time_series.h"
3+
#include "memilio/utils/logging.h"
4+
#include "memilio/compartments/simulation.h"
5+
#include "ode_sir/model.h"
6+
7+
int main()
8+
{
9+
mio::set_log_level(mio::LogLevel::debug);
10+
11+
double t0 = 0;
12+
double tmax = 50;
13+
double dt = 0.1;
14+
15+
mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt);
16+
17+
double cont_freq = 10; // see Polymod study
18+
19+
double nb_total_t0 = 10000, nb_inf_t0 = 50, nb_rec_t0 = 10;
20+
21+
const size_t num_groups = 3;
22+
mio::osir::Model model(num_groups);
23+
24+
double fact = 1.0 / num_groups;
25+
26+
auto& params = model.parameters;
27+
28+
for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); i++) {
29+
model.populations[{i, mio::osir::InfectionState::Infected}] = fact * nb_inf_t0;
30+
model.populations[{i, mio::osir::InfectionState::Recovered}] = fact * nb_rec_t0;
31+
model.populations.set_difference_from_group_total<mio::AgeGroup>({i, mio::osir::InfectionState::Susceptible},
32+
fact * nb_total_t0);
33+
34+
model.parameters.get<mio::osir::TimeInfected>()[i] = 2.0;
35+
model.parameters.get<mio::osir::TransmissionProbabilityOnContact>()[i] = 0.3;
36+
}
37+
38+
mio::ContactMatrixGroup& contact_matrix = params.get<mio::osir::ContactPatterns>();
39+
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_groups, num_groups, fact * cont_freq));
40+
contact_matrix.add_damping(Eigen::MatrixXd::Constant(num_groups, num_groups, 0.7), mio::SimulationTime(30.));
41+
42+
model.apply_constraints();
43+
44+
mio::TimeSeries<double> sir = simulate(t0, tmax, dt, model);
45+
46+
std::vector<std::string> vars = {"S", "I", "R"};
47+
printf("Number of time points :%d\n", static_cast<int>(sir.get_num_time_points()));
48+
printf("People in\n");
49+
50+
for (size_t k = 0; k < (size_t)mio::osir::InfectionState::Count; k++) {
51+
double dummy = 0;
52+
53+
for (size_t i = 0; i < (size_t)params.get_num_groups(); i++) {
54+
printf("\t %s[%d]: %.0f", vars[k].c_str(), (int)i,
55+
sir.get_last_value()[k + (size_t)mio::osir::InfectionState::Count * (int)i]);
56+
dummy += sir.get_last_value()[k + (size_t)mio::osir::InfectionState::Count * (int)i];
57+
}
58+
59+
printf("\t %s_otal: %.0f\n", vars[k].c_str(), dummy);
60+
}
61+
}

0 commit comments

Comments
 (0)