1350 Implement RSV model

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,20 +140,22 @@ 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
             "-Wno-stringop-overread;"
         )
     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)

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})

cpp/examples/ode_mseirs4.cpp

@@ -0,0 +1,94 @@
+/* 
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Henrik Zunker
+*
+* 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.
+*/
+#include "ode_mseirs4/model.h"
+#include "memilio/compartments/simulation.h"
+#include <iostream>
+
+int main()
+{
+    using FP = double;
+    mio::omseirs4::Model<FP> model;
+    auto& params = model.parameters;
+
+    // Example parameter values for day-based time unit (t in days)
+    params.get<mio::omseirs4::BaseTransmissionRate<FP>>()     = 0.4; // b0 per day
+    params.get<mio::omseirs4::SeasonalAmplitude<FP>>()        = 0.15; // b1
+    params.get<mio::omseirs4::SeasonalPhase<FP>>()            = 0.0; // phi; for phase shift use 2*pi*offsetDays/365
+    params.get<mio::omseirs4::NaturalBirthDeathRate<FP>>()    = 1.0 / (70.0 * 365.0); // mu per day
+    params.get<mio::omseirs4::LossMaternalImmunityRate<FP>>() = 1.0 / 90.0; // xi per day (~3 months)
+    params.get<mio::omseirs4::ProgressionRate<FP>>()          = 1.0 / 7.0; // sigma per day (≈7 days latent)
+    params.get<mio::omseirs4::RecoveryRate<FP>>()             = 1.0 / 14.0; // nu per day (≈14 days infectious)
+    params.get<mio::omseirs4::ImmunityWaningRate<FP>>()       = 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>(mio::omseirs4::InfectionState::MaternalImmune)}] =
+        5000.0;
+
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::E1)}] = 300.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::E2)}] = 150.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::E3)}] = 80.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::E4)}] = 70.0;
+
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::I1)}] = 200.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::I2)}] = 100.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::I3)}] = 50.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::I4)}] = 50.0;
+
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::R1)}] = 40000.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::R2)}] = 30000.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::R3)}] = 20000.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::R4)}] = 10000.0;
+
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::S2)}] = 100000.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::S3)}] = 50000.0;
+    model.populations[{mio::Index<mio::omseirs4::InfectionState>(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>(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;
+}

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
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+target_compile_options(ode_mseirs4 PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

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 <https://doi.org/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.

cpp/models/ode_mseirs4/infection_state.h

@@ -0,0 +1,52 @@
+/* 
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Henrik Zunker
+*
+* 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.
+*/
+#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

cpp/models/ode_mseirs4/model.cpp

@@ -0,0 +1,28 @@
+/* 
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Henrik Zunker
+*
+* 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.
+*/
+#include "ode_mseirs4/model.h"
+
+namespace mio
+{
+namespace omseirs4
+{
+
+} // namespace omseirs4
+} // namespace mio