diff --git a/cpp/examples/ode_seir_metapop.cpp b/cpp/examples/ode_seir_metapop.cpp index f0132ea78f..606bd4eb2e 100644 --- a/cpp/examples/ode_seir_metapop.cpp +++ b/cpp/examples/ode_seir_metapop.cpp @@ -68,12 +68,10 @@ int main() for (size_t index = 0; index < (size_t)interpolated_result.get_num_time_points(); index++) { printf("\n %f", interpolated_result.get_times()[index]); for (size_t i = 0; i < (size_t)model.parameters.get_num_regions(); i++) { - - printf("\t %.5f ", - interpolated_result.get_value(index)[(size_t)model.parameters.get_num_regions() * - ((size_t)mio::oseirmetapop::InfectionState::Count - 1) + - i] / - model.populations.get_group_total(mio::regions::Region(i)) * 100); + auto infected_idx = model.populations.get_flat_index( + {mio::regions::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Infected}); + printf("\t %.5f ", interpolated_result.get_value(index)[infected_idx] / + model.populations.get_group_total(mio::regions::Region(i)) * 100); } } printf("\n"); diff --git a/cpp/memilio/math/stepper_wrapper.h b/cpp/memilio/math/stepper_wrapper.h index 4f8eef52ef..6260531ae2 100644 --- a/cpp/memilio/math/stepper_wrapper.h +++ b/cpp/memilio/math/stepper_wrapper.h @@ -53,8 +53,9 @@ struct step_adjuster : public boost::numeric::odeint::default_step_adjuster class ControlledStepper> +template + class ControlledStepper> class ControlledStepperWrapper : public mio::OdeIntegratorCore { using Algebra = boost::numeric::odeint::vector_space_algebra; @@ -204,8 +205,9 @@ class ControlledStepperWrapper : public mio::OdeIntegratorCore * @brief This is a non-adaptive IntegratorCore. It creates and manages an instance of an explicit stepper from * boost::numeric::odeint, wrapped as mio::IntegratorCore. */ -template class ExplicitStepper> +template + class ExplicitStepper> class ExplicitStepperWrapper : public mio::OdeIntegratorCore { public: diff --git a/cpp/models/ode_seir_metapop/infection_state.h b/cpp/models/ode_seir_metapop/infection_state.h index 65a570d9ad..37a089d750 100644 --- a/cpp/models/ode_seir_metapop/infection_state.h +++ b/cpp/models/ode_seir_metapop/infection_state.h @@ -17,7 +17,6 @@ * See the License for the specific language governing permissions and * limitations under the License. */ - #ifndef SEIRMETAPOP_INFECTIONSTATE_H #define SEIRMETAPOP_INFECTIONSTATE_H diff --git a/pycode/examples/simulation/ode_seir_metapop.py b/pycode/examples/simulation/ode_seir_metapop.py new file mode 100644 index 0000000000..94bcb89a8d --- /dev/null +++ b/pycode/examples/simulation/ode_seir_metapop.py @@ -0,0 +1,83 @@ +############################################################################# +# Copyright (C) 2020-2026 MEmilio +# +# Authors: Carlotta Gerstein +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# + +import numpy as np + +import memilio.simulation as mio +import memilio.simulation.oseir_metapop as oseir_metapop +from memilio.simulation import AgeGroup, LogLevel, set_log_level +set_log_level(LogLevel.Off) + +t0 = 0 +tmax = 10 +dt = 0.1 + +model = oseir_metapop.Model(3, 1) + +for i in range(model.parameters.num_regions.get()): + # Set infection state stay times (in days) + model.populations[mio.Region(i), AgeGroup( + 0), oseir_metapop.InfectionState.Susceptible] = 10000 + +model.populations[mio.Region(0), AgeGroup( + 0), oseir_metapop.InfectionState.Infected].value += 100 +model.populations[mio.Region(0), AgeGroup( + 0), oseir_metapop.InfectionState.Susceptible].value -= 100 + +mobility_data_commuter = np.array( + [[0.4, 0.3, 0.3], [0.2, 0.7, 0.1], [0.4, 0.1, 0.5]]) +model.set_commuting_strengths(mobility_data_commuter) + +# Set contact frequency +model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones( + (3, 3)) * 2.7 + +model.parameters.TimeExposed[mio.Region(0), AgeGroup(0)] = 3. +model.parameters.TimeExposed[mio.Region(1), AgeGroup(0)] = 4. +model.parameters.TimeExposed[mio.Region(2), AgeGroup(0)] = 5. +model.parameters.TimeInfected[mio.Region(0), AgeGroup(0)] = 7. +model.parameters.TimeInfected[mio.Region(1), AgeGroup(0)] = 8. +model.parameters.TimeInfected[mio.Region(2), AgeGroup(0)] = 9. +model.parameters.TransmissionProbabilityOnContact[mio.Region( + 0), AgeGroup(0)] = 0.07333 +model.parameters.TransmissionProbabilityOnContact[mio.Region( + 1), AgeGroup(0)] = 0.07333 +model.parameters.TransmissionProbabilityOnContact[mio.Region( + 2), AgeGroup(0)] = 0.07333 + +result = oseir_metapop.simulate(t0, tmax, dt, model) +interpolated_result = oseir_metapop.interpolate_simulation_result(result) + +print("Infected individuals per Region over time [%]:") +print("".join(f"\t Region {i}" for i in range( + model.parameters.num_regions.get()))) + +for t_idx in range(interpolated_result.get_num_time_points()): + row = f" {interpolated_result.get_times()[t_idx]:.6f}" + values = interpolated_result.get_value(t_idx) + + for i in range(model.parameters.num_regions.get()): + population = model.populations.get_group_total_Region(mio.Region(i)) + infected_idx = model.populations.get_flat_index( + (mio.Region(i), AgeGroup(0), oseir_metapop.InfectionState.Infected)) + percentage = values[infected_idx] / population * 100 + row += f"\t {percentage:.5f} " + + print(row) diff --git a/pycode/memilio-simulation/CMakeLists.txt b/pycode/memilio-simulation/CMakeLists.txt index 10204b0e00..c2ec86222e 100644 --- a/pycode/memilio-simulation/CMakeLists.txt +++ b/pycode/memilio-simulation/CMakeLists.txt @@ -110,6 +110,11 @@ add_pymio_module(_simulation_oseir SOURCES memilio/simulation/bindings/models/oseir.cpp ) +add_pymio_module(_simulation_oseir_metapop + LINKED_LIBRARIES memilio ode_seir_metapop + SOURCES memilio/simulation/bindings/models/oseir_metapop.cpp +) + add_pymio_module(_simulation_oseirdb LINKED_LIBRARIES memilio ode_seirdb SOURCES memilio/simulation/bindings/models/oseirdb.cpp diff --git a/pycode/memilio-simulation/memilio/simulation/__init__.py b/pycode/memilio-simulation/memilio/simulation/__init__.py index c21aeb66ff..56e1455e66 100644 --- a/pycode/memilio-simulation/memilio/simulation/__init__.py +++ b/pycode/memilio-simulation/memilio/simulation/__init__.py @@ -39,6 +39,9 @@ def __getattr__(attr): elif attr == "oseir": import memilio.simulation.oseir as oseir return oseir + elif attr == "oseir_metapop": + import memilio.simulation.oseir_metapop as oseir_metapop + return oseir_metapop elif attr == "oseirdb": import memilio.simulation.oseirdb as oseirdb return oseirdb diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/epidemiology/populations.h b/pycode/memilio-simulation/memilio/simulation/bindings/epidemiology/populations.h index 0b0a18e729..bd330d0e67 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/epidemiology/populations.h +++ b/pycode/memilio-simulation/memilio/simulation/bindings/epidemiology/populations.h @@ -24,6 +24,7 @@ #include "utils/custom_index_array.h" #include "memilio/utils/custom_index_array.h" #include "memilio/epidemiology/populations.h" +#include "geography/regions.h" #include "pybind11/pybind11.h" diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/geography/regions.h b/pycode/memilio-simulation/memilio/simulation/bindings/geography/regions.h new file mode 100644 index 0000000000..0d09f4f69c --- /dev/null +++ b/pycode/memilio-simulation/memilio/simulation/bindings/geography/regions.h @@ -0,0 +1,38 @@ +/* +* Copyright (C) 2020-2026 MEmilio +* +* Authors: Carlotta Gerstein +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#ifndef PYMIO_REGION_H +#define PYMIO_REGION_H + +#include "memilio/geography/regions.h" + +#include "pybind11/pybind11.h" + +namespace pymio +{ + +template <> +inline std::string pretty_name() +{ + return "Region"; +} + +} // namespace pymio + +#endif //PYMIO_REGION_H diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/oseir_metapop.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/oseir_metapop.cpp new file mode 100644 index 0000000000..65e801c840 --- /dev/null +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/oseir_metapop.cpp @@ -0,0 +1,108 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +// Includes from pymio +#include "pybind_util.h" +#include "compartments/simulation.h" +#include "compartments/flow_simulation.h" +#include "compartments/compartmental_model.h" +#include "data/analyze_result.h" +#include "epidemiology/age_group.h" +#include "epidemiology/populations.h" +#include "utils/custom_index_array.h" +#include "utils/parameter_set.h" +#include "utils/index.h" + +// Includes from MEmilio +#include "models/ode_seir_metapop/model.h" +#include "models/ode_seir_metapop/infection_state.h" +#include "models/ode_seir_metapop/parameters.h" +#include "memilio/compartments/simulation.h" +#include "memilio/data/analyze_result.h" + +#include "pybind11/pybind11.h" +#include "pybind11/eigen.h" + +namespace py = pybind11; + +namespace pymio +{ +// specialization of pretty_name +template <> +inline std::string pretty_name() +{ + return "InfectionState"; +} + +} // namespace pymio + +PYBIND11_MODULE(_simulation_oseir_metapop, m) +{ + pymio::bind_interpolate_result_methods(m); + + pymio::iterable_enum(m, "InfectionState", py::module_local{}) + .value("Susceptible", mio::oseirmetapop::InfectionState::Susceptible) + .value("Exposed", mio::oseirmetapop::InfectionState::Exposed) + .value("Infected", mio::oseirmetapop::InfectionState::Infected) + .value("Recovered", mio::oseirmetapop::InfectionState::Recovered); + + pymio::bind_ParameterSet, pymio::EnablePickling::Required>( + m, "ParametersBase"); + + pymio::bind_CustomIndexArray, mio::oseirmetapop::Region, mio::AgeGroup>( + m, "RegionAgeGroupArray"); + using RegionAgeGroupPopulations = mio::Populations; + pymio::bind_Population(m, "RegionAgeGroupPopulations", mio::Tag{}); + + pymio::bind_class, pymio::EnablePickling::Never, + mio::oseirmetapop::ParametersBase>(m, "Parameters") + .def(py::init(), py::arg("num_regions"), py::arg("num_agegroups")) + .def("check_constraints", &mio::oseirmetapop::Parameters::check_constraints) + .def("apply_constraints", &mio::oseirmetapop::Parameters::apply_constraints) + .def_property_readonly("num_regions", &mio::oseirmetapop::Parameters::get_num_regions) + .def_property_readonly("num_agegroups", &mio::oseirmetapop::Parameters::get_num_agegroups); + + using Populations = + mio::Populations; + + pymio::bind_Population(m, "Populations", mio::Tag::Populations>{}); + pymio::bind_CompartmentalModel, pymio::EnablePickling::Never>(m, "ModelBase"); + + using Model = mio::oseirmetapop::Model; + pymio::bind_class>>(m, "Model") + .def(py::init(), py::arg("num_regions"), py::arg("num_agegroups")) + .def("set_commuting_strengths", py::overload_cast(&Model::set_commuting_strengths), + py::arg("commuting_strengths")) + .def("set_commuting_strengths_identity", py::overload_cast<>(&Model::set_commuting_strengths)); + + pymio::bind_Simulation>(m, "Simulation"); + pymio::bind_Flow_Simulation>(m, "FlowSimulation"); + + m.def("simulate", &mio::simulate, "Simulates an ODE SEIR metapopulation model from t0 to tmax.", + py::arg("t0"), py::arg("tmax"), py::arg("dt"), py::arg("model"), py::arg("integrator") = py::none()); + m.def("simulate_flows", &mio::simulate_flows, + "Simulates an ODE SEIR metapopulation model with flows from t0 to tmax.", py::arg("t0"), py::arg("tmax"), + py::arg("dt"), py::arg("model"), py::arg("integrator") = py::none()); + + m.attr("__version__") = "dev"; +} diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/pybind_util.h b/pycode/memilio-simulation/memilio/simulation/bindings/pybind_util.h index dc48195817..822ded1839 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/pybind_util.h +++ b/pycode/memilio-simulation/memilio/simulation/bindings/pybind_util.h @@ -255,7 +255,7 @@ auto bind_Range(pybind11::module_& m, const std::string& class_name) .def( "__iter__", [](Range& self) { - return self; + return Iterator{{self.begin(), self.end()}}; }, pybind11::keep_alive<1, 0>{}) //keep alive the Range as long as there is an iterator .def( diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/simulation.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/simulation.cpp index 5bb3793ce1..297dcadbb2 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/simulation.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/simulation.cpp @@ -59,6 +59,9 @@ PYBIND11_MODULE(_simulation, m) pymio::bind_CustomIndexArray, mio::AgeGroup>(m, "AgeGroupArray"); pymio::bind_class>(m, "AgeGroup") .def(py::init()); + pymio::bind_Index(m, "IndexRegion"); + pymio::bind_class>(m, "Region") + .def(py::init()); pymio::bind_CustomIndexArray(m, "AgeGroupSimulationDayArray"); pymio::bind_class>(m, diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/utils/index.h b/pycode/memilio-simulation/memilio/simulation/bindings/utils/index.h index cbaf6f8315..3d497cdfc2 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/utils/index.h +++ b/pycode/memilio-simulation/memilio/simulation/bindings/utils/index.h @@ -51,6 +51,7 @@ void bind_Index(pybind11::module_& m, std::string const& name) c.def(pybind11::init(), pybind11::arg("value")); c.def(pybind11::self == pybind11::self); c.def(pybind11::self != pybind11::self); + c.def("get", &mio::Index::get); bind_Index_members_if_enum(c); } diff --git a/pycode/memilio-simulation/memilio/simulation/oseir_metapop.py b/pycode/memilio-simulation/memilio/simulation/oseir_metapop.py new file mode 100644 index 0000000000..d713e8a551 --- /dev/null +++ b/pycode/memilio-simulation/memilio/simulation/oseir_metapop.py @@ -0,0 +1,25 @@ +############################################################################# +# Copyright (C) 2020-2025 MEmilio +# +# Authors: Henrik Zunker +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# + +""" +Python bindings for the MEmilio ODE SEIR metapopulation model. +""" + +from memilio.simulation._simulation_oseir_metapop import * diff --git a/pycode/memilio-simulation/tests/test_oseir_metapop.py b/pycode/memilio-simulation/tests/test_oseir_metapop.py new file mode 100644 index 0000000000..5ff6a6273a --- /dev/null +++ b/pycode/memilio-simulation/tests/test_oseir_metapop.py @@ -0,0 +1,83 @@ +############################################################################# +# Copyright (C) 2020-2026 MEmilio +# +# Authors: Henrik Zunker +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# +import unittest + +import numpy as np + +from memilio.simulation import AgeGroup, Region +from memilio.simulation.oseir_metapop import ( + FlowSimulation, + InfectionState as State, + Model, + simulate_flows, +) + + +class Test_oseir_metapop_integration(unittest.TestCase): + """Python binding tests for ODE SEIR metapopulation model.""" + + def setUp(self): + model = Model(3, 1) + A0 = AgeGroup(0) + + for region in range(3): + model.populations[Region(region), A0, State.Susceptible] = 1000.0 + + model.populations[Region(0), A0, State.Susceptible] = 990.0 + model.populations[Region(0), A0, State.Infected] = 10.0 + + model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones( + (1, 1)) * 2.7 + model.set_commuting_strengths_identity() + + self.model = model + + def test_population_after_commuting_is_bound(self): + A0 = AgeGroup(0) + population_after_commuting = self.model.parameters.PopulationAfterCommuting + + self.assertAlmostEqual( + population_after_commuting[Region(0), A0].value, 1000.0) + self.assertAlmostEqual( + population_after_commuting.get_group_total_Region(Region(0)), + 1000.0) + + def test_simulate_flows_simple(self): + flow_sim_results = simulate_flows( + t0=0.0, tmax=1.0, dt=0.1, model=self.model) + compartments = flow_sim_results[0] + flows = flow_sim_results[1] + + self.assertAlmostEqual(compartments.get_time(0), 0.0) + self.assertAlmostEqual(compartments.get_last_time(), 1.0) + self.assertEqual(len(compartments.get_last_value()), 12) + self.assertEqual(len(flows.get_last_value()), 9) + + def test_flow_simulation_simple(self): + sim = FlowSimulation(self.model, t0=0.0, dt=0.1) + sim.advance(tmax=1.0) + + self.assertAlmostEqual(sim.result.get_last_time(), 1.0) + self.assertEqual(len(sim.result.get_last_value()), 12) + self.assertEqual(len(sim.flows.get_last_value()), 9) + + +if __name__ == '__main__': + unittest.main()