953 Implement basic stochastic equation based model (#954)

- Implemented stochastic version of SIR and SIRS model, with specialized Simulations.
- Implemented Euler-Maruyama method reusing EulerIntegratorCore and a custom advance function.
- Time step dt is now updated during integration. This is currently only used by the new SDE models.
- Protected IntegratorCore for stochastic models to avoid unintended use of ODE integrators.
- Small fixes: avoid two dangling references, avoid a copy in FlowSimulation.

Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com>
Co-authored-by: Martin J. Kühn <62713180+mknaranja@users.noreply.github.com>
Co-authored-by: HenrZu <henrik.zunker@dlr.de>

Author: 4 people

cpp/CMakeLists.txt

@@ -124,6 +124,8 @@ if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/ide_seir)
     add_subdirectory(models/ode_seir)
     add_subdirectory(models/ode_sir)
+    add_subdirectory(models/sde_sir)
+    add_subdirectory(models/sde_sirs)
 endif()
 
 if(MEMILIO_BUILD_EXAMPLES)

cpp/examples/CMakeLists.txt

@@ -20,9 +20,19 @@ target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNI
 
 add_executable(ode_sir_example ode_sir.cpp)
 target_link_libraries(ode_sir_example PRIVATE memilio ode_sir)
+target_compile_options(ode_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(sde_sir_example sde_sir.cpp)
+target_link_libraries(sde_sir_example PRIVATE memilio sde_sir)
+target_compile_options(sde_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(sde_sirs_example sde_sirs.cpp)
+target_link_libraries(sde_sirs_example PRIVATE memilio sde_sirs)
+target_compile_options(sde_sirs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
 add_executable(seir_flows_example ode_seir_flows.cpp)
 target_link_libraries(seir_flows_example PRIVATE memilio ode_seir)
+target_compile_options(seir_flows_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
 add_executable(ode_secir_example ode_secir.cpp)
 target_link_libraries(ode_secir_example PRIVATE memilio ode_secir)

cpp/examples/ode_sir.cpp

@@ -18,11 +18,12 @@
 * limitations under the License.
 */
 
-#include "ode_sir/model.h"
-#include "ode_sir/infection_state.h"
-#include "ode_sir/parameters.h"
 #include "memilio/compartments/simulation.h"
+#include "memilio/math/euler.h"
 #include "memilio/utils/logging.h"
+#include "ode_sir/infection_state.h"
+#include "ode_sir/model.h"
+#include "ode_sir/parameters.h"
 
 int main()
 {

cpp/examples/sde_sir.cpp

@@ -0,0 +1,55 @@
+/* 
+* Copyright (C) 2020-2024 MEmilio
+*
+* Authors: Nils Wassmuth, Rene Schmieding, 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 "memilio/utils/logging.h"
+#include "sde_sir/model.h"
+#include "sde_sir/simulation.h"
+
+int main()
+{
+    mio::set_log_level(mio::LogLevel::debug);
+
+    double t0   = 0.;
+    double tmax = 5.;
+    double dt   = 0.1;
+
+    double total_population = 10000;
+
+    mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt);
+
+    mio::ssir::Model model;
+
+    model.populations[{mio::Index<mio::ssir::InfectionState>(mio::ssir::InfectionState::Infected)}]  = 100;
+    model.populations[{mio::Index<mio::ssir::InfectionState>(mio::ssir::InfectionState::Recovered)}] = 1000;
+    model.populations[{mio::Index<mio::ssir::InfectionState>(mio::ssir::InfectionState::Susceptible)}] =
+        total_population -
+        model.populations[{mio::Index<mio::ssir::InfectionState>(mio::ssir::InfectionState::Infected)}] -
+        model.populations[{mio::Index<mio::ssir::InfectionState>(mio::ssir::InfectionState::Recovered)}];
+    model.parameters.set<mio::ssir::TimeInfected>(10);
+    model.parameters.set<mio::ssir::TransmissionProbabilityOnContact>(1);
+    model.parameters.get<mio::ssir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
+    model.parameters.get<mio::ssir::ContactPatterns>().add_damping(0.6, mio::SimulationTime(12.5));
+
+    model.check_constraints();
+
+    auto sir = mio::ssir::simulate(t0, tmax, dt, model);
+
+    sir.print_table();
+}

cpp/examples/sde_sirs.cpp

@@ -0,0 +1,56 @@
+/* 
+* Copyright (C) 2020-2024 MEmilio
+*
+* Authors: Nils Wassmuth, Rene Schmieding, 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 "memilio/utils/logging.h"
+#include "sde_sirs/model.h"
+#include "sde_sirs/simulation.h"
+
+int main()
+{
+    mio::set_log_level(mio::LogLevel::debug);
+
+    double t0   = 0.;
+    double tmax = 5.;
+    double dt   = 0.001;
+
+    double total_population = 10000;
+
+    mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt);
+
+    mio::ssirs::Model model;
+
+    model.populations[{mio::Index<mio::ssirs::InfectionState>(mio::ssirs::InfectionState::Infected)}]  = 100;
+    model.populations[{mio::Index<mio::ssirs::InfectionState>(mio::ssirs::InfectionState::Recovered)}] = 1000;
+    model.populations[{mio::Index<mio::ssirs::InfectionState>(mio::ssirs::InfectionState::Susceptible)}] =
+        total_population -
+        model.populations[{mio::Index<mio::ssirs::InfectionState>(mio::ssirs::InfectionState::Infected)}] -
+        model.populations[{mio::Index<mio::ssirs::InfectionState>(mio::ssirs::InfectionState::Recovered)}];
+    model.parameters.set<mio::ssirs::TimeInfected>(10);
+    model.parameters.set<mio::ssirs::TimeImmune>(100);
+    model.parameters.set<mio::ssirs::TransmissionProbabilityOnContact>(1);
+    model.parameters.get<mio::ssirs::ContactPatterns>().get_baseline()(0, 0) = 20.7;
+    model.parameters.get<mio::ssirs::ContactPatterns>().add_damping(0.6, mio::SimulationTime(12.5));
+
+    model.check_constraints();
+
+    auto ssirs = mio::ssirs::simulate(t0, tmax, dt, model);
+
+    ssirs.print_table({"Susceptible", "Infected", "Recovered"});
+}

cpp/memilio/compartments/flow_simulation.h

@@ -37,7 +37,7 @@ class FlowSimulation : public Simulation<M>
 
     /**
      * @brief Set up the simulation with an ODE solver.
-     * @param[in] model An instance of a compartmental model.
+     * @param[in] model An instance of a flow model.
      * @param[in] t0 Start time.
      * @param[in] dt Initial step size of integration.
      */
@@ -60,8 +60,8 @@ class FlowSimulation : public Simulation<M>
         assert(m_flow_result.get_num_time_points() == this->get_result().get_num_time_points());
         auto result = this->get_ode_integrator().advance(
             [this](auto&& flows, auto&& t, auto&& dflows_dt) {
-                auto pop_result = this->get_result();
-                auto model      = this->get_model();
+                const auto& pop_result = this->get_result();
+                const auto& model      = this->get_model();
                 // compute current population
                 //   flows contains the accumulated outflows of each compartment for each target compartment at time t.
                 //   Using that the ODEs are linear expressions of the flows, get_derivatives can compute the total change
@@ -105,7 +105,7 @@ class FlowSimulation : public Simulation<M>
     }
     /** @} */
 
-private:
+protected:
     /**
      * @brief Computes the distribution of the Population to the InfectionState%s based on the simulated flows.
      * Uses the same method as the DerivFunction used in advance to compute the population given the flows and initial
@@ -128,18 +128,22 @@ class FlowSimulation : public Simulation<M>
     }
 
     Eigen::VectorXd m_pop; ///< pre-allocated temporary, used in right_hand_side()
+
+private:
     mio::TimeSeries<ScalarType> m_flow_result; ///< flow result of the simulation
 };
 
 /**
- * @brief simulate simulates a compartmental model and returns the same results as simulate and also the flows.
- * @param[in] t0 start time
- * @param[in] tmax end time
- * @param[in] dt initial step size of integration
- * @param[in] model: An instance of a compartmental model
- * @return a Tuple of two TimeSeries to represent the final simulation result and flows
- * @tparam Model a compartment model type
- * @tparam Sim a simulation type that can simulate the model.
+ * @brief Run a FlowSimulation of a FlowModel.
+ * @param[in] t0 Start time.
+ * @param[in] tmax End time.
+ * @param[in] dt Initial step size of integration.
+ * @param[in] model An instance of a FlowModel.
+ * @param[in] integrator Optionally override the IntegratorCore used by the FlowSimulation.
+ * @return The simulation result as two TimeSeries. The first describes the compartments at each time point,
+ *         the second gives the corresponding flows that lead from t0 to each time point.
+ * @tparam Model The particular Model derived from FlowModel to simulate.
+ * @tparam Sim A FlowSimulation that can simulate the model.
  */
 template <class Model, class Sim = FlowSimulation<Model>>
 std::vector<TimeSeries<ScalarType>> simulate_flows(double t0, double tmax, double dt, Model const& model,

cpp/memilio/compartments/simulation.h

@@ -25,7 +25,6 @@
 #include "memilio/utils/metaprogramming.h"
 #include "memilio/math/stepper_wrapper.h"
 #include "memilio/utils/time_series.h"
-#include "memilio/math/euler.h"
 
 namespace mio
 {
@@ -189,14 +188,15 @@ using is_compartment_model_simulation =
                                   is_compartment_model<typename Sim::Model>::value)>;
 
 /**
- * @brief simulate simulates a compartmental model
- * @param[in] t0 start time
- * @param[in] tmax end time
- * @param[in] dt initial step size of integration
- * @param[in] model: An instance of a compartmental model
- * @return a TimeSeries to represent the final simulation result
- * @tparam Model a compartment model type
- * @tparam Sim a simulation type that can simulate the model.
+ * @brief Run a Simulation of a CompartmentalModel.
+ * @param[in] t0 Start time.
+ * @param[in] tmax End time.
+ * @param[in] dt Initial step size of integration.
+ * @param[in] model An instance of a CompartmentalModel.
+ * @param[in] integrator Optionally override the IntegratorCore used by the Simulation.
+ * @return A TimeSeries to represent the final Simulation result
+ * @tparam Model The particular Model derived from CompartmentModel to simulate.
+ * @tparam Sim A Simulation that can simulate the model.
  */
 template <class Model, class Sim = Simulation<Model>>
 TimeSeries<ScalarType> simulate(double t0, double tmax, double dt, Model const& model,

cpp/memilio/epidemiology/dynamic_npis.h

@@ -217,7 +217,7 @@ std::vector<size_t> get_damping_indices(const DampingExpr& damping_expr, Damping
 {
     std::vector<size_t> indices;
     for (size_t i = 0; i < damping_expr.get_dampings().size(); ++i) {
-        auto& d = damping_expr.get_dampings()[i];
+        const auto d = damping_expr.get_dampings()[i];
         if (d.get_level() == lvl && d.get_type() == type && double(d.get_time()) > double(begin) &&
             double(d.get_time()) < double(end)) {
             indices.push_back(i);

cpp/memilio/io/mobility_io.h

@@ -96,7 +96,7 @@ IOResult<void> write_graph(const Graph<Model, MigrationParameters>& graph, const
     //one file for the model (parameters and population)
     for (auto inode = size_t(0); inode < graph.nodes().size(); ++inode) {
         //node
-        auto& node = graph.nodes()[inode];
+        const auto node = graph.nodes()[inode];
         BOOST_OUTCOME_TRY(js_node_model, serialize_json(node.property, ioflags));
         Json::Value js_node(Json::objectValue);
         js_node["NodeId"]  = node.id;
@@ -133,7 +133,8 @@ IOResult<void> write_graph(const Graph<Model, MigrationParameters>& graph, const
  * @param read_edges boolean value that decides whether the edges of the graph should also be read in.
  */
 template <class Model>
-IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& directory, int ioflags = IOF_None, bool read_edges = true)
+IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& directory, int ioflags = IOF_None,
+                                                       bool read_edges = true)
 {
     std::string abs_path;
     if (!file_exists(directory, abs_path)) {
@@ -160,7 +161,7 @@ IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& direct
     }
 
     //read edges; nodes must already be available for that)
-    if(read_edges){
+    if (read_edges) {
         for (auto inode = size_t(0); inode < graph.nodes().size(); ++inode) {
             //list of edges
             auto edge_filename = path_join(abs_path, "GraphEdges_node" + std::to_string(inode) + ".json");
@@ -177,7 +178,7 @@ IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& direct
                 if (end_node_idx >= graph.nodes().size()) {
                     log_error("EndNodeIndex not in range of number of graph nodes.");
                     return failure(StatusCode::OutOfRange,
-                                edge_filename + ", EndNodeIndex not in range of number of graph nodes.");
+                                   edge_filename + ", EndNodeIndex not in range of number of graph nodes.");
                 }
                 BOOST_OUTCOME_TRY(parameters, deserialize_json(e["Parameters"], Tag<MigrationParameters>{}, ioflags));
                 graph.add_edge(start_node_idx, end_node_idx, parameters);

cpp/memilio/math/integrator.cpp

@@ -19,6 +19,7 @@
 */
 #include "memilio/math/integrator.h"
 #include "memilio/utils/logging.h"
+#include <cstddef>
 
 namespace mio
 {
@@ -30,32 +31,34 @@ Eigen::Ref<Eigen::VectorXd> OdeIntegrator::advance(const DerivFunction& f, const
     assert(tmax > t0);
     assert(dt > 0);
 
-    const size_t nb_steps = (int)(ceil((tmax - t0) / dt)); // estimated number of time steps (if equidistant)
+    const size_t num_steps =
+        static_cast<size_t>(ceil((tmax - t0) / dt)); // estimated number of time steps (if equidistant)
 
-    results.reserve(results.get_num_time_points() + nb_steps);
+    results.reserve(results.get_num_time_points() + num_steps);
 
     bool step_okay = true;
 
-    double t = t0;
-    size_t i = results.get_num_time_points() - 1;
+    double dt_restore = 0; // used to restore dt if dt was decreased to reach tmax
+    double t          = t0;
+    size_t i          = results.get_num_time_points() - 1;
     while (std::abs((tmax - t) / (tmax - t0)) > 1e-10) {
         //we don't make timesteps too small as the error estimator of an adaptive integrator
         //may not be able to handle it. this is very conservative and maybe unnecessary,
         //but also unlikely to happen. may need to be reevaluated
 
-        auto dt_eff = std::min(dt, tmax - t);
+        if (dt > tmax - t) {
+            dt_restore = dt;
+            dt         = tmax - t;
+        }
         results.add_time_point();
-        step_okay &= m_core->step(f, results[i], t, dt_eff, results[i + 1]);
+        step_okay &= m_core->step(f, results[i], t, dt, results[i + 1]);
         results.get_last_time() = t;
 
         ++i;
-
-        if (std::abs((tmax - t) / (tmax - t0)) > 1e-10 || dt_eff > dt) {
-            //store dt only if it's not the last step as it is probably smaller than required for tolerances
-            //except if the step function returns a bigger step size so as to not lose efficiency
-            dt = dt_eff;
-        }
     }
+    // if dt was decreased to reach tmax in the last time iteration,
+    // we restore it as it is now probably smaller than required for tolerances
+    dt = std::max(dt, dt_restore);
 
     if (!step_okay) {
         log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");