From d03aa5b641c7dbad93e204b1ed13bfca4624a84a Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 25 Mar 2026 22:01:55 +0100 Subject: [PATCH 1/3] stage aligned approach idea implementation --- cpp/examples/CMakeLists.txt | 26 +- cpp/memilio/CMakeLists.txt | 1 + cpp/memilio/mobility/stage_aligned_mobility.h | 72 +++ cpp/models/ode_seir/CMakeLists.txt | 1 + cpp/models/ode_seir/mobility.h | 419 ++++++++++++++++++ 5 files changed, 508 insertions(+), 11 deletions(-) create mode 100644 cpp/memilio/mobility/stage_aligned_mobility.h create mode 100644 cpp/models/ode_seir/mobility.h diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index e63da90ac1..4d43d41b48 100644 --- a/cpp/examples/CMakeLists.txt +++ b/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) diff --git a/cpp/memilio/CMakeLists.txt b/cpp/memilio/CMakeLists.txt index 658733a7f8..9095dd1de3 100644 --- a/cpp/memilio/CMakeLists.txt +++ b/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 diff --git a/cpp/memilio/mobility/stage_aligned_mobility.h b/cpp/memilio/mobility/stage_aligned_mobility.h new file mode 100644 index 0000000000..a7b3c1bd0e --- /dev/null +++ b/cpp/memilio/mobility/stage_aligned_mobility.h @@ -0,0 +1,72 @@ +/* +* Copyright (C) 2020-2026 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* 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 + +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 +auto make_stage_aligned_mobility_sim(FP t0, FP dt, Graph, MobilityEdge>& graph) +{ + // Switch all node integrators to fixed-step RK4 + for (auto& node : graph.nodes()) { + node.property.get_simulation().set_integrator_core( + std::make_unique>()); + } + return make_mobility_sim(t0, dt, graph); +} + +template +auto make_stage_aligned_mobility_sim(FP t0, FP dt, Graph, MobilityEdge>&& graph) +{ + // Switch all node integrators to fixed-step RK4 + for (auto& node : graph.nodes()) { + node.property.get_simulation().set_integrator_core( + std::make_unique>()); + } + + return make_mobility_sim(t0, dt, std::move(graph)); +} +/** @} */ + +} // namespace mio + +#endif // MIO_MOBILITY_STAGE_ALIGNED_H diff --git a/cpp/models/ode_seir/CMakeLists.txt b/cpp/models/ode_seir/CMakeLists.txt index c94c6b0aa7..0d1758b6c0 100644 --- a/cpp/models/ode_seir/CMakeLists.txt +++ b/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) diff --git a/cpp/models/ode_seir/mobility.h b/cpp/models/ode_seir/mobility.h new file mode 100644 index 0000000000..31e327cb98 --- /dev/null +++ b/cpp/models/ode_seir/mobility.h @@ -0,0 +1,419 @@ +/* +* Copyright (C) 2020-2026 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* 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 OSEIR_MOBILITY_H +#define OSEIR_MOBILITY_H + +#include "ode_seir/model.h" +#include "memilio/compartments/flow_simulation.h" +#include "memilio/compartments/simulation.h" +#include "memilio/math/integrator.h" +#include "memilio/math/eigen.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/utils/time_series.h" +#include +#include +#include +#include +#include + +namespace mio +{ +namespace oseir +{ + +/** + * @brief Pre-computed SEIR rate evaluator for efficient stage-aligned mobility. + * + * Exploits the factorization f(x) = B · D(z,t) · x (cf. Zunker et al. 2026, + * https://doi.org/10.48550/arXiv.2603.11275) where the rate matrix D depends + * only on the total (aggregated) population z, not on the individual commuter + * subpopulation. + * + * @tparam FP Floating point type. + */ +template +struct SeirCoefficientEvaluator { + const Model& model; + const size_t NG; ///< Number of age groups. + const size_t NC; ///< Total number of compartments (= 4 · NG for SEIR). + + // Pre-cached constant transition rates (independent of state z) + std::vector rate_E; ///< 1 / TimeExposed per age group. + std::vector rate_I; ///< 1 / TimeInfected per age group. + std::vector beta; ///< TransmissionProbabilityOnContact per age group. + + // Pre-cached flat population indices + std::vector idxS, idxE, idxI, idxR; + + explicit SeirCoefficientEvaluator(const Model& m) + : model(m) + , NG(static_cast(m.parameters.get_num_groups())) + , NC(static_cast(m.populations.get_num_compartments())) + , rate_E(NG) + , rate_I(NG) + , beta(NG) + , idxS(NG) + , idxE(NG) + , idxI(NG) + , idxR(NG) + { + for (size_t g = 0; g < NG; ++g) { + auto ag = AgeGroup(static_cast(g)); + idxS[g] = m.populations.get_flat_index({ag, InfectionState::Susceptible}); + idxE[g] = m.populations.get_flat_index({ag, InfectionState::Exposed}); + idxI[g] = m.populations.get_flat_index({ag, InfectionState::Infected}); + idxR[g] = m.populations.get_flat_index({ag, InfectionState::Recovered}); + rate_E[g] = FP(1) / m.parameters.template get>()[ag]; + rate_I[g] = FP(1) / m.parameters.template get>()[ag]; + beta[g] = m.parameters.template get>()[ag]; + } + } + + /** + * @brief Compute force of infection lambda_i from total population z. + * + * @param[in] z Total population vector (length NC). + * @param[in] t Current time. + * @param[out] lambda Force of infection per age group (length NG). + */ + void compute_foi(const Eigen::VectorX& z, FP t, std::vector& lambda) const + { + std::fill(lambda.begin(), lambda.end(), FP(0)); + auto cm = model.parameters.template get>().get_cont_freq_mat().get_matrix_at( + SimulationTime(t)); + for (size_t j = 0; j < NG; ++j) { + const FP Nj = z[idxS[j]] + z[idxE[j]] + z[idxI[j]] + z[idxR[j]]; + const FP Ij_divNj = (Nj < Limits::zero_tolerance()) ? FP(0) : z[idxI[j]] / Nj; + for (size_t i = 0; i < NG; ++i) { + lambda[i] += cm(static_cast(i), static_cast(j)) * beta[i] * Ij_divNj; + } + } + } + + /** + * @brief Apply SEIR dynamics to any subpopulation x using pre-computed lambda. + * + * @param[in] lambda Pre-computed force of infection (from compute_foi). + * @param[in] x Population vector (totals or commuters, length NC). + * @param[out] dxdt Derivative output vector (length NC). + */ + void apply_rhs(const std::vector& lambda, const Eigen::VectorX& x, Eigen::VectorX& dxdt) const + { + dxdt.setZero(); + for (size_t g = 0; g < NG; ++g) { + const FP flow_SE = lambda[g] * x[idxS[g]]; + const FP flow_EI = rate_E[g] * x[idxE[g]]; + const FP flow_IR = rate_I[g] * x[idxI[g]]; + dxdt[idxS[g]] = -flow_SE; + dxdt[idxE[g]] = flow_SE - flow_EI; + dxdt[idxI[g]] = flow_EI - flow_IR; + dxdt[idxR[g]] = flow_IR; + } + } +}; + +/** + * @brief SEIR-specific fixed-step RK4 integrator that caches lambda at every stage. + * + * During each `step()`, the classical 4-stage RK4 is performed using + * SeirCoefficientEvaluator instead of the generic DerivFunction. + * The force-of-infection lambda_s computed at each stage is stored + * in an internal cache. + * + * After the node integration, `integrate_commuters()` reuses these cached + * values to advance commuter subpopulations without any additional compute_foi + * calls. + * + * The DerivFunction parameter in step() is ignored. The integrator uses its + * own SeirCoefficientEvaluator which must match the model attached to the + * simulation. Construct with the same model reference. + * + * @tparam FP Floating point type. + */ +template +class SeirRK4IntegratorCore : public OdeIntegratorCore +{ +public: + /// Cached lambda-vectors at the 4 RK stages. + struct StageCache { + std::array, 4> lambda; ///< lambda_s for stages s=0..3. + FP t0 = FP(0); ///< Start time of the cached step. + FP dt = FP(0); ///< Step size of the cached step. + bool valid = false; ///< True after a successful step(). + }; + + /** + * @brief Construct from an oseir::Model. + * @param model The SEIR model whose parameters define the RHS. + * Must outlive this integrator. + */ + explicit SeirRK4IntegratorCore(const Model& model) + : OdeIntegratorCore(FP{}, FP{}) // fixed-step: no dt bounds + , m_eval(model) + { + const auto NC = static_cast(m_eval.NC); + m_k1.resize(NC); + m_k2.resize(NC); + m_k3.resize(NC); + m_k4.resize(NC); + m_tmp.resize(NC); + for (auto& l : m_cache.lambda) { + l.resize(m_eval.NG); + } + } + + std::unique_ptr> clone() const override + { + return std::make_unique(*this); + } + + /** + * @brief One fixed-step explicit RK4 integration step with lambda-caching. + */ + bool step(const DerivFunction& /*f*/, Eigen::Ref> yt, FP& t, FP& dt, + Eigen::Ref> ytp1) const override + { + const FP half_dt = FP(0.5) * dt; + + // Stage 1 (t) + m_eval.compute_foi(yt, t, m_cache.lambda[0]); + m_eval.apply_rhs(m_cache.lambda[0], yt, m_k1); + + // Stage 2 (t + dt/2) + m_tmp.noalias() = yt + half_dt * m_k1; + m_eval.compute_foi(m_tmp, t + half_dt, m_cache.lambda[1]); + m_eval.apply_rhs(m_cache.lambda[1], m_tmp, m_k2); + + // Stage 3 (t + dt/2) + m_tmp.noalias() = yt + half_dt * m_k2; + m_eval.compute_foi(m_tmp, t + half_dt, m_cache.lambda[2]); + m_eval.apply_rhs(m_cache.lambda[2], m_tmp, m_k3); + + // Stage 4 (t + dt) + m_tmp.noalias() = yt + dt * m_k3; + m_eval.compute_foi(m_tmp, t + dt, m_cache.lambda[3]); + m_eval.apply_rhs(m_cache.lambda[3], m_tmp, m_k4); + + // Final update + ytp1.noalias() = yt + (dt / FP(6)) * (m_k1 + FP(2) * m_k2 + FP(2) * m_k3 + m_k4); + + // Fill cache + m_cache.t0 = t; + m_cache.dt = dt; + m_cache.valid = true; + + t += dt; + return true; + } + + /** + * @brief Integrate a commuter subpopulation using cached lambda. + * + * @param[in,out] c Commuter state vector; updated in-place to t + dt. + * @param[in] dt Step size (must match the cached step size). + */ + void integrate_commuters(Eigen::Ref::Vector> c, FP dt) const + { + assert(m_cache.valid && "integrate_commuters called without a preceding step()"); + const FP half_dt = FP(0.5) * dt; + + // Stage 1: apply cached lambda[0] to commuters + m_eval.apply_rhs(m_cache.lambda[0], c, m_k1); + + // Stage 2: apply cached lambda[1] + m_tmp.noalias() = c + half_dt * m_k1; + m_eval.apply_rhs(m_cache.lambda[1], m_tmp, m_k2); + + // Stage 3: apply cached lambda[2] + m_tmp.noalias() = c + half_dt * m_k2; + m_eval.apply_rhs(m_cache.lambda[2], m_tmp, m_k3); + + // Stage 4: apply cached lambda[3] + m_tmp.noalias() = c + dt * m_k3; + m_eval.apply_rhs(m_cache.lambda[3], m_tmp, m_k4); + + // Final update + c += (dt / FP(6)) * (m_k1 + FP(2) * m_k2 + FP(2) * m_k3 + m_k4); + } + + /// @brief Whether the lambda-cache is valid (filled by the last step()). + bool has_valid_cache() const + { + return m_cache.valid; + } + + /// @brief Access the stage cache (e.g., for testing). + const StageCache& get_cache() const + { + return m_cache; + } + + /// @brief Access the evaluator (e.g., for testing). + const SeirCoefficientEvaluator& get_evaluator() const + { + return m_eval; + } + +private: + SeirCoefficientEvaluator m_eval; + mutable StageCache m_cache; + mutable Eigen::VectorX m_k1, m_k2, m_k3, m_k4, m_tmp; +}; + +namespace detail +{ + +/** + * @brief Stage-aligned RK4 co-integration of commuter and total population. + * + * Fallback path when the node does NOT use SeirRK4IntegratorCore. + */ +template +void calculate_mobility_returns_seir(const SeirCoefficientEvaluator& eval, + Eigen::Ref::Vector> mobile_population, + Eigen::Ref::Vector> total, FP t, FP dt) +{ + const size_t NC = eval.NC; + const size_t NG = eval.NG; + + Eigen::VectorX z = total; + Eigen::VectorX c = mobile_population; + + Eigen::VectorX k1_z(NC), k2_z(NC), k3_z(NC), k4_z(NC); + Eigen::VectorX k1_c(NC), k2_c(NC), k3_c(NC), k4_c(NC); + Eigen::VectorX z_s(NC), c_s(NC); + std::vector lambda(NG); + + const FP half_dt = FP(0.5) * dt; + + // ---- Stage 1 (t) ---- + eval.compute_foi(z, t, lambda); + eval.apply_rhs(lambda, z, k1_z); + eval.apply_rhs(lambda, c, k1_c); + + // ---- Stage 2 (t + dt/2) ---- + z_s.noalias() = z + half_dt * k1_z; + c_s.noalias() = c + half_dt * k1_c; + eval.compute_foi(z_s, t + half_dt, lambda); + eval.apply_rhs(lambda, z_s, k2_z); + eval.apply_rhs(lambda, c_s, k2_c); + + // ---- Stage 3 (t + dt/2) ---- + z_s.noalias() = z + half_dt * k2_z; + c_s.noalias() = c + half_dt * k2_c; + eval.compute_foi(z_s, t + half_dt, lambda); + eval.apply_rhs(lambda, z_s, k3_z); + eval.apply_rhs(lambda, c_s, k3_c); + + // ---- Stage 4 (t + dt) ---- + z_s.noalias() = z + dt * k3_z; + c_s.noalias() = c + dt * k3_c; + eval.compute_foi(z_s, t + dt, lambda); + eval.apply_rhs(lambda, z_s, k4_z); + eval.apply_rhs(lambda, c_s, k4_c); + + // ---- Final update (only commuters) ---- + mobile_population = c + (dt / FP(6)) * (k1_c + FP(2) * k2_c + FP(2) * k3_c + k4_c); +} + +} // namespace detail + +/** + * @brief Stage-aligned mobility returns for SEIR Simulation nodes. + * + * Ideal (SeirRK4IntegratorCore on node): + * Uses cached lambda from the node's integration -> only 4× apply_rhs, + * zero compute_foi calls, regardless of how many edges return to this node. + * + * Fallback (any other integrator): + * Full co-integration with 4× compute_foi + 8× apply_rhs. + */ +template +void calculate_mobility_returns(Eigen::Ref::Vector> mobile_population, + const Simulation>& sim, + Eigen::Ref::Vector> total, FP t, FP dt) +{ + // Try ideal path: reuse cached lambda from node's SeirRK4IntegratorCore + const auto* seir_core = dynamic_cast*>(&sim.get_integrator_core()); + if (seir_core != nullptr && seir_core->has_valid_cache()) { + seir_core->integrate_commuters(mobile_population, dt); + return; + } + + // Fallback: full co-integration (e.g. when using adaptive integrator) + SeirCoefficientEvaluator eval(sim.get_model()); + detail::calculate_mobility_returns_seir(eval, mobile_population, total, t, dt); +} + +/** + * @brief Stage-aligned mobility returns for SEIR FlowSimulation nodes. + * @see Simulation overload above. + */ +template +void calculate_mobility_returns(Eigen::Ref::Vector> mobile_population, + const FlowSimulation>& sim, + Eigen::Ref::Vector> total, FP t, FP dt) +{ + // FlowSimulation: currently no caching (integrator is different), always fallback + SeirCoefficientEvaluator eval(sim.get_model()); + detail::calculate_mobility_returns_seir(eval, mobile_population, total, t, dt); +} + +/** + * @brief Create a stage-aligned SEIR mobility simulation with lambda-caching. + * + * Sets SeirRK4IntegratorCore on every node and synchronizes the node's + * internal dt with the graph dt. + * + * @param t0 Start time. + * @param dt Graph time step (= RK4 step size = mobility interval). + * @param graph Mobility graph (will be moved). + * @return A GraphSimulation ready to advance. + */ +template +auto make_seir_mobility_sim(FP t0, FP dt, + Graph>>, MobilityEdge>&& graph) +{ + for (auto& node : graph.nodes()) { + auto& sim = node.property.get_simulation(); + // Set the custom SEIR RK4 integrator with lambda-caching + sim.set_integrator_core(std::make_unique>(sim.get_model())); + // Synchronize node dt with graph dt so exactly one RK4 step is taken per graph step + sim.get_dt() = dt; + } + return make_mobility_sim(t0, dt, std::move(graph)); +} + +/// @overload For lvalue graphs. +template +auto make_seir_mobility_sim(FP t0, FP dt, Graph>>, MobilityEdge>& graph) +{ + for (auto& node : graph.nodes()) { + auto& sim = node.property.get_simulation(); + sim.set_integrator_core(std::make_unique>(sim.get_model())); + sim.get_dt() = dt; + } + return make_mobility_sim(t0, dt, graph); +} + +} // namespace oseir +} // namespace mio + +#endif // OSEIR_MOBILITY_H From e3535b37c5087fe7972170ec0f7948fbb4055e4b Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 25 Mar 2026 22:04:27 +0100 Subject: [PATCH 2/3] 2 patch example --- cpp/examples/ode_seir_graph.cpp | 107 ++++++++++++++++++++++++++++++++ 1 file changed, 107 insertions(+) create mode 100644 cpp/examples/ode_seir_graph.cpp diff --git a/cpp/examples/ode_seir_graph.cpp b/cpp/examples/ode_seir_graph.cpp new file mode 100644 index 0000000000..2d7de86687 --- /dev/null +++ b/cpp/examples/ode_seir_graph.cpp @@ -0,0 +1,107 @@ +/* +* Copyright (C) 2020-2026 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* 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 +#include + +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 model(static_cast(num_groups)); + + model.parameters.set>(5.2); + model.parameters.set>(6.0); + model.parameters.set>(0.1); + + auto& cm = model.parameters.get>(); + 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(mio::oseir::InfectionState::Count) * static_cast(num_groups); + + mio::Graph>>, mio::MobilityEdge> 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(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; +} From 7a64d501af6aed06eeed234535f0f84be4d344c0 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 25 Mar 2026 22:10:32 +0100 Subject: [PATCH 3/3] tests --- cpp/tests/CMakeLists.txt | 1 + cpp/tests/test_oseir_mobility.cpp | 300 ++++++++++++++++++++++++++++++ 2 files changed, 301 insertions(+) create mode 100644 cpp/tests/test_oseir_mobility.cpp diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 6447c6eabc..b747cdc539 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -21,6 +21,7 @@ set(TESTSOURCES test_sde_sirs.cpp test_sde_seirvv.cpp test_mobility.cpp + test_oseir_mobility.cpp test_date.cpp test_eigen_util.cpp test_odesecir_ageres.cpp diff --git a/cpp/tests/test_oseir_mobility.cpp b/cpp/tests/test_oseir_mobility.cpp new file mode 100644 index 0000000000..8197d7c666 --- /dev/null +++ b/cpp/tests/test_oseir_mobility.cpp @@ -0,0 +1,300 @@ +/* +* Copyright (C) 2020-2026 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* 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/compartments/simulation.h" +#include "memilio/compartments/flow_simulation.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/math/stepper_wrapper.h" +#include "ode_seir/model.h" +#include "ode_seir/mobility.h" +#include +#include +#include + +namespace +{ + +mio::oseir::Model create_seir_model(size_t num_agegroups, double total_pop = 10000.0) +{ + mio::oseir::Model model(static_cast(num_agegroups)); + + for (size_t g = 0; g < num_agegroups; ++g) { + auto ag = mio::AgeGroup(static_cast(g)); + model.populations[{ag, mio::oseir::InfectionState::Exposed}] = 50.0; + model.populations[{ag, mio::oseir::InfectionState::Infected}] = 30.0; + model.populations[{ag, mio::oseir::InfectionState::Recovered}] = 20.0; + model.populations[{ag, mio::oseir::InfectionState::Susceptible}] = + total_pop / num_agegroups - 50.0 - 30.0 - 20.0; + } + + model.parameters.set>(5.2); + model.parameters.set>(6.0); + model.parameters.set>(0.1); + + auto& cm = model.parameters.get>(); + cm.get_cont_freq_mat()[0].get_baseline().setConstant(2.7); + + return model; +} + +} // namespace + +TEST(TestOseirMobility, noMobilityMatchesSingleSim) +{ + const double t0 = 0.0; + const double tmax = 10.0; + const double dt = 0.5; + + auto model = create_seir_model(1); + + // Graph simulation with zero mobility coefficients + mio::Graph>>, + mio::MobilityEdge> + g; + g.add_node(0, model, t0); + g.add_node(1, model, t0); + g.add_edge(0, 1, Eigen::VectorXd::Constant(4, 0.0)); + g.add_edge(1, 0, Eigen::VectorXd::Constant(4, 0.0)); + + auto graph_sim = mio::oseir::make_seir_mobility_sim(t0, dt, std::move(g)); + graph_sim.advance(tmax); + + // Single simulation with SeirRK4IntegratorCore (same as used in graph) + auto single_sim = mio::Simulation>(model, t0); + single_sim.set_integrator_core(std::make_unique>(single_sim.get_model())); + single_sim.get_dt() = dt; + single_sim.advance(tmax); + + auto diff = (graph_sim.get_graph().nodes()[0].property.get_result().get_last_value() - + single_sim.get_result().get_last_value()) + .norm(); + EXPECT_NEAR(diff, 0.0, 1e-10); +} + +TEST(TestOseirMobility, mobilityReturnsPreservation) +{ + const double t0 = 0.0; + const double tmax = 5.0; + const double dt = 0.5; + + auto model1 = create_seir_model(1, 10000.0); + auto model2 = create_seir_model(1, 5000.0); + + mio::Graph>>, + mio::MobilityEdge> + g; + g.add_node(0, model1, t0); + g.add_node(1, model2, t0); + g.add_edge(0, 1, Eigen::VectorXd::Constant(4, 0.05)); + g.add_edge(1, 0, Eigen::VectorXd::Constant(4, 0.03)); + + auto graph_sim = mio::oseir::make_seir_mobility_sim(t0, dt, std::move(g)); + graph_sim.advance(tmax); + + double total_pop = 0.0; + for (auto& n : graph_sim.get_graph().nodes()) { + total_pop += n.property.get_result().get_last_value().sum(); + } + EXPECT_NEAR(total_pop, 15000.0, 1e-6); +} + +TEST(TestOseirMobility, seirIntegratorMatchesBoostRK4) +{ + const double t0 = 0.0; + const double tmax = 10.0; + const double dt = 0.5; + + auto model = create_seir_model(1); + + // Run 1: SeirRK4IntegratorCore + auto sim1 = mio::Simulation>(model, t0); + sim1.set_integrator_core(std::make_unique>(sim1.get_model())); + sim1.get_dt() = dt; + sim1.advance(tmax); + + // Run 2: Standard ExplicitStepperWrapper + auto sim2 = mio::Simulation>(model, t0); + sim2.set_integrator_core( + std::make_unique>()); + sim2.get_dt() = dt; + sim2.advance(tmax); + + auto diff = (sim1.get_result().get_last_value() - sim2.get_result().get_last_value()).norm(); + EXPECT_NEAR(diff, 0.0, 1e-10) << "SeirRK4IntegratorCore must match boost runge_kutta4.\n" + << "Custom: " << sim1.get_result().get_last_value().transpose() << "\n" + << "Boost: " << sim2.get_result().get_last_value().transpose(); +} + +TEST(TestOseirMobility, multiAgeGroupMobilityReturns) +{ + const double t0 = 0.0; + const double tmax = 5.0; + const double dt = 0.5; + const size_t NG = 3; + + auto model = create_seir_model(NG, 30000.0); + + mio::Graph>>, + mio::MobilityEdge> + g; + g.add_node(0, model, t0); + g.add_node(1, model, t0); + + const size_t NC = 4 * NG; + g.add_edge(0, 1, Eigen::VectorXd::Constant(NC, 0.05)); + g.add_edge(1, 0, Eigen::VectorXd::Constant(NC, 0.05)); + + auto graph_sim = mio::oseir::make_seir_mobility_sim(t0, dt, std::move(g)); + graph_sim.advance(tmax); + + double total_pop = 0.0; + for (auto& n : graph_sim.get_graph().nodes()) { + total_pop += n.property.get_result().get_last_value().sum(); + } + EXPECT_NEAR(total_pop, 60000.0, 1e-4); + + // Check non-negativity + for (auto& n : graph_sim.get_graph().nodes()) { + auto v = n.property.get_result().get_last_value(); + for (int i = 0; i < v.size(); ++i) { + EXPECT_GE(v[i], -1e-12) << "Negative compartment at index " << i; + } + } +} + +TEST(TestOseirMobility, evaluatorMatchesGetDerivatives) +{ + auto model = create_seir_model(2, 20000.0); + const double t = 3.5; + + Eigen::VectorXd z(8); + z << 8000, 500, 300, 1200, 7000, 600, 400, 2000; + + Eigen::VectorXd dydt_ref(8); + model.get_derivatives(z, z, t, dydt_ref); + + mio::oseir::SeirCoefficientEvaluator eval(model); + std::vector lambda(2); + eval.compute_foi(z, t, lambda); + Eigen::VectorXd dydt_eval(8); + eval.apply_rhs(lambda, z, dydt_eval); + + EXPECT_NEAR((dydt_ref - dydt_eval).norm(), 0.0, 1e-12) + << "SeirCoefficientEvaluator must match model.get_derivatives for totals.\n" + << "Reference: " << dydt_ref.transpose() << "\n" + << "Evaluator: " << dydt_eval.transpose(); +} + +TEST(TestOseirMobility, evaluatorMatchesSubpopDerivatives) +{ + auto model = create_seir_model(2, 20000.0); + const double t = 1.0; + + Eigen::VectorXd z(8); + z << 8000, 500, 300, 1200, 7000, 600, 400, 2000; + + Eigen::VectorXd c(8); + c << 200, 15, 10, 30, 180, 20, 12, 50; + + Eigen::VectorXd dydt_ref(8); + model.get_derivatives(z, c, t, dydt_ref); + + mio::oseir::SeirCoefficientEvaluator eval(model); + std::vector lambda(2); + eval.compute_foi(z, t, lambda); + Eigen::VectorXd dydt_eval(8); + eval.apply_rhs(lambda, c, dydt_eval); + + EXPECT_NEAR((dydt_ref - dydt_eval).norm(), 0.0, 1e-12) + << "SeirCoefficientEvaluator must match model.get_derivatives for subpopulations.\n" + << "Reference: " << dydt_ref.transpose() << "\n" + << "Evaluator: " << dydt_eval.transpose(); +} + +TEST(TestOseirMobility, cachedReturnsMatchFallback) +{ + const double t0 = 0.0; + const double dt = 0.5; + + auto model = create_seir_model(1); + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 200.0; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000.0 - 50.0 - 200.0 - 20.0; + + Eigen::VectorXd total(4); + total << 9730, 50, 200, 20; + + Eigen::VectorXd commuters(4); + commuters << 500, 3, 10, 2; + + // (a) Cached path: step() fills cache, then integrate_commuters + mio::oseir::SeirRK4IntegratorCore integrator(model); + Eigen::VectorXd ytp1(4); + double t_step = t0; + double dt_step = dt; + integrator.step(mio::DerivFunction{}, total, t_step, dt_step, ytp1); + ASSERT_TRUE(integrator.has_valid_cache()); + + Eigen::VectorXd c_cached = commuters; + integrator.integrate_commuters(c_cached, dt); + + // (b) Fallback path: full co-integration + mio::oseir::SeirCoefficientEvaluator eval(model); + Eigen::VectorXd c_fallback = commuters; + mio::oseir::detail::calculate_mobility_returns_seir(eval, c_fallback, total, t0, dt); + + EXPECT_NEAR((c_cached - c_fallback).norm(), 0.0, 1e-12) + << "Cached path must exactly match fallback co-integration.\n" + << "Cached: " << c_cached.transpose() << "\n" + << "Fallback: " << c_fallback.transpose(); +} + +TEST(TestOseirMobility, fullGraphSimNonNegative) +{ + const double t0 = 0.0; + const double tmax = 15.0; + const double dt = 0.5; + + auto model = create_seir_model(1); + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 200.0; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000.0 - 50.0 - 200.0 - 20.0; + + mio::Graph>>, + mio::MobilityEdge> + g; + g.add_node(0, model, t0); + g.add_node(1, model, t0); + g.add_edge(0, 1, Eigen::VectorXd::Constant(4, 0.1)); + g.add_edge(1, 0, Eigen::VectorXd::Constant(4, 0.1)); + + auto graph_sim = mio::oseir::make_seir_mobility_sim(t0, dt, std::move(g)); + graph_sim.advance(tmax); + + double total_pop = 0.0; + for (auto& n : graph_sim.get_graph().nodes()) { + total_pop += n.property.get_result().get_last_value().sum(); + } + EXPECT_NEAR(total_pop, 20000.0, 1e-4); + + for (auto& n : graph_sim.get_graph().nodes()) { + auto v = n.property.get_result().get_last_value(); + for (int i = 0; i < v.size(); ++i) { + EXPECT_GE(v[i], -1e-6) << "Negative compartment at index " << i; + } + } +}