Skip to content

Commit 800bdd2

Browse files
alexgit0HenrZuMaxBetzDLRmknaranja
authored
766 implementation of ode sir model (#775)
Co-authored-by: Henrik Zunker <[email protected]> Co-authored-by: MaxBetzDLR <[email protected]> Co-authored-by: Martin J. Kühn <[email protected]>
1 parent ce0a1f5 commit 800bdd2

File tree

17 files changed

+919
-0
lines changed

17 files changed

+919
-0
lines changed

‎cpp/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -122,6 +122,7 @@ if(MEMILIO_BUILD_MODELS)
122122
add_subdirectory(models/ide_secir)
123123
add_subdirectory(models/ide_seir)
124124
add_subdirectory(models/ode_seir)
125+
add_subdirectory(models/ode_sir)
125126
endif()
126127

127128
if(MEMILIO_BUILD_EXAMPLES)

‎cpp/examples/CMakeLists.txt‎

Lines changed: 4 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_sir_example ode_sir.cpp)
22+
target_link_libraries(ode_sir_example PRIVATE memilio ode_sir)
23+
24+
2125
add_executable(ode_secir_example ode_secir.cpp)
2226
target_link_libraries(ode_secir_example PRIVATE memilio ode_secir)
2327
target_compile_options(ode_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/examples/ode_sir.cpp‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
/*
2+
* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC)
3+
*
4+
* Authors: Daniel Abele, Jan Kleinert, Martin J. Kuehn
5+
*
6+
* Contact: Martin J. Kuehn <[email protected]>
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 "ode_sir/model.h"
22+
#include "ode_sir/infection_state.h"
23+
#include "ode_sir/parameters.h"
24+
#include "memilio/compartments/simulation.h"
25+
#include "memilio/utils/logging.h"
26+
27+
int main()
28+
{
29+
mio::set_log_level(mio::LogLevel::debug);
30+
31+
double t0 = 0.;
32+
double tmax = 50.;
33+
double dt = 0.1002004008016032;
34+
35+
double total_population = 1061000;
36+
37+
mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt);
38+
39+
mio::osir::Model model;
40+
41+
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Infected)}] = 1000;
42+
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Recovered)}] = 1000;
43+
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Susceptible)}] =
44+
total_population -
45+
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Infected)}] -
46+
model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Recovered)}];
47+
model.parameters.set<mio::osir::TimeInfected>(2);
48+
model.parameters.set<mio::osir::TransmissionProbabilityOnContact>(1);
49+
model.parameters.get<mio::osir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
50+
model.parameters.get<mio::osir::ContactPatterns>().add_damping(0.6, mio::SimulationTime(12.5));
51+
52+
auto integrator = std::make_shared<mio::EulerIntegratorCore>();
53+
54+
model.check_constraints();
55+
56+
auto sir = simulate(t0, tmax, dt, model, integrator);
57+
58+
bool print_to_terminal = true;
59+
60+
if (print_to_terminal) {
61+
std::vector<std::string> vars = {"S", "I", "R"};
62+
printf("\n # t");
63+
for (size_t k = 0; k < (size_t)mio::osir::InfectionState::Count; k++) {
64+
printf(" %s", vars[k].c_str());
65+
}
66+
67+
auto num_points = static_cast<size_t>(sir.get_num_time_points());
68+
for (size_t i = 0; i < num_points; i++) {
69+
printf("\n%.14f ", sir.get_time(i));
70+
Eigen::VectorXd res_j = sir.get_value(i);
71+
for (size_t j = 0; j < (size_t)mio::osir::InfectionState::Count; j++) {
72+
printf(" %.14f", res_j[j]);
73+
}
74+
}
75+
76+
Eigen::VectorXd res_j = sir.get_last_value();
77+
printf("number total: %f", res_j[0] + res_j[1] + res_j[2]);
78+
}
79+
}

‎cpp/models/ode_sir/CMakeLists.txt‎

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

‎cpp/models/ode_sir/README.md‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
2+
# ODE SIR compartment model
3+
4+
This model is a very simple ODE model with only three compartments and few parameters, mostly for demonstration of the MEmilio framework:
5+
- Susceptible, may become infected at any time
6+
- Infected, will be recovered after some time
7+
- Recovered, recovered from infectious process (dead or recovered)
8+
9+
We assume simulations over short periods of time, so that the population size can be considered constant and birth as well as (natural) mortality rates can be ignored.
10+
11+
Below is an overview of the model architecture and its compartments.
12+
13+
![SIR_model](https://github.com/DLR-SC/memilio/assets/69154294/01c9a2ae-2f5c-4bad-b7f0-34de651f2c73)
14+
| Mathematical variable | C++ variable name | Description |
15+
|---------------------------- | --------------- | -------------------------------------------------------------------------------------------------- |
16+
| $\phi$ | `ContactPatterns` | Daily contact rate / Number of daily contacts. |
17+
| $\rho$ | `TransmissionProbabilityOnContact` | Transmission risk for people located in the Susceptible compartment. |
18+
| $N$ | `populations.get_total()` | Total population. |
19+
| $T_{I}$ | `TimeInfected` | Time in days an individual stays in the Infected compartment. |
20+
21+
An example can be found in [examples/ode_sir.cpp](../../examples/ode_sir.cpp)

‎cpp/models/ode_sir/infection_state.h‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
/*
2+
* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC)
3+
*
4+
* Authors: Daniel Abele, Jan Kleinert, Martin J. Kuehn
5+
*
6+
* Contact: Martin J. Kuehn <[email protected]>
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+
#ifndef SIR_INFECTIONSTATE_H
22+
#define SIR_INFECTIONSTATE_H
23+
24+
namespace mio
25+
{
26+
namespace osir
27+
{
28+
29+
/**
30+
* @brief The InfectionState enum describes the possible
31+
* categories for the infectious state of persons
32+
*/
33+
enum class InfectionState
34+
{
35+
Susceptible,
36+
Infected,
37+
Recovered,
38+
Count
39+
};
40+
41+
} // namespace osir
42+
} // namespace mio
43+
44+
#endif // SIR_INFECTIONSTATE_H

‎cpp/models/ode_sir/model.cpp‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
/*
2+
* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC)
3+
*
4+
* Authors: Daniel Abele, Jan Kleinert, Martin J. Kuehn
5+
*
6+
* Contact: Martin J. Kuehn <[email protected]>
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 "ode_sir/model.h"
22+
23+
namespace mio
24+
{
25+
namespace osir
26+
{
27+
28+
} // namespace osir
29+
} // namespace mio

‎cpp/models/ode_sir/model.h‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
/*
2+
* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC)
3+
*
4+
* Authors: Daniel Abele, Jan Kleinert, Martin J. Kuehn
5+
*
6+
* Contact: Martin J. Kuehn <[email protected]>
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+
#ifndef SIR_MODEL_H
22+
#define SIR_MODEL_H
23+
24+
#include "memilio/compartments/compartmentalmodel.h"
25+
#include "memilio/epidemiology/populations.h"
26+
#include "memilio/epidemiology/contact_matrix.h"
27+
#include "ode_sir/infection_state.h"
28+
#include "ode_sir/parameters.h"
29+
30+
namespace mio
31+
{
32+
namespace osir
33+
{
34+
35+
/********************
36+
* define the model *
37+
********************/
38+
39+
class Model : public CompartmentalModel<InfectionState, Populations<InfectionState>, Parameters>
40+
{
41+
using Base = CompartmentalModel<InfectionState, mio::Populations<InfectionState>, Parameters>;
42+
43+
public:
44+
Model()
45+
: Base(Populations({InfectionState::Count}, 0.), ParameterSet())
46+
{
47+
}
48+
49+
void get_derivatives(Eigen::Ref<const Eigen::VectorXd> pop, Eigen::Ref<const Eigen::VectorXd> y, double t,
50+
Eigen::Ref<Eigen::VectorXd> dydt) const override
51+
{
52+
auto& params = this->parameters;
53+
double coeffStoI = params.get<ContactPatterns>().get_matrix_at(t)(0, 0) *
54+
params.get<TransmissionProbabilityOnContact>() / populations.get_total();
55+
56+
dydt[(size_t)InfectionState::Susceptible] =
57+
-coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected];
58+
dydt[(size_t)InfectionState::Infected] =
59+
coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected] -
60+
(1.0 / params.get<TimeInfected>()) * y[(size_t)InfectionState::Infected];
61+
dydt[(size_t)InfectionState::Recovered] =
62+
(1.0 / params.get<TimeInfected>()) * y[(size_t)InfectionState::Infected];
63+
}
64+
};
65+
66+
} // namespace osir
67+
} // namespace mio
68+
69+
#endif // SIR_MODEL_H

0 commit comments

Comments
 (0)