1420 Implement SEIRDB model (#1424)

Extended the SEIR model to a SEIRDB (with D: dead but not yet buried and B: buried) structure.

Additional parameters are:
a probability to recover (instead of getting "Removed" from SEIR) with counterprobability to die
a time between dying and getting buried
a probability to infect from dead but not yet buried state D


Co-authored-by: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com>

Author: HenrZu

cpp/CMakeLists.txt

@@ -183,7 +183,8 @@ if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/glct_secir)
     add_subdirectory(models/ide_secir)
     add_subdirectory(models/ide_seir)
-    add_subdirectory(models/ode_seir)
+    add_subdirectory(models/ode_seir)
+    add_subdirectory(models/ode_seirdb)
     add_subdirectory(models/ode_seair)
     add_subdirectory(models/ode_sir)
     add_subdirectory(models/sde_sir)

cpp/examples/CMakeLists.txt

@@ -18,13 +18,17 @@ add_executable(adapt_rk_example adapt_rk_test.cpp)
 target_link_libraries(adapt_rk_example PRIVATE memilio)
 target_compile_options(adapt_rk_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
-add_executable(ode_seir_example ode_seir.cpp)
-target_link_libraries(ode_seir_example PRIVATE memilio ode_seir)
-target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
-
-add_executable(ode_seir_ageres_example ode_seir_ageres.cpp)
-target_link_libraries(ode_seir_ageres_example PRIVATE memilio ode_seir)
-target_compile_options(ode_seir_ageres_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+add_executable(ode_seir_example ode_seir.cpp)
+target_link_libraries(ode_seir_example PRIVATE memilio ode_seir)
+target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(ode_seirdb_example ode_seirdb.cpp)
+target_link_libraries(ode_seirdb_example PRIVATE memilio ode_seirdb)
+target_compile_options(ode_seirdb_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(ode_seir_ageres_example ode_seir_ageres.cpp)
+target_link_libraries(ode_seir_ageres_example PRIVATE memilio ode_seir)
+target_compile_options(ode_seir_ageres_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
 add_executable(ode_sir_example ode_sir.cpp)
 target_link_libraries(ode_sir_example PRIVATE memilio ode_sir)

cpp/examples/ode_seirdb.cpp

@@ -0,0 +1,71 @@
+/*
+* 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_seirdb/model.h"
+#include "ode_seirdb/infection_state.h"
+#include "ode_seirdb/parameters.h"
+#include "memilio/compartments/simulation.h"
+#include "memilio/utils/logging.h"
+#include "memilio/utils/time_series.h"
+
+int main()
+{
+    mio::set_log_level(mio::LogLevel::debug);
+
+    ScalarType t0   = 0;
+    ScalarType tmax = 30.;
+    ScalarType dt   = 1.0;
+
+    mio::log_info("Simulating ODE SEIRDB; t={} ... {} with dt = {}.", t0, tmax, dt);
+
+    mio::oseirdb::Model<ScalarType> model(1);
+
+    ScalarType total_population                                                    = 10000;
+    model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Exposed}]   = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Infected}]  = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Recovered}] = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Dead}]      = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Buried}]    = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Susceptible}] =
+        total_population - model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Exposed}] -
+        model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Infected}] -
+        model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Recovered}] -
+        model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Dead}] -
+        model.populations[{mio::AgeGroup(0), mio::oseirdb::InfectionState::Buried}];
+
+    model.parameters.set<mio::oseirdb::TimeExposed<ScalarType>>(5.2);
+    model.parameters.set<mio::oseirdb::TimeInfected<ScalarType>>(6.0);
+    model.parameters.set<mio::oseirdb::TimeToBurial<ScalarType>>(5.5);
+    model.parameters.set<mio::oseirdb::TransmissionProbabilityOnContact<ScalarType>>(0.1);
+    model.parameters.set<mio::oseirdb::TransmissionProbabilityFromDead<ScalarType>>(0.01);
+    model.parameters.set<mio::oseirdb::ProbabilityToRecover<ScalarType>>(0.75);
+
+    mio::ContactMatrixGroup<ScalarType>& contact_matrix =
+        model.parameters.get<mio::oseirdb::ContactPatterns<ScalarType>>();
+    contact_matrix[0].get_baseline().setConstant(2.7);
+    contact_matrix[0].add_damping(0.7, mio::SimulationTime<ScalarType>(10.));
+
+    model.check_constraints();
+
+    auto result = mio::simulate<ScalarType>(t0, tmax, dt, model);
+
+    result.print_table({"S", "E", "I", "R", "D", "B"});
+    std::cout << "\nnumber total: " << result.get_last_value().sum() << "\n";
+}

cpp/models/ode_seirdb/CMakeLists.txt

@@ -0,0 +1,12 @@
+add_library(ode_seirdb
+    infection_state.h
+    model.h
+    model.cpp
+    parameters.h
+)
+target_link_libraries(ode_seirdb PUBLIC memilio)
+target_include_directories(ode_seirdb PUBLIC
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+target_compile_options(ode_seirdb PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

cpp/models/ode_seirdb/README.md

@@ -0,0 +1,7 @@
+# ODE SEIRDB compartment model
+
+This directory contains a SEIRDB (Susceptible-Exposed-Infected-Recovered-Dead-Buried) model implementation based on an ODE formulation.
+It extends the SEIR model by introducing explicit Dead and Buried compartments with configurable recovery probability, burial delay, and infectiousness of the dead.
+
+To get started, check out the [official documentation](https://memilio.readthedocs.io/en/latest/cpp/models/oseirdb.html) 
+or the [code example](../../examples/ode_seirdb.cpp).

cpp/models/ode_seirdb/infection_state.h

@@ -0,0 +1,46 @@
+/*
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Daniel Abele, Jan Kleinert, Martin J. Kuehn
+*
+* 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 SEIRDB_INFECTIONSTATE_H
+#define SEIRDB_INFECTIONSTATE_H
+
+namespace mio
+{
+namespace oseirdb
+{
+
+/**
+     * @brief The InfectionState enum describes the possible
+     * categories for the infectious state of persons
+     */
+enum class InfectionState
+{
+    Susceptible,
+    Exposed,
+    Infected,
+    Recovered,
+    Dead,
+    Buried,
+    Count
+};
+
+} // namespace oseirdb
+} // namespace mio
+
+#endif // SEIRDB_INFECTIONSTATE_H

cpp/models/ode_seirdb/model.cpp

@@ -0,0 +1,28 @@
+/*
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Daniel Abele, Jan Kleinert, Martin J. Kuehn
+*
+* 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_seirdb/model.h"
+
+namespace mio
+{
+namespace oseirdb
+{
+
+} // namespace oseirdb
+} // namespace mio

cpp/models/ode_seirdb/model.h

@@ -0,0 +1,153 @@
+/*
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Daniel Abele, Jan Kleinert, Martin J. Kuehn
+*
+* 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 SEIRDB_MODEL_H
+#define SEIRDB_MODEL_H
+
+#include "memilio/compartments/flow_model.h"
+#include "memilio/config.h"
+#include "memilio/epidemiology/age_group.h"
+#include "memilio/epidemiology/populations.h"
+#include "memilio/math/interpolation.h"
+#include "memilio/utils/time_series.h"
+#include "ode_seirdb/infection_state.h"
+#include "ode_seirdb/parameters.h"
+
+GCC_CLANG_DIAGNOSTIC(push)
+GCC_CLANG_DIAGNOSTIC(ignored "-Wshadow")
+#include <Eigen/Dense>
+GCC_CLANG_DIAGNOSTIC(pop)
+
+namespace mio
+{
+namespace oseirdb
+{
+
+/********************
+ * define the model *
+ ********************/
+
+// clang-format off
+using Flows = TypeList<Flow<InfectionState::Susceptible, InfectionState::Exposed>,
+                       Flow<InfectionState::Exposed, InfectionState::Infected>,
+                       Flow<InfectionState::Infected, InfectionState::Recovered>,
+                       Flow<InfectionState::Infected, InfectionState::Dead>,
+                       Flow<InfectionState::Dead, InfectionState::Buried>>;
+// clang-format on
+template <typename FP>
+class Model
+    : public FlowModel<FP, InfectionState, mio::Populations<FP, AgeGroup, InfectionState>, Parameters<FP>, Flows>
+{
+    using Base = FlowModel<FP, InfectionState, mio::Populations<FP, AgeGroup, InfectionState>, Parameters<FP>, Flows>;
+
+public:
+    using typename Base::ParameterSet;
+    using typename Base::Populations;
+
+    Model(const Populations& pop, const ParameterSet& params)
+        : Base(pop, params)
+    {
+    }
+
+    Model(int num_agegroups)
+        : Base(Populations({AgeGroup(num_agegroups), InfectionState::Count}), ParameterSet(AgeGroup(num_agegroups)))
+    {
+    }
+
+    void get_flows(Eigen::Ref<const Eigen::VectorX<FP>> pop, Eigen::Ref<const Eigen::VectorX<FP>> y, FP t,
+                   Eigen::Ref<Eigen::VectorX<FP>> flows) const override
+    {
+        const Index<AgeGroup> age_groups = reduce_index<Index<AgeGroup>>(this->populations.size());
+        const auto& params               = this->parameters;
+
+        for (auto i : make_index_range(age_groups)) {
+            const size_t Si = this->populations.get_flat_index({i, InfectionState::Susceptible});
+            const size_t Ei = this->populations.get_flat_index({i, InfectionState::Exposed});
+            const size_t Ii = this->populations.get_flat_index({i, InfectionState::Infected});
+            const size_t Di = this->populations.get_flat_index({i, InfectionState::Dead});
+
+            for (auto j : make_index_range(age_groups)) {
+                const size_t Sj = this->populations.get_flat_index({j, InfectionState::Susceptible});
+                const size_t Ej = this->populations.get_flat_index({j, InfectionState::Exposed});
+                const size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected});
+                const size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered});
+                const size_t Dj = this->populations.get_flat_index({j, InfectionState::Dead});
+                const size_t Bj = this->populations.get_flat_index({j, InfectionState::Buried});
+
+                const FP Nj           = pop[Sj] + pop[Ej] + pop[Ij] + pop[Rj] + pop[Dj] + pop[Bj];
+                const FP divNj        = (Nj < Limits<FP>::zero_tolerance()) ? FP(0.0) : FP(1.0 / Nj);
+                const FP contact_rate = params.template get<ContactPatterns<FP>>().get_cont_freq_mat().get_matrix_at(
+                                            SimulationTime<FP>(t))(i.get(), j.get()) *
+                                        divNj;
+                const FP coeffStoE     = contact_rate * params.template get<TransmissionProbabilityOnContact<FP>>()[i];
+                const FP coeffStoEDead = contact_rate * params.template get<TransmissionProbabilityFromDead<FP>>()[i];
+
+                flows[Base::template get_flat_flow_index<InfectionState::Susceptible, InfectionState::Exposed>(i)] +=
+                    y[Si] * (coeffStoE * pop[Ij] + coeffStoEDead * pop[Dj]);
+            }
+            flows[Base::template get_flat_flow_index<InfectionState::Exposed, InfectionState::Infected>(i)] =
+                (1.0 / params.template get<TimeExposed<FP>>()[i]) * y[Ei];
+
+            const FP inv_time_infected = 1.0 / params.template get<TimeInfected<FP>>()[i];
+            const FP prob_recover      = params.template get<ProbabilityToRecover<FP>>()[i];
+            const FP prob_die          = 1.0 - prob_recover;
+            flows[Base::template get_flat_flow_index<InfectionState::Infected, InfectionState::Recovered>(i)] =
+                inv_time_infected * prob_recover * y[Ii];
+            flows[Base::template get_flat_flow_index<InfectionState::Infected, InfectionState::Dead>(i)] =
+                inv_time_infected * prob_die * y[Ii];
+            flows[Base::template get_flat_flow_index<InfectionState::Dead, InfectionState::Buried>(i)] =
+                (1.0 / params.template get<TimeToBurial<FP>>()[i]) * y[Di];
+        }
+    }
+
+    /**
+     * serialize this.
+     * @see mio::serialize
+     */
+    template <class IOContext>
+    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 <class IOContext>
+    static IOResult<Model> deserialize(IOContext& io)
+    {
+        auto obj = io.expect_object("Model");
+        auto par = obj.expect_element("Parameters", Tag<ParameterSet>{});
+        auto pop = obj.expect_element("Populations", Tag<Populations>{});
+        return apply(
+            io,
+            [](auto&& par_, auto&& pop_) {
+                return Model{pop_, par_};
+            },
+            par, pop);
+    }
+};
+
+} // namespace oseirdb
+} // namespace mio
+
+#endif // SEIRDB_MODEL_H