diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 0a1744deab..3ce345b3fd 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -131,6 +131,7 @@ if(MEMILIO_ENABLE_WARNINGS) "-Wall;-Wextra;-Wshadow;--pedantic;") endif() endif() + # exclude some warnings we accept if(CMAKE_CXX_COMPILER_ID MATCHES "GNU") string(APPEND MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS @@ -139,13 +140,14 @@ elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang") string(APPEND MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS "-Wno-unknown-warning-option;-Wno-deprecated;-Wno-gnu-zero-variadic-macro-arguments;") endif() + # woyrkarounds for compiler bugs or overzealous optimization/analysis if(CMAKE_BUILD_TYPE STREQUAL "Release" OR CMAKE_BUILD_TYPE STREQUAL "RELEASE") - if (CMAKE_CXX_COMPILER_ID MATCHES "GNU") + if(CMAKE_CXX_COMPILER_ID MATCHES "GNU") string(APPEND MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS - "-Wno-restrict;" # https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105329 - "-Wno-maybe-uninitialized;" # https://gcc.gnu.org/bugzilla/buglist.cgi?quicksearch=may%20be%20uninitialized - "-Wno-stringop-overread;" + "-Wno-restrict;" # https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105329 + "-Wno-maybe-uninitialized;" # https://gcc.gnu.org/bugzilla/buglist.cgi?quicksearch=may%20be%20uninitialized + "-Wno-stringop-overread;" ) elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang") string(APPEND MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS @@ -153,6 +155,7 @@ if(CMAKE_BUILD_TYPE STREQUAL "Release" OR CMAKE_BUILD_TYPE STREQUAL "RELEASE") ) endif() endif() + # finalize string by setting warnings as errors if(MEMILIO_ENABLE_WARNINGS_AS_ERRORS) if(CMAKE_CXX_COMPILER_ID MATCHES "MSVC") @@ -164,8 +167,6 @@ if(MEMILIO_ENABLE_WARNINGS_AS_ERRORS) endif() endif() - - # add parts of the project include(thirdparty/CMakeLists.txt) add_subdirectory(memilio/ad) @@ -190,6 +191,7 @@ if(MEMILIO_BUILD_MODELS) add_subdirectory(models/graph_abm) add_subdirectory(models/smm) add_subdirectory(models/hybrid) + add_subdirectory(models/ode_mseirs4) endif() if(MEMILIO_BUILD_EXAMPLES) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 980164d324..39a12539e1 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -196,3 +196,7 @@ if(MEMILIO_HAS_JSONCPP) target_link_libraries(ide_initialization_example PRIVATE memilio ide_secir) target_compile_options(ide_initialization_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) endif() + +add_executable(ode_mseirs4_example ode_mseirs4.cpp) +target_link_libraries(ode_mseirs4_example PRIVATE memilio ode_mseirs4) +target_compile_options(ode_mseirs4_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/examples/ode_mseirs4.cpp b/cpp/examples/ode_mseirs4.cpp new file mode 100644 index 0000000000..82cfbdb236 --- /dev/null +++ b/cpp/examples/ode_mseirs4.cpp @@ -0,0 +1,94 @@ +/* +* 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. +*/ +#include "ode_mseirs4/model.h" +#include "memilio/compartments/simulation.h" +#include + +int main() +{ + using FP = double; + mio::omseirs4::Model model; + auto& params = model.parameters; + + // Example parameter values for day-based time unit (t in days) + params.get>() = 0.4; // b0 per day + params.get>() = 0.15; // b1 + params.get>() = 0.0; // phi; for phase shift use 2*pi*offsetDays/365 + params.get>() = 1.0 / (70.0 * 365.0); // mu per day + params.get>() = 1.0 / 90.0; // xi per day (~3 months) + params.get>() = 1.0 / 7.0; // sigma per day (≈7 days latent) + params.get>() = 1.0 / 14.0; // nu per day (≈14 days infectious) + params.get>() = 1.0 / (5.0 * 365.0); // gamma per day (5 years) + // factors default already set (0.5, 0.35, 0.25) as in the paper https://doi.org/10.1016/S0025-5564(01)00066-9. + + // Initial population (absolute counts), set all compartments explicitly. + double N = 1e6; + + // Initial populations. + model.populations[{mio::Index(mio::omseirs4::InfectionState::MaternalImmune)}] = + 5000.0; + + model.populations[{mio::Index(mio::omseirs4::InfectionState::E1)}] = 300.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::E2)}] = 150.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::E3)}] = 80.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::E4)}] = 70.0; + + model.populations[{mio::Index(mio::omseirs4::InfectionState::I1)}] = 200.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::I2)}] = 100.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::I3)}] = 50.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::I4)}] = 50.0; + + model.populations[{mio::Index(mio::omseirs4::InfectionState::R1)}] = 40000.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::R2)}] = 30000.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::R3)}] = 20000.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::R4)}] = 10000.0; + + model.populations[{mio::Index(mio::omseirs4::InfectionState::S2)}] = 100000.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::S3)}] = 50000.0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::S4)}] = 50000.0; + + // Compute S1 as residual to match N + double assigned = 5000.0 + (300.0 + 150.0 + 80.0 + 70.0) + (200.0 + 100.0 + 50.0 + 50.0) + + (40000.0 + 30000.0 + 20000.0 + 10000.0) + (100000.0 + 50000.0 + 50000.0); + double S1 = N - assigned; + if (S1 < 0) + S1 = 0; + model.populations[{mio::Index(mio::omseirs4::InfectionState::S1)}] = S1; + + model.check_constraints(); + + // simulate + double t0 = 0.0; + double tmax = 20.0; // days + double dt = 1.0; // daily output + auto result = mio::simulate(t0, tmax, dt, model); + + // print header + std::cout << "t M S1 S2 S3 S4 E1 E2 E3 E4 I1 I2 I3 I4 R1 R2 R3 R4\n"; + for (size_t i = 0; i < (size_t)result.get_num_time_points(); ++i) { + std::cout << result.get_time(i); + const auto& y = result.get_value(i); + for (size_t k = 0; k < (size_t)mio::omseirs4::InfectionState::Count; ++k) { + std::cout << ' ' << y[(Eigen::Index)k]; + } + std::cout << '\n'; + } + return 0; +} diff --git a/cpp/models/ode_mseirs4/CMakeLists.txt b/cpp/models/ode_mseirs4/CMakeLists.txt new file mode 100644 index 0000000000..c93e61e9cb --- /dev/null +++ b/cpp/models/ode_mseirs4/CMakeLists.txt @@ -0,0 +1,12 @@ +add_library(ode_mseirs4 + infection_state.h + parameters.h + model.h + model.cpp +) +target_link_libraries(ode_mseirs4 PUBLIC memilio) +target_include_directories(ode_mseirs4 PUBLIC + $ + $ +) +target_compile_options(ode_mseirs4 PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ode_mseirs4/README.md b/cpp/models/ode_mseirs4/README.md new file mode 100644 index 0000000000..866e35a8ab --- /dev/null +++ b/cpp/models/ode_mseirs4/README.md @@ -0,0 +1,18 @@ +# ODE MSEIRS4 Model + +Implementation of the MSEIRS4 model (maternal immunity + susceptible stages S1-S4 + exposed, infectious, recovered stages with four parallel infection/recovery classes). +This model is designed for Respiratory Syncytial Virus (RSV) and is based on: + +- Weber A, Weber M, Milligan P. (2001). *Modeling epidemics caused by respiratory syncytial virus (RSV).* Mathematical Biosciences 172(2): 95–113. `DOI:10.1016/S0025-5564(01)00066-9 `_ + + +## Meaning of indices 1–4 (S1–S4, E1–E4, I1–I4, R1–R4) + +- S1: fully susceptible after loss of maternal immunity (highest susceptibility). +- S2: susceptible after first infection (R1 → S2 via waning; reduced susceptibility). +- S3: susceptible after second infection (R2 → S3). +- S4: susceptible after ≥3 infections (R3 → S4 and R4 → S4; lowest susceptibility). + +Correspondingly, E_k/I_k/R_k represent exposed/infectious/recovered states of the k-th infection episode. All infectious classes (I1..I4) contribute to transmission equally in the basic formulation of the paper. + +The model includes seasonality by adjusting the transmission rate with a cosine function. diff --git a/cpp/models/ode_mseirs4/infection_state.h b/cpp/models/ode_mseirs4/infection_state.h new file mode 100644 index 0000000000..eb6aa8ddbd --- /dev/null +++ b/cpp/models/ode_mseirs4/infection_state.h @@ -0,0 +1,52 @@ +/* +* 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. +*/ +#ifndef ODE_MSEIRS4_INFECTION_STATE_H +#define ODE_MSEIRS4_INFECTION_STATE_H + +namespace mio +{ +namespace omseirs4 +{ +enum class InfectionState +{ + MaternalImmune = 0, + S1, + S2, + S3, + S4, + E1, + E2, + E3, + E4, + I1, + I2, + I3, + I4, + R1, + R2, + R3, + R4, + Count +}; + +} // namespace omseirs4 +} // namespace mio + +#endif // ODE_MSEIRS4_INFECTION_STATE_H diff --git a/cpp/models/ode_mseirs4/model.cpp b/cpp/models/ode_mseirs4/model.cpp new file mode 100644 index 0000000000..a67c446b4f --- /dev/null +++ b/cpp/models/ode_mseirs4/model.cpp @@ -0,0 +1,28 @@ +/* +* 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. +*/ +#include "ode_mseirs4/model.h" + +namespace mio +{ +namespace omseirs4 +{ + +} // namespace omseirs4 +} // namespace mio diff --git a/cpp/models/ode_mseirs4/model.h b/cpp/models/ode_mseirs4/model.h new file mode 100644 index 0000000000..07a95b3100 --- /dev/null +++ b/cpp/models/ode_mseirs4/model.h @@ -0,0 +1,158 @@ +/* +* 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. +*/ +#ifndef ODE_MSEIRS4_MODEL_H +#define ODE_MSEIRS4_MODEL_H + +#include "memilio/compartments/compartmental_model.h" +#include "memilio/epidemiology/populations.h" +#include "ode_mseirs4/infection_state.h" +#include "ode_mseirs4/parameters.h" + +#include + +namespace mio +{ +namespace omseirs4 +{ + +template +class Model : public mio::CompartmentalModel, Parameters> +{ + using Base = mio::CompartmentalModel, Parameters>; + +public: + using typename Base::ParameterSet; + using typename Base::Populations; + + Model(const Populations& pop, const ParameterSet& params) + : Base(pop, params) + { + } + + Model() + : Base(Populations({InfectionState::Count}, 0.), ParameterSet()) + { + } + + void get_derivatives(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> dydt) const override + { + auto& params = this->parameters; + const FP b0 = params.template get>(); + const FP b1 = params.template get>(); + const FP phi = params.template get>(); + const FP mu = params.template get>(); + const FP xi = params.template get>(); + const FP sigma = params.template get>(); + const FP nu = params.template get>(); + const FP gamma = params.template get>(); + const FP f2 = params.template get>(); + const FP f3 = params.template get>(); + const FP f4 = params.template get>(); + + const FP two_pi = FP(2) * std::numbers::pi_v; + const FP beta1 = b0 * (FP(1) + b1 * std::cos(two_pi * (t / FP(365.0)) + phi)); + const FP beta2 = f2 * beta1; + const FP beta3 = f3 * beta1; + const FP beta4 = f4 * beta1; + + auto idx = [&](InfectionState s) { + return static_cast(s); + }; + + FP I_total = y[idx(InfectionState::I1)] + y[idx(InfectionState::I2)] + y[idx(InfectionState::I3)] + + y[idx(InfectionState::I4)]; + FP R_total = y[idx(InfectionState::R1)] + y[idx(InfectionState::R2)] + y[idx(InfectionState::R3)] + + y[idx(InfectionState::R4)]; + FP N = pop.sum(); + + // dM + dydt[idx(InfectionState::MaternalImmune)] = mu * R_total - (xi + mu) * y[idx(InfectionState::MaternalImmune)]; + + // dS1 + dydt[idx(InfectionState::S1)] = mu * (N - R_total) + xi * y[idx(InfectionState::MaternalImmune)] - + mu * y[idx(InfectionState::S1)] - beta1 * I_total * y[idx(InfectionState::S1)]; + + // dE1..E4 + dydt[idx(InfectionState::E1)] = + beta1 * I_total * y[idx(InfectionState::S1)] - (mu + sigma) * y[idx(InfectionState::E1)]; + dydt[idx(InfectionState::E2)] = + beta2 * I_total * y[idx(InfectionState::S2)] - (mu + sigma) * y[idx(InfectionState::E2)]; + dydt[idx(InfectionState::E3)] = + beta3 * I_total * y[idx(InfectionState::S3)] - (mu + sigma) * y[idx(InfectionState::E3)]; + dydt[idx(InfectionState::E4)] = + beta4 * I_total * y[idx(InfectionState::S4)] - (mu + sigma) * y[idx(InfectionState::E4)]; + + // dI1..I4 + dydt[idx(InfectionState::I1)] = sigma * y[idx(InfectionState::E1)] - (nu + mu) * y[idx(InfectionState::I1)]; + dydt[idx(InfectionState::I2)] = sigma * y[idx(InfectionState::E2)] - (nu + mu) * y[idx(InfectionState::I2)]; + dydt[idx(InfectionState::I3)] = sigma * y[idx(InfectionState::E3)] - (nu + mu) * y[idx(InfectionState::I3)]; + dydt[idx(InfectionState::I4)] = sigma * y[idx(InfectionState::E4)] - (nu + mu) * y[idx(InfectionState::I4)]; + + // dR1..R4 + dydt[idx(InfectionState::R1)] = nu * y[idx(InfectionState::I1)] - (mu + gamma) * y[idx(InfectionState::R1)]; + dydt[idx(InfectionState::R2)] = nu * y[idx(InfectionState::I2)] - (mu + gamma) * y[idx(InfectionState::R2)]; + dydt[idx(InfectionState::R3)] = nu * y[idx(InfectionState::I3)] - (mu + gamma) * y[idx(InfectionState::R3)]; + dydt[idx(InfectionState::R4)] = nu * y[idx(InfectionState::I4)] - (mu + gamma) * y[idx(InfectionState::R4)]; + + // dS2,S3,S4 + dydt[idx(InfectionState::S2)] = gamma * y[idx(InfectionState::R1)] - mu * y[idx(InfectionState::S2)] - + beta2 * I_total * y[idx(InfectionState::S2)]; + dydt[idx(InfectionState::S3)] = gamma * y[idx(InfectionState::R2)] - mu * y[idx(InfectionState::S3)] - + beta3 * I_total * y[idx(InfectionState::S3)]; + dydt[idx(InfectionState::S4)] = gamma * (y[idx(InfectionState::R3)] + y[idx(InfectionState::R4)]) - + mu * y[idx(InfectionState::S4)] - beta4 * I_total * y[idx(InfectionState::S4)]; + } + + /** + * serialize this. + * @see mio::serialize + */ + template + void serialize(IOContext& io) const + { + auto obj = io.create_object("Model"); + obj.add_element("Parameters", this->parameters); + obj.add_element("Populations", this->populations); + } + + /** + * deserialize an object of this class. + * @see mio::deserialize + */ + template + static IOResult deserialize(IOContext& io) + { + auto obj = io.expect_object("Model"); + auto par = obj.expect_element("Parameters", Tag{}); + auto pop = obj.expect_element("Populations", Tag{}); + return apply( + io, + [](auto&& par_, auto&& pop_) { + return Model{pop_, par_}; + }, + par, pop); + } +}; + +} // namespace omseirs4 +} // namespace mio + +#endif // ODE_MSEIRS4_MODEL_H diff --git a/cpp/models/ode_mseirs4/parameters.h b/cpp/models/ode_mseirs4/parameters.h new file mode 100644 index 0000000000..540de58f64 --- /dev/null +++ b/cpp/models/ode_mseirs4/parameters.h @@ -0,0 +1,244 @@ +/* +* 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. +*/ +#ifndef ODE_MSEIRS4_PARAMETERS_H +#define ODE_MSEIRS4_PARAMETERS_H + +#include "memilio/utils/parameter_set.h" +#include "memilio/utils/uncertain_value.h" +#include "memilio/utils/logging.h" +#include + +namespace mio +{ +namespace omseirs4 +{ +template +struct BaseTransmissionRate { // b0 + using Type = UncertainValue; + static Type get_default() + { + return Type(1.0); + } + static std::string name() + { + return "BaseTransmissionRate"; + } +}; + +template +struct SeasonalAmplitude { // b1 in [0,1] + using Type = UncertainValue; + static Type get_default() + { + return Type(0.0); + } + static std::string name() + { + return "SeasonalAmplitude"; + } +}; + +template +struct SeasonalPhase { // phi in radians + using Type = UncertainValue; + static Type get_default() + { + return Type(0.0); + } + static std::string name() + { + return "SeasonalPhase"; + } +}; + +template +struct NaturalBirthDeathRate { // mu + using Type = UncertainValue; + static Type get_default() + { + return Type(0.0); + } + static std::string name() + { + return "NaturalBirthDeathRate"; + } +}; + +template +struct LossMaternalImmunityRate { // xi + using Type = UncertainValue; + static Type get_default() + { + return Type(0.0); + } + static std::string name() + { + return "LossMaternalImmunityRate"; + } +}; + +template +struct ProgressionRate { // sigma, E->I + using Type = UncertainValue; + static Type get_default() + { + return Type(1.0); + } + static std::string name() + { + return "ProgressionRate"; + } +}; + +template +struct RecoveryRate { // nu, I->R + using Type = UncertainValue; + static Type get_default() + { + return Type(1.0); + } + static std::string name() + { + return "RecoveryRate"; + } +}; + +template +struct ImmunityWaningRate { // gamma, R -> S stages + using Type = UncertainValue; + static Type get_default() + { + return Type(0.0); + } + static std::string name() + { + return "ImmunityWaningRate"; + } +}; + +template +struct Beta2Factor { // f2 multiplier for beta1 + using Type = UncertainValue; + static Type get_default() + { + return Type(0.5); + } + static std::string name() + { + return "Beta2Factor"; + } +}; + +template +struct Beta3Factor { // f3 multiplier for beta1 + using Type = UncertainValue; + static Type get_default() + { + return Type(0.35); + } + static std::string name() + { + return "Beta3Factor"; + } +}; + +template +struct Beta4Factor { // f4 multiplier for beta1 + using Type = UncertainValue; + static Type get_default() + { + return Type(0.25); + } + static std::string name() + { + return "Beta4Factor"; + } +}; + +template +using ParametersBase = + ParameterSet, SeasonalAmplitude, SeasonalPhase, NaturalBirthDeathRate, + LossMaternalImmunityRate, ProgressionRate, RecoveryRate, ImmunityWaningRate, + Beta2Factor, Beta3Factor, Beta4Factor>; + +template +class Parameters : public ParametersBase +{ +public: + Parameters() + : ParametersBase() + { + } + + bool apply_constraints() + { + bool corrected = false; + auto clamp_nonneg = [&](auto& v) { + if (v < 0) { + v = 0; + corrected = true; + } + }; + clamp_nonneg(this->template get>()); + clamp_nonneg(this->template get>()); + if (this->template get>() > 1) { + this->template get>() = 1; + corrected = true; + } + clamp_nonneg(this->template get>()); + clamp_nonneg(this->template get>()); + clamp_nonneg(this->template get>()); + clamp_nonneg(this->template get>()); + clamp_nonneg(this->template get>()); + clamp_nonneg(this->template get>()); + clamp_nonneg(this->template get>()); + clamp_nonneg(this->template get>()); + return corrected; + } + + bool check_constraints() const + { + if (this->template get>() < 0 || this->template get>() > 1) + return true; + return false; + } +private: + Parameters(ParametersBase&& base) + : ParametersBase(std::move(base)) + { + } + +public: + /** + * deserialize an object of this class. + * @see mio::deserialize + */ + template + static IOResult deserialize(IOContext& io) + { + BOOST_OUTCOME_TRY(auto&& base, ParametersBase::deserialize(io)); + return success(Parameters(std::move(base))); + } +}; + +} // namespace omseirs4 +} // namespace mio + +#endif // ODE_MSEIRS4_PARAMETERS_H diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 2994eb2b38..01db381a7e 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -87,6 +87,7 @@ set(TESTSOURCES test_graph_abm.cpp test_abstract_parameter_dist.cpp test_temporal_hybrid_model.cpp + test_odemseirs4.cpp test_uncertain_ad_compatibility.cpp ) diff --git a/cpp/tests/test_odemseirs4.cpp b/cpp/tests/test_odemseirs4.cpp new file mode 100644 index 0000000000..4621f939b6 --- /dev/null +++ b/cpp/tests/test_odemseirs4.cpp @@ -0,0 +1,207 @@ +/* +* 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. +*/ +#include "memilio/compartments/simulation.h" +#include "memilio/math/euler.h" +#include "memilio/utils/time_series.h" +#include "ode_mseirs4/model.h" +#include "ode_mseirs4/infection_state.h" +#include "ode_mseirs4/parameters.h" + +#include + +TEST(TestOdeMseirs4, simulateDefault) +{ + double t0 = 0; + double tmax = 1; + double dt = 0.1; + + mio::omseirs4::Model model; + mio::TimeSeries result = simulate(t0, tmax, dt, model); + + EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); +} + +class ModelTestOdeMseirs4 : public testing::Test +{ +public: + ModelTestOdeMseirs4() + : model() + { + } + + double t0{}; + double tmax{}; + double dt{}; + double total_population{}; + mio::omseirs4::Model model; + +protected: + void SetUp() override + { + t0 = 0.0; + tmax = 1.0; + dt = 0.5; + + total_population = 1'000'000.0; + + // Initialize parameters (per day) and disable births/deaths for conservation test + auto& params = model.parameters; + params.get>() = 0.4; // b0 + params.get>() = 0.0; // no seasonality + params.get>() = 0.0; // phase + params.get>() = 0.0; // mu = 0 + params.get>() = 1.0 / 90; // xi + params.get>() = 1.0 / 7; // sigma + params.get>() = 1.0 / 14; // nu + params.get>() = 0.0; // gamma = 0 + // beta factors use defaults + + double M = 5'000.0; + double E1 = 300.0, E2 = 150.0, E3 = 80.0, E4 = 70.0; + double I1 = 200.0, I2 = 100.0, I3 = 50.0, I4 = 50.0; + double R1 = 40'000.0, R2 = 30'000.0, R3 = 20'000.0, R4 = 10'000.0; + double S2 = 100'000.0, S3 = 50'000.0, S4 = 50'000.0; + double assigned = M + (E1 + E2 + E3 + E4) + (I1 + I2 + I3 + I4) + (R1 + R2 + R3 + R4) + (S2 + S3 + S4); + double S1 = std::max(0.0, total_population - assigned); + + using IS = mio::omseirs4::InfectionState; + model.populations[{mio::Index(IS::MaternalImmune)}] = M; + model.populations[{mio::Index(IS::E1)}] = E1; + model.populations[{mio::Index(IS::E2)}] = E2; + model.populations[{mio::Index(IS::E3)}] = E3; + model.populations[{mio::Index(IS::E4)}] = E4; + model.populations[{mio::Index(IS::I1)}] = I1; + model.populations[{mio::Index(IS::I2)}] = I2; + model.populations[{mio::Index(IS::I3)}] = I3; + model.populations[{mio::Index(IS::I4)}] = I4; + model.populations[{mio::Index(IS::R1)}] = R1; + model.populations[{mio::Index(IS::R2)}] = R2; + model.populations[{mio::Index(IS::R3)}] = R3; + model.populations[{mio::Index(IS::R4)}] = R4; + model.populations[{mio::Index(IS::S1)}] = S1; + model.populations[{mio::Index(IS::S2)}] = S2; + model.populations[{mio::Index(IS::S3)}] = S3; + model.populations[{mio::Index(IS::S4)}] = S4; + + model.check_constraints(); + } +}; + +TEST_F(ModelTestOdeMseirs4, checkPopulationConservation) +{ + auto result = mio::simulate>(t0, tmax, dt, model); + double total_pop_last = result.get_last_value().sum(); + EXPECT_NEAR(total_pop_last, total_population, 1e-7); +} + +TEST(TestOdeMseirs4, apply_constraints_parameters) +{ + mio::omseirs4::Model model; + + auto& params = model.parameters; + + // valid default first + EXPECT_EQ(params.apply_constraints(), false); + + // negative values should be clamped to 0 + params.get>() = -0.5; + params.get>() = -1.0; + params.get>() = -1.0; + params.get>() = -1.0; + params.get>() = -1.0; + params.get>() = -1.0; + params.get>() = -0.1; + params.get>() = -0.2; + params.get>() = -0.3; + params.get>() = 1.5; // will be clamped to 1 + + EXPECT_EQ(params.apply_constraints(), true); + + EXPECT_DOUBLE_EQ((double)params.get>(), 0.0); + EXPECT_DOUBLE_EQ((double)params.get>(), 0.0); + EXPECT_DOUBLE_EQ((double)params.get>(), 0.0); + EXPECT_DOUBLE_EQ((double)params.get>(), 0.0); + EXPECT_DOUBLE_EQ((double)params.get>(), 0.0); + EXPECT_DOUBLE_EQ((double)params.get>(), 0.0); + EXPECT_DOUBLE_EQ((double)params.get>(), 0.0); + EXPECT_DOUBLE_EQ((double)params.get>(), 0.0); + EXPECT_DOUBLE_EQ((double)params.get>(), 0.0); + EXPECT_DOUBLE_EQ((double)params.get>(), 1.0); +} + +TEST(TestOdeMseirs4, check_constraints_parameters) +{ + mio::omseirs4::Model model; + auto& params = model.parameters; + + // default is valid + ASSERT_EQ(params.check_constraints(), false); + + // amplitude must be in [0,1] + params.get>() = -0.1; + ASSERT_EQ(params.check_constraints(), true); + + params.get>() = 0.5; + ASSERT_EQ(params.check_constraints(), false); + + params.get>() = 1.2; + ASSERT_EQ(params.check_constraints(), true); +} + +TEST(TestOdeMseirs4, population_zero_no_nan) +{ + mio::omseirs4::Model model; + model.populations.set_total(0.0); + + auto dydt = Eigen::VectorXd((Eigen::Index)mio::omseirs4::InfectionState::Count); + dydt.setZero(); + auto y0 = model.get_initial_values(); + model.get_derivatives(y0, y0, 0.0, dydt); + + for (int i = 0; i < dydt.size(); ++i) { + EXPECT_FALSE(std::isnan(dydt[i])); + } +} + +TEST(TestOdeMseirs4, Simulation) +{ + double t0 = 0; + double tmax = 1; + double dt = 1; + + mio::omseirs4::Model model; + + using IS = mio::omseirs4::InfectionState; + model.populations[{mio::Index(IS::I1)}] = 10.0; + model.populations[{mio::Index(IS::S1)}] = 990.0; + + model.parameters.get>() = 0.2; + model.parameters.get>() = 1.0 / 5.0; + model.parameters.get>() = 1.0 / 7.0; + model.parameters.get>() = 0.0; + model.parameters.get>() = 0.0; + + model.check_constraints(); + + std::unique_ptr> integrator = std::make_unique>(); + auto sim = mio::simulate(t0, tmax, dt, model, std::move(integrator)); + + EXPECT_EQ(sim.get_num_time_points(), 2); +} diff --git a/docs/source/cpp/models/omseirs4.rst b/docs/source/cpp/models/omseirs4.rst new file mode 100644 index 0000000000..dc6be8e818 --- /dev/null +++ b/docs/source/cpp/models/omseirs4.rst @@ -0,0 +1,143 @@ +MSEIRS4 model (ODE) +=================== + +The ODE-MSEIRS4 module models a pathogen with partial and waning immunity across multiple infection episodes, +including a maternal immunity class. It is suited for pathogens with repeat infections and seasonality. + +The infection states and the transitions are visualized in the following graph. + +.. image:: https://martinkuehn.eu/research/images/ode_mseirs4.png + :alt: mseirs4_model + + +This implementation is designed for Respiratory Syncytial Virus (RSV) and is based on: + +- Weber A, Weber M, Milligan P. (2001). *Modeling epidemics caused by respiratory syncytial virus (RSV).* Mathematical Biosciences 172(2): 95–113. `DOI:10.1016/S0025-5564(01)00066-9 `_ + +Infection States +---------------- + +The model contains the following InfectionStates: + +- `MaternalImmune` (M) +- `S1`, `S2`, `S3`, `S4` (susceptible classes by infection history) +- `E1`, `E2`, `E3`, `E4` (exposed/latent for the k-th infection episode) +- `I1`, `I2`, `I3`, `I4` (infectious for the k-th infection episode) +- `R1`, `R2`, `R3`, `R4` (recovered following the k-th infection episode) + +Meaning of indices 1–4 +---------------------- + +- S1: fully susceptible after loss of maternal immunity (highest susceptibility; seasonal force β1(t)). +- S2: susceptible after first infection (R1 → S2 via waning; reduced susceptibility; β2(t)=f2·β1(t)). +- S3: susceptible after second infection (R2 → S3; β3(t)=f3·β1(t)). +- S4: susceptible after ≥3 infections (R3 → S4 and R4 → S4; lowest susceptibility; β4(t)=f4·β1(t)). + +The multipliers :math:`f_2, f_3, f_4` are dimensionless susceptibility reductions for S2–S4. +All infectious classes (I1..I4) contribute equally to transmission in the basic formulation. + +Infection State Transitions +--------------------------- + +The model is implemented as a standard CompartmentalModel. The following transitions occur: + +- Births enter M and some enter S1 +- M → S1 (loss of maternal immunity) +- S_k → E_k +- E_k → I_k +- I_k → R_k +- R1 → S2, R2 → S3, R3 → S4 and R4 → S4 (waning immunity) +- Natural deaths apply to all compartments + +Seasonality +----------- + +Time unit is days. Seasonality follows a yearly cosine: + +.. math:: + + \beta_1(t) = b_0\,\big(1 + b_1\,\cos(2\pi\,t/365 + \varphi)\big),\quad \beta_k(t) = f_k\,\beta_1(t)\ (k=2,3,4). + +Here, :math:`b_0` is the base transmission rate (per day), :math:`b_1\in[0,1]` the seasonal amplitude, and :math:`\varphi` a phase (radians). + +Parameters +---------- + +The model implements the following parameters: + +.. list-table:: + :header-rows: 1 + :widths: 30 30 40 + + * - Mathematical symbol + - C++ parameter name + - Description + * - :math:`b_0` + - ``BaseTransmissionRate`` + - Base transmission rate (per day). + * - :math:`b_1` + - ``SeasonalAmplitude`` + - Seasonal amplitude. + * - :math:`\varphi` + - ``SeasonalPhase`` + - Phase of the cosine forcing (radians). + * - :math:`\mu` + - ``NaturalBirthDeathRate`` + - Natural birth/death rate (per day). + * - :math:`\xi` + - ``LossMaternalImmunityRate`` + - Rate of losing maternal immunity M→S1 (per day). + * - :math:`\sigma` + - ``ProgressionRate`` + - Progression rate E→I (per day). + * - :math:`\nu` + - ``RecoveryRate`` + - Recovery rate I→R (per day). + * - :math:`\gamma` + - ``ImmunityWaningRate`` + - Waning rate R→S stage (per day). + * - :math:`f_2` + - ``Beta2Factor`` + - Susceptibility multiplier for S2. + * - :math:`f_3` + - ``Beta3Factor`` + - Susceptibility multiplier for S3. + * - :math:`f_4` + - ``Beta4Factor`` + - Susceptibility multiplier for S4. + +Initial conditions +------------------ + +Initial conditions are absolute counts in each InfectionState; totals may be set directly. + +Example (see `examples/ode_mseirs4.cpp `_) shows a complete initialization and simulation. + +Simulation +---------- + +Run a standard simulation via: + +.. code-block:: cpp + + double t0 = 0.0; // days + double tmax = 3650; // 10 years + double dt = 1.0; // daily step + auto timeseries = mio::simulate(t0, tmax, dt, model); + +Output +------ + +The output is a ``mio::TimeSeries`` of compartment sizes over time. Use ``print_table`` or export to CSV as needed. + +Notes +----- + +- Homogeneous mixing; no age/contact matrices in this variant. +- All rates are per day. +- As in the paper, the model keeps N approximately constant if births and deaths balance. + +References +---------- + +- Weber A, Weber M, Milligan P. (2001). *Modeling epidemics caused by respiratory syncytial virus (RSV).* Mathematical Biosciences 172(2): 95–113. `DOI:10.1016/S0025-5564(01)00066-9 `_ diff --git a/docs/source/cpp/ode.rst b/docs/source/cpp/ode.rst index 83a2de9ce4..570d3acd83 100644 --- a/docs/source/cpp/ode.rst +++ b/docs/source/cpp/ode.rst @@ -144,3 +144,4 @@ List of models models/osecir models/osecirvvs models/osecirts + models/omseirs4 diff --git a/pycode/examples/simulation/ode_omseirs4_simple.py b/pycode/examples/simulation/ode_omseirs4_simple.py new file mode 100644 index 0000000000..d2ccac13fa --- /dev/null +++ b/pycode/examples/simulation/ode_omseirs4_simple.py @@ -0,0 +1,96 @@ +############################################################################# +# 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. +############################################################################# +import argparse + +import numpy as np + +from memilio.simulation.omseirs4 import InfectionState as State +from memilio.simulation.omseirs4 import Model, simulate + +# MSEIRS4 model (maternal immunity + 4 infection levels) as described in https://doi.org/10.1016/S0025-5564(01)00066-9. + + +def run_ode_omseirs4_simulation(): + # Create model and set parameters (per day) + model = Model() + params = model.parameters + + params.BaseTransmissionRate.value = 0.4 + params.SeasonalAmplitude.value = 0.15 + params.SeasonalPhase.value = 0.0 + params.NaturalBirthDeathRate.value = 1.0 / (70.0 * 365.0) + params.LossMaternalImmunityRate.value = 1.0 / 90.0 + params.ProgressionRate.value = 1.0 / 7.0 + params.RecoveryRate.value = 1.0 / 14.0 + params.ImmunityWaningRate.value = 1.0 / (5.0 * 365.0) + params.Beta2Factor.value = 0.5 + params.Beta3Factor.value = 0.35 + params.Beta4Factor.value = 0.25 + + # Initial populations + N = 1_000_000.0 + pop = model.populations + + pop[State.MaternalImmune] = 5000.0 + + pop[State.E1] = 300.0 + pop[State.E2] = 150.0 + pop[State.E3] = 80.0 + pop[State.E4] = 70.0 + + pop[State.I1] = 200.0 + pop[State.I2] = 100.0 + pop[State.I3] = 50.0 + pop[State.I4] = 50.0 + + pop[State.R1] = 40000.0 + pop[State.R2] = 30000.0 + pop[State.R3] = 20000.0 + pop[State.R4] = 10000.0 + + pop[State.S2] = 100000.0 + pop[State.S3] = 50000.0 + pop[State.S4] = 50000.0 + + assigned = 5000.0 + (300.0 + 150.0 + 80.0 + 70.0) + (200.0 + 100.0 + 50.0 + 50.0) + \ + (40000.0 + 30000.0 + 20000.0 + 10000.0) + (100000.0 + 50000.0 + 50000.0) + s1 = max(0.0, N - assigned) + pop[State.S1] = s1 + + model.check_constraints() + + # simulate 10 days with daily time steps + t0 = 0.0 + tmax = 10.0 + dt = 1.0 + + result = simulate(t0, tmax, dt, model) + + # print last results + print("Last results:") + print(f"Time: {result.get_last_time()}") + print(f"compartments: {result.get_last_value()}") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("ode_omseirs4_simple", + description="Simple example demonstrating the ODE MSEIRS4 RSV model") + args = parser.parse_args() + run_ode_omseirs4_simulation() diff --git a/pycode/memilio-simulation/CMakeLists.txt b/pycode/memilio-simulation/CMakeLists.txt index 975579d3fd..ecececcb2c 100644 --- a/pycode/memilio-simulation/CMakeLists.txt +++ b/pycode/memilio-simulation/CMakeLists.txt @@ -42,14 +42,12 @@ function(add_pymio_module target_name) set(multiValueArgs LINKED_LIBRARIES SOURCES) cmake_parse_arguments(PYBIND11_MODULE "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) - pybind11_add_module(${target_name} MODULE ${PYBIND11_MODULE_SOURCES}) target_link_libraries(${target_name} PRIVATE ${PYBIND11_MODULE_LINKED_LIBRARIES}) target_include_directories(${target_name} PRIVATE memilio/simulation/bindings) install(TARGETS ${target_name} LIBRARY DESTINATION memilio/simulation) endfunction() - # build python extensions add_pymio_module(_simulation_abm LINKED_LIBRARIES memilio abm @@ -88,3 +86,8 @@ add_pymio_module(_simulation_osecirvvs LINKED_LIBRARIES memilio ode_secirvvs SOURCES memilio/simulation/bindings/models/osecirvvs.cpp ) + +add_pymio_module(_simulation_omseirs4 + LINKED_LIBRARIES memilio ode_mseirs4 + SOURCES memilio/simulation/bindings/models/omseirs4.cpp +) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/omseirs4.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/omseirs4.cpp new file mode 100644 index 0000000000..faf7eb78a5 --- /dev/null +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/omseirs4.cpp @@ -0,0 +1,116 @@ +/* +* 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/compartmental_model.h" +#include "epidemiology/populations.h" +#include "utils/parameter_set.h" +#include "utils/index.h" + +// Includes from MEmilio +#include "ode_mseirs4/model.h" +#include "ode_mseirs4/infection_state.h" +#include "ode_mseirs4/parameters.h" +#include "memilio/data/analyze_result.h" + +#include "pybind11/pybind11.h" +#include "Eigen/Core" +#include + +namespace py = pybind11; + +namespace pymio +{ +// specialization of pretty_name +template <> +inline std::string pretty_name() +{ + return "InfectionState"; +} +} // namespace pymio + +PYBIND11_MODULE(_simulation_omseirs4, m) +{ + // interpolation helpers + m.def("interpolate_simulation_result", + static_cast (*)(const mio::TimeSeries&, const double)>( + &mio::interpolate_simulation_result), + py::arg("ts"), py::arg("abs_tol") = 1e-14); + + m.def("interpolate_simulation_result", + static_cast (*)(const mio::TimeSeries&, const std::vector&)>( + &mio::interpolate_simulation_result), + py::arg("ts"), py::arg("interpolation_times")); + + m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results>); + + // InfectionState enum + pymio::iterable_enum(m, "InfectionState") + .value("MaternalImmune", mio::omseirs4::InfectionState::MaternalImmune) + .value("S1", mio::omseirs4::InfectionState::S1) + .value("S2", mio::omseirs4::InfectionState::S2) + .value("S3", mio::omseirs4::InfectionState::S3) + .value("S4", mio::omseirs4::InfectionState::S4) + .value("E1", mio::omseirs4::InfectionState::E1) + .value("E2", mio::omseirs4::InfectionState::E2) + .value("E3", mio::omseirs4::InfectionState::E3) + .value("E4", mio::omseirs4::InfectionState::E4) + .value("I1", mio::omseirs4::InfectionState::I1) + .value("I2", mio::omseirs4::InfectionState::I2) + .value("I3", mio::omseirs4::InfectionState::I3) + .value("I4", mio::omseirs4::InfectionState::I4) + .value("R1", mio::omseirs4::InfectionState::R1) + .value("R2", mio::omseirs4::InfectionState::R2) + .value("R3", mio::omseirs4::InfectionState::R3) + .value("R4", mio::omseirs4::InfectionState::R4); + + // Parameters + pymio::bind_ParameterSet, pymio::EnablePickling::Required>(m, + "ParametersBase"); + + pymio::bind_class, pymio::EnablePickling::Required, + mio::omseirs4::ParametersBase>(m, "Parameters") + .def(py::init<>()) + .def("check_constraints", &mio::omseirs4::Parameters::check_constraints) + .def("apply_constraints", &mio::omseirs4::Parameters::apply_constraints); + + // Populations and Model + using Populations = mio::Populations; + pymio::bind_Population(m, "Populations", mio::Tag{}); + pymio::bind_CompartmentalModel, + pymio::EnablePickling::Never>(m, "ModelBase"); + pymio::bind_class< + mio::omseirs4::Model, pymio::EnablePickling::Required, + mio::CompartmentalModel>>( + m, "Model") + .def(py::init<>()); + + // Simulation and simulate() + using Sim = mio::Simulation>; + pymio::bind_Simulation(m, "Simulation"); + + m.def("simulate", &mio::simulate>, + "Simulates an ODE MSEIRS4 model 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/omseirs4.py b/pycode/memilio-simulation/memilio/simulation/omseirs4.py new file mode 100644 index 0000000000..62511000ca --- /dev/null +++ b/pycode/memilio-simulation/memilio/simulation/omseirs4.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 ODE MSEIRS4 (RSV) model from https://doi.org/10.1016/S0025-5564(01)00066-9 +""" + +from memilio.simulation._simulation_omseirs4 import * diff --git a/pycode/memilio-simulation/memilio/simulation_test/test_omseirs4.py b/pycode/memilio-simulation/memilio/simulation_test/test_omseirs4.py new file mode 100644 index 0000000000..cc5bae3edf --- /dev/null +++ b/pycode/memilio-simulation/memilio/simulation_test/test_omseirs4.py @@ -0,0 +1,124 @@ +############################################################################# +# 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. +############################################################################# +import unittest + +from memilio.simulation.omseirs4 import InfectionState as State +from memilio.simulation.omseirs4 import Model, Simulation, simulate +import memilio.simulation as mio + + +class Test_omseirs4_integration(unittest.TestCase): + """Python binding tests for ODE MSEIRS4 model.""" + + def setUp(self): + # create a simple single-population model + model = Model() + + # parameters + model.parameters.BaseTransmissionRate.value = 0.2 + model.parameters.SeasonalAmplitude.value = 0.0 + model.parameters.SeasonalPhase.value = 0.0 + model.parameters.NaturalBirthDeathRate.value = 0.0 + model.parameters.LossMaternalImmunityRate.value = 1.0 / 90.0 + model.parameters.ProgressionRate.value = 1.0 / 7.0 + model.parameters.RecoveryRate.value = 1.0 / 14.0 + model.parameters.ImmunityWaningRate.value = 0.0 + # Beta2/3/4Factor keep defaults + + # initialize populations + N = 100000.0 + assigned = 0.0 + model.populations[State.E1] = 20.0 + model.populations[State.I1] = 10.0 + assigned += 30.0 + model.populations[State.R1] = 100.0 + model.populations[State.R2] = 50.0 + assigned += 150.0 + model.populations[State.S2] = 500.0 + model.populations[State.S3] = 200.0 + model.populations[State.S4] = 100.0 + assigned += 800.0 + model.populations[State.MaternalImmune] = 50.0 + assigned += 50.0 + model.populations[State.S1] = max(0.0, N - assigned) + + model.check_constraints() + + self.model = model + + def test_simulate_simple(self): + integrator = mio.EulerIntegratorCore() + result = simulate(t0=0.0, tmax=50.0, dt=0.5, + model=self.model, integrator=integrator) + self.assertAlmostEqual(result.get_time(0), 0.0) + self.assertAlmostEqual(result.get_time(1), 0.5) + self.assertAlmostEqual(result.get_last_time(), 50.0) + # vector length equals number of infection states + self.assertEqual(len(result.get_last_value()), + 17) # 17 infection states + + def test_simulation_simple(self): + sim = Simulation(self.model, t0=0.0, dt=1.0) + sim.advance(tmax=10.0) + res = sim.result + self.assertAlmostEqual(res.get_time(0), 0.0) + self.assertAlmostEqual(res.get_last_time(), 10.0) + self.assertEqual(len(res.get_last_value()), 17) # 17 infection states + + def test_check_and_apply_constraints_parameters(self): + model = Model() + p = model.parameters + + # default valid + self.assertEqual(p.check_constraints(), False) + + p.SeasonalAmplitude.value = -0.1 + self.assertEqual(p.check_constraints(), True) + p.SeasonalAmplitude.value = 0.5 + self.assertEqual(p.check_constraints(), False) + p.SeasonalAmplitude.value = 1.2 + self.assertEqual(p.check_constraints(), True) + + p.BaseTransmissionRate.value = -0.3 + p.ProgressionRate.value = -1.0 + p.RecoveryRate.value = -1.0 + p.ImmunityWaningRate.value = -1.0 + p.NaturalBirthDeathRate.value = -1.0 + p.LossMaternalImmunityRate.value = -1.0 + p.Beta2Factor.value = -0.1 + p.Beta3Factor.value = -0.2 + p.Beta4Factor.value = -0.3 + p.SeasonalAmplitude.value = 1.5 + changed = p.apply_constraints() + self.assertTrue(changed) + self.assertAlmostEqual(p.BaseTransmissionRate.value, 0.0) + self.assertAlmostEqual(p.ProgressionRate.value, 0.0) + self.assertAlmostEqual(p.RecoveryRate.value, 0.0) + self.assertAlmostEqual(p.ImmunityWaningRate.value, 0.0) + self.assertAlmostEqual(p.NaturalBirthDeathRate.value, 0.0) + self.assertAlmostEqual(p.LossMaternalImmunityRate.value, 0.0) + self.assertAlmostEqual(p.Beta2Factor.value, 0.0) + self.assertAlmostEqual(p.Beta3Factor.value, 0.0) + self.assertAlmostEqual(p.Beta4Factor.value, 0.0) + self.assertAlmostEqual(p.SeasonalAmplitude.value, 1.0) + + +if __name__ == '__main__': + unittest.main()