Skip to content

Commit fc0ebaf

Browse files
lenaploetzkeannawendler
lenaploetzke
and
annawendler
authored
946 Implement model based on Generalized Linear Chain Trick (#1058)
- Added a GLCT-SECIR model (generalized linear chain trick model) that allows to use phase-type distributed stay times in the compartments. - Added tests. - Moved calculate_compartments function from the LCT model to the lct_infection_state such that the GLCT model can use the function as well. Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com>
1 parent bc417ce commit fc0ebaf

File tree

17 files changed

+2024
-61
lines changed

17 files changed

+2024
-61
lines changed

‎cpp/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -146,6 +146,7 @@ if(MEMILIO_BUILD_MODELS)
146146
add_subdirectory(models/ode_secir)
147147
add_subdirectory(models/ode_secirvvs)
148148
add_subdirectory(models/lct_secir)
149+
add_subdirectory(models/glct_secir)
149150
add_subdirectory(models/ide_secir)
150151
add_subdirectory(models/ide_seir)
151152
add_subdirectory(models/ode_seir)

‎cpp/examples/CMakeLists.txt‎

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

118+
add_executable(glct_secir_example glct_secir.cpp)
119+
target_link_libraries(glct_secir_example PRIVATE memilio glct_secir)
120+
target_compile_options(glct_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
121+
118122
add_executable(history_example history.cpp)
119123
target_link_libraries(history_example PRIVATE memilio)
120124
target_compile_options(history_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/examples/glct_secir.cpp‎

Lines changed: 211 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,211 @@
1+
/*
2+
* Copyright (C) 2020-2024 MEmilio
3+
*
4+
* Authors: Lena Ploetzke
5+
*
6+
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
7+
*
8+
* Licensed under the Apache License, Version 2.0 (the "License");
9+
* you may not use this file except in compliance with the License.
10+
* You may obtain a copy of the License at
11+
*
12+
* http://www.apache.org/licenses/LICENSE-2.0
13+
*
14+
* Unless required by applicable law or agreed to in writing, software
15+
* distributed under the License is distributed on an "AS IS" BASIS,
16+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17+
* See the License for the specific language governing permissions and
18+
* limitations under the License.
19+
*/
20+
21+
#include "glct_secir/model.h"
22+
#include "glct_secir/parameters.h"
23+
#include "memilio/config.h"
24+
#include "memilio/utils/time_series.h"
25+
#include "memilio/epidemiology/uncertain_matrix.h"
26+
#include "memilio/math/eigen.h"
27+
#include "memilio/utils/logging.h"
28+
#include "memilio/compartments/simulation.h"
29+
#include "memilio/data/analyze_result.h"
30+
31+
#include <vector>
32+
33+
int main()
34+
{
35+
// Simple example to demonstrate how to run a simulation using an GLCT SECIR model.
36+
// This example recreates the lct_secir.cpp example using the GLCT model.
37+
// This means, that we use the corresponding initial numbers, parameters and numbers of subcompartments that are
38+
// equivalent to the choices in the LCT example.
39+
// Parameters, initial values and the number of subcompartments are not meant to represent a realistic scenario.
40+
// We need to double the choices of the number of subcompartments for some compartments,
41+
// as we define different strains for different transition possibilities. For both strains, the same number of
42+
// subcompartments is chosen. The transition probabilities are defined in the StartingProbabilities.
43+
constexpr size_t NumExposed = 2, NumInfectedNoSymptoms = 6, NumInfectedSymptoms = 2, NumInfectedSevere = 2,
44+
NumInfectedCritical = 10;
45+
using Model = mio::glsecir::Model<NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms, NumInfectedSevere,
46+
NumInfectedCritical>;
47+
using LctState = Model::LctState;
48+
using InfectionState = LctState::InfectionState;
49+
50+
Model model;
51+
52+
const ScalarType tmax = 10;
53+
const ScalarType t0 = 0;
54+
const ScalarType dt_init = 10; ///< Initially used step size for adaptive method.
55+
56+
// Define some epidemiological parameters needed throughput the model definition and initialization.
57+
const ScalarType timeExposed = 3.2;
58+
const ScalarType timeInfectedNoSymptoms = 2.;
59+
const ScalarType timeInfectedSymptoms = 5.8;
60+
const ScalarType timeInfectedSevere = 9.5;
61+
const ScalarType timeInfectedCritical = 7.1;
62+
const ScalarType recoveredPerInfectedNoSymptoms = 0.09;
63+
const ScalarType severePerInfectedSymptoms = 0.2;
64+
const ScalarType criticalPerSevere = 0.25;
65+
const ScalarType deathsPerCritical = 0.3;
66+
67+
// Define the initial values with the distribution of the population into subcompartments.
68+
// This method of defining the initial values using a vector of vectors is not necessary, but should remind you
69+
// how the entries of the initial value vector relate to the defined template parameters of the model or the number
70+
// of subcompartments. It is also possible to define the initial values directly.
71+
// Split the initial population defined in the lct_secir.cpp example according to the transition probabilities in
72+
// the two strains per compartment.
73+
std::vector<std::vector<ScalarType>> initial_populations = {
74+
{750},
75+
{30, 20},
76+
{20 * (1 - recoveredPerInfectedNoSymptoms), 10 * (1 - recoveredPerInfectedNoSymptoms),
77+
10 * (1 - recoveredPerInfectedNoSymptoms), 20 * recoveredPerInfectedNoSymptoms,
78+
10 * recoveredPerInfectedNoSymptoms, 10 * recoveredPerInfectedNoSymptoms},
79+
{50 * severePerInfectedSymptoms, 50 * (1 - severePerInfectedSymptoms)},
80+
{50 * criticalPerSevere, 50 * (1 - criticalPerSevere)},
81+
{10 * deathsPerCritical, 10 * deathsPerCritical, 5 * deathsPerCritical, 3 * deathsPerCritical,
82+
2 * deathsPerCritical, 10 * (1 - deathsPerCritical), 10 * (1 - deathsPerCritical), 5 * (1 - deathsPerCritical),
83+
3 * (1 - deathsPerCritical), 2 * (1 - deathsPerCritical)},
84+
{20},
85+
{10}};
86+
87+
// Assert that initial_populations has the right shape.
88+
if (initial_populations.size() != (size_t)InfectionState::Count) {
89+
mio::log_error("The number of vectors in initial_populations does not match the number of InfectionStates.");
90+
return 1;
91+
}
92+
if ((initial_populations[(size_t)InfectionState::Susceptible].size() !=
93+
LctState::get_num_subcompartments<InfectionState::Susceptible>()) ||
94+
(initial_populations[(size_t)InfectionState::Exposed].size() != NumExposed) ||
95+
(initial_populations[(size_t)InfectionState::InfectedNoSymptoms].size() != NumInfectedNoSymptoms) ||
96+
(initial_populations[(size_t)InfectionState::InfectedSymptoms].size() != NumInfectedSymptoms) ||
97+
(initial_populations[(size_t)InfectionState::InfectedSevere].size() != NumInfectedSevere) ||
98+
(initial_populations[(size_t)InfectionState::InfectedCritical].size() != NumInfectedCritical) ||
99+
(initial_populations[(size_t)InfectionState::Recovered].size() !=
100+
LctState::get_num_subcompartments<InfectionState::Recovered>()) ||
101+
(initial_populations[(size_t)InfectionState::Dead].size() !=
102+
LctState::get_num_subcompartments<InfectionState::Dead>())) {
103+
mio::log_error("The length of at least one vector in initial_populations does not match the related number of "
104+
"subcompartments.");
105+
return 1;
106+
}
107+
108+
// Transfer the initial values in initial_populations to the model.
109+
std::vector<ScalarType> flat_initial_populations;
110+
for (auto&& vec : initial_populations) {
111+
flat_initial_populations.insert(flat_initial_populations.end(), vec.begin(), vec.end());
112+
}
113+
for (size_t i = 0; i < LctState::Count; i++) {
114+
model.populations[mio::Index<LctState>(i)] = flat_initial_populations[i];
115+
}
116+
117+
// Set Parameters according to the LCT model definitions, e.g. use Erlang-distributions.
118+
// Exposed.
119+
// The get_default of the StartingProbabilities returns the first unit vector of the defined size.
120+
// It is necessary to set it although the default method is used to define the length of the vector.
121+
model.parameters.get<mio::glsecir::StartingProbabilitiesExposed>() =
122+
mio::glsecir::StartingProbabilitiesExposed().get_default(
123+
LctState::get_num_subcompartments<InfectionState::Exposed>());
124+
// The get_default function returns the TransitionMatrix that is required to have an Erlang-distributed
125+
// stay time with an average of timeExposed.
126+
model.parameters.get<mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms>() =
127+
mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms().get_default(
128+
LctState::get_num_subcompartments<InfectionState::Exposed>(), timeExposed);
129+
// This definition of the StartingProbability and the TransitionMatrix lead to an Erlang-distributed latent stage.
130+
131+
// InfectedNoSymptoms.
132+
// For InfectedNoSymptoms, two strains has to be defined, one for the Transition
133+
// InfectedNoSymptomsToInfectedSymptoms and one for InfectedNoSymptomsToRecovered.
134+
// The strains have a length of NumInfectedNoSymptoms/2. each.
135+
// The transition probability is included in the StartingProbability vector.
136+
mio::Vector<ScalarType> StartingProbabilitiesInfectedNoSymptoms =
137+
mio::Vector<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>());
138+
StartingProbabilitiesInfectedNoSymptoms[0] = 1 - recoveredPerInfectedNoSymptoms;
139+
StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)(
140+
LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>() / 2.)] = recoveredPerInfectedNoSymptoms;
141+
model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedNoSymptoms>() =
142+
StartingProbabilitiesInfectedNoSymptoms;
143+
// Define equal TransitionMatrices for the strains.
144+
// They follow the same Erlang-distribution such that we get the same result as with one strain in the LCT model.
145+
model.parameters.get<mio::glsecir::TransitionMatrixInfectedNoSymptomsToInfectedSymptoms>() =
146+
mio::glsecir::TransitionMatrixInfectedNoSymptomsToInfectedSymptoms().get_default(
147+
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>() / 2.),
148+
timeInfectedNoSymptoms);
149+
model.parameters.get<mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered>() =
150+
mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered().get_default(
151+
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>() / 2.),
152+
timeInfectedNoSymptoms);
153+
// Do the same for all compartments.
154+
// InfectedSymptoms.
155+
mio::Vector<ScalarType> StartingProbabilitiesInfectedSymptoms =
156+
mio::Vector<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>());
157+
StartingProbabilitiesInfectedSymptoms[0] = severePerInfectedSymptoms;
158+
StartingProbabilitiesInfectedSymptoms[(Eigen::Index)(
159+
LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>() / 2.)] = 1 - severePerInfectedSymptoms;
160+
model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedSymptoms>() = StartingProbabilitiesInfectedSymptoms;
161+
model.parameters.get<mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere>() =
162+
mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere().get_default(
163+
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>() / 2.), timeInfectedSymptoms);
164+
model.parameters.get<mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered>() =
165+
mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default(
166+
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>() / 2.), timeInfectedSymptoms);
167+
// InfectedSevere.
168+
mio::Vector<ScalarType> StartingProbabilitiesInfectedSevere =
169+
mio::Vector<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedSevere>());
170+
StartingProbabilitiesInfectedSevere[0] = criticalPerSevere;
171+
StartingProbabilitiesInfectedSevere[(Eigen::Index)(
172+
LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 2.)] = 1 - criticalPerSevere;
173+
model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedSevere>() = StartingProbabilitiesInfectedSevere;
174+
model.parameters.get<mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical>() =
175+
mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical().get_default(
176+
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 2.), timeInfectedSevere);
177+
model.parameters.get<mio::glsecir::TransitionMatrixInfectedSevereToRecovered>() =
178+
mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default(
179+
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 2.), timeInfectedSevere);
180+
// InfectedCritical.
181+
mio::Vector<ScalarType> StartingProbabilitiesInfectedCritical =
182+
mio::Vector<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedCritical>());
183+
StartingProbabilitiesInfectedCritical[0] = deathsPerCritical;
184+
StartingProbabilitiesInfectedCritical[(Eigen::Index)(
185+
LctState::get_num_subcompartments<InfectionState::InfectedCritical>() / 2.)] = 1 - deathsPerCritical;
186+
model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedCritical>() = StartingProbabilitiesInfectedCritical;
187+
model.parameters.get<mio::glsecir::TransitionMatrixInfectedCriticalToDead>() =
188+
mio::glsecir::TransitionMatrixInfectedCriticalToDead().get_default(
189+
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedCritical>() / 2.), timeInfectedCritical);
190+
model.parameters.get<mio::glsecir::TransitionMatrixInfectedCriticalToRecovered>() =
191+
mio::glsecir::TransitionMatrixInfectedCriticalToRecovered().get_default(
192+
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedCritical>() / 2.), timeInfectedCritical);
193+
194+
model.parameters.get<mio::glsecir::TransmissionProbabilityOnContact>() = 0.05;
195+
196+
mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::glsecir::ContactPatterns>();
197+
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10));
198+
// From SimulationTime 5, the contact pattern is reduced to 30% of the initial value.
199+
contact_matrix[0].add_damping(0.7, mio::SimulationTime(5.));
200+
201+
model.parameters.get<mio::glsecir::RelativeTransmissionNoSymptoms>() = 0.7;
202+
model.parameters.get<mio::glsecir::RiskOfInfectionFromSymptomatic>() = 0.25;
203+
204+
// Perform a simulation.
205+
mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(t0, tmax, dt_init, model);
206+
// The simulation result is divided by subcompartments as in the LCT model.
207+
// We call the function calculate_compartments to get a result according to the InfectionStates.
208+
mio::TimeSeries<ScalarType> population_no_subcompartments = model.calculate_compartments(result);
209+
auto interpolated_result = mio::interpolate_simulation_result(population_no_subcompartments, 0.1);
210+
interpolated_result.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 10, 4);
211+
}

‎cpp/examples/lct_secir.cpp‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -47,9 +47,9 @@ int main()
4747

4848
// Variable defines whether the class Initializer is used to define an initial vector from flows or whether a manually
4949
// defined initial vector is used to initialize the LCT model.
50-
bool use_initializer_flows = true;
50+
bool use_initializer_flows = false;
5151

52-
ScalarType tmax = 20;
52+
ScalarType tmax = 10;
5353

5454
// Set Parameters.
5555
model.parameters.get<mio::lsecir::TimeExposed>()[0] = 3.2;
@@ -155,10 +155,9 @@ int main()
155155

156156
// Perform a simulation.
157157
mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(0, tmax, 0.5, model);
158-
//result.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
159158
// The simulation result is divided by subcompartments.
160-
// We call the function calculate_comparttments to get a result according to the InfectionStates.
159+
// We call the function calculate_compartments to get a result according to the InfectionStates.
161160
mio::TimeSeries<ScalarType> population_no_subcompartments = model.calculate_compartments(result);
162161
auto interpolated_results = mio::interpolate_simulation_result(population_no_subcompartments);
163-
interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
162+
interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 12, 4);
164163
}

‎cpp/memilio/epidemiology/lct_infection_state.h‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,9 @@
2020
#ifndef MIO_EPI_LCT_INFECTION_STATE_H
2121
#define MIO_EPI_LCT_INFECTION_STATE_H
2222

23+
#include "memilio/config.h"
24+
#include "memilio/math/eigen.h"
25+
2326
#include <array>
2427

2528
namespace mio
@@ -78,6 +81,53 @@ class LctInfectionState
7881
return index;
7982
}
8083

84+
/**
85+
* @brief Cumulates a vector with the number of individuals in each subcompartment (with subcompartments
86+
* according to the LctInfectionState) to produce a Vector that divides the population only into the infection
87+
* states defined in InfectionStates.
88+
*
89+
* @param[in] subcompartments Vector with number of individuals in each subcompartment.
90+
* The size of the vector has to match the LctInfectionState.
91+
* @return Vector with accumulated values for the InfectionStates.
92+
*/
93+
static Vector<ScalarType> calculate_compartments(const Vector<ScalarType>& subcompartments)
94+
{
95+
assert(subcompartments.rows() == Count);
96+
97+
Vector<ScalarType> compartments((Eigen::Index)InfectionState::Count);
98+
// Use segment of the vector subcompartments of each InfectionState and sum up the values of subcompartments.
99+
compartments[(Eigen::Index)InfectionState::Susceptible] = subcompartments[0];
100+
compartments[(Eigen::Index)InfectionState::Exposed] =
101+
subcompartments
102+
.segment(get_first_index<InfectionState::Exposed>(), get_num_subcompartments<InfectionState::Exposed>())
103+
.sum();
104+
compartments[(Eigen::Index)InfectionState::InfectedNoSymptoms] =
105+
subcompartments
106+
.segment(get_first_index<InfectionState::InfectedNoSymptoms>(),
107+
get_num_subcompartments<InfectionState::InfectedNoSymptoms>())
108+
.sum();
109+
compartments[(Eigen::Index)InfectionState::InfectedSymptoms] =
110+
subcompartments
111+
.segment(get_first_index<InfectionState::InfectedSymptoms>(),
112+
get_num_subcompartments<InfectionState::InfectedSymptoms>())
113+
.sum();
114+
compartments[(Eigen::Index)InfectionState::InfectedSevere] =
115+
subcompartments
116+
.segment(get_first_index<InfectionState::InfectedSevere>(),
117+
get_num_subcompartments<InfectionState::InfectedSevere>())
118+
.sum();
119+
compartments[(Eigen::Index)InfectionState::InfectedCritical] =
120+
subcompartments
121+
.segment(get_first_index<InfectionState::InfectedCritical>(),
122+
get_num_subcompartments<InfectionState::InfectedCritical>())
123+
.sum();
124+
compartments[(Eigen::Index)InfectionState::Recovered] =
125+
subcompartments[get_first_index<InfectionState::Recovered>()];
126+
compartments[(Eigen::Index)InfectionState::Dead] = subcompartments[get_first_index<InfectionState::Dead>()];
127+
128+
return compartments;
129+
}
130+
81131
static constexpr size_t Count{(... + Ns)};
82132

83133
private:

‎cpp/models/glct_secir/CMakeLists.txt‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
add_library(glct_secir
2+
infection_state.h
3+
model.h
4+
model.cpp
5+
parameters.h
6+
)
7+
target_link_libraries(glct_secir PUBLIC memilio)
8+
target_include_directories(glct_secir PUBLIC
9+
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
10+
$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
11+
)
12+
target_compile_options(glct_secir PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

0 commit comments

Comments
 (0)