Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
74 changes: 74 additions & 0 deletions pycode/examples/simulation/sde_sir_simple.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
#############################################################################
# Copyright (C) 2020-2025 MEmilio
#
# Authors: Maximilian Betz
#
# Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
#
# 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 argparse

import numpy as np

from memilio.simulation import AgeGroup, Damping
from memilio.simulation.ssir import InfectionState as State
from memilio.simulation.ssir import (Model, simulate_stochastic)


def run_sde_sir_simulation():
"""Runs the c++ SDE SIR model"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"""Runs the c++ SDE SIR model"""
"""Runs SDE SIR model"""


# Define population of age groups
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand the comment. This is the total population, right? Does the SDE model even have age groups?

population = 10000

days = 5. # number of days to simulate
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
days = 5. # number of days to simulate
tmax = 5. # simulation time frame

dt = 0.1

# Initialize Parameters
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# Initialize Parameters
# Initialize model

model = Model()

# Compartment transition duration
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# Compartment transition duration
# Mean time in Infected compartment

model.parameters.TimeInfected.value = 10.

# # Compartment transition propabilities
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# # Compartment transition propabilities

model.parameters.TransmissionProbabilityOnContact.value = 1.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Transmission probability of 100%?


# Initial number of people in each compartment
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# Initial number of people in each compartment
# Initial number of people per compartment

model.populations[State.Infected] = 100
model.populations[State.Recovered] = 1000
model.populations.set_difference_from_total(
(State.Susceptible), population)

model.parameters.ContactPatterns.baseline = np.ones(
(1, 1)) * 2.7
model.parameters.ContactPatterns.minimum = np.zeros(
(1, 1))
model.parameters.ContactPatterns.add_damping(
Damping(coeffs=np.r_[0.6], t=12.5, level=0, type=0))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You only simulate 5 days but add a damping after 12.5 days


# Check logical constraints to parameters
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# Check logical constraints to parameters
# Check parameter constraints

model.check_constraints()

# Run Simulation
result = simulate_stochastic(0., days, dt, model)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
result = simulate_stochastic(0., days, dt, model)
result = simulate_stochastic(0., tmax, dt, model)


result.print_table(False, ["S", "I", "R"], 16, 5)


if __name__ == "__main__":
arg_parser = argparse.ArgumentParser(
'sde_sir_simple',
description='Simple example demonstrating the setup and simulation of the SDE SIR model.')
args = arg_parser.parse_args()
run_sde_sir_simulation()
80 changes: 80 additions & 0 deletions pycode/examples/simulation/sde_sirs_simple.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
#############################################################################
# Copyright (C) 2020-2025 MEmilio
#
# Authors: Maximilian Betz
#
# Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
#
# 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 argparse

import numpy as np

from memilio.simulation import AgeGroup, Damping
from memilio.simulation.ssirs import InfectionState as State
from memilio.simulation.ssirs import (
Model, simulate_stochastic, interpolate_simulation_result)


def run_sde_sirs_simulation():
"""Runs the c++ SDE SIRs model"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"""Runs the c++ SDE SIRs model"""
"""Runs SDE SIRS model"""


# Define population of age groups
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See comment above

population = 10000

days = 5. # number of days to simulate
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
days = 5. # number of days to simulate
tmax = 5. # simulation time frame

dt = 0.001

# Initialize Parameters
model = Model()

# Compartment transition duration
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please adjust the comments according to my review of SIR example

model.parameters.TimeInfected.value = 10.
model.parameters.TimeImmune.value = 100.

# # Compartment transition propabilities
model.parameters.TransmissionProbabilityOnContact.value = 1.

# Initial number of people in each compartment
model.populations[State.Infected] = 100
model.populations[State.Recovered] = 1000
model.populations.set_difference_from_total(
(State.Susceptible), population)

model.parameters.ContactPatterns.baseline = np.ones(
(1, 1)) * 20.7
model.parameters.ContactPatterns.minimum = np.zeros(
(1, 1))
model.parameters.ContactPatterns.add_damping(
Damping(coeffs=np.r_[0.6], t=12.5, level=0, type=0))

# Check logical constraints to parameters
model.check_constraints()

# Run Simulation
result = simulate_stochastic(0., days, dt, model)

# Interpolate results
result = interpolate_simulation_result(result)

# Print results
result.print_table(False, ["Susceptible", "Infected", "Recovered"], 16, 5)


if __name__ == "__main__":
arg_parser = argparse.ArgumentParser(
'sde_sirs_simple',
description='Simple example demonstrating the setup and simulation of the SDE SIRS model.')
args = arg_parser.parse_args()
run_sde_sirs_simulation()
10 changes: 10 additions & 0 deletions pycode/memilio-simulation/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -88,3 +88,13 @@ add_pymio_module(_simulation_osecirvvs
LINKED_LIBRARIES memilio ode_secirvvs
SOURCES memilio/simulation/bindings/models/osecirvvs.cpp
)

add_pymio_module(_simulation_ssir
LINKED_LIBRARIES memilio sde_sir
SOURCES memilio/simulation/bindings/models/ssir.cpp
)

add_pymio_module(_simulation_ssirs
LINKED_LIBRARIES memilio sde_sirs
SOURCES memilio/simulation/bindings/models/ssirs.cpp
)
7 changes: 7 additions & 0 deletions pycode/memilio-simulation/memilio/simulation/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,13 @@ def __getattr__(attr):
elif attr == "osecirvvs":
import memilio.simulation.osecirvvs as osecirvvs
return osecirvvs
elif attr == "ssir":
import memilio.simulation.ssir as ssir
return ssir

elif attr == "ssirs":
import memilio.simulation.ssirs as ssirs
return ssirs

raise AttributeError("module {!r} has no attribute "
"{!r}".format(__name__, attr))
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
/*
* Copyright (C) 2020-2025 MEmilio
*
* Authors: Maximilian Betz
*
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
*
* 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 "utils/index.h"
#include "compartments/compartmental_model.h"
#include "utils/custom_index_array.h"
#include "utils/parameter_set.h"
#include "epidemiology/populations.h"

//Includes from MEmilio
#include "sde_sir/model.h"
#include "sde_sir/infection_state.h"
#include "memilio/compartments/stochastic_simulation.h"
#include "memilio/compartments/stochastic_model.h"
#include "memilio/data/analyze_result.h"

#include "pybind11/pybind11.h"
#include <pybind11/stl.h>
#include <pybind11/functional.h>

namespace py = pybind11;

namespace pymio
{
//specialization of pretty_name
template <>
inline std::string pretty_name<mio::ssir::InfectionState>()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is that used for?

{
return "InfectionState";
}

} // namespace pymio

PYBIND11_MODULE(_simulation_ssir, m)
{
m.def("interpolate_simulation_result",
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
&mio::interpolate_simulation_result),
py::arg("ts"), py::arg("abs_tol") = 1e-14);

m.def("interpolate_simulation_result",
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
&mio::interpolate_simulation_result),
py::arg("ts"), py::arg("interpolation_times"));

m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);

pymio::iterable_enum<mio::ssir::InfectionState>(m, "InfectionState")
.value("Susceptible", mio::ssir::InfectionState::Susceptible)
.value("Infected", mio::ssir::InfectionState::Infected)
.value("Recovered", mio::ssir::InfectionState::Recovered);

pymio::bind_ParameterSet<mio::ssir::ParametersBase, pymio::EnablePickling::Never>(m, "ParametersBase");

pymio::bind_class<mio::ssir::Parameters, pymio::EnablePickling::Required, mio::ssir::ParametersBase>(m,
"Parameters")
.def(py::init<>())
.def("check_constraints", &mio::ssir::Parameters::check_constraints)
.def("apply_constraints", &mio::ssir::Parameters::apply_constraints);

using Populations = mio::Populations<double, mio::ssir::InfectionState>;
pymio::bind_Population(m, "Populations", mio::Tag<mio::ssir::Model::Populations>{});
pymio::bind_CompartmentalModel<mio::ssir::InfectionState, Populations, mio::ssir::Parameters,
pymio::EnablePickling::Never>(m, "ModelBase");
pymio::bind_class<mio::ssir::Model, pymio::EnablePickling::Never,
mio::CompartmentalModel<double, mio::ssir::InfectionState, Populations, mio::ssir::Parameters>>(
m, "Model")
.def(py::init<>());

m.def(
"simulate_stochastic",
[](double t0, double tmax, double dt, mio::ssir::Model const& model) {
return mio::simulate_stochastic<double, mio::ssir::Model>(t0, tmax, dt, model);
},
"Simulates an SDE SIR model from t0 to tmax.", py::arg("t0"), py::arg("tmax"), py::arg("dt"), py::arg("model"));

m.attr("__version__") = "dev";
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
/*
* Copyright (C) 2020-2025 MEmilio
*
* Authors: Maximilian Betz
*
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
*
* 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 "utils/index.h"
#include "compartments/compartmental_model.h"
#include "utils/custom_index_array.h"
#include "utils/parameter_set.h"
#include "epidemiology/populations.h"

//Includes from MEmilio
#include "sde_sirs/model.h"
#include "sde_sirs/infection_state.h"
#include "memilio/compartments/stochastic_simulation.h"
#include "memilio/compartments/stochastic_model.h"
#include "memilio/data/analyze_result.h"

#include "pybind11/pybind11.h"
#include <pybind11/stl.h>
#include <pybind11/functional.h>

namespace py = pybind11;

namespace pymio
{
//specialization of pretty_name
template <>
inline std::string pretty_name<mio::ssirs::InfectionState>()
{
return "InfectionState";
}

} // namespace pymio

PYBIND11_MODULE(_simulation_ssirs, m)
{
m.def("interpolate_simulation_result",
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
&mio::interpolate_simulation_result),
py::arg("ts"), py::arg("abs_tol") = 1e-14);

m.def("interpolate_simulation_result",
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
&mio::interpolate_simulation_result),
py::arg("ts"), py::arg("interpolation_times"));

m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);

pymio::iterable_enum<mio::ssirs::InfectionState>(m, "InfectionState")
.value("Susceptible", mio::ssirs::InfectionState::Susceptible)
.value("Infected", mio::ssirs::InfectionState::Infected)
.value("Recovered", mio::ssirs::InfectionState::Recovered);

pymio::bind_ParameterSet<mio::ssirs::ParametersBase, pymio::EnablePickling::Never>(m, "ParametersBase");

pymio::bind_class<mio::ssirs::Parameters, pymio::EnablePickling::Required, mio::ssirs::ParametersBase>(m,
"Parameters")
.def(py::init<>())
.def("check_constraints", &mio::ssirs::Parameters::check_constraints)
.def("apply_constraints", &mio::ssirs::Parameters::apply_constraints);

using Populations = mio::Populations<double, mio::ssirs::InfectionState>;
pymio::bind_Population(m, "Populations", mio::Tag<mio::ssirs::Model::Populations>{});
pymio::bind_CompartmentalModel<mio::ssirs::InfectionState, Populations, mio::ssirs::Parameters,
pymio::EnablePickling::Never>(m, "ModelBase");
pymio::bind_class<mio::ssirs::Model, pymio::EnablePickling::Never,
mio::CompartmentalModel<double, mio::ssirs::InfectionState, Populations, mio::ssirs::Parameters>>(
m, "Model")
.def(py::init<>());

m.def(
"simulate_stochastic",
[](double t0, double tmax, double dt, mio::ssirs::Model const& model) {
return mio::simulate_stochastic<double, mio::ssirs::Model>(t0, tmax, dt, model);
},
"Simulates an SDE SIR model from t0 to tmax.", py::arg("t0"), py::arg("tmax"), py::arg("dt"), py::arg("model"));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"Simulates an SDE SIR model from t0 to tmax.", py::arg("t0"), py::arg("tmax"), py::arg("dt"), py::arg("model"));
"Simulates an SDE SIRS model from t0 to tmax.", py::arg("t0"), py::arg("tmax"), py::arg("dt"), py::arg("model"));


m.attr("__version__") = "dev";
}
Loading