From 1c1e730e391065349593ede065fd4da791998784 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 20 Mar 2026 10:40:15 +0100 Subject: [PATCH 1/9] new functionality + tests --- cpp/examples/CMakeLists.txt | 26 +- .../ode_secir_traveltime_mobility.cpp | 160 ++++ cpp/memilio/CMakeLists.txt | 2 + .../metapopulation_mobility_traveltime.cpp | 24 + .../metapopulation_mobility_traveltime.h | 734 ++++++++++++++++++ cpp/tests/CMakeLists.txt | 1 + cpp/tests/test_mobility_traveltime.cpp | 208 +++++ 7 files changed, 1144 insertions(+), 11 deletions(-) create mode 100644 cpp/examples/ode_secir_traveltime_mobility.cpp create mode 100644 cpp/memilio/mobility/metapopulation_mobility_traveltime.cpp create mode 100644 cpp/memilio/mobility/metapopulation_mobility_traveltime.h create mode 100644 cpp/tests/test_mobility_traveltime.cpp diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index e63da90ac1..b76e0e7692 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -18,17 +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_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_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) @@ -108,6 +108,10 @@ add_executable(ode_secir_graph_example ode_secir_graph.cpp) target_link_libraries(ode_secir_graph_example PRIVATE memilio ode_secir) target_compile_options(ode_secir_graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ode_secir_traveltime_example ode_secir_traveltime_mobility.cpp) +target_link_libraries(ode_secir_traveltime_example PRIVATE memilio ode_secir) +target_compile_options(ode_secir_traveltime_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(graph_stochastic_mobility_example graph_stochastic_mobility.cpp) target_link_libraries(graph_stochastic_mobility_example PRIVATE memilio ode_secir) target_compile_options(graph_stochastic_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/examples/ode_secir_traveltime_mobility.cpp b/cpp/examples/ode_secir_traveltime_mobility.cpp new file mode 100644 index 0000000000..095329d1b2 --- /dev/null +++ b/cpp/examples/ode_secir_traveltime_mobility.cpp @@ -0,0 +1,160 @@ +/* +* 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. +*/ + +/** + * Example for the travel-time-aware metapopulation mobility model with three patches. + * + * Only A has initial infections. Individuals from A commute to B via C, spending time in C's mobility node during transit. + * Due to the time spent in the mobility node, C residents get infected early on, and B receives infections both from + * arriving A commuters and from C residents who commute directly to B. + * + * This example demonstrates the travel-time-aware Graph-ODE model from: + * + * H. Zunker et al., "Novel travel time aware metapopulation models ...", (2025), doi.org/10.1371/journal.pcbi.1012630 + */ + +#include "memilio/config.h" +#include "memilio/mobility/metapopulation_mobility_traveltime.h" +#include "ode_secir/model.h" +#include "ode_secir/infection_state.h" +#include "ode_secir/parameters.h" +#include "memilio/compartments/simulation.h" + +#include +#include +#include + +using FP = double; +using Model = mio::osecir::Model; +using Sim = mio::osecir::Simulation; + +static Model build_model(double N, double I0) +{ + Model model(1 /*age groups*/); + + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}] = N - I0; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = I0; + + model.parameters.set>(0); + model.parameters.set>(0.0); + + model.parameters.get>()[mio::AgeGroup(0)] = 3.2; + model.parameters.get>()[mio::AgeGroup(0)] = 2.0; + model.parameters.get>()[mio::AgeGroup(0)] = 5.8; + model.parameters.get>()[mio::AgeGroup(0)] = 9.5; + model.parameters.get>()[mio::AgeGroup(0)] = 7.1; + + model.parameters.get>()[mio::AgeGroup(0)] = 0.08; + model.parameters.get>()[mio::AgeGroup(0)] = 0.7; + model.parameters.get>()[mio::AgeGroup(0)] = 0.09; + model.parameters.get>()[mio::AgeGroup(0)] = 0.25; + model.parameters.get>()[mio::AgeGroup(0)] = 0.45; + model.parameters.get>() = 35; + model.parameters.get>()[mio::AgeGroup(0)] = 0.2; + model.parameters.get>()[mio::AgeGroup(0)] = 0.25; + model.parameters.get>()[mio::AgeGroup(0)] = 0.3; + + auto& cm = model.parameters.get>(); + cm.get_cont_freq_mat()[0].get_baseline() = Eigen::MatrixX::Constant(1, 1, 10.0); + + return model; +} + +static void print_row(const std::string& label, const auto& node_prop) +{ + const auto& last = node_prop.local_sim.get_result().get_last_value(); + const int iS = static_cast(mio::osecir::InfectionState::Susceptible); + const int iE = static_cast(mio::osecir::InfectionState::Exposed); + const int iINS = static_cast(mio::osecir::InfectionState::InfectedNoSymptoms); + std::cout << std::left << std::setw(8) << label << " S=" << std::setw(10) << static_cast(last[iS]) + << " E=" << std::setw(8) << static_cast(last[iE]) << " INS=" << std::setw(8) + << static_cast(last[iINS]) << " Total=" << static_cast(last.sum()) << "\n"; +} + +int main() +{ + mio::set_log_level(mio::LogLevel::warn); + + const FP t0 = 0.0; + const FP tmax = 30.0; + const FP dt = 0.5; + + // -- Population sizes ------------------------------------------------- + const double N_A = 500'000.0; + const double N_B = 300'000.0; + const double N_C = 80'000.0; + const double I0_A = N_A * 0.001; // 0.1 % infected at start + const double I0_B = 0.0; + const double I0_C = 0.0; + + // -- Mobility parameters ----------------------------------------------- + // A -> B trip passes through C + // Choosing 1 h per patch means A traveler spend 1 h in C's mobility node + const FP tt_patch_AC_CB = 1.0 / 24.0; // 1 h per patch + const FP tt_direct_BC = 0.5 / 24.0; // 30 min for direct B <-> C travel + const FP stay_dur_AB = 8.0 / 24.0; + const FP stay_dur_C = 6.0 / 24.0; + + const FP mob_coeff_AB = 0.05; // 5 % of A commute to B + const FP mob_coeff_BC = 0.08; // 8 % of B/C commute to each other + + const int nc = static_cast(mio::osecir::InfectionState::Count); + const Eigen::VectorX coeffs_AB = Eigen::VectorX::Constant(nc, mob_coeff_AB); + const Eigen::VectorX coeffs_BC = Eigen::VectorX::Constant(nc, mob_coeff_BC); + + // Nodes (index 0=A, 1=B, 2=C) + mio::Graph, mio::TravelTimeEdge> graph; + graph.add_node(0, build_model(N_A, I0_A), t0, stay_dur_AB); + graph.add_node(1, build_model(N_B, I0_B), t0, stay_dur_AB); + graph.add_node(2, build_model(N_C, I0_C), t0, stay_dur_C); + + // A -> B via C + const FP tt_AB = 2.0 * tt_patch_AC_CB; + const std::vector path_AB = {0, 2, 1}; // A -> C -> B + graph.add_edge(0, 1, mio::MobilityParameters(coeffs_AB), tt_AB, path_AB); + + // B <-> C + graph.add_edge(1, 2, mio::MobilityParameters(coeffs_BC), tt_direct_BC, std::vector{1, 2}); + graph.add_edge(2, 1, mio::MobilityParameters(coeffs_BC), tt_direct_BC, std::vector{2, 1}); + + // -- Run simulation ---------------------------------------------------- + auto sim = mio::make_traveltime_sim(t0, dt, std::move(graph)); + sim.advance(tmax); + + // -- Results ----------------------------------------------------------- + std::cout << "Patch A starts with " << I0_A << " infected. A traveler from A travels to B via C.\n" + << "C is transit-only for A->B; C also has its own B<->C travelers.\n" + << "B and C start healthy and infections reach them through transit mixing.\n\n"; + + std::cout << std::string(62, '-') << "\n"; + std::cout << "\n Results after " << tmax << " days \n\n"; + print_row("A", sim.get_graph().nodes()[0].property); + print_row("B", sim.get_graph().nodes()[1].property); + print_row("C", sim.get_graph().nodes()[2].property); + std::cout << std::string(62, '-') << "\n"; + + std::cout << "\nNote: C imports infections early because A commuters\n" + << " spend ~1 h inside C's mobility node,\n" + << " infecting susceptible there. This is not captured by a\n" + << " standard metapopulation model without travel time,\n" + << " where transit is instantaneous.\n"; + + return 0; +} diff --git a/cpp/memilio/CMakeLists.txt b/cpp/memilio/CMakeLists.txt index 658733a7f8..c0efd3be11 100644 --- a/cpp/memilio/CMakeLists.txt +++ b/cpp/memilio/CMakeLists.txt @@ -77,6 +77,8 @@ add_library(memilio mobility/metapopulation_mobility_instant.cpp mobility/metapopulation_mobility_stochastic.h mobility/metapopulation_mobility_stochastic.cpp + mobility/metapopulation_mobility_traveltime.h + mobility/metapopulation_mobility_traveltime.cpp mobility/graph_simulation.h mobility/graph_simulation.cpp mobility/graph.h diff --git a/cpp/memilio/mobility/metapopulation_mobility_traveltime.cpp b/cpp/memilio/mobility/metapopulation_mobility_traveltime.cpp new file mode 100644 index 0000000000..79a59b5de0 --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.cpp @@ -0,0 +1,24 @@ +/* +* 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/mobility/metapopulation_mobility_traveltime.h" + +namespace mio +{ +} // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h new file mode 100644 index 0000000000..dc15d61082 --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h @@ -0,0 +1,734 @@ +/* +* 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_TRAVELTIME_H +#define MIO_MOBILITY_TRAVELTIME_H + +#include "memilio/mobility/graph.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/math/eigen.h" +#include "memilio/math/euler.h" +#include "memilio/math/floating_point.h" +#include "memilio/utils/logging.h" +#include "memilio/utils/time_series.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace mio +{ + +/** + * @brief Property for a graph node in the travel-time-aware mobility model. + * + * Each node holds two SimulationNode instances: + * - `local_sim`: Standard local ODE simulation + * - `mobility_sim`: ODE simulation representing individuals currently in transit + * through this node's mobility model + * + * Additionally, `stay_duration` defines how long inbound travelers remain at + * this node before returning home. On a day-scale of [0,1], this is the fraction + * of the day spent at the destination. + * + * @tparam FP Floating-point type (e.g., double). + * @tparam Sim Simulation type (e.g., mio::osecir::Simulation). + */ +template +struct TravelTimeNodeProperty { + + /** + * @brief Constructor if both nodes share the same model and t0. + * + * @param model The epidemiological model. + * @param t0 Start time. + * @param stay Stay duration in [0, 1). + */ + TravelTimeNodeProperty(const typename Sim::Model& model, FP t0, FP stay) + : local_sim(model, t0) + , mobility_sim(model, t0) + , stay_duration(stay) + { + } + + /** + * @brief Construct from two simulation arguments and a stay duration. + * + * @param local_args Arguments forwarded to local_sim. + * @param mobility_args Arguments forwarded to mobility_sim. + * @param stay Stay duration in [0, 1). + */ + template + TravelTimeNodeProperty(std::tuple local_args, std::tuple mobility_args, FP stay) + : local_sim(std::make_from_tuple>(std::move(local_args))) + , mobility_sim(std::make_from_tuple>(std::move(mobility_args))) + , stay_duration(stay) + { + } + + SimulationNode local_sim; ///< Local population dynamics. + SimulationNode mobility_sim; ///< Transit population dynamics. + FP stay_duration; ///< Fraction of day spent at destination, in [0, 1). +}; + +/** + * @brief Edge property for the travel-time-aware mobility model. + * + * Extends the standard MobilityParameters with: + * - `travel_time`: Total travel time for one direction (as fraction of a day). + * - `path`: Ordered list of node indices traversed during the trip + * (includes origin and destination, i.e., [from, r1, r2, ..., to]). + * + * The travel time is divided equally among all units in the path. + * + * @tparam FP Floating-point type. + */ +template +struct TravelTimeEdge { + + /** + * @brief Construct from mobility parameters, travel time, and path. + * @param params Mobility coefficients. + * @param tt One-way travel time (fraction of day). + * @param route Ordered list of node indices (origin, intermediates, destination). + */ + TravelTimeEdge(const MobilityParameters& params, FP tt, std::vector route) + : parameters(params) + , travel_time(tt) + , path(std::move(route)) + , mobile_population(params.get_coefficients().get_shape().rows()) + , return_times(0) + { + } + + /** + * @brief Construct from coefficient vector, travel time, and path. + * @param coeffs Mobility coefficient vector. + * @param tt One-way travel time (fraction of day). + * @param route Ordered list of node indices. + */ + TravelTimeEdge(const Eigen::VectorX& coeffs, FP tt, std::vector route) + : parameters(coeffs) + , travel_time(tt) + , path(std::move(route)) + , mobile_population(coeffs.rows()) + , return_times(0) + { + } + + MobilityParameters parameters; ///< Mobility coefficients and NPIs. + FP travel_time; ///< One-way travel time as fraction of day. + std::vector path; ///< Node indices along the route (origin -> destination). + + // Commuter state tracking (one entry per commuter group, for advance function). + TimeSeries mobile_population; ///< Compartment counts for people currently in transit. + TimeSeries return_times; ///< Times at which the mobile population returns. +}; + +/** + * @brief Correct small negative compartment values caused by the auxiliary Euler + * step to estimate the infection states of travelers. + * + * The vector is partitioned into `num_age_groups` equally sized blocks. + * Any value below `tolerance` is zeroed out and its deficit is added to the + * largest compartment in the same age-group block. This preserves the total + * population within each age group. + * + * @tparam FP Floating-point type. + * @param vec Compartment vector to correct (in-place). + * @param num_age_groups Number of age groups. + * @param tolerance Threshold below which a value is considered negative. + * @param max_iter Maximum number of correction iterations. + */ + +template +void correct_negative_compartments(Eigen::Ref> vec, size_t num_age_groups, FP tolerance = -1e-7, + size_t max_iter = 100) +{ + const size_t n_comparts = static_cast(vec.size()) / num_age_groups; + for (size_t iter = 0; iter < max_iter; ++iter) { + if (vec.minCoeff() >= tolerance) { + return; + } + Eigen::Index min_idx; + const FP min_val = vec.minCoeff(&min_idx); + const auto grp = static_cast(min_idx) / n_comparts; + const auto beg = static_cast(grp * n_comparts); + const auto end = beg + static_cast(n_comparts); + Eigen::Index max_idx; + vec.segment(beg, end - beg).maxCoeff(&max_idx); + max_idx += beg; + vec(max_idx) += min_val; + vec(min_idx) = FP{0}; + } + if (vec.minCoeff() < tolerance) { + log_error("correct_negative_compartments: could not correct all negative values after {} iterations.", + max_iter); + } +} + +/** + * @brief Daily mobility schedule computed once from the graph topology. + * + * The day is discretised into `n_steps` equal time slices of width `1/n_steps`. + * For every edge the schedule records: + * - `node_at_step[edge][step]`: which node index holds the travelers at that step. + * - `in_mobility[edge][step]`: true, if the travelers are in a mobility model in that step. + * - `first_mobility_step[edge]`: first step at which mobility starts. + * + * Derived per-node schedules tell the integrator when it must stop and + * synchronise with the mobility exchange: + * - `local_breakpoints[node]`: steps at which local_sim must be synchronized. + * - `mobility_breakpoints[node]`: steps at which mobility_sim must be synchronized. + */ +struct TravelTimeSchedule { + std::vector> node_at_step; ///< [edge][step] -> node index + std::vector> in_mobility; ///< [edge][step] -> mobility flag + std::vector first_mobility_step; ///< [edge] -> first step with mobility + std::vector> local_breakpoints; ///< [node] -> sorted step indices + std::vector> mobility_breakpoints; ///< [node] -> sorted step indices + /// For each step: which edges exchange travelers + std::vector> edges_at_step; + /// For each step: which nodes advance their local_sim + std::vector> local_nodes_at_step; + /// For each step: which nodes advance their mobility_sim + std::vector> mobility_nodes_at_step; +}; + +/** + * @brief Compute the TravelTimeSchedule for a given graph. + * + * @tparam FP Floating-point type. + * @tparam NodeProp TravelTimeNodeProperty type. + * @tparam EdgeProp TravelTimeEdge type. + * @param graph The simulation graph. + * @param n_steps Number of time steps per day (default: 100). + * @param eps Floating-point tolerance for time comparisons. + * @return Fully populated TravelTimeSchedule. + */ +template +TravelTimeSchedule compute_traveltime_schedule(const Graph& graph, size_t n_steps = 100, + double eps = 1e-10) +{ + const size_t n_edges = graph.edges().size(); + const size_t n_nodes = graph.nodes().size(); + const double dt_step = 1.0 / static_cast(n_steps); + + TravelTimeSchedule sched; + sched.node_at_step.resize(n_edges); + sched.in_mobility.resize(n_edges); + sched.first_mobility_step.resize(n_edges, n_steps); + + // -- 1. Build per-edge schedules -------------------------------------- + for (size_t ei = 0; ei < n_edges; ++ei) { + const auto& e = graph.edges()[ei]; + const auto& edge = e.property; + const size_t dest_node = e.end_node_idx; + const size_t src_node = e.start_node_idx; + const FP stay = graph.nodes()[dest_node].property.stay_duration; + const size_t n_patches = edge.path.size(); // includes start and end nodes + // Travel time per patch (split evenly). + const double tt_per_patch = + std::max(dt_step, std::round(edge.travel_time / static_cast(n_patches) / dt_step) * dt_step); + const double total_travel = tt_per_patch * static_cast(n_patches); + + // Start of outbound travel: everything before is local. + const double t_depart = std::max(0.0, 1.0 - 2.0 * total_travel - static_cast(stay)); + const size_t step_depart = static_cast((t_depart + eps) / dt_step); + // Arrival at destination. + const size_t steps_travel = static_cast((total_travel + eps) / dt_step); + const size_t step_arrive = step_depart + steps_travel; + // Departure from destination for return trip. + const size_t step_leave_dst = n_steps - steps_travel; + + sched.node_at_step[ei].resize(n_steps, src_node); + sched.in_mobility[ei].resize(n_steps, false); + + // Outbound travel: step_depart … step_arrive-1 + sched.first_mobility_step[ei] = step_depart; + size_t current_step = step_depart; + for (size_t patch = 0; patch < n_patches; ++patch) { + const size_t patch_steps = static_cast((tt_per_patch + eps) / dt_step); + for (size_t s = 0; s < patch_steps && (current_step + s) < n_steps; ++s) { + sched.node_at_step[ei][current_step + s] = edge.path[patch]; + sched.in_mobility[ei][current_step + s] = true; + } + current_step += patch_steps; + } + + // Stay at destination: step_arrive … step_leave_dst-1 + for (size_t s = step_arrive; s < step_leave_dst && s < n_steps; ++s) { + sched.node_at_step[ei][s] = dest_node; + sched.in_mobility[ei][s] = false; + } + + // Return travel: step_leave_dst … n_steps-1 (reversed path) + current_step = step_leave_dst; + for (size_t patch = n_patches; patch > 0 && current_step < n_steps; --patch) { + const size_t patch_steps = static_cast((tt_per_patch + eps) / dt_step); + for (size_t s = 0; s < patch_steps && (current_step + s) < n_steps; ++s) { + sched.node_at_step[ei][current_step + s] = edge.path[patch - 1]; + sched.in_mobility[ei][current_step + s] = true; + } + current_step += patch_steps; + } + } + + // -- 2. Derive per-node breakpoint schedules --------------------------- + sched.local_breakpoints.resize(n_nodes); + sched.mobility_breakpoints.resize(n_nodes); + + for (size_t ni = 0; ni < n_nodes; ++ni) { + sched.local_breakpoints[ni].push_back(0); + sched.mobility_breakpoints[ni].push_back(0); + + auto prev_local_edges = std::vector{}; + auto prev_mob_edges = std::vector{}; + for (size_t t = 1; t < n_steps; ++t) { + std::vector local_e, mob_e; + for (size_t ei = 0; ei < n_edges; ++ei) { + if (sched.node_at_step[ei][t] == ni) { + if (sched.in_mobility[ei][t]) + mob_e.push_back(ei); + else + local_e.push_back(ei); + } + } + if (local_e != prev_local_edges) { + sched.local_breakpoints[ni].push_back(t); + prev_local_edges = local_e; + } + if (mob_e != prev_mob_edges) { + sched.mobility_breakpoints[ni].push_back(t); + prev_mob_edges = mob_e; + } + } + // Ensure final step is always present. + if (sched.local_breakpoints[ni].back() != n_steps - 1) + sched.local_breakpoints[ni].push_back(n_steps - 1); + if (sched.mobility_breakpoints[ni].back() != n_steps - 1) + sched.mobility_breakpoints[ni].push_back(n_steps - 1); + } + + // -- 3. Build per-step lookup tables ---------------------------------- + sched.edges_at_step.resize(n_steps); + sched.local_nodes_at_step.resize(n_steps); + sched.mobility_nodes_at_step.resize(n_steps); + + for (size_t t = 0; t < n_steps; ++t) { + // Edges active at this step (first mobility step or node/mode change). + for (size_t ei = 0; ei < n_edges; ++ei) { + if (t < sched.first_mobility_step[ei]) + continue; + bool add = (t == sched.first_mobility_step[ei]); + if (!add && t > 0) { + add = (sched.node_at_step[ei][t] != sched.node_at_step[ei][t - 1]) || + (sched.in_mobility[ei][t] != sched.in_mobility[ei][t - 1]); + } + if (add) + sched.edges_at_step[t].push_back(ei); + } + // Nodes whose local_sim breakpoint falls at t. + for (size_t ni = 0; ni < n_nodes; ++ni) { + if (std::binary_search(sched.local_breakpoints[ni].begin(), sched.local_breakpoints[ni].end(), t)) + sched.local_nodes_at_step[t].push_back(ni); + if (std::binary_search(sched.mobility_breakpoints[ni].begin(), sched.mobility_breakpoints[ni].end(), t)) + sched.mobility_nodes_at_step[t].push_back(ni); + } + } + + return sched; +} + +/** + * @brief Graph simulation with travel-time-aware mobility. + * + * This class implements the travel-time-aware metapopulation model from: + * H. Zunker et al., "Novel travel time aware metapopulation models ...", + * (2024), https://doi.org/10.1371/journal.pcbi.1012630 + * + * Key differences to `metapopulation_mobility_instant.h`: + * - Realistic travel and stay durations + * - Transmission on travel (even with other travelers from different origins/destinations) + * - Explicit trip paths with intermediate nodes (e.g., for transit hubs) + * - More complex day schedule with multiple synchronisation points + * + * Day structure (normalised to [0, 1]): + * 0 --- local --> t_depart --- outbound transit --> t_arrive --- stay --> t_return --- inbound transit --> 1 + * + * The day is discretised into `n_steps`. The schedule is precomputed + * once by compute_traveltime_schedule before the simulation starts. + * + * + * @tparam FP Floating-point type (e.g., double). + * @tparam GraphT Graph type (Graph, TravelTimeEdge>). + */ +template +class GraphSimulationTravelTime +{ +public: + using Graph = GraphT; + using NodeProp = typename GraphT::NodeProperty; + using EdgeProp = typename GraphT::EdgeProperty; + + /** + * @brief Construct the simulation. + * + * @param t0 Simulation start time. + * @param dt Time step. + * @param graph The simulation graph. + * @param n_steps Number of sub-steps per day (default 100). + */ + GraphSimulationTravelTime(FP t0, FP dt, GraphT graph, size_t n_steps = 100) + : m_t(t0) + , m_dt(dt) + , m_graph(std::move(graph)) + , m_n_steps(n_steps) + , m_schedule(compute_traveltime_schedule(m_graph, n_steps)) + { + // Reset mobility_sim compartments to zero as transit nodes always start empty. + for (auto& n : m_graph.nodes()) { + n.property.mobility_sim.get_result().get_last_value().setZero(); + } + } + + FP get_t() const + { + return m_t; + } + + GraphT& get_graph() + { + return m_graph; + } + + const GraphT& get_graph() const + { + return m_graph; + } + + /** + * @brief Advance the simulation up to t_max. + * + * Processes each day in three phases per sub-step: + * 1. Move commuters along active edges (init or transfer to next node). + * 2. Advance local nodes up to the next breakpoint. + * 3. Update commuter infection states in the transit nodes. + * + * At the end of each day, transit nodes are emptied (they should already be empty up to a certain tolerance) + * and for longer simulations (more than interpolation_threshold (default 20) days) results are interpolated to restrict memory usage to one time point per day. + * + * @param t_max Maximum simulation time. + * @param interpolation_threshold Time after which results are interpolated to daily resolution. + */ + void advance(FP t_max, FP interpolation_threshold = FP{20}) + { + const FP sub_dt = FP{1} / static_cast(m_n_steps); + + // Earliest commuter departure across all edges. + FP t_first_departure = compute_first_departure(); + + while (m_t < t_max - FP{1e-10}) { + if (m_t + t_first_departure > t_max) { + // If t_max is reached before the first commuter departure of the next day, + // no mobility exchange occurs. Only advance local sims and stop. + for (auto& n : m_graph.nodes()) { + n.property.local_sim.get_simulation().advance(t_max); + } + } + m_t = t_max; + break; + } + + for (size_t step = 0; step < m_n_steps; ++step) { + const FP t_sub = m_t + static_cast(step) * sub_dt; + + // Phase 1: Move commuters along edges. + move_commuters(step, t_sub); + + // Phase 2: Advance local simulations. + advance_local_sims(step, t_sub, sub_dt); + + // Phase 3: Update infection states of commuters in transit. + update_commuter_states(step, t_sub, sub_dt); + } + + // End-of-day: return remaining transit individuals, reset mobility nodes. + end_of_day(m_t + FP{1}); + + m_t += FP{1}; + + // Interpolate to daily resolution after interpolation_threshold days to limit memory usage. + if (m_t > interpolation_threshold) { + interpolate_results(); + } + } +} + +private : FP m_t; +FP m_dt; +GraphT m_graph; +size_t m_n_steps; +TravelTimeSchedule m_schedule; + +FP compute_first_departure() const +{ + FP earliest = FP{1}; + for (size_t ei = 0; ei < m_graph.edges().size(); ++ei) { + const FP t_dep = static_cast(m_schedule.first_mobility_step[ei]) / m_n_steps; + earliest = std::min(earliest, t_dep); + } + return earliest; +} + +void move_commuters(size_t step, FP t_sub) +{ + for (size_t ei : m_schedule.edges_at_step[step]) { + if (step == m_schedule.first_mobility_step[ei]) { + // First mobility step: extract commuters from local_sim of origin. + init_mobility(ei, t_sub); + } + else { + // Subsequent steps: move commuters from previous to current transit node. + transfer_commuters(ei, step, t_sub); + } + } +} + +/// Extract commuters from local_sim at the origin node and place them in the +/// first transit node's mobility_sim. +void init_mobility(size_t ei, FP t_sub) +{ + auto& e = m_graph.edges()[ei]; + auto& edge = e.property; + auto& node_from = m_graph.nodes()[e.start_node_idx].property.local_sim; + + const auto& last_val = node_from.get_result().get_last_value(); + const auto& coeffs = edge.parameters.get_coefficients().get_matrix_at(SimulationTime(t_sub)); + const auto factors = get_mobility_factors(node_from, t_sub, last_val); + + Eigen::VectorX travel_subgroup = (last_val.array() * coeffs.array() * factors.array()).matrix(); + + // Clamp negatives + travel_subgroup = travel_subgroup.cwiseMax(FP{0}); + + // Optional: test commuters + test_commuters(node_from, travel_subgroup, t_sub); + + edge.mobile_population.add_time_point(t_sub, travel_subgroup); + edge.return_times.add_time_point(t_sub); + + // Move from origin local_sim. + Eigen::Ref> origin_val = node_from.get_result().get_last_value(); + origin_val -= travel_subgroup; + + // Place in first transit node's mobility_sim. + const size_t first_transit = edge.path.front(); + auto& mob_sim = m_graph.nodes()[first_transit].property.mobility_sim; + mob_sim.get_result().get_last_value() += travel_subgroup; +} + +/// Move commuters from the previous node to the current node in the schedule. +void transfer_commuters(size_t ei, size_t step, FP /*t_sub*/) +{ + auto& e = m_graph.edges()[ei]; + auto& edge = e.property; + + if (edge.mobile_population.get_num_time_points() == 0) + return; + + const size_t prev_node = m_schedule.node_at_step[ei][step - 1]; + const size_t curr_node = m_schedule.node_at_step[ei][step]; + const bool prev_mob = m_schedule.in_mobility[ei][step - 1]; + const bool curr_mob = m_schedule.in_mobility[ei][step]; + + if (prev_node == curr_node && prev_mob == curr_mob) + return; + + Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); + size_t num_age = static_cast( + m_graph.nodes()[prev_node].property.local_sim.get_simulation().get_model().parameters.get_num_groups()); + correct_negative_compartments(travel_subgroup, num_age); + + // Remove from previous location. + auto& sim_from = + prev_mob ? m_graph.nodes()[prev_node].property.mobility_sim : m_graph.nodes()[prev_node].property.local_sim; + sim_from.get_result().get_last_value() -= travel_subgroup; + + // Add to new location. + auto& sim_to = + curr_mob ? m_graph.nodes()[curr_node].property.mobility_sim : m_graph.nodes()[curr_node].property.local_sim; + sim_to.get_result().get_last_value() += travel_subgroup; + + correct_negative_compartments(sim_from.get_result().get_last_value(), num_age); + correct_negative_compartments(sim_to.get_result().get_last_value(), num_age); +} + +void advance_local_sims(size_t step, FP t_sub, FP sub_dt) +{ + for (size_t ni : m_schedule.local_nodes_at_step[step]) { + // Determine the next breakpoint to find dt for this node. + const FP node_dt = get_local_dt(ni, step, sub_dt); + m_graph.nodes()[ni].property.local_sim.get_simulation().advance(t_sub + node_dt); + } +} + +FP get_local_dt(size_t ni, size_t step, FP sub_dt) const +{ + const auto& bps = m_schedule.local_breakpoints[ni]; + const auto it = std::lower_bound(bps.begin(), bps.end(), step); + if (it == bps.end() || std::next(it) == bps.end()) { + // Last breakpoint: advance to end of day. + return (static_cast(m_n_steps) - static_cast(step)) * sub_dt; + } + return static_cast(*std::next(it) - step) * sub_dt; +} + +void update_commuter_states(size_t step, FP t_sub, FP sub_dt) +{ + for (size_t ei : m_schedule.edges_at_step[step]) { + auto& e = m_graph.edges()[ei]; + auto& edge = e.property; + + if (edge.mobile_population.get_num_time_points() == 0) + continue; + if (!m_schedule.in_mobility[ei][step]) + continue; // If not in transit, skip. + + const size_t curr_node = m_schedule.node_at_step[ei][step]; + auto& mob_node = m_graph.nodes()[curr_node].property.mobility_sim; + auto& model = mob_node.get_simulation().get_model(); + + // Auxiliary Euler step for commuter sub-population state estimation + FP mob_dt = get_mobility_dt(ei, step, sub_dt); + Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); + const Eigen::VectorX total_mob = mob_node.get_result().get_last_value(); + + Eigen::VectorX y0 = travel_subgroup.eval(); + Eigen::VectorX y1 = Eigen::VectorX::Zero(y0.size()); + FP t_euler = t_sub; + auto deriv_fn = [&](Eigen::Ref> y, FP /*t*/, Eigen::Ref> dydt) { + model.get_derivatives(total_mob, y, t_sub, dydt); + }; + DerivFunction f = deriv_fn; + EulerIntegratorCore().step(f, y0, t_euler, mob_dt, y1); + travel_subgroup = y1; + + // The auxiliary Euler heuristic can is prone to overshooting, especially when a subpopulation share is + // high and the dynamics are fast. To prevent negative compartment values, we apply a correction step. + size_t num_age = static_cast(model.parameters.get_num_groups()); + correct_negative_compartments(travel_subgroup, num_age); + } +} + +FP get_mobility_dt(size_t ei, size_t step, FP sub_dt) const +{ + const size_t curr_node = m_schedule.node_at_step[ei][step]; + const auto& bps = m_schedule.mobility_breakpoints[curr_node]; + const auto it = std::lower_bound(bps.begin(), bps.end(), step); + if (it == bps.end() || std::next(it) == bps.end()) { + return (static_cast(m_n_steps) - static_cast(step)) * sub_dt; + } + return static_cast(*std::next(it) - step) * sub_dt; +} + +void end_of_day(FP t_end) +{ + // Return all still-mobile commuters to their home local_sim. + for (size_t ei = 0; ei < m_graph.edges().size(); ++ei) { + auto& e = m_graph.edges()[ei]; + auto& edge = e.property; + + if (edge.mobile_population.get_num_time_points() == 0) + continue; + + // Final node in schedule is always the origin (home). + const size_t home_node = e.start_node_idx; + const size_t last_step = m_n_steps - 1; + const bool last_mob = m_schedule.in_mobility[ei][last_step]; + auto& sim_last = last_mob ? m_graph.nodes()[m_schedule.node_at_step[ei][last_step]].property.mobility_sim + : m_graph.nodes()[m_schedule.node_at_step[ei][last_step]].property.local_sim; + + Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); + size_t num_age = static_cast( + m_graph.nodes()[home_node].property.local_sim.get_simulation().get_model().parameters.get_num_groups()); + correct_negative_compartments(travel_subgroup, num_age); + + sim_last.get_result().get_last_value() -= travel_subgroup; + m_graph.nodes()[home_node].property.local_sim.get_result().get_last_value() += travel_subgroup; + + // Clear edge commuter tracking. + for (Eigen::Index i = edge.mobile_population.get_num_time_points() - 1; i >= 0; --i) { + edge.mobile_population.remove_time_point(i); + edge.return_times.remove_time_point(i); + } + } + + // Reset mobility_sim to zero: transit nodes start empty each day. + for (auto& n : m_graph.nodes()) { + n.property.mobility_sim.get_result().get_last_value().setZero(); + } + + mio::unused(t_end); +} + +/// Keep only one time point per day in local_sim results (daily resolution). +void interpolate_results() +{ + for (auto& n : m_graph.nodes()) { + auto& res = n.property.local_sim.get_simulation().get_result(); + if (res.get_num_time_points() > 2) { + const FP t_last = res.get_last_time(); + const auto v_last = res.get_last_value().eval(); + while (res.get_num_time_points() > 1) { + res.remove_last_time_point(); + } + res.add_time_point(t_last, v_last); + } + } +} +}; // namespace mio + +/** + * @brief Create a travel-time-aware graph simulation. + * + * @tparam FP Floating-point type. + * @tparam NodeProp TravelTimeNodeProperty specialisation. + * @tparam EdgeProp TravelTimeEdge specialisation. + * @param t0 Simulation start time. + * @param dt Day-level time step. + * @param graph Simulation graph. + * @param n_steps Number of day sub-steps (default 100). + * @return GraphSimulationTravelTime instance. + */ +template +GraphSimulationTravelTime> +make_traveltime_sim(FP t0, FP dt, Graph graph, size_t n_steps = 100) +{ + return GraphSimulationTravelTime>(t0, dt, std::move(graph), n_steps); +} + +} // namespace mio + +#endif // MIO_MOBILITY_TRAVELTIME_H diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 6447c6eabc..eae6e056f5 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_mobility_traveltime.cpp test_date.cpp test_eigen_util.cpp test_odesecir_ageres.cpp diff --git a/cpp/tests/test_mobility_traveltime.cpp b/cpp/tests/test_mobility_traveltime.cpp new file mode 100644 index 0000000000..1ddfbf9ca9 --- /dev/null +++ b/cpp/tests/test_mobility_traveltime.cpp @@ -0,0 +1,208 @@ +/* +* 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/mobility/metapopulation_mobility_instant.h" +#include "memilio/mobility/metapopulation_mobility_traveltime.h" +#include "memilio/compartments/simulation.h" +#include "memilio/math/euler.h" +#include "ode_seir/model.h" +#include "ode_seir/infection_state.h" + +#include "gtest/gtest.h" +#include +#include + +using FP = double; +using Model = mio::oseir::Model; +using Sim = mio::Simulation; + +static Model make_seir_model(double N, double E0) +{ + Model model(1); + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = static_cast(N - E0); + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = static_cast(E0); + model.populations.set_total(static_cast(N)); + model.parameters.get>().get_cont_freq_mat()[0].get_baseline().setConstant(10.0); + model.parameters.set>(0.1); + model.parameters.set>(4.0); + model.parameters.set>(10.0); + return model; +} + +static mio::Graph, mio::TravelTimeEdge> +make_two_node_graph(double N0, double E0, double N1, double E1, FP travel_time, FP stay, FP mob_coeff) +{ + const int n_comp = static_cast(mio::oseir::InfectionState::Count); + Eigen::VectorX coeffs = Eigen::VectorX::Constant(n_comp, mob_coeff); + + mio::Graph, mio::TravelTimeEdge> g; + g.add_node(0, make_seir_model(N0, E0), FP{0}, stay); + g.add_node(1, make_seir_model(N1, E1), FP{0}, stay); + g.add_edge(0, 1, mio::MobilityParameters(coeffs), travel_time, std::vector{0, 1}); + g.add_edge(1, 0, mio::MobilityParameters(coeffs), travel_time, std::vector{1, 0}); + return g; +} + +TEST(TravelTimeMobility, ScheduleBasicStructure) +{ + const FP travel_time = 0.05; + const FP stay_dur = 0.4; + const size_t n_steps = 100; + + auto graph = make_two_node_graph(1000.0, 10.0, 800.0, 0.0, travel_time, stay_dur, 0.1); + auto sched = mio::compute_traveltime_schedule(graph, n_steps); + + // Schedule must have one entry per edge. + ASSERT_EQ(sched.node_at_step.size(), graph.edges().size()); + ASSERT_EQ(sched.in_mobility.size(), graph.edges().size()); + + // Each edge schedule must have n_steps entries. + for (size_t ei = 0; ei < graph.edges().size(); ++ei) { + ASSERT_EQ(sched.node_at_step[ei].size(), n_steps); + ASSERT_EQ(sched.in_mobility[ei].size(), n_steps); + } + + // Breakpoints must be non-empty for every node. + ASSERT_EQ(sched.local_breakpoints.size(), graph.nodes().size()); + for (size_t ni = 0; ni < graph.nodes().size(); ++ni) { + EXPECT_FALSE(sched.local_breakpoints[ni].empty()) << "local_breakpoints empty for node " << ni; + } + + // first_mobility_step must be valid indices. + for (size_t ei = 0; ei < graph.edges().size(); ++ei) { + EXPECT_LT(sched.first_mobility_step[ei], n_steps) << "first_mobility_step out of range for edge " << ei; + } +} + +TEST(TravelTimeMobility, PopulationConservation) +{ + const FP t0 = 0.0; + const FP tmax = 3.0; + const FP travel_time = 0.04; + const FP stay_dur = 0.3; + const double N0 = 5000.0; + const double N1 = 3000.0; + + auto graph = make_two_node_graph(N0, 50.0, N1, 0.0, travel_time, stay_dur, 0.1); + auto sim = mio::make_traveltime_sim(t0, FP{1.0}, std::move(graph)); + sim.advance(tmax); + + const FP total_local0 = sim.get_graph().nodes()[0].property.local_sim.get_result().get_last_value().sum(); + const FP total_local1 = sim.get_graph().nodes()[1].property.local_sim.get_result().get_last_value().sum(); + const FP total_mob0 = sim.get_graph().nodes()[0].property.mobility_sim.get_result().get_last_value().sum(); + const FP total_mob1 = sim.get_graph().nodes()[1].property.mobility_sim.get_result().get_last_value().sum(); + + // At the end of complete days, all commuters must be home. + EXPECT_NEAR(total_mob0, 0.0, 5.0) << "Mobility node 0 not empty at end of simulation"; + EXPECT_NEAR(total_mob1, 0.0, 5.0) << "Mobility node 1 not empty at end of simulation"; + + // Total system population is conserved (SEIR: no deaths). + EXPECT_NEAR(total_local0 + total_local1, N0 + N1, 5.0) << "Total population not conserved"; +} + +TEST(TravelTimeMobility, ZeroMobilityCoefficient) +{ + const FP t0 = 0.0; + const FP tmax = 5.0; + + // Reference plain single-node simulation. + Model model0 = make_seir_model(2000.0, 20.0); + auto single = mio::Simulation(model0, t0); + single.advance(tmax); + const Eigen::VectorX ref = single.get_result().get_last_value().eval(); + + // TravelTime graph with zero mobility coefficient. + auto graph = make_two_node_graph(2000.0, 20.0, 1500.0, 0.0, 0.05, 0.4, 0.0); + auto sim = mio::make_traveltime_sim(t0, FP{1.0}, std::move(graph)); + sim.advance(tmax); + + const auto& v_tt = sim.get_graph().nodes()[0].property.local_sim.get_result().get_last_value(); + + // With zero mobility, node 0 should evolve as the isolated simulation. + EXPECT_NEAR((v_tt - ref).norm(), 0.0, 1.0) << "Zero-mobility node deviates from isolated simulation"; +} + +TEST(TravelTimeMobility, LargeTravelTimeDoesNotCrash) +{ + const FP t0 = 0.0; + + // Nearly all-day travel: 49% + 1% stay + 49% return = 99% of day in travel. + auto graph = make_two_node_graph(1000.0, 10.0, 800.0, 0.0, 0.49, 0.01, 0.05); + + EXPECT_NO_THROW({ + auto sim = mio::make_traveltime_sim(t0, FP{1.0}, std::move(graph)); + sim.advance(FP{2.0}); + }); +} + +TEST(TravelTimeMobility, CorrectNegativeCompartments) +{ + // 2 age groups, 4 compartments each == size 8. + // Use values below the default tolerance (-1e-7) to trigger correction. + Eigen::VectorX v(8); + v << 100.0, 50.0, -1e-5, 30.0, // group 0: negative below tolerance at index 2 + 200.0, 80.0, 40.0, -2e-5; // group 1: negative below tolerance at index 7 + + const FP sum0_before = v.head(4).sum(); + const FP sum1_before = v.tail(4).sum(); + + mio::correct_negative_compartments(v, 2); + + for (int i = 0; i < v.size(); ++i) { + EXPECT_GE(v[i], FP{0}) << "Negative value remains at index " << i; + } + + // Population within each age group is conserved + EXPECT_NEAR(v.head(4).sum(), sum0_before, 1e-12) << "Population not conserved in group 0"; + EXPECT_NEAR(v.tail(4).sum(), sum1_before, 1e-12) << "Population not conserved in group 1"; + + // Values that are within tolerance (-1e-7) should NOT be corrected + Eigen::VectorX v2(4); + v2 << 100.0, 50.0, -1e-8, 30.0; + mio::correct_negative_compartments(v2, 1); + EXPECT_LT(v2[2], FP{0}) << "Value above tolerance should not have been corrected"; +} + +TEST(TravelTimeMobility, ScheduleNodeAssignmentTwoNode) +{ + const FP travel_time = 0.1; + const FP stay_dur = 0.3; + const size_t n_steps = 100; + + auto graph = make_two_node_graph(1000.0, 0.0, 800.0, 0.0, travel_time, stay_dur, 0.1); + auto sched = mio::compute_traveltime_schedule(graph, n_steps); + + // Edge 0->1 is edge index 0 (first edge added). + const size_t ei = 0; + ASSERT_LT(ei, sched.node_at_step.size()); + + const size_t first = sched.first_mobility_step[ei]; + EXPECT_LT(first, n_steps); + + // During the stay phase, in_mobility must be false and node must be dest (1). + bool found_stay = false; + for (size_t s = first; s < n_steps; ++s) { + if (!sched.in_mobility[ei][s]) { + found_stay = true; + EXPECT_EQ(sched.node_at_step[ei][s], size_t{1}) << "Expected destination node during stay at step " << s; + } + } + EXPECT_TRUE(found_stay) << "No stay phase found in edge schedule"; +} From a442043b4df6c98a9295ed208592f33c4fd4a541 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 20 Mar 2026 12:06:13 +0100 Subject: [PATCH 2/9] rm wrong } --- .../metapopulation_mobility_traveltime.h | 402 +++++++++--------- 1 file changed, 201 insertions(+), 201 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h index dc15d61082..5922bb6de5 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h @@ -456,259 +456,259 @@ class GraphSimulationTravelTime for (auto& n : m_graph.nodes()) { n.property.local_sim.get_simulation().advance(t_max); } + m_t = t_max; + break; } - m_t = t_max; - break; - } - for (size_t step = 0; step < m_n_steps; ++step) { - const FP t_sub = m_t + static_cast(step) * sub_dt; + for (size_t step = 0; step < m_n_steps; ++step) { + const FP t_sub = m_t + static_cast(step) * sub_dt; - // Phase 1: Move commuters along edges. - move_commuters(step, t_sub); + // Phase 1: Move commuters along edges. + move_commuters(step, t_sub); - // Phase 2: Advance local simulations. - advance_local_sims(step, t_sub, sub_dt); + // Phase 2: Advance local simulations. + advance_local_sims(step, t_sub, sub_dt); - // Phase 3: Update infection states of commuters in transit. - update_commuter_states(step, t_sub, sub_dt); - } + // Phase 3: Update infection states of commuters in transit. + update_commuter_states(step, t_sub, sub_dt); + } - // End-of-day: return remaining transit individuals, reset mobility nodes. - end_of_day(m_t + FP{1}); + // End-of-day: return remaining transit individuals, reset mobility nodes. + end_of_day(m_t + FP{1}); - m_t += FP{1}; + m_t += FP{1}; - // Interpolate to daily resolution after interpolation_threshold days to limit memory usage. - if (m_t > interpolation_threshold) { - interpolate_results(); + // Interpolate to daily resolution after interpolation_threshold days to limit memory usage. + if (m_t > interpolation_threshold) { + interpolate_results(); + } } } -} -private : FP m_t; -FP m_dt; -GraphT m_graph; -size_t m_n_steps; -TravelTimeSchedule m_schedule; +private: + FP m_t; + FP m_dt; + GraphT m_graph; + size_t m_n_steps; + TravelTimeSchedule m_schedule; -FP compute_first_departure() const -{ - FP earliest = FP{1}; - for (size_t ei = 0; ei < m_graph.edges().size(); ++ei) { - const FP t_dep = static_cast(m_schedule.first_mobility_step[ei]) / m_n_steps; - earliest = std::min(earliest, t_dep); + FP compute_first_departure() const + { + FP earliest = FP{1}; + for (size_t ei = 0; ei < m_graph.edges().size(); ++ei) { + const FP t_dep = static_cast(m_schedule.first_mobility_step[ei]) / m_n_steps; + earliest = std::min(earliest, t_dep); + } + return earliest; } - return earliest; -} -void move_commuters(size_t step, FP t_sub) -{ - for (size_t ei : m_schedule.edges_at_step[step]) { - if (step == m_schedule.first_mobility_step[ei]) { - // First mobility step: extract commuters from local_sim of origin. - init_mobility(ei, t_sub); - } - else { - // Subsequent steps: move commuters from previous to current transit node. - transfer_commuters(ei, step, t_sub); + void move_commuters(size_t step, FP t_sub) + { + for (size_t ei : m_schedule.edges_at_step[step]) { + if (step == m_schedule.first_mobility_step[ei]) { + // First mobility step: extract commuters from local_sim of origin. + init_mobility(ei, t_sub); + } + else { + // Subsequent steps: move commuters from previous to current transit node. + transfer_commuters(ei, step, t_sub); + } } } -} - -/// Extract commuters from local_sim at the origin node and place them in the -/// first transit node's mobility_sim. -void init_mobility(size_t ei, FP t_sub) -{ - auto& e = m_graph.edges()[ei]; - auto& edge = e.property; - auto& node_from = m_graph.nodes()[e.start_node_idx].property.local_sim; - const auto& last_val = node_from.get_result().get_last_value(); - const auto& coeffs = edge.parameters.get_coefficients().get_matrix_at(SimulationTime(t_sub)); - const auto factors = get_mobility_factors(node_from, t_sub, last_val); + /// Extract commuters from local_sim at the origin node and place them in the + /// first transit node's mobility_sim. + void init_mobility(size_t ei, FP t_sub) + { + auto& e = m_graph.edges()[ei]; + auto& edge = e.property; + auto& node_from = m_graph.nodes()[e.start_node_idx].property.local_sim; - Eigen::VectorX travel_subgroup = (last_val.array() * coeffs.array() * factors.array()).matrix(); + const auto& last_val = node_from.get_result().get_last_value(); + const auto& coeffs = edge.parameters.get_coefficients().get_matrix_at(SimulationTime(t_sub)); + const auto factors = get_mobility_factors(node_from, t_sub, last_val); - // Clamp negatives - travel_subgroup = travel_subgroup.cwiseMax(FP{0}); + Eigen::VectorX travel_subgroup = (last_val.array() * coeffs.array() * factors.array()).matrix(); - // Optional: test commuters - test_commuters(node_from, travel_subgroup, t_sub); + // Clamp negatives + travel_subgroup = travel_subgroup.cwiseMax(FP{0}); - edge.mobile_population.add_time_point(t_sub, travel_subgroup); - edge.return_times.add_time_point(t_sub); + // Optional: test commuters + test_commuters(node_from, travel_subgroup, t_sub); - // Move from origin local_sim. - Eigen::Ref> origin_val = node_from.get_result().get_last_value(); - origin_val -= travel_subgroup; + edge.mobile_population.add_time_point(t_sub, travel_subgroup); + edge.return_times.add_time_point(t_sub); - // Place in first transit node's mobility_sim. - const size_t first_transit = edge.path.front(); - auto& mob_sim = m_graph.nodes()[first_transit].property.mobility_sim; - mob_sim.get_result().get_last_value() += travel_subgroup; -} + // Move from origin local_sim. + Eigen::Ref> origin_val = node_from.get_result().get_last_value(); + origin_val -= travel_subgroup; -/// Move commuters from the previous node to the current node in the schedule. -void transfer_commuters(size_t ei, size_t step, FP /*t_sub*/) -{ - auto& e = m_graph.edges()[ei]; - auto& edge = e.property; + // Place in first transit node's mobility_sim. + const size_t first_transit = edge.path.front(); + auto& mob_sim = m_graph.nodes()[first_transit].property.mobility_sim; + mob_sim.get_result().get_last_value() += travel_subgroup; + } - if (edge.mobile_population.get_num_time_points() == 0) - return; + /// Move commuters from the previous node to the current node in the schedule. + void transfer_commuters(size_t ei, size_t step, FP /*t_sub*/) + { + auto& e = m_graph.edges()[ei]; + auto& edge = e.property; - const size_t prev_node = m_schedule.node_at_step[ei][step - 1]; - const size_t curr_node = m_schedule.node_at_step[ei][step]; - const bool prev_mob = m_schedule.in_mobility[ei][step - 1]; - const bool curr_mob = m_schedule.in_mobility[ei][step]; + if (edge.mobile_population.get_num_time_points() == 0) + return; - if (prev_node == curr_node && prev_mob == curr_mob) - return; + const size_t prev_node = m_schedule.node_at_step[ei][step - 1]; + const size_t curr_node = m_schedule.node_at_step[ei][step]; + const bool prev_mob = m_schedule.in_mobility[ei][step - 1]; + const bool curr_mob = m_schedule.in_mobility[ei][step]; - Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); - size_t num_age = static_cast( - m_graph.nodes()[prev_node].property.local_sim.get_simulation().get_model().parameters.get_num_groups()); - correct_negative_compartments(travel_subgroup, num_age); + if (prev_node == curr_node && prev_mob == curr_mob) + return; - // Remove from previous location. - auto& sim_from = - prev_mob ? m_graph.nodes()[prev_node].property.mobility_sim : m_graph.nodes()[prev_node].property.local_sim; - sim_from.get_result().get_last_value() -= travel_subgroup; + Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); + size_t num_age = static_cast( + m_graph.nodes()[prev_node].property.local_sim.get_simulation().get_model().parameters.get_num_groups()); + correct_negative_compartments(travel_subgroup, num_age); - // Add to new location. - auto& sim_to = - curr_mob ? m_graph.nodes()[curr_node].property.mobility_sim : m_graph.nodes()[curr_node].property.local_sim; - sim_to.get_result().get_last_value() += travel_subgroup; + // Remove from previous location. + auto& sim_from = + prev_mob ? m_graph.nodes()[prev_node].property.mobility_sim : m_graph.nodes()[prev_node].property.local_sim; + sim_from.get_result().get_last_value() -= travel_subgroup; - correct_negative_compartments(sim_from.get_result().get_last_value(), num_age); - correct_negative_compartments(sim_to.get_result().get_last_value(), num_age); -} + // Add to new location. + auto& sim_to = + curr_mob ? m_graph.nodes()[curr_node].property.mobility_sim : m_graph.nodes()[curr_node].property.local_sim; + sim_to.get_result().get_last_value() += travel_subgroup; -void advance_local_sims(size_t step, FP t_sub, FP sub_dt) -{ - for (size_t ni : m_schedule.local_nodes_at_step[step]) { - // Determine the next breakpoint to find dt for this node. - const FP node_dt = get_local_dt(ni, step, sub_dt); - m_graph.nodes()[ni].property.local_sim.get_simulation().advance(t_sub + node_dt); + correct_negative_compartments(sim_from.get_result().get_last_value(), num_age); + correct_negative_compartments(sim_to.get_result().get_last_value(), num_age); } -} -FP get_local_dt(size_t ni, size_t step, FP sub_dt) const -{ - const auto& bps = m_schedule.local_breakpoints[ni]; - const auto it = std::lower_bound(bps.begin(), bps.end(), step); - if (it == bps.end() || std::next(it) == bps.end()) { - // Last breakpoint: advance to end of day. - return (static_cast(m_n_steps) - static_cast(step)) * sub_dt; + void advance_local_sims(size_t step, FP t_sub, FP sub_dt) + { + for (size_t ni : m_schedule.local_nodes_at_step[step]) { + // Determine the next breakpoint to find dt for this node. + const FP node_dt = get_local_dt(ni, step, sub_dt); + m_graph.nodes()[ni].property.local_sim.get_simulation().advance(t_sub + node_dt); + } } - return static_cast(*std::next(it) - step) * sub_dt; -} -void update_commuter_states(size_t step, FP t_sub, FP sub_dt) -{ - for (size_t ei : m_schedule.edges_at_step[step]) { - auto& e = m_graph.edges()[ei]; - auto& edge = e.property; - - if (edge.mobile_population.get_num_time_points() == 0) - continue; - if (!m_schedule.in_mobility[ei][step]) - continue; // If not in transit, skip. + FP get_local_dt(size_t ni, size_t step, FP sub_dt) const + { + const auto& bps = m_schedule.local_breakpoints[ni]; + const auto it = std::lower_bound(bps.begin(), bps.end(), step); + if (it == bps.end() || std::next(it) == bps.end()) { + // Last breakpoint: advance to end of day. + return (static_cast(m_n_steps) - static_cast(step)) * sub_dt; + } + return static_cast(*std::next(it) - step) * sub_dt; + } - const size_t curr_node = m_schedule.node_at_step[ei][step]; - auto& mob_node = m_graph.nodes()[curr_node].property.mobility_sim; - auto& model = mob_node.get_simulation().get_model(); + void update_commuter_states(size_t step, FP t_sub, FP sub_dt) + { + for (size_t ei : m_schedule.edges_at_step[step]) { + auto& e = m_graph.edges()[ei]; + auto& edge = e.property; - // Auxiliary Euler step for commuter sub-population state estimation - FP mob_dt = get_mobility_dt(ei, step, sub_dt); - Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); - const Eigen::VectorX total_mob = mob_node.get_result().get_last_value(); - - Eigen::VectorX y0 = travel_subgroup.eval(); - Eigen::VectorX y1 = Eigen::VectorX::Zero(y0.size()); - FP t_euler = t_sub; - auto deriv_fn = [&](Eigen::Ref> y, FP /*t*/, Eigen::Ref> dydt) { - model.get_derivatives(total_mob, y, t_sub, dydt); - }; - DerivFunction f = deriv_fn; - EulerIntegratorCore().step(f, y0, t_euler, mob_dt, y1); - travel_subgroup = y1; - - // The auxiliary Euler heuristic can is prone to overshooting, especially when a subpopulation share is - // high and the dynamics are fast. To prevent negative compartment values, we apply a correction step. - size_t num_age = static_cast(model.parameters.get_num_groups()); - correct_negative_compartments(travel_subgroup, num_age); + if (edge.mobile_population.get_num_time_points() == 0) + continue; + if (!m_schedule.in_mobility[ei][step]) + continue; // If not in transit, skip. + + const size_t curr_node = m_schedule.node_at_step[ei][step]; + auto& mob_node = m_graph.nodes()[curr_node].property.mobility_sim; + auto& model = mob_node.get_simulation().get_model(); + + // Auxiliary Euler step for commuter sub-population state estimation + FP mob_dt = get_mobility_dt(ei, step, sub_dt); + Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); + const Eigen::VectorX total_mob = mob_node.get_result().get_last_value(); + + Eigen::VectorX y0 = travel_subgroup.eval(); + Eigen::VectorX y1 = Eigen::VectorX::Zero(y0.size()); + FP t_euler = t_sub; + auto deriv_fn = [&](Eigen::Ref> y, FP /*t*/, Eigen::Ref> dydt) { + model.get_derivatives(total_mob, y, t_sub, dydt); + }; + DerivFunction f = deriv_fn; + EulerIntegratorCore().step(f, y0, t_euler, mob_dt, y1); + travel_subgroup = y1; + + // The auxiliary Euler heuristic can is prone to overshooting, especially when a subpopulation share is + // high and the dynamics are fast. To prevent negative compartment values, we apply a correction step. + size_t num_age = static_cast(model.parameters.get_num_groups()); + correct_negative_compartments(travel_subgroup, num_age); + } } -} -FP get_mobility_dt(size_t ei, size_t step, FP sub_dt) const -{ - const size_t curr_node = m_schedule.node_at_step[ei][step]; - const auto& bps = m_schedule.mobility_breakpoints[curr_node]; - const auto it = std::lower_bound(bps.begin(), bps.end(), step); - if (it == bps.end() || std::next(it) == bps.end()) { - return (static_cast(m_n_steps) - static_cast(step)) * sub_dt; + FP get_mobility_dt(size_t ei, size_t step, FP sub_dt) const + { + const size_t curr_node = m_schedule.node_at_step[ei][step]; + const auto& bps = m_schedule.mobility_breakpoints[curr_node]; + const auto it = std::lower_bound(bps.begin(), bps.end(), step); + if (it == bps.end() || std::next(it) == bps.end()) { + return (static_cast(m_n_steps) - static_cast(step)) * sub_dt; + } + return static_cast(*std::next(it) - step) * sub_dt; } - return static_cast(*std::next(it) - step) * sub_dt; -} - -void end_of_day(FP t_end) -{ - // Return all still-mobile commuters to their home local_sim. - for (size_t ei = 0; ei < m_graph.edges().size(); ++ei) { - auto& e = m_graph.edges()[ei]; - auto& edge = e.property; - - if (edge.mobile_population.get_num_time_points() == 0) - continue; - // Final node in schedule is always the origin (home). - const size_t home_node = e.start_node_idx; - const size_t last_step = m_n_steps - 1; - const bool last_mob = m_schedule.in_mobility[ei][last_step]; - auto& sim_last = last_mob ? m_graph.nodes()[m_schedule.node_at_step[ei][last_step]].property.mobility_sim - : m_graph.nodes()[m_schedule.node_at_step[ei][last_step]].property.local_sim; + void end_of_day(FP t_end) + { + // Return all still-mobile commuters to their home local_sim. + for (size_t ei = 0; ei < m_graph.edges().size(); ++ei) { + auto& e = m_graph.edges()[ei]; + auto& edge = e.property; - Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); - size_t num_age = static_cast( - m_graph.nodes()[home_node].property.local_sim.get_simulation().get_model().parameters.get_num_groups()); - correct_negative_compartments(travel_subgroup, num_age); + if (edge.mobile_population.get_num_time_points() == 0) + continue; - sim_last.get_result().get_last_value() -= travel_subgroup; - m_graph.nodes()[home_node].property.local_sim.get_result().get_last_value() += travel_subgroup; + // Final node in schedule is always the origin (home). + const size_t home_node = e.start_node_idx; + const size_t last_step = m_n_steps - 1; + const bool last_mob = m_schedule.in_mobility[ei][last_step]; + auto& sim_last = last_mob ? m_graph.nodes()[m_schedule.node_at_step[ei][last_step]].property.mobility_sim + : m_graph.nodes()[m_schedule.node_at_step[ei][last_step]].property.local_sim; + + Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); + size_t num_age = static_cast( + m_graph.nodes()[home_node].property.local_sim.get_simulation().get_model().parameters.get_num_groups()); + correct_negative_compartments(travel_subgroup, num_age); + + sim_last.get_result().get_last_value() -= travel_subgroup; + m_graph.nodes()[home_node].property.local_sim.get_result().get_last_value() += travel_subgroup; + + // Clear edge commuter tracking. + for (Eigen::Index i = edge.mobile_population.get_num_time_points() - 1; i >= 0; --i) { + edge.mobile_population.remove_time_point(i); + edge.return_times.remove_time_point(i); + } + } - // Clear edge commuter tracking. - for (Eigen::Index i = edge.mobile_population.get_num_time_points() - 1; i >= 0; --i) { - edge.mobile_population.remove_time_point(i); - edge.return_times.remove_time_point(i); + // Reset mobility_sim to zero: transit nodes start empty each day. + for (auto& n : m_graph.nodes()) { + n.property.mobility_sim.get_result().get_last_value().setZero(); } - } - // Reset mobility_sim to zero: transit nodes start empty each day. - for (auto& n : m_graph.nodes()) { - n.property.mobility_sim.get_result().get_last_value().setZero(); + mio::unused(t_end); } - mio::unused(t_end); -} - -/// Keep only one time point per day in local_sim results (daily resolution). -void interpolate_results() -{ - for (auto& n : m_graph.nodes()) { - auto& res = n.property.local_sim.get_simulation().get_result(); - if (res.get_num_time_points() > 2) { - const FP t_last = res.get_last_time(); - const auto v_last = res.get_last_value().eval(); - while (res.get_num_time_points() > 1) { - res.remove_last_time_point(); + /// Keep only one time point per day in local_sim results (daily resolution). + void interpolate_results() + { + for (auto& n : m_graph.nodes()) { + auto& res = n.property.local_sim.get_simulation().get_result(); + if (res.get_num_time_points() > 2) { + const FP t_last = res.get_last_time(); + const auto v_last = res.get_last_value().eval(); + while (res.get_num_time_points() > 1) { + res.remove_last_time_point(); + } + res.add_time_point(t_last, v_last); } - res.add_time_point(t_last, v_last); } } -} -}; // namespace mio +}; // class GraphSimulationTravelTime /** * @brief Create a travel-time-aware graph simulation. From 1cd6c740768a4e4e8cf54c26a11acebbedea6f5b Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 20 Mar 2026 12:16:38 +0100 Subject: [PATCH 3/9] better example, all share transit node C --- .../ode_secir_traveltime_mobility.cpp | 32 +++++++++++++------ 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/cpp/examples/ode_secir_traveltime_mobility.cpp b/cpp/examples/ode_secir_traveltime_mobility.cpp index 095329d1b2..9d43b4d473 100644 --- a/cpp/examples/ode_secir_traveltime_mobility.cpp +++ b/cpp/examples/ode_secir_traveltime_mobility.cpp @@ -105,10 +105,20 @@ int main() const double I0_C = 0.0; // -- Mobility parameters ----------------------------------------------- - // A -> B trip passes through C - // Choosing 1 h per patch means A traveler spend 1 h in C's mobility node - const FP tt_patch_AC_CB = 1.0 / 24.0; // 1 h per patch - const FP tt_direct_BC = 0.5 / 24.0; // 30 min for direct B <-> C travel + // A -> B trip passes through C (path = {A, C, B}, n_patches = 3). + // travel_time_AB = 2 * tt_patch_AC_CB; with n_steps=100 this yields tt_per_patch = 0.03 (3 sub-steps/patch). + // step_depart_AB = 48 -> A's commuters are in C's mobility_sim at sub-steps 51–53 (outbound) + // and at sub-steps 94–96 (return). + // + // B <-> C (path = {B,C} or {C,B}, n_patches = 2). + // tt_direct_BC = 2/24 -> tt_per_patch = 0.04 (4 sub-steps/patch). + // C->B: step_depart = 50 -> C->B commuters are in C's mobility_sim at sub-steps 50–53 + // -> overlaps with A's outbound transit (51–53) in C's mobility_sim. + // B->C return: in C's mobility_sim at sub-steps 92–95 + // -> overlaps with A's return transit (94–96) in C's mobility_sim. + // All three groups therefore share C's mobility_sim and can infect each other during transit. + const FP tt_patch_AC_CB = 1.0 / 24.0; // 1 h per patch (A <-> C) + const FP tt_direct_BC = 2.0 / 24.0; // 2 h one-way B <-> C (tt_per_patch=0.04 -> schedule overlap in C) const FP stay_dur_AB = 8.0 / 24.0; const FP stay_dur_C = 6.0 / 24.0; @@ -150,11 +160,15 @@ int main() print_row("C", sim.get_graph().nodes()[2].property); std::cout << std::string(62, '-') << "\n"; - std::cout << "\nNote: C imports infections early because A commuters\n" - << " spend ~1 h inside C's mobility node,\n" - << " infecting susceptible there. This is not captured by a\n" - << " standard metapopulation model without travel time,\n" - << " where transit is instantaneous.\n"; + std::cout << "\nSchedule overlap (n_steps=100, sub-step = 0.01 day):\n" + << " A outbound through C : sub-steps 51-53\n" + << " C->B outbound in C : sub-steps 50-53 -> overlap 51-53\n" + << " A return through C : sub-steps 94-96\n" + << " B->C return in C : sub-steps 92-95 -> overlap 94-95\n" + << "\nAll three groups share C's mobility_sim during these steps,\n" + << "so infections can spread between A-transit commuters and C<->B commuters.\n" + << "This effect is not captured by standard metapopulation models\n" + << "without explicit travel time.\n"; return 0; } From 2dfea6cd8352ebc3b4459550b55dc4ef85ae1abc Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 20 Mar 2026 12:32:39 +0100 Subject: [PATCH 4/9] gcc 13 bug with dangling reference --- .../metapopulation_mobility_traveltime.h | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h index 5922bb6de5..fbbaef5471 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h @@ -241,8 +241,9 @@ TravelTimeSchedule compute_traveltime_schedule(const Graph& sched.first_mobility_step.resize(n_edges, n_steps); // -- 1. Build per-edge schedules -------------------------------------- + auto all_edges = graph.edges(); for (size_t ei = 0; ei < n_edges; ++ei) { - const auto& e = graph.edges()[ei]; + const auto& e = all_edges[ei]; const auto& edge = e.property; const size_t dest_node = e.end_node_idx; const size_t src_node = e.start_node_idx; @@ -520,7 +521,8 @@ class GraphSimulationTravelTime /// first transit node's mobility_sim. void init_mobility(size_t ei, FP t_sub) { - auto& e = m_graph.edges()[ei]; + auto all_edges = m_graph.edges(); + auto& e = all_edges[ei]; auto& edge = e.property; auto& node_from = m_graph.nodes()[e.start_node_idx].property.local_sim; @@ -552,8 +554,9 @@ class GraphSimulationTravelTime /// Move commuters from the previous node to the current node in the schedule. void transfer_commuters(size_t ei, size_t step, FP /*t_sub*/) { - auto& e = m_graph.edges()[ei]; - auto& edge = e.property; + auto all_edges = m_graph.edges(); + auto& e = all_edges[ei]; + auto& edge = e.property; if (edge.mobile_population.get_num_time_points() == 0) return; @@ -607,8 +610,9 @@ class GraphSimulationTravelTime void update_commuter_states(size_t step, FP t_sub, FP sub_dt) { + auto all_edges = m_graph.edges(); for (size_t ei : m_schedule.edges_at_step[step]) { - auto& e = m_graph.edges()[ei]; + auto& e = all_edges[ei]; auto& edge = e.property; if (edge.mobile_population.get_num_time_points() == 0) @@ -656,8 +660,9 @@ class GraphSimulationTravelTime void end_of_day(FP t_end) { // Return all still-mobile commuters to their home local_sim. - for (size_t ei = 0; ei < m_graph.edges().size(); ++ei) { - auto& e = m_graph.edges()[ei]; + auto all_edges = m_graph.edges(); + for (size_t ei = 0; ei < all_edges.size(); ++ei) { + auto& e = all_edges[ei]; auto& edge = e.property; if (edge.mobile_population.get_num_time_points() == 0) From 6bc733bcf44dfe845a587f90a84d4074558a9b9c Mon Sep 17 00:00:00 2001 From: Henrik Zunker <69154294+HenrZu@users.noreply.github.com> Date: Fri, 27 Mar 2026 10:22:28 +0100 Subject: [PATCH 5/9] [ci skip] Apply suggestions from code review Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com> --- cpp/memilio/mobility/metapopulation_mobility_traveltime.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h index fbbaef5471..40d810eaad 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h @@ -17,8 +17,8 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#ifndef MIO_MOBILITY_TRAVELTIME_H -#define MIO_MOBILITY_TRAVELTIME_H +#ifndef MIO_MOBILITY_METAPOPULATION_MOBILITY_TRAVELTIME_H +#define MIO_MOBILITY_METAPOPULATION_MOBILITY_TRAVELTIME_H #include "memilio/mobility/graph.h" #include "memilio/mobility/graph_simulation.h" @@ -736,4 +736,4 @@ make_traveltime_sim(FP t0, FP dt, Graph graph, size_t n_step } // namespace mio -#endif // MIO_MOBILITY_TRAVELTIME_H +#endif // MIO_MOBILITY_METAPOPULATION_MOBILITY_TRAVELTIME_H From 95c2ca47137695fd847668a017ee9020b21fae21 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 27 Mar 2026 11:26:06 +0100 Subject: [PATCH 6/9] [ci skip] traveltime_h review suggestions --- .../metapopulation_mobility_traveltime.h | 99 +++++++++---------- 1 file changed, 44 insertions(+), 55 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h index 40d810eaad..7e0257715c 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h @@ -26,6 +26,7 @@ #include "memilio/math/eigen.h" #include "memilio/math/euler.h" #include "memilio/math/floating_point.h" +#include "memilio/math/math_utils.h" #include "memilio/utils/logging.h" #include "memilio/utils/time_series.h" @@ -150,41 +151,25 @@ struct TravelTimeEdge { * @brief Correct small negative compartment values caused by the auxiliary Euler * step to estimate the infection states of travelers. * - * The vector is partitioned into `num_age_groups` equally sized blocks. - * Any value below `tolerance` is zeroed out and its deficit is added to the - * largest compartment in the same age-group block. This preserves the total - * population within each age group. + * The vector is partitioned into `num_age_groups` equally sized blocks. For each + * block, mio::map_to_nonnegative is applied, which redistributes negative values + * proportionally across the positive entries while preserving the nonnegative sum. * * @tparam FP Floating-point type. * @param vec Compartment vector to correct (in-place). - * @param num_age_groups Number of age groups. - * @param tolerance Threshold below which a value is considered negative. - * @param max_iter Maximum number of correction iterations. + * @param num_age_groups Number of age groups */ - template -void correct_negative_compartments(Eigen::Ref> vec, size_t num_age_groups, FP tolerance = -1e-7, - size_t max_iter = 100) +void correct_negative_compartments(Eigen::Ref> vec, size_t num_age_groups) { - const size_t n_comparts = static_cast(vec.size()) / num_age_groups; - for (size_t iter = 0; iter < max_iter; ++iter) { - if (vec.minCoeff() >= tolerance) { - return; + const Eigen::Index n_comparts = static_cast(vec.size()) / static_cast(num_age_groups); + for (size_t grp = 0; grp < num_age_groups; ++grp) { + auto slice = vec.segment(static_cast(grp) * n_comparts, n_comparts); + auto result = map_to_nonnegative(slice); + if (!result) { + log_error("correct_negative_compartments: could not map age group {} to nonnegative values: {}", grp, + result.error().message()); } - Eigen::Index min_idx; - const FP min_val = vec.minCoeff(&min_idx); - const auto grp = static_cast(min_idx) / n_comparts; - const auto beg = static_cast(grp * n_comparts); - const auto end = beg + static_cast(n_comparts); - Eigen::Index max_idx; - vec.segment(beg, end - beg).maxCoeff(&max_idx); - max_idx += beg; - vec(max_idx) += min_val; - vec(min_idx) = FP{0}; - } - if (vec.minCoeff() < tolerance) { - log_error("correct_negative_compartments: could not correct all negative values after {} iterations.", - max_iter); } } @@ -223,13 +208,15 @@ struct TravelTimeSchedule { * @tparam NodeProp TravelTimeNodeProperty type. * @tparam EdgeProp TravelTimeEdge type. * @param graph The simulation graph. - * @param n_steps Number of time steps per day (default: 100). + * @param n_steps Number of time steps per day (default: 120). The value 120 corresponds to a resolution of + * 12 minutes per step and is divisible by many common travel-time fractions (e.g. 0.05, 0.1, + * 0.25), so typical commute durations map directly to exact integer step counts without rounding. * @param eps Floating-point tolerance for time comparisons. * @return Fully populated TravelTimeSchedule. */ -template -TravelTimeSchedule compute_traveltime_schedule(const Graph& graph, size_t n_steps = 100, - double eps = 1e-10) +template +TravelTimeSchedule compute_traveltime_schedule(const Graph, TravelTimeEdge>& graph, + size_t n_steps = 120, FP eps = Limits::zero_tolerance()) { const size_t n_edges = graph.edges().size(); const size_t n_nodes = graph.nodes().size(); @@ -240,7 +227,7 @@ TravelTimeSchedule compute_traveltime_schedule(const Graph& sched.in_mobility.resize(n_edges); sched.first_mobility_step.resize(n_edges, n_steps); - // -- 1. Build per-edge schedules -------------------------------------- + // 1. Build per-edge schedules auto all_edges = graph.edges(); for (size_t ei = 0; ei < n_edges; ++ei) { const auto& e = all_edges[ei]; @@ -296,7 +283,7 @@ TravelTimeSchedule compute_traveltime_schedule(const Graph& } } - // -- 2. Derive per-node breakpoint schedules --------------------------- + // 2. Derive per-node breakpoint schedules sched.local_breakpoints.resize(n_nodes); sched.mobility_breakpoints.resize(n_nodes); @@ -332,7 +319,7 @@ TravelTimeSchedule compute_traveltime_schedule(const Graph& sched.mobility_breakpoints[ni].push_back(n_steps - 1); } - // -- 3. Build per-step lookup tables ---------------------------------- + // 3. Build per-step lookup tables sched.edges_at_step.resize(n_steps); sched.local_nodes_at_step.resize(n_steps); sched.mobility_nodes_at_step.resize(n_steps); @@ -401,12 +388,12 @@ class GraphSimulationTravelTime * @param graph The simulation graph. * @param n_steps Number of sub-steps per day (default 100). */ - GraphSimulationTravelTime(FP t0, FP dt, GraphT graph, size_t n_steps = 100) + GraphSimulationTravelTime(FP t0, FP dt, GraphT graph, size_t n_steps = 120) : m_t(t0) , m_dt(dt) , m_graph(std::move(graph)) , m_n_steps(n_steps) - , m_schedule(compute_traveltime_schedule(m_graph, n_steps)) + , m_schedule(compute_traveltime_schedule(m_graph, n_steps)) { // Reset mobility_sim compartments to zero as transit nodes always start empty. for (auto& n : m_graph.nodes()) { @@ -438,20 +425,22 @@ class GraphSimulationTravelTime * 3. Update commuter infection states in the transit nodes. * * At the end of each day, transit nodes are emptied (they should already be empty up to a certain tolerance) - * and for longer simulations (more than interpolation_threshold (default 20) days) results are interpolated to restrict memory usage to one time point per day. + * and for longer simulations (more than interpolation_threshold (default 20) days) results are reduced to + * one time point per day to limit memory usage. * * @param t_max Maximum simulation time. - * @param interpolation_threshold Time after which results are interpolated to daily resolution. + * @param interpolation_threshold Time after which results are reduced to daily resolution. */ void advance(FP t_max, FP interpolation_threshold = FP{20}) { const FP sub_dt = FP{1} / static_cast(m_n_steps); // Earliest commuter departure across all edges. - FP t_first_departure = compute_first_departure(); + const FP t_first_departure = compute_first_departure(); - while (m_t < t_max - FP{1e-10}) { - if (m_t + t_first_departure > t_max) { + const auto n_days = static_cast(t_max - m_t); + for (size_t day = 0; day < n_days; ++day) { + if (m_t + t_first_departure >= t_max) { // If t_max is reached before the first commuter departure of the next day, // no mobility exchange occurs. Only advance local sims and stop. for (auto& n : m_graph.nodes()) { @@ -475,13 +464,13 @@ class GraphSimulationTravelTime } // End-of-day: return remaining transit individuals, reset mobility nodes. - end_of_day(m_t + FP{1}); + end_of_day(); m_t += FP{1}; - // Interpolate to daily resolution after interpolation_threshold days to limit memory usage. + // Reduce to daily resolution after interpolation_threshold days to limit memory usage. if (m_t > interpolation_threshold) { - interpolate_results(); + reduce_to_daily_resolution(); } } } @@ -512,7 +501,7 @@ class GraphSimulationTravelTime } else { // Subsequent steps: move commuters from previous to current transit node. - transfer_commuters(ei, step, t_sub); + transfer_commuters(ei, step); } } } @@ -552,7 +541,7 @@ class GraphSimulationTravelTime } /// Move commuters from the previous node to the current node in the schedule. - void transfer_commuters(size_t ei, size_t step, FP /*t_sub*/) + void transfer_commuters(size_t ei, size_t step) { auto all_edges = m_graph.edges(); auto& e = all_edges[ei]; @@ -657,7 +646,7 @@ class GraphSimulationTravelTime return static_cast(*std::next(it) - step) * sub_dt; } - void end_of_day(FP t_end) + void end_of_day() { // Return all still-mobile commuters to their home local_sim. auto all_edges = m_graph.edges(); @@ -694,12 +683,10 @@ class GraphSimulationTravelTime for (auto& n : m_graph.nodes()) { n.property.mobility_sim.get_result().get_last_value().setZero(); } - - mio::unused(t_end); } /// Keep only one time point per day in local_sim results (daily resolution). - void interpolate_results() + void reduce_to_daily_resolution() { for (auto& n : m_graph.nodes()) { auto& res = n.property.local_sim.get_simulation().get_result(); @@ -727,11 +714,13 @@ class GraphSimulationTravelTime * @param n_steps Number of day sub-steps (default 100). * @return GraphSimulationTravelTime instance. */ -template -GraphSimulationTravelTime> -make_traveltime_sim(FP t0, FP dt, Graph graph, size_t n_steps = 100) +template +GraphSimulationTravelTime, TravelTimeEdge>> +make_traveltime_sim(FP t0, FP dt, Graph, TravelTimeEdge> graph, + size_t n_steps = 120) { - return GraphSimulationTravelTime>(t0, dt, std::move(graph), n_steps); + return GraphSimulationTravelTime, TravelTimeEdge>>( + t0, dt, std::move(graph), n_steps); } } // namespace mio From a6df67e2886ed7fa4195e1da46222302b3db2000 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 27 Mar 2026 13:03:53 +0100 Subject: [PATCH 7/9] day to normalized time period, class for scheduler in own file, adjust tests --- .../metapopulation_mobility_traveltime.h | 232 ++---------- cpp/memilio/mobility/traveltime_schedule.h | 350 ++++++++++++++++++ cpp/tests/test_mobility_traveltime.cpp | 43 +-- 3 files changed, 397 insertions(+), 228 deletions(-) create mode 100644 cpp/memilio/mobility/traveltime_schedule.h diff --git a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h index 7e0257715c..552d83f036 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h @@ -23,6 +23,7 @@ #include "memilio/mobility/graph.h" #include "memilio/mobility/graph_simulation.h" #include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/mobility/traveltime_schedule.h" #include "memilio/math/eigen.h" #include "memilio/math/euler.h" #include "memilio/math/floating_point.h" @@ -173,182 +174,6 @@ void correct_negative_compartments(Eigen::Ref> vec, size_t nu } } -/** - * @brief Daily mobility schedule computed once from the graph topology. - * - * The day is discretised into `n_steps` equal time slices of width `1/n_steps`. - * For every edge the schedule records: - * - `node_at_step[edge][step]`: which node index holds the travelers at that step. - * - `in_mobility[edge][step]`: true, if the travelers are in a mobility model in that step. - * - `first_mobility_step[edge]`: first step at which mobility starts. - * - * Derived per-node schedules tell the integrator when it must stop and - * synchronise with the mobility exchange: - * - `local_breakpoints[node]`: steps at which local_sim must be synchronized. - * - `mobility_breakpoints[node]`: steps at which mobility_sim must be synchronized. - */ -struct TravelTimeSchedule { - std::vector> node_at_step; ///< [edge][step] -> node index - std::vector> in_mobility; ///< [edge][step] -> mobility flag - std::vector first_mobility_step; ///< [edge] -> first step with mobility - std::vector> local_breakpoints; ///< [node] -> sorted step indices - std::vector> mobility_breakpoints; ///< [node] -> sorted step indices - /// For each step: which edges exchange travelers - std::vector> edges_at_step; - /// For each step: which nodes advance their local_sim - std::vector> local_nodes_at_step; - /// For each step: which nodes advance their mobility_sim - std::vector> mobility_nodes_at_step; -}; - -/** - * @brief Compute the TravelTimeSchedule for a given graph. - * - * @tparam FP Floating-point type. - * @tparam NodeProp TravelTimeNodeProperty type. - * @tparam EdgeProp TravelTimeEdge type. - * @param graph The simulation graph. - * @param n_steps Number of time steps per day (default: 120). The value 120 corresponds to a resolution of - * 12 minutes per step and is divisible by many common travel-time fractions (e.g. 0.05, 0.1, - * 0.25), so typical commute durations map directly to exact integer step counts without rounding. - * @param eps Floating-point tolerance for time comparisons. - * @return Fully populated TravelTimeSchedule. - */ -template -TravelTimeSchedule compute_traveltime_schedule(const Graph, TravelTimeEdge>& graph, - size_t n_steps = 120, FP eps = Limits::zero_tolerance()) -{ - const size_t n_edges = graph.edges().size(); - const size_t n_nodes = graph.nodes().size(); - const double dt_step = 1.0 / static_cast(n_steps); - - TravelTimeSchedule sched; - sched.node_at_step.resize(n_edges); - sched.in_mobility.resize(n_edges); - sched.first_mobility_step.resize(n_edges, n_steps); - - // 1. Build per-edge schedules - auto all_edges = graph.edges(); - for (size_t ei = 0; ei < n_edges; ++ei) { - const auto& e = all_edges[ei]; - const auto& edge = e.property; - const size_t dest_node = e.end_node_idx; - const size_t src_node = e.start_node_idx; - const FP stay = graph.nodes()[dest_node].property.stay_duration; - const size_t n_patches = edge.path.size(); // includes start and end nodes - // Travel time per patch (split evenly). - const double tt_per_patch = - std::max(dt_step, std::round(edge.travel_time / static_cast(n_patches) / dt_step) * dt_step); - const double total_travel = tt_per_patch * static_cast(n_patches); - - // Start of outbound travel: everything before is local. - const double t_depart = std::max(0.0, 1.0 - 2.0 * total_travel - static_cast(stay)); - const size_t step_depart = static_cast((t_depart + eps) / dt_step); - // Arrival at destination. - const size_t steps_travel = static_cast((total_travel + eps) / dt_step); - const size_t step_arrive = step_depart + steps_travel; - // Departure from destination for return trip. - const size_t step_leave_dst = n_steps - steps_travel; - - sched.node_at_step[ei].resize(n_steps, src_node); - sched.in_mobility[ei].resize(n_steps, false); - - // Outbound travel: step_depart … step_arrive-1 - sched.first_mobility_step[ei] = step_depart; - size_t current_step = step_depart; - for (size_t patch = 0; patch < n_patches; ++patch) { - const size_t patch_steps = static_cast((tt_per_patch + eps) / dt_step); - for (size_t s = 0; s < patch_steps && (current_step + s) < n_steps; ++s) { - sched.node_at_step[ei][current_step + s] = edge.path[patch]; - sched.in_mobility[ei][current_step + s] = true; - } - current_step += patch_steps; - } - - // Stay at destination: step_arrive … step_leave_dst-1 - for (size_t s = step_arrive; s < step_leave_dst && s < n_steps; ++s) { - sched.node_at_step[ei][s] = dest_node; - sched.in_mobility[ei][s] = false; - } - - // Return travel: step_leave_dst … n_steps-1 (reversed path) - current_step = step_leave_dst; - for (size_t patch = n_patches; patch > 0 && current_step < n_steps; --patch) { - const size_t patch_steps = static_cast((tt_per_patch + eps) / dt_step); - for (size_t s = 0; s < patch_steps && (current_step + s) < n_steps; ++s) { - sched.node_at_step[ei][current_step + s] = edge.path[patch - 1]; - sched.in_mobility[ei][current_step + s] = true; - } - current_step += patch_steps; - } - } - - // 2. Derive per-node breakpoint schedules - sched.local_breakpoints.resize(n_nodes); - sched.mobility_breakpoints.resize(n_nodes); - - for (size_t ni = 0; ni < n_nodes; ++ni) { - sched.local_breakpoints[ni].push_back(0); - sched.mobility_breakpoints[ni].push_back(0); - - auto prev_local_edges = std::vector{}; - auto prev_mob_edges = std::vector{}; - for (size_t t = 1; t < n_steps; ++t) { - std::vector local_e, mob_e; - for (size_t ei = 0; ei < n_edges; ++ei) { - if (sched.node_at_step[ei][t] == ni) { - if (sched.in_mobility[ei][t]) - mob_e.push_back(ei); - else - local_e.push_back(ei); - } - } - if (local_e != prev_local_edges) { - sched.local_breakpoints[ni].push_back(t); - prev_local_edges = local_e; - } - if (mob_e != prev_mob_edges) { - sched.mobility_breakpoints[ni].push_back(t); - prev_mob_edges = mob_e; - } - } - // Ensure final step is always present. - if (sched.local_breakpoints[ni].back() != n_steps - 1) - sched.local_breakpoints[ni].push_back(n_steps - 1); - if (sched.mobility_breakpoints[ni].back() != n_steps - 1) - sched.mobility_breakpoints[ni].push_back(n_steps - 1); - } - - // 3. Build per-step lookup tables - sched.edges_at_step.resize(n_steps); - sched.local_nodes_at_step.resize(n_steps); - sched.mobility_nodes_at_step.resize(n_steps); - - for (size_t t = 0; t < n_steps; ++t) { - // Edges active at this step (first mobility step or node/mode change). - for (size_t ei = 0; ei < n_edges; ++ei) { - if (t < sched.first_mobility_step[ei]) - continue; - bool add = (t == sched.first_mobility_step[ei]); - if (!add && t > 0) { - add = (sched.node_at_step[ei][t] != sched.node_at_step[ei][t - 1]) || - (sched.in_mobility[ei][t] != sched.in_mobility[ei][t - 1]); - } - if (add) - sched.edges_at_step[t].push_back(ei); - } - // Nodes whose local_sim breakpoint falls at t. - for (size_t ni = 0; ni < n_nodes; ++ni) { - if (std::binary_search(sched.local_breakpoints[ni].begin(), sched.local_breakpoints[ni].end(), t)) - sched.local_nodes_at_step[t].push_back(ni); - if (std::binary_search(sched.mobility_breakpoints[ni].begin(), sched.mobility_breakpoints[ni].end(), t)) - sched.mobility_nodes_at_step[t].push_back(ni); - } - } - - return sched; -} - /** * @brief Graph simulation with travel-time-aware mobility. * @@ -360,14 +185,13 @@ TravelTimeSchedule compute_traveltime_schedule(const Graph t_depart --- outbound transit --> t_arrive --- stay --> t_return --- inbound transit --> 1 * - * The day is discretised into `n_steps`. The schedule is precomputed - * once by compute_traveltime_schedule before the simulation starts. - * + * The period is discretised into `n_steps` sub-steps. The schedule is precomputed + * once by the TravelTimeSchedule constructor before the simulation starts. * * @tparam FP Floating-point type (e.g., double). * @tparam GraphT Graph type (Graph, TravelTimeEdge>). @@ -384,16 +208,16 @@ class GraphSimulationTravelTime * @brief Construct the simulation. * * @param t0 Simulation start time. - * @param dt Time step. + * @param dt Time step (length of one normalized time period in simulation units). * @param graph The simulation graph. - * @param n_steps Number of sub-steps per day (default 100). + * @param n_steps Number of sub-steps per normalized time period (default 120). */ GraphSimulationTravelTime(FP t0, FP dt, GraphT graph, size_t n_steps = 120) : m_t(t0) , m_dt(dt) , m_graph(std::move(graph)) , m_n_steps(n_steps) - , m_schedule(compute_traveltime_schedule(m_graph, n_steps)) + , m_schedule(m_graph, n_steps, Limits::zero_tolerance()) { // Reset mobility_sim compartments to zero as transit nodes always start empty. for (auto& n : m_graph.nodes()) { @@ -486,7 +310,7 @@ class GraphSimulationTravelTime { FP earliest = FP{1}; for (size_t ei = 0; ei < m_graph.edges().size(); ++ei) { - const FP t_dep = static_cast(m_schedule.first_mobility_step[ei]) / m_n_steps; + const FP t_dep = static_cast(m_schedule.first_mobility_step(ei)) / m_n_steps; earliest = std::min(earliest, t_dep); } return earliest; @@ -494,8 +318,8 @@ class GraphSimulationTravelTime void move_commuters(size_t step, FP t_sub) { - for (size_t ei : m_schedule.edges_at_step[step]) { - if (step == m_schedule.first_mobility_step[ei]) { + for (size_t ei : m_schedule.edges_at(step)) { + if (step == m_schedule.first_mobility_step(ei)) { // First mobility step: extract commuters from local_sim of origin. init_mobility(ei, t_sub); } @@ -550,10 +374,10 @@ class GraphSimulationTravelTime if (edge.mobile_population.get_num_time_points() == 0) return; - const size_t prev_node = m_schedule.node_at_step[ei][step - 1]; - const size_t curr_node = m_schedule.node_at_step[ei][step]; - const bool prev_mob = m_schedule.in_mobility[ei][step - 1]; - const bool curr_mob = m_schedule.in_mobility[ei][step]; + const size_t prev_node = m_schedule.node_at(ei, step - 1); + const size_t curr_node = m_schedule.node_at(ei, step); + const bool prev_mob = m_schedule.in_mobility_at(ei, step - 1); + const bool curr_mob = m_schedule.in_mobility_at(ei, step); if (prev_node == curr_node && prev_mob == curr_mob) return; @@ -579,7 +403,7 @@ class GraphSimulationTravelTime void advance_local_sims(size_t step, FP t_sub, FP sub_dt) { - for (size_t ni : m_schedule.local_nodes_at_step[step]) { + for (size_t ni : m_schedule.local_nodes_at(step)) { // Determine the next breakpoint to find dt for this node. const FP node_dt = get_local_dt(ni, step, sub_dt); m_graph.nodes()[ni].property.local_sim.get_simulation().advance(t_sub + node_dt); @@ -588,10 +412,10 @@ class GraphSimulationTravelTime FP get_local_dt(size_t ni, size_t step, FP sub_dt) const { - const auto& bps = m_schedule.local_breakpoints[ni]; + const auto& bps = m_schedule.local_breakpoints(ni); const auto it = std::lower_bound(bps.begin(), bps.end(), step); if (it == bps.end() || std::next(it) == bps.end()) { - // Last breakpoint: advance to end of day. + // Last breakpoint: advance to end of normalized time period. return (static_cast(m_n_steps) - static_cast(step)) * sub_dt; } return static_cast(*std::next(it) - step) * sub_dt; @@ -600,16 +424,16 @@ class GraphSimulationTravelTime void update_commuter_states(size_t step, FP t_sub, FP sub_dt) { auto all_edges = m_graph.edges(); - for (size_t ei : m_schedule.edges_at_step[step]) { + for (size_t ei : m_schedule.edges_at(step)) { auto& e = all_edges[ei]; auto& edge = e.property; if (edge.mobile_population.get_num_time_points() == 0) continue; - if (!m_schedule.in_mobility[ei][step]) + if (!m_schedule.in_mobility_at(ei, step)) continue; // If not in transit, skip. - const size_t curr_node = m_schedule.node_at_step[ei][step]; + const size_t curr_node = m_schedule.node_at(ei, step); auto& mob_node = m_graph.nodes()[curr_node].property.mobility_sim; auto& model = mob_node.get_simulation().get_model(); @@ -637,8 +461,8 @@ class GraphSimulationTravelTime FP get_mobility_dt(size_t ei, size_t step, FP sub_dt) const { - const size_t curr_node = m_schedule.node_at_step[ei][step]; - const auto& bps = m_schedule.mobility_breakpoints[curr_node]; + const size_t curr_node = m_schedule.node_at(ei, step); + const auto& bps = m_schedule.mobility_breakpoints(curr_node); const auto it = std::lower_bound(bps.begin(), bps.end(), step); if (it == bps.end() || std::next(it) == bps.end()) { return (static_cast(m_n_steps) - static_cast(step)) * sub_dt; @@ -660,9 +484,9 @@ class GraphSimulationTravelTime // Final node in schedule is always the origin (home). const size_t home_node = e.start_node_idx; const size_t last_step = m_n_steps - 1; - const bool last_mob = m_schedule.in_mobility[ei][last_step]; - auto& sim_last = last_mob ? m_graph.nodes()[m_schedule.node_at_step[ei][last_step]].property.mobility_sim - : m_graph.nodes()[m_schedule.node_at_step[ei][last_step]].property.local_sim; + const bool last_mob = m_schedule.in_mobility_at(ei, last_step); + auto& sim_last = last_mob ? m_graph.nodes()[m_schedule.node_at(ei, last_step)].property.mobility_sim + : m_graph.nodes()[m_schedule.node_at(ei, last_step)].property.local_sim; Eigen::Ref> travel_subgroup = edge.mobile_population.get_last_value(); size_t num_age = static_cast( @@ -709,9 +533,9 @@ class GraphSimulationTravelTime * @tparam NodeProp TravelTimeNodeProperty specialisation. * @tparam EdgeProp TravelTimeEdge specialisation. * @param t0 Simulation start time. - * @param dt Day-level time step. + * @param dt Time step (length of one normalized time period in simulation units). * @param graph Simulation graph. - * @param n_steps Number of day sub-steps (default 100). + * @param n_steps Number of sub-steps per normalized time period (default 120). * @return GraphSimulationTravelTime instance. */ template diff --git a/cpp/memilio/mobility/traveltime_schedule.h b/cpp/memilio/mobility/traveltime_schedule.h new file mode 100644 index 0000000000..1dd54786ca --- /dev/null +++ b/cpp/memilio/mobility/traveltime_schedule.h @@ -0,0 +1,350 @@ +/* +* 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_TRAVELTIME_SCHEDULE_H +#define MIO_MOBILITY_TRAVELTIME_SCHEDULE_H + +#include "memilio/math/floating_point.h" + +#include +#include +#include +#include +#include + +namespace mio +{ + +/** + * @brief Precomputed mobility schedule for the travel-time-aware metapopulation model. + * + * Discretises one normalized time period [0, 1) into n_steps equal slices of width + * `1/n_steps`, and records for every graph edge and every sub-step: + * - which node index currently holds the edge's mobile sub-population, and + * - whether that sub-population is in a `mobility_sim` (in transit) or a + * `local_sim` (at home or at the destination). + * + * Derived per-node breakpoint lists tell the ODE integrator when it must stop and + * synchronise state with the commuter exchange: + * - `local_breakpoints(ni)`: sub-steps at which `local_sim` of node ni must synchronise. + * - `mobility_breakpoints(ni)`: sub-steps at which `mobility_sim` of node ni must synchronise. + * + * Per-step tables to accelerate the inner advance loop: + * - `edges_at(step)`: edges whose commuters change node or mode at this step. + * - `local_nodes_at(step)`: nodes whose `local_sim` breakpoint falls at this step. + * - `mobility_nodes_at(step)`: nodes whose `mobility_sim` breakpoint falls at this step. + * + * The schedule is computed once and then reused throughout the simulation. + * Use the constructor to build it and all data is read-only afterwards. + */ +class TravelTimeSchedule +{ +public: + /** + * @brief Compute the schedule from a graph. + * + * @tparam FP Floating-point type of the graph's node/edge properties. + * @tparam GraphT Graph type, expected to be + * `Graph, TravelTimeEdge>`. + * @param graph The simulation graph. + * @param n_steps Number of sub-steps per normalized time period (default: 120). + * 120 corresponds to 12-minute intervals when the period is one day + * and is evenly divisible by common commute-time fractions such as + * 0.05, 0.1, and 0.25, so typical trip durations map to exact integer + * step counts without rounding. + * @param eps Floating-point tolerance used when converting continuous time + * fractions to step indices to avoid off-by-one errors due to + * rounding. + */ + template + TravelTimeSchedule(const GraphT& graph, size_t n_steps = 120, FP eps = Limits::zero_tolerance()); + + /** + * @brief Node index at which edge mobile sub-population are at sub-step step. + * + * Before mobility starts (step < first_mobility_step(ei)) the edge's commuters are + * implicitly at the origin. Calling this function for those steps still returns the + * source node. + */ + size_t node_at(size_t ei, size_t step) const + { + assert(ei < m_n_edges && step < m_n_steps); + return m_node_at_step[ei * m_n_steps + step]; + } + + /** + * @brief Whether edge commuters are in a `mobility_sim` at sub-step step. + * + * Returns `true` during in-transit phases (outbound and return travel) and `false` + * during the stay at the destination or while still at the origin. + */ + bool in_mobility_at(size_t ei, size_t step) const + { + assert(ei < m_n_edges && step < m_n_steps); + return m_in_mobility[ei * m_n_steps + step]; + } + + /** + * @brief First sub-step at which edge commuters begin their outbound trip. + * + * All sub-steps before this index belong to the local phase at the origin. + * Equals n_steps if the edge is never active (e.g. zero travel time or zero + * mobility coefficient). + */ + size_t first_mobility_step(size_t ei) const + { + assert(ei < m_n_edges); + return m_first_mobility_step[ei]; + } + + /** + * @brief Sorted list of sub-steps at which node local_sim must synchronize. + * + * A breakpoint is inserted whenever the set of edges whose commuters are present in + * the local population of this node changes. The integrator must stop at each + * breakpoint to add or remove commuters before continuing. + */ + const std::vector& local_breakpoints(size_t ni) const + { + assert(ni < m_n_nodes); + return m_local_breakpoints[ni]; + } + + /** + * @brief Sorted list of sub-steps at which node mobility_sim must synchronize. + * + * Analogous to `local_breakpoints` but for the transit population simulation. + */ + const std::vector& mobility_breakpoints(size_t ni) const + { + assert(ni < m_n_nodes); + return m_mobility_breakpoints[ni]; + } + + /** + * @brief Edge indices whose commuters change their current node or mobility mode at + * sub-step step (including the first departure step). + * + * Used by the advance loop to efficiently process only the edges that require a + * commuter transfer at each step. + */ + const std::vector& edges_at(size_t step) const + { + assert(step < m_n_steps); + return m_edges_at_step[step]; + } + + /** + * @brief Indices of nodes whose `local_sim` breakpoint falls at sub-step step. + * + * Allows the advance loop to advance only the affected nodes at each sub-step. + */ + const std::vector& local_nodes_at(size_t step) const + { + assert(step < m_n_steps); + return m_local_nodes_at_step[step]; + } + + /** + * @brief Indices of nodes whose `mobility_sim` breakpoint falls at sub-step step. + */ + const std::vector& mobility_nodes_at(size_t step) const + { + assert(step < m_n_steps); + return m_mobility_nodes_at_step[step]; + } + + size_t num_steps() const + { + return m_n_steps; + } + size_t num_edges() const + { + return m_n_edges; + } + size_t num_nodes() const + { + return m_n_nodes; + } + +private: + size_t m_n_edges; + size_t m_n_nodes; + size_t m_n_steps; + + /// Flat per-edge/per-step node index: `m_node_at_step[ei * n_steps + step]` returns the node index at which + /// the edge's commuters are located at this step. + std::vector m_node_at_step; + + /// Flat per-edge/per-step in-transit flag: `m_in_mobility[ei * n_steps + step]` returns whether + /// the edge's commuters are in transit at this step. + std::vector m_in_mobility; + + std::vector m_first_mobility_step; ///< first step with commuter movement + + std::vector> m_local_breakpoints; ///< sorted step indices + std::vector> m_mobility_breakpoints; ///< sorted step indices + + std::vector> m_edges_at_step; ///< edge indices at each step + std::vector> m_local_nodes_at_step; ///< node indices at each step + std::vector> m_mobility_nodes_at_step; ///< node indices at each step +}; + +template +TravelTimeSchedule::TravelTimeSchedule(const GraphT& graph, size_t n_steps, FP eps) + : m_n_edges(graph.edges().size()) + , m_n_nodes(graph.nodes().size()) + , m_n_steps(n_steps) + , m_node_at_step(m_n_edges * n_steps, 0) + , m_in_mobility(m_n_edges * n_steps, false) + , m_first_mobility_step(m_n_edges, n_steps) + , m_local_breakpoints(m_n_nodes) + , m_mobility_breakpoints(m_n_nodes) + , m_edges_at_step(n_steps) + , m_local_nodes_at_step(n_steps) + , m_mobility_nodes_at_step(n_steps) +{ + const double dt_step = 1.0 / static_cast(n_steps); + const double eps_d = static_cast(eps); + + auto all_edges = graph.edges(); + + // 1. Build per-edge schedules: initialise default node to source, then overwrite the mobility and stay phases + for (size_t ei = 0; ei < m_n_edges; ++ei) { + const size_t src = all_edges[ei].start_node_idx; + std::fill(m_node_at_step.begin() + static_cast(ei * n_steps), + m_node_at_step.begin() + static_cast((ei + 1) * n_steps), src); + } + + for (size_t ei = 0; ei < m_n_edges; ++ei) { + const auto& e = all_edges[ei]; + const auto& edge = e.property; + const size_t dest_node = e.end_node_idx; + const size_t n_patches = edge.path.size(); // includes origin and destination nodes + + // Travel time per path patch, snapped to the nearest step boundary. + const double tt_per_patch = std::max( + dt_step, + std::round(static_cast(edge.travel_time) / static_cast(n_patches) / dt_step) * dt_step); + const double total_travel = tt_per_patch * static_cast(n_patches); + + // The normalized time period is partitioned as: + // [0, t_depart) local at origin + // [t_depart, t_arrive) outbound transit + // [t_arrive, t_leave_dst) stay at destination + // [t_leave_dst, 1) return transit + const double t_depart = std::max(0.0, 1.0 - 2.0 * total_travel - + static_cast(graph.nodes()[dest_node].property.stay_duration)); + + const size_t step_depart = static_cast((t_depart + eps_d) / dt_step); + const size_t steps_travel = static_cast((total_travel + eps_d) / dt_step); + const size_t step_arrive = step_depart + steps_travel; + const size_t step_leave_dst = n_steps - steps_travel; + + m_first_mobility_step[ei] = step_depart; + + // Outbound travel: step_depart … step_arrive-1 + size_t current_step = step_depart; + for (size_t patch = 0; patch < n_patches; ++patch) { + const size_t patch_steps = static_cast((tt_per_patch + eps_d) / dt_step); + for (size_t s = 0; s < patch_steps && (current_step + s) < n_steps; ++s) { + m_node_at_step[ei * n_steps + current_step + s] = edge.path[patch]; + m_in_mobility[ei * n_steps + current_step + s] = true; + } + current_step += patch_steps; + } + + // Stay at destination: step_arrive … step_leave_dst-1 + for (size_t s = step_arrive; s < step_leave_dst && s < n_steps; ++s) { + m_node_at_step[ei * n_steps + s] = dest_node; + m_in_mobility[ei * n_steps + s] = false; + } + + // Return travel: step_leave_dst … n_steps-1 (reversed path) + current_step = step_leave_dst; + for (size_t patch = n_patches; patch > 0 && current_step < n_steps; --patch) { + const size_t patch_steps = static_cast((tt_per_patch + eps_d) / dt_step); + for (size_t s = 0; s < patch_steps && (current_step + s) < n_steps; ++s) { + m_node_at_step[ei * n_steps + current_step + s] = edge.path[patch - 1]; + m_in_mobility[ei * n_steps + current_step + s] = true; + } + current_step += patch_steps; + } + } + + // 2. Derive per-node breakpoint schedules + for (size_t ni = 0; ni < m_n_nodes; ++ni) { + m_local_breakpoints[ni].push_back(0); + m_mobility_breakpoints[ni].push_back(0); + + auto prev_local_edges = std::vector{}; + auto prev_mob_edges = std::vector{}; + for (size_t t = 1; t < n_steps; ++t) { + std::vector local_e, mob_e; + for (size_t ei = 0; ei < m_n_edges; ++ei) { + if (m_node_at_step[ei * n_steps + t] == ni) { + if (m_in_mobility[ei * n_steps + t]) + mob_e.push_back(ei); + else + local_e.push_back(ei); + } + } + if (local_e != prev_local_edges) { + m_local_breakpoints[ni].push_back(t); + prev_local_edges = local_e; + } + if (mob_e != prev_mob_edges) { + m_mobility_breakpoints[ni].push_back(t); + prev_mob_edges = mob_e; + } + } + // Ensure the final sub-step is always present as a breakpoint. + if (m_local_breakpoints[ni].back() != n_steps - 1) + m_local_breakpoints[ni].push_back(n_steps - 1); + if (m_mobility_breakpoints[ni].back() != n_steps - 1) + m_mobility_breakpoints[ni].push_back(n_steps - 1); + } + + // 3. Build per-step lookup tables + for (size_t t = 0; t < n_steps; ++t) { + // Edges whose commuters start moving or change node/mode at this step. + for (size_t ei = 0; ei < m_n_edges; ++ei) { + if (t < m_first_mobility_step[ei]) + continue; + bool add = (t == m_first_mobility_step[ei]); + if (!add && t > 0) { + add = (m_node_at_step[ei * n_steps + t] != m_node_at_step[ei * n_steps + t - 1]) || + (m_in_mobility[ei * n_steps + t] != m_in_mobility[ei * n_steps + t - 1]); + } + if (add) + m_edges_at_step[t].push_back(ei); + } + // Nodes with a breakpoint at this step. + for (size_t ni = 0; ni < m_n_nodes; ++ni) { + if (std::binary_search(m_local_breakpoints[ni].begin(), m_local_breakpoints[ni].end(), t)) + m_local_nodes_at_step[t].push_back(ni); + if (std::binary_search(m_mobility_breakpoints[ni].begin(), m_mobility_breakpoints[ni].end(), t)) + m_mobility_nodes_at_step[t].push_back(ni); + } + } +} + +} // namespace mio + +#endif // MIO_MOBILITY_TRAVELTIME_SCHEDULE_H diff --git a/cpp/tests/test_mobility_traveltime.cpp b/cpp/tests/test_mobility_traveltime.cpp index 1ddfbf9ca9..2bb8847744 100644 --- a/cpp/tests/test_mobility_traveltime.cpp +++ b/cpp/tests/test_mobility_traveltime.cpp @@ -20,12 +20,15 @@ #include "memilio/mobility/metapopulation_mobility_instant.h" #include "memilio/mobility/metapopulation_mobility_traveltime.h" +#include "memilio/mobility/traveltime_schedule.h" #include "memilio/compartments/simulation.h" #include "memilio/math/euler.h" +#include "memilio/math/floating_point.h" #include "ode_seir/model.h" #include "ode_seir/infection_state.h" #include "gtest/gtest.h" +#include "utils.h" #include #include @@ -67,32 +70,29 @@ TEST(TravelTimeMobility, ScheduleBasicStructure) const size_t n_steps = 100; auto graph = make_two_node_graph(1000.0, 10.0, 800.0, 0.0, travel_time, stay_dur, 0.1); - auto sched = mio::compute_traveltime_schedule(graph, n_steps); + auto sched = mio::TravelTimeSchedule(graph, n_steps, mio::Limits::zero_tolerance()); // Schedule must have one entry per edge. - ASSERT_EQ(sched.node_at_step.size(), graph.edges().size()); - ASSERT_EQ(sched.in_mobility.size(), graph.edges().size()); + ASSERT_EQ(sched.num_edges(), graph.edges().size()); // Each edge schedule must have n_steps entries. - for (size_t ei = 0; ei < graph.edges().size(); ++ei) { - ASSERT_EQ(sched.node_at_step[ei].size(), n_steps); - ASSERT_EQ(sched.in_mobility[ei].size(), n_steps); - } + ASSERT_EQ(sched.num_steps(), n_steps); // Breakpoints must be non-empty for every node. - ASSERT_EQ(sched.local_breakpoints.size(), graph.nodes().size()); + ASSERT_EQ(sched.num_nodes(), graph.nodes().size()); for (size_t ni = 0; ni < graph.nodes().size(); ++ni) { - EXPECT_FALSE(sched.local_breakpoints[ni].empty()) << "local_breakpoints empty for node " << ni; + EXPECT_FALSE(sched.local_breakpoints(ni).empty()) << "local_breakpoints empty for node " << ni; } // first_mobility_step must be valid indices. for (size_t ei = 0; ei < graph.edges().size(); ++ei) { - EXPECT_LT(sched.first_mobility_step[ei], n_steps) << "first_mobility_step out of range for edge " << ei; + EXPECT_LT(sched.first_mobility_step(ei), n_steps) << "first_mobility_step out of range for edge " << ei; } } TEST(TravelTimeMobility, PopulationConservation) { + mio::LogLevelOverride llo(mio::LogLevel::off); const FP t0 = 0.0; const FP tmax = 3.0; const FP travel_time = 0.04; @@ -141,6 +141,7 @@ TEST(TravelTimeMobility, ZeroMobilityCoefficient) TEST(TravelTimeMobility, LargeTravelTimeDoesNotCrash) { + mio::LogLevelOverride llo(mio::LogLevel::off); // suppress expected errors from correct_negative_compartments const FP t0 = 0.0; // Nearly all-day travel: 49% + 1% stay + 49% return = 99% of day in travel. @@ -155,10 +156,10 @@ TEST(TravelTimeMobility, LargeTravelTimeDoesNotCrash) TEST(TravelTimeMobility, CorrectNegativeCompartments) { // 2 age groups, 4 compartments each == size 8. - // Use values below the default tolerance (-1e-7) to trigger correction. + // Negative values are corrected via map_to_nonnegative per age group slice. Eigen::VectorX v(8); - v << 100.0, 50.0, -1e-5, 30.0, // group 0: negative below tolerance at index 2 - 200.0, 80.0, 40.0, -2e-5; // group 1: negative below tolerance at index 7 + v << 100.0, 50.0, -1e-5, 30.0, // group 0: negative at index 2 + 200.0, 80.0, 40.0, -2e-5; // group 1: negative at index 7 const FP sum0_before = v.head(4).sum(); const FP sum1_before = v.tail(4).sum(); @@ -172,12 +173,6 @@ TEST(TravelTimeMobility, CorrectNegativeCompartments) // Population within each age group is conserved EXPECT_NEAR(v.head(4).sum(), sum0_before, 1e-12) << "Population not conserved in group 0"; EXPECT_NEAR(v.tail(4).sum(), sum1_before, 1e-12) << "Population not conserved in group 1"; - - // Values that are within tolerance (-1e-7) should NOT be corrected - Eigen::VectorX v2(4); - v2 << 100.0, 50.0, -1e-8, 30.0; - mio::correct_negative_compartments(v2, 1); - EXPECT_LT(v2[2], FP{0}) << "Value above tolerance should not have been corrected"; } TEST(TravelTimeMobility, ScheduleNodeAssignmentTwoNode) @@ -187,21 +182,21 @@ TEST(TravelTimeMobility, ScheduleNodeAssignmentTwoNode) const size_t n_steps = 100; auto graph = make_two_node_graph(1000.0, 0.0, 800.0, 0.0, travel_time, stay_dur, 0.1); - auto sched = mio::compute_traveltime_schedule(graph, n_steps); + auto sched = mio::TravelTimeSchedule(graph, n_steps, mio::Limits::zero_tolerance()); // Edge 0->1 is edge index 0 (first edge added). const size_t ei = 0; - ASSERT_LT(ei, sched.node_at_step.size()); + ASSERT_LT(ei, sched.num_edges()); - const size_t first = sched.first_mobility_step[ei]; + const size_t first = sched.first_mobility_step(ei); EXPECT_LT(first, n_steps); // During the stay phase, in_mobility must be false and node must be dest (1). bool found_stay = false; for (size_t s = first; s < n_steps; ++s) { - if (!sched.in_mobility[ei][s]) { + if (!sched.in_mobility_at(ei, s)) { found_stay = true; - EXPECT_EQ(sched.node_at_step[ei][s], size_t{1}) << "Expected destination node during stay at step " << s; + EXPECT_EQ(sched.node_at(ei, s), size_t{1}) << "Expected destination node during stay at step " << s; } } EXPECT_TRUE(found_stay) << "No stay phase found in edge schedule"; From 0203630cde3d139bd960c68a99ef286c472f9d9e Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 27 Mar 2026 13:12:30 +0100 Subject: [PATCH 8/9] [ci skip] precompute first_departure_time from step --- .../metapopulation_mobility_traveltime.h | 3 +-- cpp/memilio/mobility/traveltime_schedule.h | 21 ++++++++++++++++--- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h index 552d83f036..808d2ae055 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h @@ -310,8 +310,7 @@ class GraphSimulationTravelTime { FP earliest = FP{1}; for (size_t ei = 0; ei < m_graph.edges().size(); ++ei) { - const FP t_dep = static_cast(m_schedule.first_mobility_step(ei)) / m_n_steps; - earliest = std::min(earliest, t_dep); + earliest = std::min(earliest, static_cast(m_schedule.first_departure_time(ei))); } return earliest; } diff --git a/cpp/memilio/mobility/traveltime_schedule.h b/cpp/memilio/mobility/traveltime_schedule.h index 1dd54786ca..f86dc8eb17 100644 --- a/cpp/memilio/mobility/traveltime_schedule.h +++ b/cpp/memilio/mobility/traveltime_schedule.h @@ -101,7 +101,7 @@ class TravelTimeSchedule } /** - * @brief First sub-step at which edge commuters begin their outbound trip. + * @brief First sub-step index at which edge commuters begin their outbound trip. * * All sub-steps before this index belong to the local phase at the origin. * Equals n_steps if the edge is never active (e.g. zero travel time or zero @@ -113,6 +113,18 @@ class TravelTimeSchedule return m_first_mobility_step[ei]; } + /** + * @brief Precomputed departure time for edge ei as a fraction of the normalized + * time period, i.e. `first_mobility_step(ei) / n_steps`. + * + * Stored explicitly to avoid the division in the inner advance loop. + */ + double first_departure_time(size_t ei) const + { + assert(ei < m_n_edges); + return m_first_departure_time[ei]; + } + /** * @brief Sorted list of sub-steps at which node local_sim must synchronize. * @@ -196,7 +208,8 @@ class TravelTimeSchedule /// the edge's commuters are in transit at this step. std::vector m_in_mobility; - std::vector m_first_mobility_step; ///< first step with commuter movement + std::vector m_first_mobility_step; ///< first step index with commuter movement + std::vector m_first_departure_time; ///< precomputed departure time = first_mobility_step / n_steps std::vector> m_local_breakpoints; ///< sorted step indices std::vector> m_mobility_breakpoints; ///< sorted step indices @@ -214,6 +227,7 @@ TravelTimeSchedule::TravelTimeSchedule(const GraphT& graph, size_t n_steps, FP e , m_node_at_step(m_n_edges * n_steps, 0) , m_in_mobility(m_n_edges * n_steps, false) , m_first_mobility_step(m_n_edges, n_steps) + , m_first_departure_time(m_n_edges, 1.0) , m_local_breakpoints(m_n_nodes) , m_mobility_breakpoints(m_n_nodes) , m_edges_at_step(n_steps) @@ -257,7 +271,8 @@ TravelTimeSchedule::TravelTimeSchedule(const GraphT& graph, size_t n_steps, FP e const size_t step_arrive = step_depart + steps_travel; const size_t step_leave_dst = n_steps - steps_travel; - m_first_mobility_step[ei] = step_depart; + m_first_mobility_step[ei] = step_depart; + m_first_departure_time[ei] = static_cast(step_depart) * dt_step; // Outbound travel: step_depart … step_arrive-1 size_t current_step = step_depart; From ec279fbf35085a2cd0a829aff59199d6e08b2a6a Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 27 Mar 2026 14:36:07 +0100 Subject: [PATCH 9/9] Reverted correct_negative_compartments to the iterative max->min redistribution approach --- .../ode_secir_traveltime_mobility.cpp | 30 +++++++------- .../metapopulation_mobility_traveltime.h | 40 +++++++++++++------ 2 files changed, 42 insertions(+), 28 deletions(-) diff --git a/cpp/examples/ode_secir_traveltime_mobility.cpp b/cpp/examples/ode_secir_traveltime_mobility.cpp index 9d43b4d473..8514b47682 100644 --- a/cpp/examples/ode_secir_traveltime_mobility.cpp +++ b/cpp/examples/ode_secir_traveltime_mobility.cpp @@ -106,19 +106,22 @@ int main() // -- Mobility parameters ----------------------------------------------- // A -> B trip passes through C (path = {A, C, B}, n_patches = 3). - // travel_time_AB = 2 * tt_patch_AC_CB; with n_steps=100 this yields tt_per_patch = 0.03 (3 sub-steps/patch). - // step_depart_AB = 48 -> A's commuters are in C's mobility_sim at sub-steps 51–53 (outbound) - // and at sub-steps 94–96 (return). + // travel_time_AB = 2 * tt_patch_AC_CB; with n_steps=120 this yields tt_per_patch = 3 sub-steps/patch. + // step_depart_AB = 62 -> A's commuters are in C's mobility_sim at outbound sub-steps 65–67 + // and at return sub-steps 114–116. // // B <-> C (path = {B,C} or {C,B}, n_patches = 2). - // tt_direct_BC = 2/24 -> tt_per_patch = 0.04 (4 sub-steps/patch). - // C->B: step_depart = 50 -> C->B commuters are in C's mobility_sim at sub-steps 50–53 - // -> overlaps with A's outbound transit (51–53) in C's mobility_sim. - // B->C return: in C's mobility_sim at sub-steps 92–95 - // -> overlaps with A's return transit (94–96) in C's mobility_sim. + // tt_direct_BC = 2/24 -> tt_per_patch = 5 sub-steps/patch. + // C->B: step_depart = 60 -> C->B commuters in C's mobility_sim at outbound sub-steps 60–64, + // and return sub-steps 115–119. + // B->C: step_depart = 70 -> B->C commuters in C's mobility_sim at outbound sub-steps 75–79, + // and return sub-steps 110–114. + // During the return phase, all three groups share C's mobility_sim: + // A->B return (114–116), B->C return (110–114) at sub-step 114 + // A->B return (114–116), C->B return (115–119) at sub-steps 115–116 // All three groups therefore share C's mobility_sim and can infect each other during transit. const FP tt_patch_AC_CB = 1.0 / 24.0; // 1 h per patch (A <-> C) - const FP tt_direct_BC = 2.0 / 24.0; // 2 h one-way B <-> C (tt_per_patch=0.04 -> schedule overlap in C) + const FP tt_direct_BC = 2.0 / 24.0; // 2 h one-way B <-> C (5 sub-steps/patch -> schedule overlap in C) const FP stay_dur_AB = 8.0 / 24.0; const FP stay_dur_C = 6.0 / 24.0; @@ -160,13 +163,8 @@ int main() print_row("C", sim.get_graph().nodes()[2].property); std::cout << std::string(62, '-') << "\n"; - std::cout << "\nSchedule overlap (n_steps=100, sub-step = 0.01 day):\n" - << " A outbound through C : sub-steps 51-53\n" - << " C->B outbound in C : sub-steps 50-53 -> overlap 51-53\n" - << " A return through C : sub-steps 94-96\n" - << " B->C return in C : sub-steps 92-95 -> overlap 94-95\n" - << "\nAll three groups share C's mobility_sim during these steps,\n" - << "so infections can spread between A-transit commuters and C<->B commuters.\n" + std::cout << "\nDuring the return phase, A-transit commuters and C<->B commuters\n" + << "share C's mobility_sim and can infect each other during transit.\n" << "This effect is not captured by standard metapopulation models\n" << "without explicit travel time.\n"; diff --git a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h index 808d2ae055..75d5957a4f 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_traveltime.h +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h @@ -153,23 +153,39 @@ struct TravelTimeEdge { * step to estimate the infection states of travelers. * * The vector is partitioned into `num_age_groups` equally sized blocks. For each - * block, mio::map_to_nonnegative is applied, which redistributes negative values - * proportionally across the positive entries while preserving the nonnegative sum. + * block, negative values are iteratively redistributed to the largest positive + * compartment in the same block until all values are within `tolerance` of zero, + * preserving the total population sum exactly. * - * @tparam FP Floating-point type. - * @param vec Compartment vector to correct (in-place). - * @param num_age_groups Number of age groups + * @tparam FP Floating-point type. + * @param vec Compartment vector to correct (in-place). + * @param num_age_groups Number of age groups. + * @param tolerance Values above this threshold are considered nonnegative. + * @param max_iter Maximum redistribution iterations per age group. */ template -void correct_negative_compartments(Eigen::Ref> vec, size_t num_age_groups) +void correct_negative_compartments(Eigen::Ref> vec, size_t num_age_groups, + FP tolerance = static_cast(-1e-7), size_t max_iter = 100) { - const Eigen::Index n_comparts = static_cast(vec.size()) / static_cast(num_age_groups); + const size_t n_comparts = vec.size() / num_age_groups; for (size_t grp = 0; grp < num_age_groups; ++grp) { - auto slice = vec.segment(static_cast(grp) * n_comparts, n_comparts); - auto result = map_to_nonnegative(slice); - if (!result) { - log_error("correct_negative_compartments: could not map age group {} to nonnegative values: {}", grp, - result.error().message()); + const auto beg = static_cast(grp * n_comparts); + const auto len = static_cast(n_comparts); + auto slice = vec.segment(beg, len); + for (size_t iter = 0; iter < max_iter; ++iter) { + Eigen::Index min_idx; + if (slice.minCoeff(&min_idx) >= tolerance) { + break; + } + Eigen::Index max_idx; + slice.maxCoeff(&max_idx); + slice(max_idx) += slice(min_idx); + slice(min_idx) = FP{0}; + } + if (slice.minCoeff() < tolerance) { + log_error("correct_negative_compartments: could not correct all negative values in age group {} " + "after {} iterations.", + grp, max_iter); } } }