Skip to content

Commit f5e415a

Browse files
HenrZumknaranja
HenrZu
and
mknaranja
authored
945 Allow contact increase for simulation of larger events (#975)
- Remove the assertion for negative Damping coefficients to allow modeling of contact increases for the simulation of larger events - Add example and software tests for contact changes in ODE SECIR model Co-authored-by: Martin Kühn <[email protected]>
1 parent b1b8f37 commit f5e415a

8 files changed

Lines changed: 215 additions & 12 deletions

File tree

‎cpp/examples/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,10 @@ add_executable(ode_secir_example ode_secir.cpp)
3838
target_link_libraries(ode_secir_example PRIVATE memilio ode_secir)
3939
target_compile_options(ode_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
4040

41+
add_executable(ode_secir_contact_changes ode_secir_contact_changes.cpp)
42+
target_link_libraries(ode_secir_contact_changes PRIVATE memilio ode_secir)
43+
target_compile_options(ode_secir_contact_changes PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
44+
4145
add_executable(ode_secirvvs_example ode_secirvvs.cpp)
4246
target_link_libraries(ode_secirvvs_example PRIVATE memilio ode_secirvvs)
4347
target_compile_options(ode_secirvvs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/examples/ode_secir_contact_changes.cpp‎

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
/*
2+
* Copyright (C) 2020-2024 MEmilio
3+
*
4+
* Authors: Daniel Abele, 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+
#include "ode_secir/model.h"
21+
#include "memilio/compartments/simulation.h"
22+
#include "memilio/utils/logging.h"
23+
#include "memilio/math/euler.h"
24+
25+
/*
26+
* This example demonstrates how to realize contact behavior changes with any of our ordinary differential
27+
* equation-based models. This example can thus easily be adapted for other models like osecirvvs, oseir, osir etc.
28+
* We print out the flows (i.e., new transmissions, infections, hospitalizations etc. per time point.) As we use an
29+
* fixed-time step Explicit Euler, we can compare them.
30+
*/
31+
int main()
32+
{
33+
mio::set_log_level(mio::LogLevel::warn);
34+
35+
double t0 = 0;
36+
double dt = 0.1;
37+
double tmax = 1;
38+
39+
double cont_freq = 10;
40+
41+
double nb_total_t0 = 1000, nb_inf_t0 = 10;
42+
43+
auto integrator = std::make_shared<mio::EulerIntegratorCore>();
44+
45+
// default model run to be compared against
46+
mio::osecir::Model model_a(1);
47+
const auto indx_flow_SE =
48+
model_a.get_flat_flow_index<mio::osecir::InfectionState::Susceptible, mio::osecir::InfectionState::Exposed>(
49+
{mio::AgeGroup(0)});
50+
51+
model_a.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0;
52+
model_a.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
53+
nb_total_t0);
54+
mio::ContactMatrixGroup& contact_matrix_a = model_a.parameters.get<mio::osecir::ContactPatterns>();
55+
contact_matrix_a[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
56+
// set probability of transmission and risk of infection to 1.
57+
model_a.parameters.get<mio::osecir::TransmissionProbabilityOnContact>() = 1.0;
58+
model_a.parameters.get<mio::osecir::RiskOfInfectionFromSymptomatic>() = 1.0;
59+
auto result_a = simulate_flows(t0, tmax, dt, model_a, integrator);
60+
result_a[1].print_table({"S->E", "E->I_NS", "I_NS->I_Sy", "I_NS->R", "I_NSC->I_SyC", "I_NSC->R", "I_Sy->I_Sev",
61+
"I_Sy->R", "I_SyC->I_Sev", "I_SyC->R", "I_Sev->I_Crit", "I_Sev->R", "I_Sev->D",
62+
"I_Crit->D", "I_Crit->R"},
63+
4, 4);
64+
std::cout << "With default contacts, the number of new transmissions (flow from S->E) in first time step is: "
65+
<< result_a[1].get_value(1)[indx_flow_SE] << ".\n";
66+
67+
// The contacts are halfed: reduced transmission through damping with value 0.5
68+
mio::osecir::Model model_b{model_a};
69+
model_b.populations.set_total(nb_total_t0);
70+
model_b.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0;
71+
model_b.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
72+
nb_total_t0);
73+
mio::ContactMatrixGroup& contact_matrix_b = model_b.parameters.get<mio::osecir::ContactPatterns>();
74+
contact_matrix_b[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
75+
contact_matrix_b[0].add_damping(0.5, mio::SimulationTime(0.)); // contact reduction happens here!
76+
auto result_b = simulate_flows(t0, tmax, dt, model_b, integrator);
77+
result_b[1].print_table({"S->E", "E->I_NS", "I_NS->I_Sy", "I_NS->R", "I_NSC->I_SyC", "I_NSC->R", "I_Sy->I_Sev",
78+
"I_Sy->R", "I_SyC->I_Sev", "I_SyC->R", "I_Sev->I_Crit", "I_Sev->R", "I_Sev->D",
79+
"I_Crit->D", "I_Crit->R"},
80+
4, 4);
81+
std::cout << "With contacts reduced to a half of the original example, the number of new transmissions (flow from "
82+
"S->E) in first time step is: "
83+
<< result_b[1].get_value(1)[indx_flow_SE] << ".\n";
84+
85+
// No contacts at all: no transmission through damping with value 1.
86+
mio::osecir::Model model_c{model_a};
87+
model_c.populations.set_total(nb_total_t0);
88+
model_c.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0;
89+
model_c.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
90+
nb_total_t0);
91+
mio::ContactMatrixGroup& contact_matrix_c = model_c.parameters.get<mio::osecir::ContactPatterns>();
92+
contact_matrix_c[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
93+
contact_matrix_c[0].add_damping(1., mio::SimulationTime(0.)); // contact reduction happens here!
94+
auto result_c = simulate_flows(t0, tmax, dt, model_c, integrator);
95+
result_c[1].print_table({"S->E", "E->I_NS", "I_NS->I_Sy", "I_NS->R", "I_NSC->I_SyC", "I_NSC->R", "I_Sy->I_Sev",
96+
"I_Sy->R", "I_SyC->I_Sev", "I_SyC->R", "I_Sev->I_Crit", "I_Sev->R", "I_Sev->D",
97+
"I_Crit->D", "I_Crit->R"},
98+
4, 4);
99+
std::cout
100+
<< "With contacts reduced to zero, the number of new transmissions (flow from S->E) in first time step is: "
101+
<< result_c[1].get_value(1)[indx_flow_SE] << ".\n";
102+
103+
// The contacts are doubled: increased transmission through damping with value -1.
104+
mio::osecir::Model model_d{model_a};
105+
model_d.populations.set_total(nb_total_t0);
106+
model_d.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0;
107+
model_d.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
108+
nb_total_t0);
109+
mio::ContactMatrixGroup& contact_matrix_d = model_d.parameters.get<mio::osecir::ContactPatterns>();
110+
contact_matrix_d[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
111+
contact_matrix_d[0].add_damping(-1., mio::SimulationTime(0.)); // contact increase happens here!
112+
auto result_d = simulate_flows(t0, tmax, dt, model_d, integrator);
113+
result_d[1].print_table({"S->E", "E->I_NS", "I_NS->I_Sy", "I_NS->R", "I_NSC->I_SyC", "I_NSC->R", "I_Sy->I_Sev",
114+
"I_Sy->R", "I_SyC->I_Sev", "I_SyC->R", "I_Sev->I_Crit", "I_Sev->R", "I_Sev->D",
115+
"I_Crit->D", "I_Crit->R"},
116+
4, 4);
117+
std::cout << "With contacts doubled, the number of new transmissions (flow from S->E) in first time step is: "
118+
<< result_d[1].get_value(1)[indx_flow_SE] << "\n";
119+
}

‎cpp/examples/ode_seir.cpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ int main()
2929

3030
double t0 = 0;
3131
double tmax = 50.;
32-
double dt = 0.1;
32+
double dt = 1.0;
3333

3434
mio::log_info("Simulating ODE SEIR; t={} ... {} with dt = {}.", t0, tmax, dt);
3535

@@ -50,7 +50,9 @@ int main()
5050
model.parameters.set<mio::oseir::TimeInfected>(2);
5151

5252
model.parameters.get<mio::oseir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
53-
model.parameters.get<mio::oseir::ContactPatterns>().add_damping(0.6, mio::SimulationTime(12.5));
53+
54+
// contacts increase by 100% after 12.5 days
55+
model.parameters.get<mio::oseir::ContactPatterns>().add_damping(-1., mio::SimulationTime(12.5));
5456

5557
model.check_constraints();
5658

‎cpp/memilio/epidemiology/damping.h‎

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -95,8 +95,7 @@ class Damping : public std::tuple<typename S::Matrix, DampingLevel, DampingType,
9595
Damping(const Eigen::MatrixBase<ME>& m, DampingLevel level, DampingType type, SimulationTime t)
9696
: Base(m, level, type, t)
9797
{
98-
assert((get_coeffs().array() >= 0.).all() && (get_coeffs().array() <= 1.).all() &&
99-
"damping coefficient out of range");
98+
assert((get_coeffs().array() <= 1.).all() && "damping coefficient out of range");
10099
}
101100

102101
/**
@@ -352,7 +351,8 @@ class Dampings
352351
*/
353352
auto get_matrix_at(SimulationTime t) const
354353
{
355-
assert(!m_accumulated_dampings_cached.empty() && "Cache is not current. Did you disable the automatic cache update?");
354+
assert(!m_accumulated_dampings_cached.empty() &&
355+
"Cache is not current. Did you disable the automatic cache update?");
356356
auto ub =
357357
std::upper_bound(m_accumulated_dampings_cached.begin(), m_accumulated_dampings_cached.end(),
358358
std::make_tuple(t), [](auto&& tup1, auto&& tup2) {
@@ -555,8 +555,7 @@ void Dampings<D>::update_cache()
555555
//update active damping
556556
update_active_dampings(damping, active_by_type, sum_by_level);
557557
auto combined_damping = inclusive_exclusive_sum(sum_by_level);
558-
assert((combined_damping.array() <= 1).all() && (combined_damping.array() >= 0).all() &&
559-
"unexpected error, accumulated damping out of range.");
558+
assert((combined_damping.array() <= 1).all() && "unexpected error, accumulated damping out of range.");
560559
if (floating_point_equal(double(get<SimulationTime>(damping)),
561560
double(get<SimulationTime>(m_accumulated_dampings_cached.back())), 1e-15, 1e-15)) {
562561
std::get<Matrix>(m_accumulated_dampings_cached.back()) = combined_damping;
@@ -604,13 +603,13 @@ void Dampings<S>::update_active_dampings(
604603
//replace active of the same type and level
605604
auto& active_same_type = *iter_active_same_type;
606605
//find active with the same level
607-
auto& sum_same_level = *std::find_if(sum_by_level.begin(), sum_by_level.end(), [&damping](auto& sum) {
606+
auto& sum_same_level = *std::find_if(sum_by_level.begin(), sum_by_level.end(), [&damping](auto& sum) {
608607
return get<DampingLevel>(sum) == get<DampingLevel>(damping);
609608
});
610609
//remove active with the same type and level and add new one
611610
get<MatrixIdx>(sum_same_level) += get<MatrixIdx>(damping) - get<MatrixIdx>(active_same_type).get();
612611
//avoid negative values due to rounding error if e.g. a previous damping is lifted
613-
get<MatrixIdx>(sum_same_level) = get<MatrixIdx>(sum_same_level).cwiseMax(0.);
612+
get<MatrixIdx>(sum_same_level) = get<MatrixIdx>(sum_same_level).cwiseMax(0.);
614613
get<MatrixIdx>(active_same_type) = get<MatrixIdx>(damping);
615614
}
616615
else {

‎cpp/models/ode_secir/parameters.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -221,7 +221,7 @@ struct RiskOfInfectionFromSymptomatic {
221221
using Type = CustomIndexArray<UncertainValue, AgeGroup>;
222222
static Type get_default(AgeGroup size)
223223
{
224-
return Type(size, 0.);
224+
return Type(size, 1.);
225225
}
226226
static std::string name()
227227
{

‎cpp/models/ode_secirvvs/parameters.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -281,7 +281,7 @@ struct RiskOfInfectionFromSymptomatic {
281281
using Type = CustomIndexArray<UncertainValue, AgeGroup>;
282282
static Type get_default(AgeGroup size)
283283
{
284-
return Type(size, 0.);
284+
return Type(size, 1.);
285285
}
286286
static std::string name()
287287
{

‎cpp/tests/test_damping.cpp‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,19 @@ TEST(TestDampings, dampingOfSameType)
8080
EXPECT_THAT(print_wrap(dampings.get_matrix_at(1e5)), MatrixNear(D2));
8181
}
8282

83+
TEST(TestDampings, dampingsNegative)
84+
{
85+
mio::Dampings<mio::Damping<mio::RectMatrixShape>> dampings(2, 2);
86+
auto D1 = -4.25;
87+
auto D2 = (Eigen::MatrixXd(2, 2) << 0.25, -0.5, 0.75, -1).finished();
88+
dampings.add(D1, mio::DampingLevel(7), mio::DampingType(3), mio::SimulationTime(0.5));
89+
dampings.add(D2, mio::DampingLevel(13), mio::DampingType(3), mio::SimulationTime(2.0));
90+
EXPECT_EQ(print_wrap(dampings.get_matrix_at(-1e5)), print_wrap(Eigen::MatrixXd::Zero(2, 2)));
91+
EXPECT_EQ(print_wrap(dampings.get_matrix_at(-0.5)), print_wrap(Eigen::MatrixXd::Zero(2, 2)));
92+
EXPECT_THAT(print_wrap(dampings.get_matrix_at(0.5 + 1e-32)), MatrixNear(Eigen::MatrixXd::Constant(2, 2, D1)));
93+
EXPECT_THAT(print_wrap(dampings.get_matrix_at(1e5)), MatrixNear((D1 + D2.array() - D1 * D2.array()).matrix()));
94+
}
95+
8396
TEST(TestDampings, dampingsCombined)
8497
{
8598
mio::Dampings<mio::Damping<mio::SquareMatrixShape>> dampings(2);

‎cpp/tests/test_odesecir.cpp‎

Lines changed: 67 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -555,6 +555,72 @@ TEST(TestOdeSecir, testSettersAndGetters)
555555
EXPECT_EQ(vec[21], model.parameters.get<mio::osecir::Seasonality>());
556556
}
557557

558+
TEST(TestOdeSecir, testDamping)
559+
{
560+
// Test functionality of dampings
561+
// (initially only implemented contact reductions but now also allow contact increases).
562+
// Contact matrices with dampings are cosine-smoothed in decline/increase along one day to be C1 differentiable.
563+
// If EulerIntegratorCore with dt=1 is used, we jump across this smoothing so that we can express the relationship
564+
// between old and new transmission directly, only including damping, contact, and transmission probability values.
565+
double t0 = 0;
566+
double dt = 1;
567+
double tmax = dt;
568+
569+
double cont_freq = 10;
570+
571+
double nb_total_t0 = 1000, nb_inf_t0 = 10;
572+
573+
auto integrator = std::make_shared<mio::EulerIntegratorCore>();
574+
575+
// default model run to be compared against
576+
mio::osecir::Model model_a(1);
577+
model_a.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0;
578+
model_a.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
579+
nb_total_t0);
580+
mio::ContactMatrixGroup& contact_matrix_a = model_a.parameters.get<mio::osecir::ContactPatterns>();
581+
contact_matrix_a[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
582+
// set probability of transmission and risk of infection to 1.
583+
model_a.parameters.get<mio::osecir::TransmissionProbabilityOnContact>() = 1.0;
584+
model_a.parameters.get<mio::osecir::RiskOfInfectionFromSymptomatic>() = 1.0;
585+
auto result_a = simulate_flows(t0, tmax, dt, model_a, integrator);
586+
587+
// reduced transmission
588+
mio::osecir::Model model_b{model_a};
589+
model_b.populations.set_total(nb_total_t0);
590+
model_b.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0;
591+
model_b.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
592+
nb_total_t0);
593+
mio::ContactMatrixGroup& contact_matrix_b = model_b.parameters.get<mio::osecir::ContactPatterns>();
594+
contact_matrix_b[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
595+
contact_matrix_b[0].add_damping(0.5, mio::SimulationTime(0.));
596+
auto result_b = simulate_flows(t0, tmax, dt, model_b, integrator);
597+
EXPECT_EQ(2 * result_b[1].get_last_value()[0], result_a[1].get_last_value()[0]);
598+
599+
// no transmission
600+
mio::osecir::Model model_c{model_a};
601+
model_c.populations.set_total(nb_total_t0);
602+
model_c.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0;
603+
model_c.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
604+
nb_total_t0);
605+
mio::ContactMatrixGroup& contact_matrix_c = model_c.parameters.get<mio::osecir::ContactPatterns>();
606+
contact_matrix_c[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
607+
contact_matrix_c[0].add_damping(1., mio::SimulationTime(0.));
608+
auto result_c = simulate_flows(t0, tmax, dt, model_c, integrator);
609+
EXPECT_EQ(result_c[1].get_last_value()[0], 0.0);
610+
611+
// increased transmission to a factor of two (by +1)
612+
mio::osecir::Model model_d{model_a};
613+
model_d.populations.set_total(nb_total_t0);
614+
model_d.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0;
615+
model_d.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
616+
nb_total_t0);
617+
mio::ContactMatrixGroup& contact_matrix_d = model_d.parameters.get<mio::osecir::ContactPatterns>();
618+
contact_matrix_d[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
619+
contact_matrix_d[0].add_damping(-1., mio::SimulationTime(0.));
620+
auto result_d = simulate_flows(t0, tmax, dt, model_d, integrator);
621+
EXPECT_EQ(2 * result_a[1].get_last_value()[0], result_d[1].get_last_value()[0]);
622+
}
623+
558624
TEST(TestOdeSecir, testModelConstraints)
559625
{
560626
mio::set_log_level(
@@ -614,6 +680,7 @@ TEST(TestOdeSecir, testModelConstraints)
614680

615681
mio::TimeSeries<double> secihurd_interp = mio::interpolate_simulation_result(secihurd);
616682

683+
// Tests that infection numbers are higher in Winter season
617684
model.parameters.set<mio::osecir::StartDay>(100);
618685
model.parameters.set<mio::osecir::Seasonality>(0.5);
619686

@@ -634,7 +701,6 @@ TEST(TestOdeSecir, testModelConstraints)
634701
}
635702

636703
// temporary test for random variables
637-
638704
set_params_distributions_normal(model, t0, tmax, 0.2);
639705
model.parameters.set<mio::osecir::Seasonality>(mio::UncertainValue(0.0));
640706
model.parameters.set<mio::osecir::ICUCapacity>(mio::UncertainValue(8000));

0 commit comments

Comments
 (0)