Skip to content

Commit ce75c58

Browse files
lenaploetzkeannawendlermknaranja
authored
958 implement a check for the results of the initialization (#959)
- Optimization of IDE SECIR initialization to avoid overwriting of beforehand set values. - Check and return of initialization could be conducted properly. - Added unit tests to check_constraints() function to increase code coverage. Co-authored-by: annawendler <[email protected]> Co-authored-by: Martin J. Kühn <[email protected]>
1 parent 04179df commit ce75c58

4 files changed

Lines changed: 199 additions & 86 deletions

File tree

‎cpp/README.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ The MEmilio C++ library contains the implementation of the epidemiological model
44

55
Directory structure:
66
- memilio: framework for developing epidemiological models with, e.g., interregional mobility implementations, nonpharmaceutical interventions (NPIs), and mathematical, programming, and IO utilities.
7-
- models: implementation of concrete models (ODE and ABM)
7+
- models: implementation of concrete models (ODE, IDE, LCT and ABM)
88
- simulations: simulation applications that were used to generate the scenarios and data for publications
99
- examples: small applications that help with using the framework and models
1010
- tests: unit tests for framework and models.

‎cpp/models/ide_secir/model.cpp‎

Lines changed: 53 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -66,14 +66,51 @@ void Model::initialize(ScalarType dt)
6666
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] -
6767
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
6868
}
69-
else {
69+
else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] > 1e-12) {
70+
// Take initialized value for Susceptibles if value can't be calculated via the standard formula.
71+
m_initialization_method = 2;
72+
// Calculate other compartment sizes for t=0.
73+
other_compartments_current_timestep(dt);
7074

75+
// R; need an initial value for R, therefore do not calculate via compute_recovered()
76+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] =
77+
m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] -
78+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
79+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
80+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
81+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
82+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
83+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
84+
}
85+
else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] > 1e-12) {
86+
// If value for Recovered is initialized and standard method is not applicable (i.e., no value for total infections
87+
// or suceptibles given directly), calculate Susceptibles via other compartments.
88+
// The calculation of other compartments' values is not dependent on Susceptibles at time 0, i.e., S(0), but only on the transitions of the past.
89+
m_initialization_method = 3;
90+
other_compartments_current_timestep(dt);
91+
92+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
93+
m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
94+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
95+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
96+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
97+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
98+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] -
99+
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
100+
}
101+
else {
71102
// compute Susceptibles at time 0 and m_forceofinfection at time -dt as initial values for discretization scheme
72103
// use m_forceofinfection at -dt to be consistent with further calculations of S (see compute_susceptibles()),
73104
// where also the value of m_forceofinfection for the previous timestep is used
74105
update_forceofinfection(dt, true);
75106
if (m_forceofinfection > 1e-12) {
76-
m_initialization_method = 2;
107+
m_initialization_method = 4;
108+
109+
/* Attention: With an inappropriate combination of parameters and initial conditions, it can happen that S
110+
becomes greater than N when this formula is used. In this case, at least one compartment is initialized
111+
with a number less than zero, so that a log_error is thrown.
112+
However, this initialization method is consistent with the numerical solver of the model equations,
113+
so it may sometimes make sense to use this method. */
77114
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
78115
m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] /
79116
(dt * m_forceofinfection);
@@ -91,46 +128,27 @@ void Model::initialize(ScalarType dt)
91128
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
92129
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
93130
}
94-
else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] > 1e-12) {
95-
//take initialized value for Susceptibles if value can't be calculated via the standard formula
96-
m_initialization_method = 3;
97-
//calculate other compartment sizes for t=0
98-
other_compartments_current_timestep(dt);
99-
100-
//R; need an initial value for R, therefore do not calculate via compute_recovered()
101-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] =
102-
m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] -
103-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
104-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
105-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
106-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
107-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
108-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
109-
}
110-
else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] > 1e-12) {
111-
//if value for Recovered is initialized and standard method is not applicable, calculate Susceptibles via other compartments
112-
//determining other compartment sizes is not dependent of Susceptibles(0), just of the transitions of the past.
113-
//calculate other compartment sizes for t=0
114-
m_initialization_method = 4;
115-
other_compartments_current_timestep(dt);
116-
117-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
118-
m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
119-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
120-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
121-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
122-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
123-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] -
124-
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
125-
}
126131
else {
127132
m_initialization_method = -1;
128133
log_error("Error occured while initializing compartments: Force of infection is evaluated to 0 and neither "
129134
"Susceptibles nor Recovered or total confirmed cases for time 0 were set. One of them should be "
130135
"larger 0.");
131136
}
132137
}
133-
// compute m_forceofinfection at time 0 needed for further simulation
138+
139+
// Verify that the initialization produces an appropriate result.
140+
// Another check would be if the sum of the compartments is equal to N, but in all initialization methods one of the compartments is initialized via N - the others.
141+
// This also means that if a compartment is greater than N, we will always have one or more compartments less than zero.
142+
// Check if all compartments are non negative.
143+
for (int i = 0; i < (int)InfectionState::Count; i++) {
144+
if (m_populations[0][i] < 0) {
145+
m_initialization_method = -2;
146+
log_error("Initialization failed. One or more initial values for populations are less than zero. It may "
147+
"help to use a different method for initialization.");
148+
}
149+
}
150+
151+
// Compute m_forceofinfection at time 0 needed for further simulation.
134152
update_forceofinfection(dt);
135153
}
136154

‎cpp/models/ide_secir/model.h‎

Lines changed: 27 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -53,23 +53,15 @@ class Model
5353
const ParameterSet& Parameterset_init = ParameterSet());
5454

5555
/**
56-
* @brief Checks constraints on model parameters.
56+
* @brief Checks constraints on model parameters and initial data.
57+
* @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false.
5758
*/
58-
void check_constraints(ScalarType dt) const
59+
bool check_constraints(ScalarType dt) const
5960
{
60-
if (!(m_populations.get_num_time_points() > 0)) {
61-
log_error("Model construction failed. No initial time point for populations.");
62-
}
63-
64-
for (int i = 0; i < (int)InfectionState::Count; i++) {
65-
if (m_populations[0][i] < 0) {
66-
log_error("Initialization failed. Initial values for populations are less than zero.");
67-
}
68-
}
69-
7061
if (!((int)m_transitions.get_num_elements() == (int)InfectionTransition::Count)) {
71-
log_error(
72-
"Initialization failed. Number of elements in transition vector does not match the required number.");
62+
log_error("A variable given for model construction is not valid. Number of elements in transition vector "
63+
"does not match the required number.");
64+
return true;
7365
}
7466

7567
ScalarType support_max = std::max(
@@ -95,17 +87,25 @@ class Model
9587
if (m_transitions.get_num_time_points() < (Eigen::Index)std::ceil(support_max / dt)) {
9688
log_error(
9789
"Initialization failed. Not enough time points for transitions given before start of simulation.");
90+
return true;
9891
}
9992

100-
parameters.check_constraints();
93+
return parameters.check_constraints();
10194
}
10295

10396
/**
104-
* @brief Calculate the number of individuals in each compartment for time 0.
105-
*
106-
* Initial transitions are used to calculate the initial compartment sizes.
107-
* @param[in] dt Time discretization step size.
108-
*/
97+
* @brief Initializes the model and sets compartment values at time zero.
98+
*
99+
* The initialization method is selected automatically based on the different values that need to be set beforehand.
100+
* Infection compartments are always computed through historic flow.
101+
* Initialization methods for susceptibles and recovered are tested in the following order:
102+
* 1.) If a positive number for the total number of confirmed cases is set, recovered is set according to that value and #Susceptible%s are derived.
103+
* 2.) If #Susceptible%s are set, recovered will be derived.
104+
* 3.) If #Recovered are set directly, #Susceptible%s are derived.
105+
* 4.) If none of the above is set with positive value, the force of infection is used as in Messina et al (2021) to set the #Susceptible%s.
106+
*
107+
* @param[in] dt Time discretization step size.
108+
*/
109109
void initialize(ScalarType dt);
110110

111111
/**
@@ -208,14 +208,15 @@ class Model
208208
}
209209

210210
/**
211-
* @brief Specifies a number associated with the method used for initialization.
211+
* @brief Returns the number associated with the method selected automatically for initialization.
212212
*
213-
* @returns 0 if the initialization method has not yet been selected,
213+
* @returns 0 if the model has not yet been initialized and the method has not yet been selected,
214214
* 1 if the method using the total number of confirmed cases at time 0 is used,
215-
* 2 if the force of infection method is used,
216-
* 3 if the initialization is calculated using a prior set value for S,
217-
* 4 if the initialization is calculated using a prior set value for R and
218-
* -1 if the initialization was not possible using any of the methods.
215+
* 2 if the initialization is calculated using a beforehand set value for S,
216+
* 3 if the initialization is calculated using a beforehand set value for R,
217+
* 4 if the force of infection method is used,
218+
* -1 if the initialization was not possible using any of the methods and
219+
* -2 if the initialization was possible using one of the provided methods but the result is not appropriate.
219220
*/
220221
int get_initialization_method()
221222
{

0 commit comments

Comments
 (0)