Skip to content

Commit f59d955

Browse files
MaxBetzDLRHenrZu
andauthored
300 python bindings for new vaccination model (#906)
Co-authored-by: Henrik Zunker <[email protected]>
1 parent 8ddb863 commit f59d955

14 files changed

Lines changed: 924 additions & 39 deletions

File tree

cpp/examples/ode_secirvvs.cpp

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -66,11 +66,16 @@ int main()
6666

6767
model.parameters.get<mio::osecirvvs::ICUCapacity>() = 100;
6868
model.parameters.get<mio::osecirvvs::TestAndTraceCapacity>() = 0.0143;
69-
model.parameters.get<mio::osecirvvs::DailyFirstVaccination>().resize(mio::SimulationDay(size_t(1000)));
70-
model.parameters.get<mio::osecirvvs::DailyFirstVaccination>().array().setConstant(5);
71-
model.parameters.get<mio::osecirvvs::DailyFullVaccination>().resize(mio::SimulationDay(size_t(1000)));
72-
model.parameters.get<mio::osecirvvs::DailyFullVaccination>().array().setConstant(3);
73-
69+
const size_t daily_vaccinations = 10;
70+
model.parameters.get<mio::osecirvvs::DailyFirstVaccination>().resize(mio::SimulationDay((size_t)tmax + 1));
71+
model.parameters.get<mio::osecirvvs::DailyFullVaccination>().resize(mio::SimulationDay((size_t)tmax + 1));
72+
for (size_t i = 0; i < tmax + 1; ++i) {
73+
auto num_vaccinations = static_cast<double>(i * daily_vaccinations);
74+
model.parameters.get<mio::osecirvvs::DailyFirstVaccination>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
75+
num_vaccinations;
76+
model.parameters.get<mio::osecirvvs::DailyFullVaccination>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
77+
num_vaccinations;
78+
}
7479
auto& contacts = model.parameters.get<mio::osecirvvs::ContactPatterns>();
7580
auto& contact_matrix = contacts.get_cont_freq_mat();
7681
contact_matrix[0].get_baseline().setConstant(0.5);

cpp/models/ode_secirvvs/parameters_io.h

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1225,6 +1225,24 @@ export_input_data_county_timeseries(std::vector<Model>&& model, const std::strin
12251225

12261226
return success();
12271227
}
1228+
#else
1229+
template <class Model>
1230+
IOResult<void> export_input_data_county_timeseries(std::vector<Model>&, const std::string&, std::vector<int> const&,
1231+
Date, const std::vector<double>&, double, int, const std::string&,
1232+
const std::string&, const std::string&)
1233+
{
1234+
mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data.");
1235+
return success();
1236+
}
1237+
1238+
template <class Model>
1239+
IOResult<void> export_input_data_county_timeseries(std::vector<Model>&&, const std::string&, std::vector<int> const&,
1240+
Date, const std::vector<double>&, double, int, const std::string&,
1241+
const std::string&, const std::string&, bool, const std::string&)
1242+
{
1243+
mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data.");
1244+
return success();
1245+
}
12281246

12291247
#endif //MEMILIO_HAS_HDF5
12301248

pycode/examples/simulation/2020_npis_sarscov2_wildtype_germany.py

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -177,8 +177,8 @@ def set_contact_matrices(self, model):
177177
minimum_file = os.path.join(
178178
self.data_dir, "contacts", "minimum_" + location + ".txt")
179179
contact_matrices[i] = mio.ContactMatrix(
180-
mio.secir.read_mobility_plain(baseline_file),
181-
mio.secir.read_mobility_plain(minimum_file)
180+
mio.read_mobility_plain(baseline_file),
181+
mio.read_mobility_plain(minimum_file)
182182
)
183183
model.parameters.ContactPatterns.cont_freq_mat = contact_matrices
184184

@@ -203,7 +203,7 @@ def set_npis(self, params, end_date):
203203

204204
typ_home = Intervention.Home.value
205205
typ_school = Intervention.SchoolClosure.value
206-
typ_home = Intervention.HomeOffice.value
206+
typ_homeoffice = Intervention.HomeOffice.value
207207
typ_gathering = Intervention.GatheringBanFacilitiesClosure.value
208208
typ_distance = Intervention.PhysicalDistanceAndMasks.value
209209
typ_senior = Intervention.SeniorAwareness.value
@@ -231,7 +231,7 @@ def school_closure(t, min, max):
231231

232232
def home_office(t, min, max):
233233
return damping_helper(
234-
t, min, max, lvl_main, typ_home, [loc_work])
234+
t, min, max, lvl_main, typ_homeoffice, [loc_work])
235235

236236
def social_events(t, min, max):
237237
return damping_helper(
Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
#############################################################################
2+
# Copyright (C) 2020-2024 MEmilio
3+
#
4+
# Authors: Martin J. Kuehn, Maximilian Betz
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+
import argparse
21+
import os
22+
from datetime import date, datetime
23+
24+
import matplotlib.pyplot as plt
25+
import numpy as np
26+
import pandas as pd
27+
28+
from memilio.simulation import AgeGroup, Damping, SimulationDay
29+
from memilio.simulation.osecirvvs import InfectionState
30+
from memilio.simulation.osecirvvs import Model, simulate
31+
32+
33+
def run_secirvvs_simulation(show_plot=True):
34+
"""
35+
Runs the c++ SECIRVVS model using a single age group
36+
and plots the results
37+
"""
38+
39+
t0 = 0
40+
tmax = 30 # number of days to simulate
41+
dt = 0.1
42+
num_groups = 1
43+
44+
# Initialize Parameters
45+
model = Model(num_groups)
46+
47+
# set parameters
48+
for i in range(num_groups):
49+
# Initial number of peaople in each compartment
50+
model.populations[AgeGroup(i), InfectionState.ExposedNaive] = 10
51+
model.populations[AgeGroup(
52+
i), InfectionState.ExposedImprovedImmunity] = 11
53+
model.populations[AgeGroup(
54+
i), InfectionState.ExposedPartialImmunity] = 12
55+
model.populations[AgeGroup(
56+
i), InfectionState.InfectedNoSymptomsNaive] = 13
57+
model.populations[AgeGroup(
58+
i), InfectionState.InfectedNoSymptomsNaiveConfirmed] = 13
59+
model.populations[AgeGroup(
60+
i), InfectionState.InfectedNoSymptomsPartialImmunity] = 14
61+
model.populations[AgeGroup(
62+
i), InfectionState.InfectedNoSymptomsPartialImmunityConfirmed] = 14
63+
model.populations[AgeGroup(
64+
i), InfectionState.InfectedNoSymptomsImprovedImmunity] = 15
65+
model.populations[AgeGroup(
66+
i), InfectionState.InfectedNoSymptomsImprovedImmunityConfirmed] = 15
67+
model.populations[AgeGroup(
68+
i), InfectionState.InfectedSymptomsNaive] = 5
69+
model.populations[AgeGroup(
70+
i), InfectionState.InfectedSymptomsNaiveConfirmed] = 5
71+
model.populations[AgeGroup(
72+
i), InfectionState.InfectedSymptomsPartialImmunity] = 6
73+
model.populations[AgeGroup(
74+
i), InfectionState.InfectedSymptomsPartialImmunityConfirmed] = 6
75+
model.populations[AgeGroup(
76+
i), InfectionState.InfectedSymptomsImprovedImmunity] = 7
77+
model.populations[AgeGroup(
78+
i), InfectionState.InfectedSymptomsImprovedImmunityConfirmed] = 7
79+
model.populations[AgeGroup(i), InfectionState.InfectedSevereNaive] = 8
80+
model.populations[AgeGroup(
81+
i), InfectionState.InfectedSevereImprovedImmunity] = 1
82+
model.populations[AgeGroup(
83+
i), InfectionState.InfectedSeverePartialImmunity] = 2
84+
model.populations[AgeGroup(
85+
i), InfectionState.InfectedCriticalNaive] = 3
86+
model.populations[AgeGroup(
87+
i), InfectionState.InfectedCriticalPartialImmunity] = 4
88+
model.populations[AgeGroup(
89+
i), InfectionState.InfectedCriticalImprovedImmunity] = 5
90+
model.populations[AgeGroup(
91+
i), InfectionState.SusceptibleImprovedImmunity] = 6
92+
model.populations[AgeGroup(
93+
i), InfectionState.SusceptiblePartialImmunity] = 7
94+
model.populations[AgeGroup(i), InfectionState.DeadNaive] = 0
95+
model.populations[AgeGroup(i), InfectionState.DeadPartialImmunity] = 0
96+
model.populations[AgeGroup(i), InfectionState.DeadImprovedImmunity] = 0
97+
model.populations.set_difference_from_group_total_AgeGroup(
98+
(AgeGroup(i), InfectionState.SusceptibleNaive), 1000)
99+
100+
model.parameters.ICUCapacity.value = 100
101+
model.parameters.TestAndTraceCapacity.value = 0.0143
102+
model.parameters.DailyFirstVaccination.resize_SimulationDay(
103+
SimulationDay(tmax + 1))
104+
model.parameters.DailyFullVaccination.resize_SimulationDay(
105+
SimulationDay(tmax + 1))
106+
daily_vaccinations = 10
107+
for i, num_vaccinations in enumerate(range(0, daily_vaccinations * (tmax + 1), daily_vaccinations)):
108+
model.parameters.DailyFirstVaccination[AgeGroup(
109+
0), SimulationDay(i)] = num_vaccinations
110+
model.parameters.DailyFullVaccination[AgeGroup(
111+
0), SimulationDay(i)] = num_vaccinations
112+
113+
# contact patterns
114+
baseline = np.ones((num_groups, num_groups)) * 0.5
115+
np.fill_diagonal(baseline, 5.0)
116+
model.parameters.ContactPatterns.cont_freq_mat[0].baseline = baseline
117+
model.parameters.ContactPatterns.cont_freq_mat.add_damping(Damping(
118+
coeffs=np.ones((num_groups, num_groups)) * 0.3, t=5.0, level=0, type=0))
119+
120+
# times
121+
model.parameters.TimeInfectedSymptoms[AgeGroup(0)] = 7
122+
model.parameters.TimeInfectedSevere[AgeGroup(0)] = 6
123+
model.parameters.TimeInfectedCritical[AgeGroup(0)] = 7
124+
125+
# probabilities
126+
model.parameters.TransmissionProbabilityOnContact[AgeGroup(0)] = 0.15
127+
model.parameters.RelativeTransmissionNoSymptoms[AgeGroup(0)] = 0.5
128+
# The precise value between Risk* (situation under control) and MaxRisk* (situation not under control)
129+
# depends on incidence and test and trace capacity
130+
model.parameters.RiskOfInfectionFromSymptomatic[AgeGroup(0)] = 0.0
131+
model.parameters.MaxRiskOfInfectionFromSymptomatic[AgeGroup(0)] = 0.4
132+
model.parameters.RecoveredPerInfectedNoSymptoms[AgeGroup(0)] = 0.2
133+
model.parameters.SeverePerInfectedSymptoms[AgeGroup(0)] = 0.1
134+
model.parameters.CriticalPerSevere[AgeGroup(0)] = 0.1
135+
model.parameters.DeathsPerCritical[AgeGroup(0)] = 0.1
136+
137+
model.parameters.ReducExposedPartialImmunity[AgeGroup(0)] = 0.8
138+
model.parameters.ReducExposedImprovedImmunity[AgeGroup(0)] = 0.331
139+
model.parameters.ReducInfectedSymptomsPartialImmunity[AgeGroup(0)] = 0.65
140+
model.parameters.ReducInfectedSymptomsImprovedImmunity[AgeGroup(0)] = 0.243
141+
model.parameters.ReducInfectedSevereCriticalDeadPartialImmunity[AgeGroup(
142+
0)] = 0.1
143+
model.parameters.ReducInfectedSevereCriticalDeadImprovedImmunity[AgeGroup(
144+
0)] = 0.091
145+
model.parameters.ReducTimeInfectedMild[AgeGroup(0)] = 0.9
146+
147+
model.parameters.Seasonality.value = 0.2
148+
149+
model.apply_constraints()
150+
151+
# Run Simulation
152+
result = simulate(t0, tmax, dt, model)
153+
154+
# # interpolate results
155+
# result = interpolate_simulation_result(result)
156+
157+
print(result.get_last_value())
158+
159+
160+
if __name__ == "__main__":
161+
arg_parser = argparse.ArgumentParser(
162+
'osecirvvs_simple',
163+
description='Simple example demonstrating the setup and simulation of the SECIRVVS model with a single age group.')
164+
args = arg_parser.parse_args()
165+
run_secirvvs_simulation(**args.__dict__)

pycode/memilio-simulation/CMakeLists.txt

Lines changed: 18 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -33,13 +33,6 @@ endif()
3333
add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/../../cpp ${CMAKE_CURRENT_BINARY_DIR}/cpp EXCLUDE_FROM_ALL)
3434

3535
# build python extensions
36-
pybind11_add_module(_simulation_secir MODULE
37-
memilio/simulation/secir.cpp
38-
)
39-
target_link_libraries(_simulation_secir PRIVATE memilio ode_secir)
40-
target_include_directories(_simulation_secir PRIVATE memilio/simulation)
41-
install(TARGETS _simulation_secir LIBRARY DESTINATION memilio)
42-
4336
pybind11_add_module(_simulation_abm MODULE
4437
memilio/simulation/abm.cpp
4538
)
@@ -62,17 +55,30 @@ target_link_libraries(_simulation PRIVATE memilio)
6255
target_include_directories(_simulation PRIVATE memilio/simulation)
6356
install(TARGETS _simulation LIBRARY DESTINATION memilio)
6457

58+
pybind11_add_module(_simulation_osir MODULE
59+
memilio/simulation/osir.cpp
60+
)
61+
target_link_libraries(_simulation_osir PRIVATE memilio ode_sir)
62+
target_include_directories(_simulation_osir PRIVATE memilio/simulation)
63+
install(TARGETS _simulation_osir LIBRARY DESTINATION memilio)
64+
6565
pybind11_add_module(_simulation_oseir MODULE
6666
memilio/simulation/oseir.cpp
6767
)
6868
target_link_libraries(_simulation_oseir PRIVATE memilio ode_seir)
6969
target_include_directories(_simulation_oseir PRIVATE memilio/simulation)
7070
install(TARGETS _simulation_oseir LIBRARY DESTINATION memilio)
7171

72-
pybind11_add_module(_simulation_osir MODULE
73-
memilio/simulation/osir.cpp
72+
pybind11_add_module(_simulation_secir MODULE
73+
memilio/simulation/secir.cpp
7474
)
75-
target_link_libraries(_simulation_osir PRIVATE memilio ode_sir)
76-
target_include_directories(_simulation_osir PRIVATE memilio/simulation)
77-
install(TARGETS _simulation_osir LIBRARY DESTINATION memilio)
75+
target_link_libraries(_simulation_secir PRIVATE memilio ode_secir)
76+
target_include_directories(_simulation_secir PRIVATE memilio/simulation)
77+
install(TARGETS _simulation_secir LIBRARY DESTINATION memilio)
7878

79+
pybind11_add_module(_simulation_osecirvvs MODULE
80+
memilio/simulation/osecirvvs.cpp
81+
)
82+
target_link_libraries(_simulation_osecirvvs PRIVATE memilio ode_secirvvs)
83+
target_include_directories(_simulation_osecirvvs PRIVATE memilio/simulation)
84+
install(TARGETS _simulation_osecirvvs LIBRARY DESTINATION memilio)

pycode/memilio-simulation/memilio/simulation/io/mobility_io.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ void bind_write_graph(pybind11::module_& m)
3434
{
3535
m.def("write_graph", [&](const mio::Graph<Model, mio::MigrationParameters>& graph, const std::string& directory) {
3636
int ioflags = mio::IOF_None;
37-
mio::write_graph<Model>(graph, directory, ioflags);
37+
auto ioresult = mio::write_graph<Model>(graph, directory, ioflags);
3838
});
3939
}
4040

pycode/memilio-simulation/memilio/simulation/io/result_io.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@ void bind_save_results(pybind11::module_& m)
3939
const std::vector<std::vector<Model>>& ensemble_params, const std::vector<int>& county_ids,
4040
const std::string& result_dir, bool save_single_runs, bool save_percentiles) {
4141
boost::filesystem::path dir(result_dir);
42-
mio::save_results<Model>(ensemble_results, ensemble_params, county_ids, dir, save_single_runs,
42+
auto ioresult = mio::save_results<Model>(ensemble_results, ensemble_params, county_ids, dir, save_single_runs,
4343
save_percentiles);
4444
return NULL;
4545
});

0 commit comments

Comments
 (0)