410 Make IDE model compatible with graph (#1203)

annawendler · reneSchm · web-flow · commit 339076246600 · 2025-04-07T16:34:57.000+02:00
Add graph simulation that doesn't require mobility so that IDE model can be used within graph

Co-authored-by: reneSchm &lt;49305466+reneSchm@users.noreply.github.com&gt;
diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt
@@ -116,6 +116,10 @@ add_executable(ide_secir_ageres_example ide_secir_ageres.cpp)
 target_link_libraries(ide_secir_ageres_example PRIVATE memilio ide_secir)
 target_compile_options(ide_secir_ageres_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(ide_secir_graph_example ide_secir_graph.cpp)
+target_link_libraries(ide_secir_graph_example PRIVATE memilio ide_secir)
+target_compile_options(ide_secir_graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 add_executable(lct_secir_example lct_secir.cpp)
 target_link_libraries(lct_secir_example PRIVATE memilio lct_secir)
 target_compile_options(lct_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
diff --git a/cpp/examples/ide_secir_graph.cpp b/cpp/examples/ide_secir_graph.cpp
@@ -0,0 +1,103 @@
+/* 
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Anna Wendler
+*
+* 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 "ide_secir/infection_state.h"
+#include "ide_secir/model.h"
+#include "ide_secir/simulation.h"
+#include "memilio/config.h"
+#include "memilio/mobility/graph.h"
+#include "memilio/mobility/metapopulation_mobility_instant.h"
+
+int main()
+{
+    // This is an example for running the IDE-SECIR model in a graph. It demonstrates a simple setup for a graph without
+    // mobility that consists of two nodes and no edges.
+
+    using Vec = mio::TimeSeries<ScalarType>::Vector;
+
+    const ScalarType t0   = 0.;
+    const ScalarType tmax = 10.;
+
+    // Set time step size of IDE solver. Note that the time step size will be constant throughout the simulation. 
+    const ScalarType dt_ide_solver = 1.;
+
+    size_t num_agegroups = 1;
+
+    mio::CustomIndexArray<ScalarType, mio::AgeGroup> N_init =
+        mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 10000.);
+    mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths_init =
+        mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 13.10462213);
+
+    // Create TimeSeries with num_transitions * num_agegroups elements where initial transitions needed for simulation
+    // will be stored. We require values for the transitions for a sufficient number of time points before the start of
+    // the simulation to initialize our model.  
+    size_t num_transitions = (size_t)mio::isecir::InfectionTransition::Count;
+    mio::TimeSeries<ScalarType> transitions_init(num_transitions * num_agegroups);
+
+    // Define vector of transitions that will be added to the time points of the TimeSeries transitions_init.
+    Vec vec_init(num_transitions * num_agegroups);
+    vec_init[(size_t)mio::isecir::InfectionTransition::SusceptibleToExposed]                 = 25.;
+    vec_init[(size_t)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms]          = 15.;
+    vec_init[(size_t)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 8.;
+    vec_init[(size_t)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered]        = 4.;
+    vec_init[(size_t)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere]     = 1.;
+    vec_init[(size_t)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered]          = 4.;
+    vec_init[(size_t)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical]     = 1.;
+    vec_init[(size_t)mio::isecir::InfectionTransition::InfectedSevereToRecovered]            = 1.;
+    vec_init[(size_t)mio::isecir::InfectionTransition::InfectedCriticalToDead]               = 1.;
+    vec_init[(size_t)mio::isecir::InfectionTransition::InfectedCriticalToRecovered]          = 1.;
+    // Multiply vec_init with dt_ide_solver so that within a time interval of length 1, always the above number of 
+    // individuals are transitioning from one compartment to another, irrespective of the chosen time step size. 
+    vec_init = vec_init * dt_ide_solver;
+
+    // In this example, we use the default transition distributions. For these distributions, setting the initial time 
+    // point of the TimeSeries transitions_init at time -10 will give us a sufficient number of time points before t0=0.
+    // For more information on this, we refer to the documentation of TransitionDistributions in
+    // models/ide_secir/parameters.h.
+    transitions_init.add_time_point(-10, vec_init);
+    // Add further time points with distance dt_ide_solver until time t0.
+    while (transitions_init.get_last_time() < t0 - dt_ide_solver / 2.) {
+        transitions_init.add_time_point(transitions_init.get_last_time() + dt_ide_solver, vec_init);
+    }
+
+    // Initialize IDE model that will be used in in each node in the graph below. 
+    mio::isecir::Model model(std::move(transitions_init), N_init, deaths_init, num_agegroups);
+    model.check_constraints(dt_ide_solver);
+
+    // Set up graph with two nodes and no edges. To each node, we pass an id, the above constructed IDE model as well 
+    // as the time step size that will be used by the numerical solver. For simplicity, we use the same model in
+    // both nodes.
+    mio::Graph<mio::SimulationNode<mio::isecir::Simulation>, mio::MobilityEdge<ScalarType>> g;
+    g.add_node(1001, model, dt_ide_solver);
+    g.add_node(1002, model, dt_ide_solver);
+
+    // We use make_no_mobilty_sim to create a GraphSimulation for the model graph we defined above. This function allows 
+    // us to create a simulation without having to define mobility between nodes. Any edges would be effectively 
+    // ignored by the simulation.
+    auto sim = mio::make_no_mobility_sim(t0, std::move(g));
+    // Run the simulation. This advances all graph nodes independently. 
+    sim.advance(tmax);
+
+    // Print results of first node.
+    auto& node_0  = sim.get_graph().nodes()[0];
+    auto result_0 = node_0.property.get_result();
+    result_0.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
+
+    return 0;
+}
diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h
@@ -403,13 +403,13 @@ class MobilityEdge
 
 private:
     MobilityParameters<FP> m_parameters;
-    TimeSeries<double> m_mobile_population;
-    TimeSeries<double> m_return_times;
+    TimeSeries<FP> m_mobile_population;
+    TimeSeries<FP> m_return_times;
     bool m_return_mobile_population;
-    double m_t_last_dynamic_npi_check               = -std::numeric_limits<double>::infinity();
-    std::pair<double, SimulationTime> m_dynamic_npi = {-std::numeric_limits<double>::max(), SimulationTime(0)};
+    FP m_t_last_dynamic_npi_check               = -std::numeric_limits<FP>::infinity();
+    std::pair<FP, SimulationTime> m_dynamic_npi = {-std::numeric_limits<FP>::max(), SimulationTime(0)};
     std::vector<std::vector<size_t>> m_saved_compartment_indices; // groups of indices from compartments to save
-    TimeSeries<double> m_mobility_results; // save results from edges + entry for the total number of commuters
+    TimeSeries<FP> m_mobility_results; // save results from edges + entry for the total number of commuters
 
     /**
      * @brief Computes a condensed version of `m_mobile_population` and stores it in `m_mobility_results`.
@@ -419,11 +419,11 @@ class MobilityEdge
      *
      * @param[in] t The current time.
      */
-    void add_mobility_result_time_point(const double t);
+    void add_mobility_result_time_point(const FP t);
 };
 
 template <typename FP>
-void MobilityEdge<FP>::add_mobility_result_time_point(const double t)
+void MobilityEdge<FP>::add_mobility_result_time_point(const FP t)
 {
     const size_t save_indices_size = this->m_saved_compartment_indices.size();
     if (save_indices_size > 0) {
@@ -435,7 +435,7 @@ void MobilityEdge<FP>::add_mobility_result_time_point(const double t)
         std::transform(this->m_saved_compartment_indices.begin(), this->m_saved_compartment_indices.end(),
                        condensed_values.data(), [&last_value](const auto& indices) {
                            return std::accumulate(indices.begin(), indices.end(), 0.0,
-                                                  [&last_value](double sum, auto i) {
+                                                  [&last_value](FP sum, auto i) {
                                                       return sum + last_value[i];
                                                   });
                        });
@@ -563,7 +563,7 @@ template <class Sim>
 void MobilityEdge<FP>::apply_mobility(FP t, FP dt, SimulationNode<Sim>& node_from, SimulationNode<Sim>& node_to)
 {
     //check dynamic npis
-    if (m_t_last_dynamic_npi_check == -std::numeric_limits<double>::infinity()) {
+    if (m_t_last_dynamic_npi_check == -std::numeric_limits<FP>::infinity()) {
         m_t_last_dynamic_npi_check = node_from.get_t0();
     }
 
@@ -574,8 +574,8 @@ void MobilityEdge<FP>::apply_mobility(FP t, FP dt, SimulationNode<Sim>& node_fro
         auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel);
         if (exceeded_threshold != dyn_npis.get_thresholds().end() &&
             (exceeded_threshold->first > m_dynamic_npi.first ||
-             t > double(m_dynamic_npi.second))) { //old NPI was weaker or is expired
-            auto t_end    = SimulationTime(t + double(dyn_npis.get_duration()));
+             t > FP(m_dynamic_npi.second))) { //old NPI was weaker or is expired
+            auto t_end    = SimulationTime(t + FP(dyn_npis.get_duration()));
             m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end);
             implement_dynamic_npis(
                 m_parameters.get_coefficients(), exceeded_threshold->second, SimulationTime(t), t_end, [this](auto& g) {
@@ -673,8 +673,8 @@ void apply_mobility(FP t, FP dt, MobilityEdge<FP>& mobilityEdge, SimulationNode<
  */
 template <typename FP, class Sim>
 GraphSimulation<Graph<SimulationNode<Sim>, MobilityEdge<FP>>, FP, FP,
-                void (*)(double, double, mio::MobilityEdge<>&, mio::SimulationNode<Sim>&, mio::SimulationNode<Sim>&),
-                void (*)(double, double, mio::SimulationNode<Sim>&)>
+                void (*)(FP, FP, mio::MobilityEdge<FP>&, mio::SimulationNode<Sim>&, mio::SimulationNode<Sim>&),
+                void (*)(FP, FP, mio::SimulationNode<Sim>&)>
 make_mobility_sim(FP t0, FP dt, const Graph<SimulationNode<Sim>, MobilityEdge<FP>>& graph)
 {
     return make_graph_sim(t0, dt, graph, static_cast<void (*)(FP, FP, SimulationNode<Sim>&)>(&advance_model<Sim>),
@@ -684,8 +684,8 @@ make_mobility_sim(FP t0, FP dt, const Graph<SimulationNode<Sim>, MobilityEdge<FP
 
 template <typename FP, class Sim>
 GraphSimulation<Graph<SimulationNode<Sim>, MobilityEdge<FP>>, FP, FP,
-                void (*)(double, double, mio::MobilityEdge<>&, mio::SimulationNode<Sim>&, mio::SimulationNode<Sim>&),
-                void (*)(double, double, mio::SimulationNode<Sim>&)>
+                void (*)(FP, FP, mio::MobilityEdge<FP>&, mio::SimulationNode<Sim>&, mio::SimulationNode<Sim>&),
+                void (*)(FP, FP, mio::SimulationNode<Sim>&)>
 make_mobility_sim(FP t0, FP dt, Graph<SimulationNode<Sim>, MobilityEdge<FP>>&& graph)
 {
     return make_graph_sim(t0, dt, std::move(graph),
@@ -696,6 +696,40 @@ make_mobility_sim(FP t0, FP dt, Graph<SimulationNode<Sim>, MobilityEdge<FP>>&& g
 
 /** @} */
 
+/**
+ * Create a graph simulation without mobility.
+ * 
+ * Note that we set the time step of the graph simulation to infinity since we do not require any exchange between the 
+ * nodes. Hence, in each node, the simulation runs until tmax when advancing the simulation without interruption.
+ *
+ * @param t0 Start time of the simulation.
+ * @param graph Set up for graph-based simulation.
+ * @{
+ */
+ template <typename FP, class Sim>
+ auto make_no_mobility_sim(FP t0, Graph<SimulationNode<Sim>, MobilityEdge<FP>>& graph)
+ {
+     using GraphSim =
+         GraphSimulation<Graph<SimulationNode<Sim>, MobilityEdge<FP>>, FP, FP,
+                         void (*)(FP, FP, mio::MobilityEdge<FP>&, mio::SimulationNode<Sim>&, mio::SimulationNode<Sim>&),
+                         void (*)(FP, FP, mio::SimulationNode<Sim>&)>;
+     return GraphSim(t0, std::numeric_limits<FP>::infinity(), graph, &advance_model<Sim>,
+                     [](FP, FP, MobilityEdge<FP>&, SimulationNode<Sim>&, SimulationNode<Sim>&) {});
+ }
+
+template <typename FP, class Sim>
+auto make_no_mobility_sim(FP t0, Graph<SimulationNode<Sim>, MobilityEdge<FP>>&& graph)
+{
+    using GraphSim =
+        GraphSimulation<Graph<SimulationNode<Sim>, MobilityEdge<FP>>, FP, FP,
+                        void (*)(FP, FP, mio::MobilityEdge<FP>&, mio::SimulationNode<Sim>&, mio::SimulationNode<Sim>&),
+                        void (*)(FP, FP, mio::SimulationNode<Sim>&)>;
+    return GraphSim(t0, std::numeric_limits<FP>::infinity(), std::move(graph), &advance_model<Sim>,
+                    [](FP, FP, MobilityEdge<FP>&, SimulationNode<Sim>&, SimulationNode<Sim>&) {});
+}
+
+/** @} */
+
 } // namespace mio
 
 #endif //METAPOPULATION_MOBILITY_INSTANT_H
diff --git a/cpp/models/ide_secir/parameters.h b/cpp/models/ide_secir/parameters.h
@@ -47,6 +47,10 @@ namespace isecir
 /**
  * @brief Transition distribution for each transition in #InfectionTransition.
  *
+ * For each transition, the corresponding transition distribution can be chosen independently. 
+ * The choice of distributions determines how many initial time points are required to initialize the model, see
+ * get_global_support_max() in models/ide_secir/model.h.
+ *
  * As a default we use SmootherCosine functions for all transitions with m_parameter=2.
  */
 struct TransitionDistributions {
diff --git a/cpp/models/ide_secir/simulation.cpp b/cpp/models/ide_secir/simulation.cpp
@@ -32,9 +32,7 @@ void Simulation::advance(ScalarType tmax)
 {
     mio::log_info("Simulating IDE-SECIR from t0 = {} until tmax = {} with dt = {}.",
                   m_model->transitions.get_last_time(), tmax, m_dt);
-    m_model->set_transitiondistributions_support_max(m_dt);
-    m_model->set_transitiondistributions_derivative(m_dt);
-    m_model->set_transitiondistributions_in_forceofinfection(m_dt);
+
     m_model->initial_compute_compartments(m_dt);
 
     // For every time step:
diff --git a/cpp/models/ide_secir/simulation.h b/cpp/models/ide_secir/simulation.h
@@ -41,12 +41,16 @@ class Simulation
     /**
      * @brief setup the Simulation for an IDE model.
      * @param[in] model An instance of the IDE model.
-     * @param[in] dt Step size of numerical solver.
+     * @param[in] dt Step size of numerical solver. Throughout the simulation, the step size will be constant. 
      */
     Simulation(Model const& model, ScalarType dt = 0.1)
         : m_model(std::make_unique<Model>(model))
         , m_dt(dt)
     {
+        assert(m_dt > 0);
+        m_model->set_transitiondistributions_support_max(m_dt);
+        m_model->set_transitiondistributions_derivative(m_dt);
+        m_model->set_transitiondistributions_in_forceofinfection(m_dt);
     }
 
     /** 
diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp
