1229 stage aligned scheme implementation

cpp/examples/CMakeLists.txt

@@ -18,17 +18,21 @@ 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_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_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_graph_example ode_seir_graph.cpp)
+target_link_libraries(ode_seir_graph_example PRIVATE memilio ode_seir)
+target_compile_options(ode_seir_graph_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_seir_graph.cpp

@@ -0,0 +1,107 @@
+/*
+* Copyright (C) 2020-2026 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.
+*/
+
+/**
+ * Two-patch SEIR metapopulation with stage-aligned mobility returns.
+ *
+ * The stage-aligned approach does:
+ *   1. A SEIR-specific RK4 integrator (SeirRK4IntegratorCore) that caches the
+ *      force-of-infection at each stage during node integration.
+ *   2. Reuse the cached values to compute mobility returns with only apply_rhs calls - no
+ *      additional compute_foi work per edge.
+ *
+ */
+
+#include "ode_seir/model.h"
+#include "ode_seir/mobility.h"
+#include "memilio/mobility/metapopulation_mobility_instant.h"
+#include "memilio/compartments/simulation.h"
+
+#include <iostream>
+#include <iomanip>
+
+int main()
+{
+    using FP = double;
+
+    mio::set_log_level(mio::LogLevel::warn);
+
+    const FP t0   = 0.0;
+    const FP tmax = 30.0;
+    const FP dt   = 0.5;
+
+    // SEIR model with single age group
+    const size_t num_groups = 1;
+    mio::oseir::Model<FP> model(static_cast<int>(num_groups));
+
+    model.parameters.set<mio::oseir::TimeExposed<FP>>(5.2);
+    model.parameters.set<mio::oseir::TimeInfected<FP>>(6.0);
+    model.parameters.set<mio::oseir::TransmissionProbabilityOnContact<FP>>(0.1);
+
+    auto& cm = model.parameters.get<mio::oseir::ContactPatterns<FP>>();
+    cm.get_cont_freq_mat()[0].get_baseline().setConstant(2.7);
+
+    // Patch 1
+    auto model1                                                                     = model;
+    model1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9700;
+    model1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}]     = 200;
+    model1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}]    = 80;
+    model1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]   = 20;
+
+    // Patch 2
+    auto model2                                                                     = model;
+    model2.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9990;
+    model2.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}]     = 5;
+    model2.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}]    = 3;
+    model2.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]   = 2;
+
+    // Build mobility graph: 10% of each compartment commutes daily between the two patches
+    const auto NC =
+        static_cast<Eigen::Index>(mio::oseir::InfectionState::Count) * static_cast<Eigen::Index>(num_groups);
+
+    mio::Graph<mio::SimulationNode<FP, mio::Simulation<FP, mio::oseir::Model<FP>>>, mio::MobilityEdge<FP>> g;
+    g.add_node(1001, model1, t0);
+    g.add_node(1002, model2, t0);
+    g.add_edge(0, 1, Eigen::VectorXd::Constant(NC, 0.1));
+    g.add_edge(1, 0, Eigen::VectorXd::Constant(NC, 0.1));
+
+    // make_seir_mobility_sim sets SeirRK4IntegratorCore on every node as integrator,
+    // synchronizes node dt with graph dt, and returns a GraphSimulation.
+    auto sim = mio::oseir::make_seir_mobility_sim<FP>(t0, dt, std::move(g));
+
+    sim.advance(tmax);
+
+    // Print results
+    std::cout << std::fixed << std::setprecision(2);
+    std::cout << "\n=== Stage-Aligned SEIR Graph Simulation Results ===\n\n";
+
+    for (size_t n = 0; n < sim.get_graph().nodes().size(); ++n) {
+        auto& result = sim.get_graph().nodes()[n].property.get_result();
+        auto last    = result.get_last_value();
+        std::cout << "Patch " << n << " at t=" << result.get_last_time() << ":\n"
+                  << "  S = " << last[0] << "\n"
+                  << "  E = " << last[1] << "\n"
+                  << "  I = " << last[2] << "\n"
+                  << "  R = " << last[3] << "\n"
+                  << "  Total = " << last.sum() << "\n\n";
+    }
+
+    return 0;
+}

cpp/memilio/CMakeLists.txt

@@ -79,6 +79,7 @@ add_library(memilio
     mobility/metapopulation_mobility_stochastic.cpp
     mobility/graph_simulation.h
     mobility/graph_simulation.cpp
+    mobility/stage_aligned_mobility.h
     mobility/graph.h
     mobility/graph.cpp
     mobility/graph_builder.h

cpp/memilio/mobility/stage_aligned_mobility.h

@@ -0,0 +1,72 @@
+/*
+* Copyright (C) 2020-2026 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 MIO_MOBILITY_STAGE_ALIGNED_H
+#define MIO_MOBILITY_STAGE_ALIGNED_H
+
+#include "memilio/mobility/metapopulation_mobility_instant.h"
+#include "memilio/math/stepper_wrapper.h"
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+
+namespace mio
+{
+
+/**
+ * @brief Create a stage-aligned mobility simulation.
+ *
+ * This simulation is replacing the default adaptive integrator (Cash-Karp 5(4))
+ * in each node with a fixed-step explicit RK4 whose step size matches the graph dt.
+ * Combined with a model-specific calculate_mobility_returns overload that does
+ * RK4 co-integration (e.g., from ode_seir/mobility.h). For more details, see:
+ * 
+ * Zunker et al. 2026, https://doi.org/10.48550/arXiv.2603.11275
+ *
+ * @param t0    Start time of the simulation.
+ * @param dt    Time step between mobility events (also used as the RK4 step size).
+ * @param graph The mobility graph with SimulationNodes and MobilityEdges.
+ * @return A GraphSimulation configured for stage-aligned mobility updates.
+ * @{
+ */
+template <typename FP, class Sim>
+auto make_stage_aligned_mobility_sim(FP t0, FP dt, Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>& graph)
+{
+    // Switch all node integrators to fixed-step RK4
+    for (auto& node : graph.nodes()) {
+        node.property.get_simulation().set_integrator_core(
+            std::make_unique<ExplicitStepperWrapper<FP, boost::numeric::odeint::runge_kutta4>>());
+    }
+    return make_mobility_sim<FP>(t0, dt, graph);
+}
+
+template <typename FP, class Sim>
+auto make_stage_aligned_mobility_sim(FP t0, FP dt, Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>&& graph)
+{
+    // Switch all node integrators to fixed-step RK4
+    for (auto& node : graph.nodes()) {
+        node.property.get_simulation().set_integrator_core(
+            std::make_unique<ExplicitStepperWrapper<FP, boost::numeric::odeint::runge_kutta4>>());
+    }
+
+    return make_mobility_sim<FP>(t0, dt, std::move(graph));
+}
+/** @} */
+
+} // namespace mio
+
+#endif // MIO_MOBILITY_STAGE_ALIGNED_H

cpp/models/ode_seir/CMakeLists.txt

@@ -2,6 +2,7 @@ add_library(ode_seir
     infection_state.h
     model.h
     model.cpp
+    mobility.h
     parameters.h
 )
 target_link_libraries(ode_seir PUBLIC memilio)