Skip to content

Commit c78bafd

Browse files
hatritannawendler
hatrit
and
annawendler
authored
379 Implement age resolution for IDE model (#1083)
- Added the possibility to use age groups in the IDE-SECIR model. - Parameters and member functions of the model were adapted accordingly. - Tests were added. Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com>
1 parent c12ef2e commit c78bafd

16 files changed

+1458
-789
lines changed

‎cpp/examples/CMakeLists.txt‎

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

118+
add_executable(ide_secir_ageres_example ide_secir_ageres.cpp)
119+
target_link_libraries(ide_secir_ageres_example PRIVATE memilio ide_secir)
120+
target_compile_options(ide_secir_ageres_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
121+
118122
add_executable(lct_secir_example lct_secir.cpp)
119123
target_link_libraries(lct_secir_example PRIVATE memilio lct_secir)
120124
target_compile_options(lct_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/examples/ide_initialization.cpp‎

Lines changed: 13 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626
#include "memilio/utils/time_series.h"
2727
#include "memilio/utils/date.h"
2828
#include "memilio/math/eigen.h"
29+
#include <cstddef>
2930
#include <string>
3031
#include <vector>
3132
#include <iostream>
@@ -60,21 +61,26 @@ int main(int argc, char** argv)
6061
// The default parameters of the IDE-SECIR model are used, so that the simulation results are not realistic and are for demonstration purpose only.
6162

6263
// Initialize model.
63-
ScalarType total_population = 80 * 1e6;
64-
ScalarType deaths = 0; // The number of deaths will be overwritten if real data is used for initialization.
65-
ScalarType dt = 0.5;
64+
size_t num_agegroups = 1;
65+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> total_population =
66+
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 80 * 1e6);
67+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths = mio::CustomIndexArray<ScalarType, mio::AgeGroup>(
68+
mio::AgeGroup(num_agegroups),
69+
0.); // The number of deaths will be overwritten if real data is used for initialization.
70+
ScalarType dt = 0.5;
6671
mio::isecir::Model model(mio::TimeSeries<ScalarType>((int)mio::isecir::InfectionTransition::Count),
67-
total_population, deaths);
72+
total_population, deaths, num_agegroups);
6873

6974
// Check provided parameters.
7075
std::string filename = setup(argc, argv);
7176
if (filename.empty()) {
7277
std::cout << "You did not provide a valid filename. A default initialization is used." << std::endl;
7378

7479
using Vec = mio::TimeSeries<ScalarType>::Vector;
75-
mio::TimeSeries<ScalarType> init((int)mio::isecir::InfectionTransition::Count);
76-
init.add_time_point<Eigen::VectorXd>(-7., Vec::Constant((int)mio::isecir::InfectionTransition::Count, 1. * dt));
77-
while (init.get_last_time() < -dt / 2) {
80+
mio::TimeSeries<ScalarType> init(num_agegroups * (size_t)mio::isecir::InfectionTransition::Count);
81+
init.add_time_point<Eigen::VectorXd>(-7.,
82+
Vec::Constant((size_t)mio::isecir::InfectionTransition::Count, 1. * dt));
83+
while (init.get_last_time() < -dt / 2.) {
7884
init.add_time_point(init.get_last_time() + dt,
7985
Vec::Constant((int)mio::isecir::InfectionTransition::Count, 1. * dt));
8086
}

‎cpp/examples/ide_secir.cpp‎

Lines changed: 35 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,9 @@
2222
#include "ide_secir/infection_state.h"
2323
#include "ide_secir/simulation.h"
2424
#include "memilio/config.h"
25+
#include "memilio/epidemiology/age_group.h"
2526
#include "memilio/math/eigen.h"
27+
#include "memilio/utils/custom_index_array.h"
2628
#include "memilio/utils/time_series.h"
2729
#include "memilio/epidemiology/uncertain_matrix.h"
2830
#include "memilio/epidemiology/state_age_function.h"
@@ -32,18 +34,23 @@ int main()
3234
{
3335
using Vec = mio::TimeSeries<ScalarType>::Vector;
3436

35-
ScalarType tmax = 10;
36-
ScalarType N = 10000;
37-
ScalarType deaths = 13.10462213;
38-
ScalarType dt = 1e-2;
37+
size_t num_agegroups = 1;
38+
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.;
3945

4046
int num_transitions = (int)mio::isecir::InfectionTransition::Count;
4147

42-
// Create TimeSeries with num_transitions elements where transitions needed for simulation will be stored.
43-
mio::TimeSeries<ScalarType> init(num_transitions);
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);
4451

4552
// Add time points for initialization of transitions.
46-
Vec vec_init(num_transitions);
53+
Vec vec_init(num_transitions * num_agegroups);
4754
vec_init[(int)mio::isecir::InfectionTransition::SusceptibleToExposed] = 25.0;
4855
vec_init[(int)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 15.0;
4956
vec_init[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 8.0;
@@ -54,7 +61,8 @@ int main()
5461
vec_init[(int)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = 1.0;
5562
vec_init[(int)mio::isecir::InfectionTransition::InfectedCriticalToDead] = 1.0;
5663
vec_init[(int)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = 1.0;
57-
vec_init = vec_init * dt;
64+
65+
vec_init = vec_init * dt;
5866
// Add initial time point to time series.
5967
init.add_time_point(-10, vec_init);
6068
// Add further time points until time 0.
@@ -63,11 +71,13 @@ int main()
6371
}
6472

6573
// Initialize model.
66-
mio::isecir::Model model(std::move(init), N, deaths);
74+
mio::isecir::Model model(std::move(init), N, deaths, num_agegroups);
6775

68-
// Uncomment one of these lines to use a different method to initialize the model using the TimeSeries init.
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.
6978
// model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 1000;
70-
// model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0;
79+
// Initialization method with Recovered.
80+
// model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0;
7181

7282
// Set working parameters.
7383
mio::SmootherCosine smoothcos(2.0);
@@ -77,23 +87,26 @@ int main()
7787
vec_delaydistrib[(int)mio::isecir::InfectionTransition::SusceptibleToExposed].set_distribution_parameter(-1.);
7888
vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]
7989
.set_distribution_parameter(4.0);
80-
model.parameters.set<mio::isecir::TransitionDistributions>(vec_delaydistrib);
90+
91+
model.parameters.get<mio::isecir::TransitionDistributions>()[mio::AgeGroup(0)] = vec_delaydistrib;
8192

8293
std::vector<ScalarType> vec_prob(num_transitions, 0.5);
8394
// The following probabilities must be 1, as there is no other way to go.
8495
vec_prob[Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)] = 1;
8596
vec_prob[Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)] = 1;
86-
model.parameters.set<mio::isecir::TransitionProbabilities>(vec_prob);
97+
model.parameters.get<mio::isecir::TransitionProbabilities>()[mio::AgeGroup(0)] = vec_prob;
8798

88-
mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, 1);
89-
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10.));
99+
mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, num_agegroups);
100+
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, 10.));
90101
model.parameters.get<mio::isecir::ContactPatterns>() = mio::UncertainContactMatrix(contact_matrix);
91102

92103
mio::ExponentialSurvivalFunction exponential(0.5);
93104
mio::StateAgeFunctionWrapper prob(exponential);
94-
model.parameters.set<mio::isecir::TransmissionProbabilityOnContact>(prob);
95-
model.parameters.set<mio::isecir::RelativeTransmissionNoSymptoms>(prob);
96-
model.parameters.set<mio::isecir::RiskOfInfectionFromSymptomatic>(prob);
105+
106+
model.parameters.get<mio::isecir::TransmissionProbabilityOnContact>()[mio::AgeGroup(0)] = prob;
107+
model.parameters.get<mio::isecir::RelativeTransmissionNoSymptoms>()[mio::AgeGroup(0)] = prob;
108+
model.parameters.get<mio::isecir::RiskOfInfectionFromSymptomatic>()[mio::AgeGroup(0)] = prob;
109+
97110
model.parameters.set<mio::isecir::Seasonality>(0.1);
98111
// Start the simulation on the 40th day of a year (i.e. in February).
99112
model.parameters.set<mio::isecir::StartDay>(40);
@@ -104,8 +117,10 @@ int main()
104117
mio::isecir::Simulation sim(model, dt);
105118
sim.advance(tmax);
106119

107-
auto interpolated_results = mio::interpolate_simulation_result(sim.get_result(), dt / 2);
120+
auto interpolated_results = mio::interpolate_simulation_result(sim.get_result(), dt / 2.);
121+
108122
interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
109123
// Uncomment this line to print the transitions.
110-
// 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"}, 16, 8);
124+
// sim.get_transitions().print_table(
125+
// {"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);
111126
}

‎cpp/examples/ide_secir_ageres.cpp‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
2+
#include "ide_secir/model.h"
3+
#include "ide_secir/infection_state.h"
4+
#include "ide_secir/simulation.h"
5+
#include "memilio/config.h"
6+
#include "memilio/epidemiology/age_group.h"
7+
#include "memilio/math/eigen.h"
8+
#include "memilio/utils/custom_index_array.h"
9+
#include "memilio/utils/time_series.h"
10+
#include "memilio/epidemiology/uncertain_matrix.h"
11+
#include "memilio/epidemiology/state_age_function.h"
12+
#include "memilio/data/analyze_result.h"
13+
14+
int main()
15+
{
16+
using Vec = mio::TimeSeries<ScalarType>::Vector;
17+
18+
size_t num_agegroups = 2;
19+
20+
ScalarType tmax = 5;
21+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> N =
22+
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 10000.);
23+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths =
24+
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 6.);
25+
ScalarType dt = 1.;
26+
27+
int num_transitions = (int)mio::isecir::InfectionTransition::Count;
28+
29+
// Create TimeSeries with num_transitions * num_agegroups elements where transitions needed for simulation will be
30+
// stored.
31+
mio::TimeSeries<ScalarType> init(num_transitions * num_agegroups);
32+
33+
// Add time points for initialization of transitions.
34+
Vec vec_init(num_transitions * num_agegroups);
35+
// Values for the Infectiontransitions are the same for all AgeGroups.
36+
for (size_t group = 0; group < num_agegroups; ++group) {
37+
vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::SusceptibleToExposed] = 25.0;
38+
vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 15.0;
39+
vec_init[group * num_transitions +
40+
(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 8.0;
41+
vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered] = 4.0;
42+
vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] =
43+
1.0;
44+
vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered] = 4.0;
45+
vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] =
46+
1.0;
47+
vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = 1.0;
48+
vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedCriticalToDead] = 1.0;
49+
vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = 1.0;
50+
}
51+
52+
// Add initial time point to time series.
53+
init.add_time_point(-10, vec_init);
54+
// Add further time points until time 0.
55+
while (init.get_last_time() < -dt / 2) {
56+
init.add_time_point(init.get_last_time() + dt, vec_init);
57+
}
58+
59+
// Initialize model.
60+
mio::isecir::Model model(std::move(init), N, deaths, num_agegroups);
61+
62+
// Uncomment these lines to use a different method to initialize the model using the TimeSeries init.
63+
// Initialization method with Susceptibles.
64+
// model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 1000;
65+
// model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Count +
66+
// (Eigen::Index)mio::isecir::InfectionState::Susceptible] = 1000;
67+
// Initialization method with Recovered.
68+
// model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0;
69+
// model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Count +
70+
// (Eigen::Index)mio::isecir::InfectionState::Recovered] = 0;
71+
72+
// Set working parameters.
73+
// First AgeGroup for Transition Distributions.
74+
mio::SmootherCosine smoothcos1(2.0);
75+
mio::StateAgeFunctionWrapper delaydistribution1(smoothcos1);
76+
std::vector<mio::StateAgeFunctionWrapper> vec_delaydistrib1(num_transitions, delaydistribution1);
77+
// TransitionDistribution is not used for SusceptibleToExposed. Therefore, the parameter can be set to any value.
78+
vec_delaydistrib1[(int)mio::isecir::InfectionTransition::SusceptibleToExposed].set_distribution_parameter(-1.);
79+
80+
model.parameters.get<mio::isecir::TransitionDistributions>()[mio::AgeGroup(0)] = vec_delaydistrib1;
81+
82+
//Second AgeGroup for Transition Distributions.
83+
mio::SmootherCosine smoothcos2(3.0);
84+
mio::StateAgeFunctionWrapper delaydistribution2(smoothcos2);
85+
std::vector<mio::StateAgeFunctionWrapper> vec_delaydistrib2(num_transitions, delaydistribution2);
86+
// TransitionDistribution is not used for SusceptibleToExposed. Therefore, the parameter can be set to any value.
87+
vec_delaydistrib2[(int)mio::isecir::InfectionTransition::SusceptibleToExposed].set_distribution_parameter(-1.);
88+
89+
model.parameters.get<mio::isecir::TransitionDistributions>()[mio::AgeGroup(1)] = vec_delaydistrib2;
90+
91+
std::vector<ScalarType> vec_prob(num_transitions, 0.5);
92+
// The following probabilities must be 1, as there is no other way to go.
93+
vec_prob[Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)] = 1;
94+
vec_prob[Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)] = 1;
95+
for (mio::AgeGroup group = mio::AgeGroup(0); group < mio::AgeGroup(num_agegroups); ++group) {
96+
model.parameters.get<mio::isecir::TransitionProbabilities>()[group] = vec_prob;
97+
}
98+
99+
mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, static_cast<Eigen::Index>(num_agegroups));
100+
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, 10.));
101+
model.parameters.get<mio::isecir::ContactPatterns>() = mio::UncertainContactMatrix(contact_matrix);
102+
103+
mio::ExponentialSurvivalFunction exponential(0.5);
104+
mio::StateAgeFunctionWrapper prob(exponential);
105+
for (mio::AgeGroup group = mio::AgeGroup(0); group < mio::AgeGroup(num_agegroups); ++group) {
106+
model.parameters.get<mio::isecir::TransmissionProbabilityOnContact>()[group] = prob;
107+
model.parameters.get<mio::isecir::RelativeTransmissionNoSymptoms>()[group] = prob;
108+
model.parameters.get<mio::isecir::RiskOfInfectionFromSymptomatic>()[group] = prob;
109+
}
110+
model.parameters.set<mio::isecir::Seasonality>(0.1);
111+
// Start the simulation on the 40th day of a year (i.e. in February).
112+
model.parameters.set<mio::isecir::StartDay>(40);
113+
114+
model.check_constraints(dt);
115+
116+
// Carry out simulation.
117+
mio::isecir::Simulation sim(model, dt);
118+
sim.advance(tmax);
119+
120+
auto interpolated_results = mio::interpolate_simulation_result(sim.get_result(), dt / 2.);
121+
122+
interpolated_results.print_table(
123+
{"S1", "E1", "C1", "I1", "H1", "U1", "R1", "D1 ", "S2", "E2", "C2", "I2", "H2", "U2", "R2", "D2 "}, 16, 8);
124+
// Uncomment this line to print the transitions.
125+
// 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",
126+
// "H->R 1", "U->D 1", "U->R 1", "S->E 2", "E->C 2", "C->I 2", "C->R 2",
127+
// "I->H 2", "I->R 2", "H->U 2", "H->R 2", "U->D 2", "U->R 2"},
128+
// 16, 8);
129+
}

‎cpp/memilio/epidemiology/state_age_function.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222
#define STATEAGEFUNCTION_H
2323

2424
#include "memilio/config.h"
25+
#include "memilio/utils/compiler_diagnostics.h"
2526
#include "memilio/utils/parameter_set.h"
2627
#include "memilio/math/smoother.h"
2728
#include "memilio/math/floating_point.h"
@@ -811,6 +812,10 @@ struct ErlangDensity : public StateAgeFunction {
811812
*/
812813
struct StateAgeFunctionWrapper {
813814

815+
StateAgeFunctionWrapper()
816+
: m_function(mio::SmootherCosine(1.0).clone())
817+
{
818+
}
814819
/**
815820
* @brief Constructs a new StateAgeFunctionWrapper object
816821
*

0 commit comments

Comments
 (0)