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..8514b47682 --- /dev/null +++ b/cpp/examples/ode_secir_traveltime_mobility.cpp @@ -0,0 +1,172 @@ +/* +* 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 (path = {A, C, B}, n_patches = 3). + // 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 = 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 (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; + + 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 << "\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"; + + 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..75d5957a4f --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_traveltime.h @@ -0,0 +1,567 @@ +/* +* 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_METAPOPULATION_MOBILITY_TRAVELTIME_H +#define MIO_MOBILITY_METAPOPULATION_MOBILITY_TRAVELTIME_H + +#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" +#include "memilio/math/math_utils.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. For each + * 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. + * @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, + FP tolerance = static_cast(-1e-7), size_t max_iter = 100) +{ + const size_t n_comparts = vec.size() / num_age_groups; + for (size_t grp = 0; grp < num_age_groups; ++grp) { + 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); + } + } +} + +/** + * @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 schedule with multiple synchronisation points per normalized time period + * + * Structure of one normalized time period [0, 1): + * 0 --- local --> t_depart --- outbound transit --> t_arrive --- stay --> t_return --- inbound transit --> 1 + * + * 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>). + */ +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 (length of one normalized time period in simulation units). + * @param graph The simulation graph. + * @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(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()) { + 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 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 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. + const FP t_first_departure = compute_first_departure(); + + 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()) { + 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}; + + // Reduce to daily resolution after interpolation_threshold days to limit memory usage. + if (m_t > interpolation_threshold) { + reduce_to_daily_resolution(); + } + } + } + +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) { + earliest = std::min(earliest, static_cast(m_schedule.first_departure_time(ei))); + } + return earliest; + } + + void move_commuters(size_t step, FP t_sub) + { + 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); + } + else { + // Subsequent steps: move commuters from previous to current transit node. + transfer_commuters(ei, step); + } + } + } + + /// 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 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; + + 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) + { + 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; + + 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; + + 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)) { + // 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 normalized time period. + 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) + { + auto all_edges = m_graph.edges(); + 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_at(ei, step)) + continue; // If not in transit, skip. + + 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(); + + // 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(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() + { + // Return all still-mobile commuters to their home local_sim. + 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) + 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_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( + 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(); + } + } + + /// Keep only one time point per day in local_sim results (daily resolution). + void reduce_to_daily_resolution() + { + 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); + } + } + } +}; // class GraphSimulationTravelTime + +/** + * @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 Time step (length of one normalized time period in simulation units). + * @param graph Simulation graph. + * @param n_steps Number of sub-steps per normalized time period (default 120). + * @return GraphSimulationTravelTime instance. + */ +template +GraphSimulationTravelTime, TravelTimeEdge>> +make_traveltime_sim(FP t0, FP dt, Graph, TravelTimeEdge> graph, + size_t n_steps = 120) +{ + return GraphSimulationTravelTime, TravelTimeEdge>>( + t0, dt, std::move(graph), n_steps); +} + +} // namespace mio + +#endif // MIO_MOBILITY_METAPOPULATION_MOBILITY_TRAVELTIME_H diff --git a/cpp/memilio/mobility/traveltime_schedule.h b/cpp/memilio/mobility/traveltime_schedule.h new file mode 100644 index 0000000000..f86dc8eb17 --- /dev/null +++ b/cpp/memilio/mobility/traveltime_schedule.h @@ -0,0 +1,365 @@ +/* +* 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 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 + * mobility coefficient). + */ + size_t first_mobility_step(size_t ei) const + { + assert(ei < m_n_edges); + 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. + * + * 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 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 + + 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_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) + , 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; + m_first_departure_time[ei] = static_cast(step_depart) * dt_step; + + // 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/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..2bb8847744 --- /dev/null +++ b/cpp/tests/test_mobility_traveltime.cpp @@ -0,0 +1,203 @@ +/* +* 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/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 + +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::TravelTimeSchedule(graph, n_steps, mio::Limits::zero_tolerance()); + + // Schedule must have one entry per edge. + ASSERT_EQ(sched.num_edges(), graph.edges().size()); + + // Each edge schedule must have n_steps entries. + ASSERT_EQ(sched.num_steps(), n_steps); + + // Breakpoints must be non-empty for every node. + 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; + } + + // 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) +{ + mio::LogLevelOverride llo(mio::LogLevel::off); + 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) +{ + 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. + 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. + // 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 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(); + + 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"; +} + +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::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.num_edges()); + + 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_at(ei, s)) { + found_stay = true; + 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"; +}