Skip to content

Commit b131549

Browse files
lenaploetzkereneSchm
lenaploetzke
and
reneSchm
authored
914 make the lct secir model a derived class of flowmodel (#1055)
- make LCT model a derived class of CompartmentalModel - use simulation function of CompartmentalModel instead of separate function - LctState works with size_t instead of int Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com>
1 parent 88e0439 commit b131549

File tree

9 files changed

+320
-591
lines changed

9 files changed

+320
-591
lines changed

‎cpp/examples/lct_secir.cpp‎

Lines changed: 37 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -19,90 +19,67 @@
1919
*/
2020

2121
#include "lct_secir/model.h"
22-
#include "lct_secir/infection_state.h"
23-
#include "lct_secir/simulation.h"
2422
#include "memilio/config.h"
2523
#include "memilio/utils/time_series.h"
2624
#include "memilio/epidemiology/uncertain_matrix.h"
2725
#include "memilio/math/eigen.h"
2826
#include "memilio/utils/logging.h"
27+
#include "memilio/compartments/simulation.h"
28+
#include "memilio/data/analyze_result.h"
2929

3030
#include <vector>
3131

3232
int main()
3333
{
3434
// Simple example to demonstrate how to run a simulation using an LCT SECIR model.
3535
// Parameters, initial values and the number of subcompartments are not meant to represent a realistic scenario.
36+
constexpr size_t NumExposed = 2, NumInfectedNoSymptoms = 3, NumInfectedSymptoms = 1, NumInfectedSevere = 1,
37+
NumInfectedCritical = 5;
38+
using Model = mio::lsecir::Model<NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms, NumInfectedSevere,
39+
NumInfectedCritical>;
40+
using LctState = Model::LctState;
41+
using InfectionState = LctState::InfectionState;
3642

37-
using Model = mio::lsecir::Model<2, 3, 1, 1, 5>;
38-
using LctState = Model::LctState;
43+
Model model;
3944

4045
ScalarType tmax = 20;
4146

42-
// Define the initial value vector init with the distribution of the population into subcompartments.
43-
// This method of defining the vector using a vector of vectors is a bit of overhead, but should remind you how
44-
// the entries of the initial value vector relate to the defined template parameters of the model or the number of subcompartments.
45-
// It is also possible to define the initial value vector directly.
47+
// Define the initial values with the distribution of the population into subcompartments.
48+
// This method of defining the initial values using a vector of vectors is not necessary, but should remind you
49+
// how the entries of the initial value vector relate to the defined template parameters of the model or the number
50+
// of subcompartments. It is also possible to define the initial values directly.
4651
std::vector<std::vector<ScalarType>> initial_populations = {{750}, {30, 20}, {20, 10, 10}, {50},
4752
{50}, {10, 10, 5, 3, 2}, {20}, {10}};
4853

4954
// Assert that initial_populations has the right shape.
50-
if (initial_populations.size() != (size_t)LctState::InfectionState::Count) {
55+
if (initial_populations.size() != (size_t)InfectionState::Count) {
5156
mio::log_error("The number of vectors in initial_populations does not match the number of InfectionStates.");
5257
return 1;
5358
}
54-
if ((initial_populations[(int)LctState::InfectionState::Susceptible].size() !=
55-
(size_t)LctState::get_num_subcompartments<LctState::InfectionState::Susceptible>()) ||
56-
(initial_populations[(int)LctState::InfectionState::Exposed].size() !=
57-
(size_t)LctState::get_num_subcompartments<LctState::InfectionState::Exposed>()) ||
58-
(initial_populations[(int)LctState::InfectionState::InfectedNoSymptoms].size() !=
59-
(size_t)LctState::get_num_subcompartments<LctState::InfectionState::InfectedNoSymptoms>()) ||
60-
(initial_populations[(int)LctState::InfectionState::InfectedSymptoms].size() !=
61-
(size_t)LctState::get_num_subcompartments<LctState::InfectionState::InfectedSymptoms>()) ||
62-
(initial_populations[(int)LctState::InfectionState::InfectedSevere].size() !=
63-
(size_t)LctState::get_num_subcompartments<LctState::InfectionState::InfectedSevere>()) ||
64-
(initial_populations[(int)LctState::InfectionState::InfectedCritical].size() !=
65-
(size_t)LctState::get_num_subcompartments<LctState::InfectionState::InfectedCritical>()) ||
66-
(initial_populations[(int)LctState::InfectionState::Recovered].size() !=
67-
(size_t)LctState::get_num_subcompartments<LctState::InfectionState::Recovered>()) ||
68-
(initial_populations[(int)LctState::InfectionState::Dead].size() !=
69-
(size_t)LctState::get_num_subcompartments<LctState::InfectionState::Dead>())) {
59+
if ((initial_populations[(size_t)InfectionState::Susceptible].size() !=
60+
LctState::get_num_subcompartments<InfectionState::Susceptible>()) ||
61+
(initial_populations[(size_t)InfectionState::Exposed].size() != NumExposed) ||
62+
(initial_populations[(size_t)InfectionState::InfectedNoSymptoms].size() != NumInfectedNoSymptoms) ||
63+
(initial_populations[(size_t)InfectionState::InfectedSymptoms].size() != NumInfectedSymptoms) ||
64+
(initial_populations[(size_t)InfectionState::InfectedSevere].size() != NumInfectedSevere) ||
65+
(initial_populations[(size_t)InfectionState::InfectedCritical].size() != NumInfectedCritical) ||
66+
(initial_populations[(size_t)InfectionState::Recovered].size() !=
67+
LctState::get_num_subcompartments<InfectionState::Recovered>()) ||
68+
(initial_populations[(size_t)InfectionState::Dead].size() !=
69+
LctState::get_num_subcompartments<InfectionState::Dead>())) {
7070
mio::log_error("The length of at least one vector in initial_populations does not match the related number of "
7171
"subcompartments.");
7272
return 1;
7373
}
7474

75-
// Transfer the initial values in initial_populations to the vector init.
76-
Eigen::VectorXd init = Eigen::VectorXd::Zero(LctState::Count);
77-
init[LctState::get_first_index<LctState::InfectionState::Susceptible>()] =
78-
initial_populations[(int)LctState::InfectionState::Susceptible][0];
79-
for (int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::Exposed>(); i++) {
80-
init[LctState::get_first_index<LctState::InfectionState::Exposed>() + i] =
81-
initial_populations[(int)LctState::InfectionState::Exposed][i];
75+
// Transfer the initial values in initial_populations to the model.
76+
std::vector<ScalarType> flat_initial_populations;
77+
for (auto&& vec : initial_populations) {
78+
flat_initial_populations.insert(flat_initial_populations.end(), vec.begin(), vec.end());
8279
}
83-
for (int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedNoSymptoms>(); i++) {
84-
init[LctState::get_first_index<LctState::InfectionState::InfectedNoSymptoms>() + i] =
85-
initial_populations[(int)LctState::InfectionState::InfectedNoSymptoms][i];
80+
for (size_t i = 0; i < LctState::Count; i++) {
81+
model.populations[mio::Index<LctState>(i)] = flat_initial_populations[i];
8682
}
87-
for (int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedSymptoms>(); i++) {
88-
init[LctState::get_first_index<LctState::InfectionState::InfectedSymptoms>() + i] =
89-
initial_populations[(int)LctState::InfectionState::InfectedSymptoms][i];
90-
}
91-
for (int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedSevere>(); i++) {
92-
init[LctState::get_first_index<LctState::InfectionState::InfectedSevere>() + i] =
93-
initial_populations[(int)LctState::InfectionState::InfectedSevere][i];
94-
}
95-
for (int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedCritical>(); i++) {
96-
init[LctState::get_first_index<LctState::InfectionState::InfectedCritical>() + i] =
97-
initial_populations[(int)LctState::InfectionState::InfectedCritical][i];
98-
}
99-
init[LctState::get_first_index<LctState::InfectionState::Recovered>()] =
100-
initial_populations[(int)LctState::InfectionState::Recovered][0];
101-
init[LctState::get_first_index<LctState::InfectionState::Dead>()] =
102-
initial_populations[(int)LctState::InfectionState::Dead][0];
103-
104-
// Initialize model.
105-
Model model(std::move(init));
10683

10784
// Set Parameters.
10885
model.parameters.get<mio::lsecir::TimeExposed>() = 3.2;
@@ -127,8 +104,10 @@ int main()
127104
model.parameters.set<mio::lsecir::DeathsPerCritical>(0.3);
128105

129106
// Perform a simulation.
130-
mio::TimeSeries<ScalarType> result = mio::lsecir::simulate(0, tmax, 0.5, model);
131-
// Calculate the distribution in the InfectionState%s without subcompartments of the result and print it.
132-
mio::TimeSeries<ScalarType> population_no_subcompartments = model.calculate_populations(result);
133-
population_no_subcompartments.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
107+
mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(0, tmax, 0.5, model);
108+
// The simulation result is divided by subcompartments.
109+
// We call the function calculate_comparttments to get a result according to the InfectionStates.
110+
mio::TimeSeries<ScalarType> population_no_subcompartments = model.calculate_compartments(result);
111+
auto interpolated_results = mio::interpolate_simulation_result(population_no_subcompartments);
112+
interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
134113
}

‎cpp/memilio/compartments/compartmentalmodel.h‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -88,18 +88,18 @@ struct CompartmentalModel {
8888
{
8989
}
9090

91-
CompartmentalModel(const CompartmentalModel&) = default;
92-
CompartmentalModel(CompartmentalModel&&) = default;
91+
CompartmentalModel(const CompartmentalModel&) = default;
92+
CompartmentalModel(CompartmentalModel&&) = default;
9393
CompartmentalModel& operator=(const CompartmentalModel&) = default;
94-
CompartmentalModel& operator=(CompartmentalModel&&) = default;
95-
virtual ~CompartmentalModel() = default;
94+
CompartmentalModel& operator=(CompartmentalModel&&) = default;
95+
virtual ~CompartmentalModel() = default;
9696

9797
//REMARK: Not pure virtual for easier java/python bindings
9898
virtual void get_derivatives(Eigen::Ref<const Vector<FP>>, Eigen::Ref<const Vector<FP>> /*y*/, FP /*t*/,
9999
Eigen::Ref<Vector<FP>> /*dydt*/) const {};
100100

101101
/**
102-
* @brief eval_right_hand_side evaulates the right-hand-side f of the ODE dydt = f(y, t)
102+
* @brief eval_right_hand_side evaluates the right-hand-side f of the ODE dydt = f(y, t)
103103
*
104104
* The heart of the compartmental model is a first order ODE dydt = f(y,t), where y is a flat
105105
* representation of all the compartmental populations at time t. This function evaluates the

‎cpp/memilio/epidemiology/lct_infection_state.h‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -31,28 +31,28 @@ namespace mio
3131
* @tparam Ns Number of subcompartments for each infection state defined in InfectionState.
3232
* The number of given template arguments must be equal to the entry Count from InfectionState.
3333
*/
34-
template <class InfectionStates, int... Ns>
34+
template <class InfectionStates, size_t... Ns>
3535
class LctInfectionState
3636
{
3737
public:
3838
using InfectionState = InfectionStates;
3939
static_assert((size_t)InfectionState::Count == sizeof...(Ns),
40-
"The number of integers provided as template parameters must be "
40+
"The number of the size_t's provided as template parameters must be "
4141
"the same as the entry Count of InfectionState.");
4242

4343
static_assert(((Ns > 0) && ...), "The number of subcompartments must be at least 1.");
4444

4545
/**
4646
* @brief Gets the number of subcompartments in an infection state.
4747
*
48-
* @tparam State: Infection state for which the number of subcompartments should be returned.
48+
* @tparam State Infection state for which the number of subcompartments should be returned.
4949
* @return Number of subcompartments for State. Returned value is always at least one.
5050
*/
5151
template <InfectionState State>
52-
static constexpr int get_num_subcompartments()
52+
static constexpr size_t get_num_subcompartments()
5353
{
5454
static_assert(State < InfectionState::Count, "State must be a a valid InfectionState.");
55-
return m_subcompartment_numbers[(int)State];
55+
return m_subcompartment_numbers[(size_t)State];
5656
}
5757

5858
/**
@@ -65,20 +65,20 @@ class LctInfectionState
6565
* Returned value is always non-negative.
6666
*/
6767
template <InfectionState State>
68-
static constexpr int get_first_index()
68+
static constexpr size_t get_first_index()
6969
{
7070
static_assert(State < InfectionState::Count, "State must be a a valid InfectionState.");
71-
int index = 0;
72-
for (int i = 0; i < (int)(State); i++) {
71+
size_t index = 0;
72+
for (size_t i = 0; i < (size_t)(State); i++) {
7373
index = index + m_subcompartment_numbers[i];
7474
}
7575
return index;
7676
}
7777

78-
static constexpr int Count{(... + Ns)};
78+
static constexpr size_t Count{(... + Ns)};
7979

8080
private:
81-
static constexpr const std::array<int, sizeof...(Ns)> m_subcompartment_numbers{
81+
static constexpr const std::array<size_t, sizeof...(Ns)> m_subcompartment_numbers{
8282
Ns...}; ///< Vector which defines the number of subcompartments for each infection state of InfectionState.
8383
};
8484

‎cpp/models/lct_secir/CMakeLists.txt‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,6 @@ add_library(lct_secir
22
infection_state.h
33
model.h
44
model.cpp
5-
simulation.h
65
parameters.h
76
parameters_io.h
87
)

0 commit comments

Comments
 (0)