diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index e7c5cb99f4..55f5d500e9 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -76,6 +76,10 @@ add_executable(ode_secirvvs_example ode_secirvvs.cpp) target_link_libraries(ode_secirvvs_example PRIVATE memilio ode_secirvvs) target_compile_options(ode_secirvvs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ode_seir_detailed_mobility ode_seir_detailed_mobility.cpp) +target_link_libraries(ode_seir_detailed_mobility PRIVATE memilio ode_seir) +target_compile_options(ode_seir_detailed_mobility PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(ode_secirts_example ode_secirts.cpp) target_link_libraries(ode_secirts_example PRIVATE memilio ode_secirts) target_compile_options(ode_secirts_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/examples/ode_seir_detailed_mobility.cpp b/cpp/examples/ode_seir_detailed_mobility.cpp new file mode 100644 index 0000000000..d757fda332 --- /dev/null +++ b/cpp/examples/ode_seir_detailed_mobility.cpp @@ -0,0 +1,101 @@ +/* +* Copyright (C) 2020-2024 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 "ode_seir/model.h" +#include "ode_seir/infection_state.h" +#include "ode_seir/parameters.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" +#include "memilio/compartments/simulation.h" +#include "memilio/compartments/flow_simulation.h" +#include "memilio/data/analyze_result.h" + +int main() +{ + const auto t0 = 0.; + const auto tmax = 10.; + const auto dt = 0.5; //time step of migration, daily migration every second step + + mio::oseir::Model<> model(1); + + // set population + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000; + + // set transition times + model.parameters.set>(1); + model.parameters.set>(1); + + // set contact matrix + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(2.7); + + //two mostly identical groups + auto model_group1 = model; + auto model_group2 = model; + + //some contact restrictions in group 1 + mio::ContactMatrixGroup& contact_matrix1 = + model_group1.parameters.get>().get_cont_freq_mat(); + contact_matrix1[0].add_damping(0.5, mio::SimulationTime(5)); + + //infection starts in group 1 + model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9990; + model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 10; + + auto sim1 = mio::FlowSimulation>(model_group1, t0, dt); + auto sim2 = mio::FlowSimulation>(model_group2, t0, dt); + double stay_time_1 = 0.3; + double stay_time_2 = 0.3; + // mio::ExtendedGraph>> g; + mio::ExtendedGraph>> g; + g.add_node(1, sim1, sim2, stay_time_1); + g.add_node(2, sim1, sim2, stay_time_2); + + double traveltime = 0.1; + std::vector path1 = {0, 1}; + std::vector path2 = {1, 0}; + g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01), traveltime, path1); + g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01), traveltime, path2); + + auto sim = mio::make_mobility_sim(t0, dt, std::move(g)); + + sim.advance(tmax); + // results node 1 + std::cout << "Results node 1" << std::endl; + auto interpolated_sim1 = + mio::interpolate_simulation_result(sim.get_graph().nodes()[0].property.base_sim.get_result()); + interpolated_sim1.print_table({"S", "E", "I", "R"}); + + // results node 1 mobility_sim + std::cout << "Mobility results node 1" << std::endl; + auto interpolated_sim1_mobility = sim.get_graph().nodes()[0].property.mobility_sim.get_result(); + interpolated_sim1_mobility.print_table({"S", "E", "I", "R"}); + + // results node 2 + std::cout << "Results node 2" << std::endl; + auto interpolated_sim2 = sim.get_graph().nodes()[1].property.base_sim.get_result(); + interpolated_sim2.print_table({"S", "E", "I", "R"}); + + // results node 2 mobility_sim + std::cout << "Mobility results node 2" << std::endl; + auto interpolated_sim2_mobility = sim.get_graph().nodes()[1].property.mobility_sim.get_result(); + interpolated_sim2_mobility.print_table({"S", "E", "I", "R"}); + + return 0; +} diff --git a/cpp/memilio/CMakeLists.txt b/cpp/memilio/CMakeLists.txt index b62d972e9a..f50dd6122a 100644 --- a/cpp/memilio/CMakeLists.txt +++ b/cpp/memilio/CMakeLists.txt @@ -66,6 +66,8 @@ add_library(memilio mobility/metapopulation_mobility_instant.cpp mobility/metapopulation_mobility_stochastic.h mobility/metapopulation_mobility_stochastic.cpp + mobility/metapopulation_mobility_detailed.cpp + mobility/metapopulation_mobility_detailed.h mobility/graph_simulation.h mobility/graph_simulation.cpp mobility/graph.h diff --git a/cpp/memilio/compartments/flow_model.h b/cpp/memilio/compartments/flow_model.h index 4825474769..4cfd820263 100644 --- a/cpp/memilio/compartments/flow_model.h +++ b/cpp/memilio/compartments/flow_model.h @@ -55,8 +55,8 @@ using filtered_tuple_t = decltype(filter_tuple(std::declval() // Remove all occurrences of OmittedTag from the types in an Index = IndexTemplate. template class IndexTemplate, class Index> -using filtered_index_t = decltype(as_index( - std::declval()))>>())); +using filtered_index_t = decltype( + as_index(std::declval()))>>())); } //namespace details @@ -203,6 +203,26 @@ class FlowModel : public CompartmentalModel return index_of_type_v, Flows>; } + /** + * @brief Returns the current flow values. + * + * @return A constant reference to an Eigen::VectorX containing the current flow values. + */ + Eigen::VectorX& get_flow_values() const + { + return m_flow_values; + } + + /** + * @brief Sets the flow values. + * + * @param flows A constant reference to an Eigen::VectorX containing flow values. + */ + void set_flow_values(const Eigen::VectorX flows) + { + m_flow_values = flows; + } + private: mutable Eigen::VectorX m_flow_values; ///< Cache to avoid allocation in get_derivatives (using get_flows). diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 1a78664ecc..746c95915f 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -35,6 +35,35 @@ #include #include +#include + +// Check if a type has a member called 'stay_duration' +template +struct has_stay_duration : std::false_type { +}; + +template +struct has_stay_duration().stay_duration)>> : std::true_type { +}; + +// Check if a type has a member called 'travel_time' +template +struct has_travel_time : std::false_type { +}; + +template +struct has_travel_time().travel_time)>> : std::true_type { +}; + +// Check if a type has a member called 'path' +template +struct has_path : std::false_type { +}; + +template +struct has_path().path)>> : std::true_type { +}; + namespace mio { @@ -42,25 +71,15 @@ namespace mio * Class that performs multiple simulation runs with randomly sampled parameters. * Can simulate mobility graphs with one simulation in each node or single simulations. * @tparam S type of simulation that runs in one node of the graph. + * @tparam ParametersGraph stores the parameters of the simulation. This is the input of ParameterStudies. + * @tparam SimulationGraph stores simulations and their results of each run. This is the output of ParameterStudies for each run. */ -template +template >, + class SimulationGraph = Graph, MobilityEdge>> class ParameterStudy { public: - /** - * The type of simulation of a single node of the graph. - */ using Simulation = S; - /** - * The Graph type that stores the parametes of the simulation. - * This is the input of ParameterStudies. - */ - using ParametersGraph = mio::Graph>; - /** - * The Graph type that stores simulations and their results of each run. - * This is the output of ParameterStudies for each run. - */ - using SimulationGraph = mio::Graph, mio::MobilityEdge>; /** * create study for graph of compartment models. @@ -336,18 +355,84 @@ class ParameterStudy } private: - //sample parameters and create simulation + /** + * @brief Adds a node to the graph using the overload based on the node's property. + * + * This function checks whether the node property type has a 'stay_duration' member. + * If so, it uses the specialized version that extracts additional simulation parameters. + * Otherwise, we just use the simulation time (m_t0) and integration step size (m_dt_integration). + * + * @tparam GraphType The type of the graph. + * @tparam NodeType The type of the node. + * @param graph The graph to which the node will be added. + * @param node The node to add. + */ + template + void add_node(GraphType& graph, const NodeType& node) + { + // Deduce the property type of the node. + using PropertyType = std::decay_t; + + // If the property type has a member called 'stay_duration', use the specialized version. + if constexpr (has_stay_duration::value) { + graph.add_node(node.id, node.property.base_sim, node.property.mobility_sim, node.property.stay_duration); + } + else { + graph.add_node(node.id, node.property, m_t0, m_dt_integration); + } + } + + /** + * @brief Adds an edge to the graph using the appropriate overload based on the edge's property. + * + * This function checks whether the edge property type has both 'travel_time' and 'path' members. + * If so, it uses the specialized version that uses additional parameters. + * Otherwise, we just use the edge property. + * + * @tparam GraphType The type of the graph. + * @tparam EdgeType The type of the edge. + * @param graph The graph to which the edge will be added. + * @param edge The edge to add. + */ + template + void add_edge(GraphType& graph, const EdgeType& edge) + { + // Deduce the property type of the edge. + using PropertyType = std::decay_t; + + // Check if the property type has both 'travel_time' and 'path' members. + if constexpr (has_travel_time::value && has_path::value) { + graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property.get_parameters(), + edge.property.travel_time, edge.property.path); + } + else { + graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property); + } + } + + /** + * @brief Creates a mobility simulation based on a sampled graph. + * @tparam SampleGraphFunction The type of the function that samples the graph. + * @param sample_graph A function that takes the internal graph and returns a sampled graph. + * @return A mobility simulation created from the sampled simulation graph. + */ template - mio::GraphSimulation create_sampled_simulation(SampleGraphFunction sample_graph) + auto create_sampled_simulation(SampleGraphFunction sample_graph) { + // Create an empty simulation graph. SimulationGraph sim_graph; + // Sample the original graph using the provided sampling function. auto sampled_graph = sample_graph(m_graph); + + // Iterate over each node in the sampled graph and add it to the simulation graph. for (auto&& node : sampled_graph.nodes()) { - sim_graph.add_node(node.id, node.property, m_t0, m_dt_integration); + add_node(sim_graph, node); } + + // Iterate over each edge in the sampled graph and add it to the simulation graph. for (auto&& edge : sampled_graph.edges()) { - sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property); + add_edge(sim_graph, edge); } return make_mobility_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); diff --git a/cpp/memilio/compartments/simulation.h b/cpp/memilio/compartments/simulation.h index a72bee0386..10a787b6b4 100644 --- a/cpp/memilio/compartments/simulation.h +++ b/cpp/memilio/compartments/simulation.h @@ -60,6 +60,19 @@ class Simulation { } + /** + * @brief Copy constructor + * @param[in] other Another simulation to copy + */ + Simulation(const Simulation& other) + : m_integratorCore(other.m_integratorCore) + , m_model(std::make_unique(*other.m_model)) + , m_integrator(m_integratorCore) + , m_result(other.m_result) + , m_dt(other.m_dt) + { + } + /** * @brief Set the integrator core used in the simulation. * @param[in] integrator A shared pointer to an object derived from IntegratorCore. diff --git a/cpp/memilio/io/mobility_io.cpp b/cpp/memilio/io/mobility_io.cpp index 5786aa949b..1dbdc111e0 100644 --- a/cpp/memilio/io/mobility_io.cpp +++ b/cpp/memilio/io/mobility_io.cpp @@ -158,6 +158,101 @@ IOResult read_mobility_plain(const std::string& filename) return success(mobility); } +IOResult read_duration_stay(const std::string& filename) +{ + // count the number of lines in the file + BOOST_OUTCOME_TRY(auto&& num_lines, count_lines(filename)); + + // if the file is empty, return an empty matrix + if (num_lines == 0) { + return success(Eigen::MatrixXd(0, 0)); + } + + // open the file + std::fstream file; + file.open(filename, std::ios::in); + + // if the file could not be opened, return an error + if (!file.is_open()) { + return failure(StatusCode::FileNotFound, filename); + } + + Eigen::VectorXd duration(num_lines); + + try { + std::string tp; + int linenumber = 0; + // read the file line by line + while (getline(file, tp)) { + auto line = split(tp, ' '); + Eigen::Index i = static_cast(linenumber); + duration(i) = std::stod(line[0]); + linenumber++; + } + } + catch (std::runtime_error& ex) { + return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what()); + } + + return success(duration); +} + +IOResult>>> read_path_mobility(const std::string& filename) +{ + BOOST_OUTCOME_TRY(auto&& num_lines, count_lines(filename)); + + if (num_lines == 0) { + std::vector>> arr(0, std::vector>(0)); + return success(arr); + } + + std::fstream file; + file.open(filename, std::ios::in); + if (!file.is_open()) { + return failure(StatusCode::FileNotFound, filename); + } + + const int num_nodes = static_cast(std::sqrt(num_lines)); + std::vector>> arr(num_nodes, std::vector>(num_nodes)); + + try { + std::string tp; + while (getline(file, tp)) { + auto line = split(tp, ' '); + auto indx_x = std::stoi(line[0]); + auto indx_y = std::stoi(line[1]); + if (indx_x != indx_y) { + auto path = std::accumulate(line.begin() + 2, line.end(), std::string("")); + + // string -> vector of integers + std::vector path_vec; + + // Remove the square brackets and \r + path = path.substr(1, path.size() - 3); + std::stringstream ss(path); + std::string token; + + // get numbers and save them in path_vec + while (std::getline(ss, token, ',')) { + path_vec.push_back(std::stoi(token)); + } + + for (int number : path_vec) { + arr[indx_x][indx_y].push_back(number); + } + } + else { + arr[indx_x][indx_y].push_back(static_cast(indx_x)); + } + } + } + catch (std::runtime_error& ex) { + return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what()); + } + + return success(arr); +} + #ifdef MEMILIO_HAS_HDF5 IOResult save_edges(const std::vector>>& ensemble_edges, const std::vector>& pairs_edges, const fs::path& result_dir, diff --git a/cpp/memilio/io/mobility_io.h b/cpp/memilio/io/mobility_io.h index 2caf2b48e4..61a45dbace 100644 --- a/cpp/memilio/io/mobility_io.h +++ b/cpp/memilio/io/mobility_io.h @@ -17,8 +17,8 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#ifndef READ_TWITTER_H -#define READ_TWITTER_H +#ifndef MIO_MOBILITY_IO_H +#define MIO_MOBILITY_IO_H #include "memilio/io/json_serializer.h" #include "memilio/mobility/graph.h" @@ -58,6 +58,22 @@ IOResult read_mobility_formatted(const std::string& filename); */ IOResult read_mobility_plain(const std::string& filename); +/** + * @brief Reads txt file containing the duration of stay in each county. + Writes it into a Eigen vector of size N, where N is the number of local entites. + * @param filename name of file to be read + * @return IOResult containing the duration of stay in each local entity + */ +IOResult read_duration_stay(const std::string& filename); + +/** + * @brief For each edge we have the path from the start node to the end node. This functions reads the file and returns the path for each edge. + * + * @param filename Filename of the file containing the paths. + * @return IOResult>>> contains the paths for each edge. + */ +IOResult>>> read_path_mobility(const std::string& filename); + #ifdef MEMILIO_HAS_JSONCPP /** @@ -214,4 +230,4 @@ IOResult save_edges(const std::vector>>& en } // namespace mio -#endif // READ_TWITTER_H +#endif // MIO_MOBILITY_IO_H diff --git a/cpp/memilio/math/floating_point.h b/cpp/memilio/math/floating_point.h index 5ac7b81ae8..4a08ccecc7 100644 --- a/cpp/memilio/math/floating_point.h +++ b/cpp/memilio/math/floating_point.h @@ -136,6 +136,20 @@ bool floating_point_greater_equal(T v1, T v2, T abs_tol = 0, T rel_tol = std::nu return !floating_point_less(v1, v2, abs_tol, rel_tol); } +/** + * @brief Rounds a value to the nearest nth decimal place. + * + * @param x The double value to be rounded. + * @param n The number of decimal places we want to round to. + * @return The rounded double value to n decimal digits. + */ +template +T round_nth_decimal(T x, size_t n) +{ + const T factor = std::pow(10.0, n); + return std::round(x * factor) / factor; +} + } // namespace mio #endif // MIO_MATH_FLOATING_POINT_H diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.cpp b/cpp/memilio/mobility/metapopulation_mobility_detailed.cpp new file mode 100644 index 0000000000..cb50b08115 --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.cpp @@ -0,0 +1,25 @@ +/* +* Copyright (C) 2020-2024 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_detailed.h" + +namespace mio +{ + +} // namespace mio \ No newline at end of file diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h new file mode 100644 index 0000000000..120ec159a0 --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -0,0 +1,1089 @@ +/* +* Copyright (C) 2020-2024 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 METAPOPULATION_MOBILITY_DETAILED_H +#define METAPOPULATION_MOBILITY_DETAILED_H + +#include "memilio/compartments/parameter_studies.h" +#include "memilio/epidemiology/simulation_day.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/utils/logging.h" +#include "memilio/utils/time_series.h" +#include "memilio/math/eigen.h" +#include "memilio/math/eigen_util.h" +#include "memilio/utils/metaprogramming.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/math/euler.h" +#include "memilio/epidemiology/contact_matrix.h" +#include "memilio/epidemiology/dynamic_npis.h" +#include "memilio/compartments/simulation.h" +#include "memilio/utils/date.h" +#include "memilio/mobility/graph.h" +#include "memilio/io/mobility_io.h" +#include "memilio/math/floating_point.h" + +#include "boost/filesystem.hpp" + +#include +#include +#include + +namespace mio +{ + +/** + * @brief Aggregates two simulation models and a stay duration for a node. + * + * @tparam Sim Type representing the simulation or model. + */ +template +struct ExtendedNodeProperty { + Sim base_sim; ///< The base simulation model representing local dynamics. + Sim mobility_sim; ///< The simulation model used for mobility-related processes. + double stay_duration; ///< Duration individuals stay in the node. + + /** + * @brief Constructs an ExtendedNodeProperty with the given simulations and stay duration. + * @param sim1 The base simulation model. + * @param sim2 The mobility simulation model. + * @param t The stay duration. + */ + ExtendedNodeProperty(Sim sim1, Sim sim2, double t) + : base_sim(sim1) + , mobility_sim(sim2) + , stay_duration(t) + { + } +}; + +/** + * @brief Extends the basic MobilityEdge with travel time and a detailed travel path. + * + * @tparam FP Floating-point type used (default is double). + */ +template +class ExtendedMobilityEdge : public MobilityEdge +{ +public: + double travel_time; ///< The travel time along this edge. + std::vector path; ///< A vector representing the travel path (e.g., sequence of node IDs). + + /** + * @brief Constructs an ExtendedMobilityEdge using mobility parameters. + * + * @param params The mobility parameters used to initialize the base MobilityEdge. + * @param tt The travel time. + * @param p A vector representing the travel path. + */ + ExtendedMobilityEdge(const MobilityParameters& params, double tt, std::vector p) + : MobilityEdge(params) + , travel_time(tt) + , path(p) + { + } + + /** + * @brief Constructs an ExtendedMobilityEdge using a vector of coefficients. + * + * @param coeffs A vector of mobility coefficients. + * @param tt The travel time. + * @param p A vector representing the travel path. + */ + ExtendedMobilityEdge(const Eigen::VectorXd& coeffs, double tt, std::vector p) + : MobilityEdge(coeffs) + , travel_time(tt) + , path(p) + { + } + + /** + * @brief Returns a reference to the mobile population that migrated along this edge. + * + * @return A reference to the migrated mobile population. + */ + auto& get_migrated() + { + return this->m_mobile_population; + } + + /** + * @brief Returns a reference to the return times associated with this edge. + * + * @return A reference to the vector containing the return times. + */ + auto& get_return_times() + { + return this->m_return_times; + } + + /** + * @brief Returns a const reference to the mobility parameters associated with this edge. + * + * @return A const reference to the mobility parameters. + */ + auto& get_parameters() const + { + return this->m_parameters; + } +}; + +template +using ExtendedGraph = Graph, ExtendedMobilityEdge>; + +// Default implementation when get_mobility_factors is not defined for Sim +template ::value, void*> = nullptr> +auto get_mobility_factors(const Sim& /*sim*/, double /*t*/, const Eigen::Ref& y) +{ + return Eigen::VectorXd::Ones(y.rows()); +} + +/** + * @brief Checks for negative values in an Eigen vector and corrects them. + * + * The vector is assumed to be partitioned into groups corresponding to different age groups. + * Each group has a size of (vec.size() / num_age_groups). If any element is below the specified + * tolerance, it is set to zero and its (negative) value is added to the maximum element in the same group. + * This adjustment is repeated until no element is below the tolerance or the maximum iteration count is exceeded. + * + * @tparam FP Floating-point type. + * @param vec An Eigen vector reference containing values to be checked and corrected. + * @param num_age_groups The number of age groups used to determine the group size. + * @param tolerance The threshold for negative values. Elements below this value are corrected (default is -1e-7). + * @param max_iterations Maximum number of iterations allowed for corrections (default is 100). + */ +template +void check_negative_values_vec(Eigen::Ref> vec, const size_t num_age_groups, FP tolerance = -1e-7, + const size_t max_iterations = 100) +{ + // Determine the number of compartments per age group. + const size_t num_comparts = vec.size() / num_age_groups; + size_t iteration_count = 0; + + // Loop until no negative values (below tolerance) remain in the vector. + while (vec.minCoeff() < tolerance) { + Eigen::Index min_index; + // Get the minimum value in the vector and its index. + const FP min_value = vec.minCoeff(&min_index); + + // Determine the current age group based on the index. + const auto curr_age_group = min_index / num_comparts; + // Compute the range of indices for the current age group. + const auto indx_begin = curr_age_group * num_comparts; + const auto indx_end = indx_begin + num_comparts; + + // Find the index of the maximum value within the same age group. + auto max_group_indx = indx_begin; + for (auto i = indx_begin; i < indx_end; ++i) { + if (vec(i) > vec(max_group_indx)) { + max_group_indx = i; + } + } + + // Correct the negative value by setting it to zero and + // adding its (negative) value to the maximum value of the group. + vec(min_index) = 0; + vec(max_group_indx) += min_value; + + // Check if the iteration count exceeds the allowed maximum. + if (iteration_count > max_iterations) { + log_error("Number of iterations exceeded in check_negative_values_vec."); + std::exit(1); + } + iteration_count++; + } +} + +/** + * @brief Finds the time index corresponding to a given time point in the simulation. + * + * This function searches for a time point (t) in the simulation's results. + * It compares the provided time with each time point stored in the simulation (starting from the latest) + * until it finds one that differs from t by less than 1e-10. If the time point is not found and + * create_new_tp is true, a new time point is added to the results and flows. Otherwise, an error is logged. + * + * @tparam FP Floating-point type used. + * @tparam Sim Simulation type that has access to results and flows. + * @param simulation The simulation containing the results and flows. + * @param t The time point to search for. + * @param create_new_tp Flag indicating whether to create a new time point if t is not found. + * @return Eigen::Index The index corresponding to the time point. + */ +template +Eigen::Index find_time_index(Sim& simulation, FP t, bool create_new_tp) +{ + // Get references to simulation results and flows. + auto& results = simulation.get_result(); + auto& flows = simulation.get_flows(); + + // Ensure that the number of time points in results and flows match. + if (results.get_num_time_points() != flows.get_num_time_points()) { + log_error("Number of time points in results (" + std::to_string(results.get_num_time_points()) + + ") and flows (" + std::to_string(flows.get_num_time_points()) + ") do not match in find_time_index."); + } + + // Start from the last time point index. + Eigen::Index t_indx = results.get_num_time_points() - 1; + // Iterate backwards until a time point is found that matches t within a small tolerance. + for (; t_indx >= 0; --t_indx) { + if (std::abs(results.get_time(t_indx) - t) < 1e-10) { + break; + } + } + + // If no matching time point is found and create_new_tp is false, log an error. + if (t_indx < 0 && !create_new_tp) { + log_error("Time point " + std::to_string(t) + " not found in find_time_index. Latest time point is " + + std::to_string(results.get_last_time())); + } + + // If t is greater than the last time point and create_new_tp is enabled, create a new time point. + if (t_indx < 0 && results.get_last_time() < t && create_new_tp) { + // Initialize a new result vector with zeros. + Eigen::VectorXd results_vec = Eigen::VectorXd::Zero(results.get_last_value().size()); + results.add_time_point(t, results_vec); + // For flows, copy the last value (since flows are accumulated). + flows.add_time_point(t, flows.get_last_value()); + t_indx = results.get_num_time_points() - 1; + } + + return t_indx; +} + +/** + * @brief Provides mobility-related functions for simulation models. + + * @tparam FP Floating-point type used. + * @tparam Sim Simulation type that has access to results, flows, and model. + */ +template +class MobilityFunctions +{ +public: + /** + * @brief Initializes mobility for a given edge at time t. + * + * This function finds the starting time point in the source simulation, calculates the number + * of commuters to be moved based on the simulation results and mobility coefficients, and adds + * a corresponding time point to the return times of the edge. It then moves the commuters from the + * source simulation to the target simulation. + * + * @param t The current time. + * @param edge The extended mobility edge to be initialized. + * @param sim_source The source simulation from which commuters are taken. + * @param sim_target The target simulation to which commuters are added. + */ + void init_mobility(FP t, ExtendedMobilityEdge& edge, Sim& sim_source, Sim& sim_target) + { + // Find the time index in the source simulation without creating a new time point. + const auto t_indx_start_mobility_sim_from = find_time_index(sim_source, t, false); + + // Initialize the number of commuters at the start of mobility by computing + // a weighted product of the result, mobility coefficients and mobility factors. + auto results_from = sim_source.get_result(); + edge.get_migrated().add_time_point( + t, (results_from.get_value(t_indx_start_mobility_sim_from).array() * + edge.get_parameters().get_coefficients().get_matrix_at(t).array() * + get_mobility_factors(sim_source, t, results_from.get_value(t_indx_start_mobility_sim_from)).array()) + .matrix()); + // Save the time point for the return times. + edge.get_return_times().add_time_point(t); + + // Move commuters to the target simulation. + // If the target simulation does not have the time point, create one. + const auto t_indx_start_mobility_sim_to = find_time_index(sim_target, t, true); + sim_target.get_result().get_value(t_indx_start_mobility_sim_to) += edge.get_migrated().get_last_value(); + sim_source.get_result().get_last_value() -= edge.get_migrated().get_last_value(); + } + + /** + * @brief Moves the migrated population from one simulation to another. + * + * This function transfers commuters from the source simulation to the target simulation, + * updating the respective result vectors. It also calls check_negative_values_vec to ensure + * that no negative values remain in the population vectors after moving. + * + * @param t The current time. + * @param edge The extended mobility edge holding the migrated population. + * @param sim_source The simulation from which commuters are removed. + * @param sim_target The simulation to which commuters are added. + */ + void move_migrated(FP t, ExtendedMobilityEdge& edge, Sim& sim_source, Sim& sim_target) + { + // Determine the number of age groups from the simulation's model parameters. + const size_t num_age_groups = static_cast(sim_source.get_model().parameters.get_num_groups()); + // Correct any negative values in the migrated population vector. + check_negative_values_vec(edge.get_migrated().get_last_value(), num_age_groups); + + // Find the arrival time index in the target simulation, creating a new time point if needed. + const auto t_indx_sim_to_arrival = find_time_index(sim_target, t, true); + // Remove the migrated population from the source simulation. + sim_source.get_result().get_last_value() -= edge.get_migrated().get_last_value(); + // Add the migrated population to the target simulation. + sim_target.get_result().get_value(t_indx_sim_to_arrival) += edge.get_migrated().get_last_value(); + + // Re-check both simulations for negative values after the move. + check_negative_values_vec(sim_source.get_result().get_last_value(), num_age_groups); + check_negative_values_vec(sim_target.get_result().get_value(t_indx_sim_to_arrival), num_age_groups); + } + + /** + * @brief Updates the status of commuters (migrated population) in the simulation. + * + * @param t The current time. + * @param dt The time increment for updating commuters. + * @param edge The extended mobility edge holding the migrated population. + * @param sim The simulation whose mobility model is being updated. + * @param is_mobility_model Flag indicating whether the simulation is a mobility model. + */ + void update_commuters(FP t, FP dt, ExtendedMobilityEdge& edge, Sim& sim, bool is_mobility_model) + { + // Find or create the time point for the current time in the simulation. + const auto t_indx_start_mobility_sim = find_time_index(sim, t, true); + // Initialize flows vector with zeros. + Eigen::VectorXd flows = Eigen::VectorXd::Zero(sim.get_flows().get_last_value().size()); + // Update the status of migrated commuters and update flows accordingly. + update_status_migrated(edge.get_migrated().get_last_value(), sim, + sim.get_result().get_value(t_indx_start_mobility_sim), t, dt, flows); + + // If the simulation holds a mobility model, update it by creating a new time point if necessary. + if (is_mobility_model) { + const auto t_indx_mobility_model = find_time_index(sim, t + dt, true); + if (t_indx_mobility_model != sim.get_result().get_num_time_points() - 1) { + log_error("Time point " + std::to_string(t + dt) + + " not the latest in update_commuters. Latest time point is " + + std::to_string(sim.get_result().get_last_time())); + } + sim.get_result().get_value(t_indx_mobility_model) += edge.get_migrated().get_last_value(); + sim.get_flows().get_value(t_indx_mobility_model) += flows; + } + } + + /** + * @brief Deletes all time points from the migrated and return times of an edge. + * + * @param edge The extended mobility edge whose time points are to be removed. + */ + void delete_migrated(ExtendedMobilityEdge& edge) + { + // Remove time points in reverse order to avoid index shifting issues. + for (Eigen::Index i = edge.get_return_times().get_num_time_points() - 1; i >= 0; --i) { + edge.get_migrated().remove_time_point(i); + edge.get_return_times().remove_time_point(i); + } + } +}; +/** + * @brief Manages the computation of schedules for a graph's edges and nodes. + */ +class ScheduleManager +{ +public: + /** + * @brief Holds all schedule-related vectors. + * Exlanation of the vectors: + * - "schedule_edges": For each edge, this vector contains the node index where an individual is located at each timestep for a day. + * - "mobility_schedule_edges": For each edge and each timestep, a boolean flag indicating if the individual is in a mobility model. + * - "mobility_int_schedule": For each node, a vector of timesteps at which the mobility model needs to be synchronized. + * - "local_int_schedule": For each node, a vector of timesteps at which the local model needs to be synchronized. + * - "edges_mobility": For each timestep, a list of edge indices where mobility is happening. + * - "nodes_mobility": For each timestep, a list of node indices where the local model is exchanging individuals. + * - "nodes_mobility_m": For each timestep, a list of node indices where the mobility model is exchanging individuals + * - "first_mobility": For each edge, the first timestep at which mobility occurs. + */ + struct Schedule { + std::vector> schedule_edges; + std::vector> mobility_schedule_edges; + std::vector> mobility_int_schedule; + std::vector> local_int_schedule; + std::vector> edges_mobility; + std::vector> nodes_mobility; + std::vector> nodes_mobility_m; + std::vector first_mobility; + }; + + /** + * @brief Construct a new ScheduleManager object. + * @param n_t Total number of timesteps. + * @param eps A epsilon as tolerance. + */ + ScheduleManager(size_t n_t, double eps = 1e-10) + : n_timesteps(n_t) + , epsilon(eps) + { + } + + /** + * @brief Compute the complete schedule for the given graph. + * + * This function calculates both the edge schedule and the node schedule. + * @tparam Graph Type of the graph. + * @param graph The graph object containing nodes and edges. + * @return Schedule The computed schedule. + */ + template + Schedule compute_schedule(const Graph& graph) + { + Schedule schedule; + calculate_edge_schedule(graph, schedule); + calculate_node_schedule(graph, schedule); + return schedule; + } + +private: + size_t n_timesteps; ///< Total number of timesteps. + double epsilon; ///< Tolerance + + /** + * @brief Fills a range within a vector with a specified value. + * @tparam T The type of elements in the vector. + * @param vec The vector to fill. + * @param start The starting index. + * @param end The ending index. + * @param value The value to fill with. + */ + template + void fill_vector_range(std::vector& vec, size_t start, size_t end, const T& value) + { + if (start < vec.size() && end <= vec.size() && start < end) { + std::fill(vec.begin() + start, vec.begin() + end, value); + } + } + + /** + * @brief Get indices of edges while iterating over all edges in the provided schedule and selects those edges + * for which the node index at timestep t is equal to node_Idx, and the corresponding mobility flag matches the mobility parameter. + * + * @param schedule The precomputed schedule containing the node positions and mobility flags for each edge. + * @param node_Idx The index of the node to check against the edge schedules. + * @param t The timestep at which to evaluate the edge positions and mobility states. + * @param mobility The mobility flag to match; set to true to select edges in mobility mode, false otherwise. + * @return std::vector A vector containing the indices of edges that satisfy both conditions. + */ + std::vector find_edges_at_time(const Schedule& schedule, size_t node_idx, size_t t, bool mobility) const + { + std::vector edges; + // Loop over all edges in the schedule. + for (size_t edge_idx = 0; edge_idx < schedule.schedule_edges.size(); ++edge_idx) { + // Find edges located at node_idx at timestep t with the specified mobility flag. + if (schedule.schedule_edges[edge_idx][t] == node_idx && + schedule.mobility_schedule_edges[edge_idx][t] == mobility) { + edges.push_back(edge_idx); + } + } + return edges; + } + + /** + * @brief Calculates the schedule for all edges in the graph. + * + * @tparam Graph Type of the graph. + * @param graph The graph object. + * @param schedule The schedule object to populate. + */ + template + void calculate_edge_schedule(const Graph& graph, Schedule& schedule) + { + schedule.schedule_edges.reserve(graph.edges().size()); + schedule.mobility_schedule_edges.reserve(graph.edges().size()); + + for (auto& e : graph.edges()) { + std::vector edge_schedule(n_timesteps, 0); + std::vector mobility_flag(n_timesteps, false); + + // Calculate travel time per region (minimum 0.01) + const double travel_time_per_region = + std::max(0.01, round_nth_decimal(e.property.travel_time / e.property.path.size(), 2)); + + // Calculate start time for mobility (minimum 0.01) + const double start_mobility = + std::max(0.01, round_nth_decimal(1 - 2 * travel_time_per_region * e.property.path.size() - + graph.nodes()[e.end_node_idx].property.stay_duration, + 2)); + + // Calculate arrival time at destination + const double arrival_time = start_mobility + travel_time_per_region * e.property.path.size(); + + // Convert times to indices + const size_t start_index = static_cast((start_mobility + epsilon) * 100); + const size_t arrival_index = static_cast((arrival_time + epsilon) * 100); + const size_t end_index = n_timesteps - (arrival_index - start_index); + + // Fill schedule before mobility with start node index. + fill_vector_range(edge_schedule, 0, start_index, e.start_node_idx); + + // Fill mobility flag during mobility period. + fill_vector_range(mobility_flag, start_index, arrival_index, true); + + // Fill schedule during mobility along the given path. + size_t current_index = start_index; + for (size_t county : e.property.path) { + size_t next_index = current_index + static_cast((travel_time_per_region + epsilon) * 100); + fill_vector_range(edge_schedule, current_index, next_index, county); + current_index = next_index; + } + + // Fill remaining schedule after mobility with destination node. + fill_vector_range(edge_schedule, current_index, end_index, e.property.path.back()); + + // Fill return mobility period after staying. + fill_vector_range(mobility_flag, end_index, n_timesteps, true); + + // Reverse the mobility for the return trip. + auto first_true = std::find(mobility_flag.begin(), mobility_flag.end(), true); + auto last_true = std::find(mobility_flag.rbegin(), mobility_flag.rend(), true); + if (first_true != mobility_flag.end() && last_true != mobility_flag.rend()) { + size_t first_index_found = std::distance(mobility_flag.begin(), first_true); + size_t mobility_duration = arrival_index - start_index; + std::vector reversed_path(edge_schedule.begin() + first_index_found, + edge_schedule.begin() + first_index_found + mobility_duration); + std::reverse(reversed_path.begin(), reversed_path.end()); + std::copy(reversed_path.begin(), reversed_path.end(), edge_schedule.end() - mobility_duration); + + // Store the computed schedules. + schedule.schedule_edges.push_back(std::move(edge_schedule)); + schedule.mobility_schedule_edges.push_back(std::move(mobility_flag)); + } + else { + log_error("Error in creating schedule for an edge."); + } + } + } + + /** + * @brief Calculates the node schedule based on edge schedules. + * + * @tparam Graph Type of the graph. + * @param graph The graph object. + * @param schedule The schedule object to populate. + */ + template + void calculate_node_schedule(const Graph& graph, Schedule& schedule) + { + // Compute integration schedules per node. + for (size_t node_idx = 0; node_idx < graph.nodes().size(); ++node_idx) { + // Local integration schedule. + std::vector local_schedule{0}; // Always start at t=0. + auto current_edges = find_edges_at_time(schedule, node_idx, 0, false); + for (size_t t = 1; t < n_timesteps; ++t) { + auto new_edges = find_edges_at_time(schedule, node_idx, t, false); + if (new_edges.size() != current_edges.size() || + !std::equal(current_edges.begin(), current_edges.end(), new_edges.begin())) { + local_schedule.push_back(t); + current_edges = new_edges; + } + } + if (local_schedule.back() != n_timesteps - 1) { + local_schedule.push_back(n_timesteps - 1); + } + schedule.local_int_schedule.push_back(local_schedule); + + // Mobility integration schedule. + std::vector mobility_schedule; + auto current_mob_edges = find_edges_at_time(schedule, node_idx, 0, true); + for (size_t t = 1; t < n_timesteps; ++t) { + auto new_mob_edges = find_edges_at_time(schedule, node_idx, t, true); + if (new_mob_edges.size() != current_mob_edges.size() || + !std::equal(current_mob_edges.begin(), current_mob_edges.end(), new_mob_edges.begin())) { + mobility_schedule.push_back(t); + current_mob_edges = new_mob_edges; + } + } + schedule.mobility_int_schedule.push_back(mobility_schedule); + } + + // Determine the first timestep with mobility for each edge. + schedule.first_mobility.reserve(graph.edges().size()); + for (size_t edge_idx = 0; edge_idx < graph.edges().size(); ++edge_idx) { + size_t t = 0; + for (; t < schedule.mobility_schedule_edges[edge_idx].size(); ++t) { + if (schedule.mobility_schedule_edges[edge_idx][t]) { + break; + } + } + schedule.first_mobility.push_back(t); + } + + // Initialize mobility-related vectors per timestep. + schedule.edges_mobility.reserve(n_timesteps); + schedule.nodes_mobility.reserve(n_timesteps); + schedule.nodes_mobility_m.reserve(n_timesteps); + + // At t = 0: collect edges with mobility starting at 0. + std::vector initial_edges; + for (size_t edge_idx = 0; edge_idx < graph.edges().size(); ++edge_idx) { + if (schedule.first_mobility[edge_idx] == 0) { + initial_edges.push_back(edge_idx); + } + } + schedule.edges_mobility.push_back(std::move(initial_edges)); + + // At t = 0: initialize with all nodes. + std::vector initial_nodes(graph.nodes().size()); + std::iota(initial_nodes.begin(), initial_nodes.end(), 0); + schedule.nodes_mobility.push_back(std::move(initial_nodes)); + + // At t = 0: nodes with mobility activity (if present in mobility_int_schedule). + std::vector initial_nodes_mob_m; + for (size_t node_idx = 0; node_idx < graph.nodes().size(); ++node_idx) { + if (std::binary_search(schedule.mobility_int_schedule[node_idx].begin(), + schedule.mobility_int_schedule[node_idx].end(), 0)) { + initial_nodes_mob_m.push_back(node_idx); + } + } + schedule.nodes_mobility_m.push_back(std::move(initial_nodes_mob_m)); + + // For each subsequent timestep, update mobility edge and node lists. + for (size_t t = 1; t < n_timesteps; ++t) { + // Identify mobility-active edges at timestep t. + std::vector edges_at_t; + for (size_t edge_idx = 0; edge_idx < graph.edges().size(); ++edge_idx) { + size_t current_node = schedule.schedule_edges[edge_idx][t]; + if (t >= schedule.first_mobility[edge_idx]) { + if (schedule.mobility_schedule_edges[edge_idx][t] && + std::binary_search(schedule.mobility_int_schedule[current_node].begin(), + schedule.mobility_int_schedule[current_node].end(), t)) { + edges_at_t.push_back(edge_idx); + } + else if (!schedule.mobility_schedule_edges[edge_idx][t] && + std::binary_search(schedule.local_int_schedule[current_node].begin(), + schedule.local_int_schedule[current_node].end(), t)) { + edges_at_t.push_back(edge_idx); + } + } + } + schedule.edges_mobility.push_back(edges_at_t); + + // Identify nodes with local and mobility integration at timestep t. + std::vector nodes_local; + std::vector nodes_mob_m; + for (size_t node_idx = 0; node_idx < graph.nodes().size(); ++node_idx) { + if (std::binary_search(schedule.local_int_schedule[node_idx].begin(), + schedule.local_int_schedule[node_idx].end(), t)) { + nodes_local.push_back(node_idx); + } + if (std::binary_search(schedule.mobility_int_schedule[node_idx].begin(), + schedule.mobility_int_schedule[node_idx].end(), t)) { + nodes_mob_m.push_back(node_idx); + } + } + schedule.nodes_mobility.push_back(nodes_local); + schedule.nodes_mobility_m.push_back(nodes_mob_m); + } + } +}; + +/** + * @brief Extended graph simulation class that adds mobility models. + * + * This class extends the base graph simulation by incorporating mobility-specific functions. + * It computes schedules for edge and node mobility and advances the simulation in discrete time steps. + * + * @tparam Graph Type of the graph. + * @tparam MobilityFunctions Type providing mobility-related functions. + */ +template +class GraphSimulationExtended + : public GraphSimulationBase, + std::function> +{ +public: + using node_function = std::function; + using edge_function = + std::function; + + /** + * @brief Construct a new GraphSimulationExtended object. + * + * @param t0 Initial simulation time. + * @param dt Time step size. + * @param g Reference to the graph. + * @param node_func Function to update nodes. + * @param modes Mobility functions object. + */ + GraphSimulationExtended(double t0, double dt, Graph& g, const node_function& node_func, MobilityFunctions modes) + : GraphSimulationBase, + std::function>(t0, dt, std::move(g), + node_func, {}) + , m_mobility_functions(modes) + { + ScheduleManager schedule_manager(100); // Assume 100 timesteps. + schedules = schedule_manager.compute_schedule(this->m_graph); + } + + /** + * @brief Construct a new GraphSimulationExtended object. + * + * @param t0 Initial simulation time. + * @param dt Time step size. + * @param g Graph to be moved. + * @param node_func Function to update nodes. + * @param modes Mobility functions object. + */ + GraphSimulationExtended(double t0, double dt, Graph&& g, const node_function& node_func, MobilityFunctions modes) + : GraphSimulationBase, + std::function>( + t0, dt, std::forward(g), node_func, {}) + , m_mobility_functions(modes) + { + ScheduleManager schedule_manager(100); // Assume 100 timesteps. + schedules = schedule_manager.compute_schedule(this->m_graph); + } + + /** + * @brief Advances the simulation until time t_max. + * + * The simulation proceeds in three steps per time interval: + * 1. Moving commuters between nodes. + * 2. Advancing the local node integrators. + * 3. Updating the status of commuters and mobility models. + * + * At the end of each day, populations in mobility nodes are reset and, for long simulations, the time series + * are interpolated to save memory. + * + * @param t_max Maximum simulation time. + */ + void advance(double t_max = 1.0) + { + // Compute earliest mobility start time among all edges. + ScalarType dt_first_mobility = std::accumulate( + this->m_graph.edges().begin(), this->m_graph.edges().end(), std::numeric_limits::max(), + [&](ScalarType current_min, const auto& e) { + auto travel_time_per_region = round_nth_decimal(e.property.travel_time / e.property.path.size(), 2); + travel_time_per_region = (travel_time_per_region < 0.01) ? 0.01 : travel_time_per_region; + auto start_mobility = + round_nth_decimal(1 - 2 * (travel_time_per_region * e.property.path.size()) - + this->m_graph.nodes()[e.end_node_idx].property.stay_duration, + 2); + if (start_mobility < 0.) + start_mobility = 0.; + return std::min(current_min, start_mobility); + }); + + // Reset mobility simulation results before starting. + for (auto& n : this->m_graph.nodes()) { + n.property.mobility_sim.get_result().get_last_value().setZero(); + } + + const auto min_dt = 0.01; + double t_begin = this->m_t - 1.0; + + // Main simulation loop. + while (this->m_t < t_max - epsilon) { + t_begin += 1; + if (this->m_t + dt_first_mobility > t_max) { + dt_first_mobility = t_max - this->m_t; + for (auto& n : this->m_graph.nodes()) { + n.property.base_sim.advance(this->m_t + dt_first_mobility); + } + break; + } + + size_t index_schedule = 0; + while (t_begin + 1 > this->m_t + 1e-10) { + // Step 1: Move commuters. + move_commuters(index_schedule); + // Step 2: Advance local nodes. + advance_local_nodes(index_schedule); + // Step 3: Update commuter statuses. + update_status_commuters(index_schedule); + + // Handle last timestep actions. + handle_last_time_step(index_schedule); + + index_schedule++; + this->m_t += min_dt; + } + // At day's end, reset compartments for mobility nodes. + for (auto& n : this->m_graph.nodes()) { + n.property.mobility_sim.get_result().get_last_value().setZero(); + } + // For simulations beyond t > 20, interpolate results to save memory. + if (this->m_t > 20) { + for (auto& n : this->m_graph.nodes()) { + // Local simulation interpolation. + auto& result_local = n.property.base_sim.get_result(); + auto interpolated_result = interpolate_simulation_result(result_local); + auto& flows_local = n.property.base_sim.get_flows(); + auto interpolated_flows = interpolate_simulation_result(flows_local); + n.property.base_sim.get_result() = interpolated_result; + n.property.base_sim.get_flows() = interpolated_flows; + + // Mobility simulation interpolation. + auto& result_mob = n.property.mobility_sim.get_result(); + interpolated_result = interpolate_simulation_result(result_mob); + auto& flows_mob = n.property.mobility_sim.get_flows(); + interpolated_flows = interpolate_simulation_result(flows_mob); + n.property.mobility_sim.get_result() = interpolated_result; + n.property.mobility_sim.get_flows() = interpolated_flows; + } + } + } + } + +private: + MobilityFunctions m_mobility_functions; ///< Mobility functions used to update simulation state. + ScheduleManager::Schedule schedules; ///< Precomputed schedules for mobility. + const double epsilon = 1e-10; ///< Tolerance + + /** + * @brief Calculates the next time step for an edge given its current schedule index. + * + * @param edge_index Index of the edge. + * @param index_schedule Current schedule index. + * @return ScalarType Next time step value. + * + * @throw runtime_error if the schedule index is not found. + */ + ScalarType calculate_next_dt(size_t edge_index, size_t index_schedule) + { + auto current_node_idx = schedules.schedule_edges[edge_index][index_schedule]; + bool in_mobility_node = schedules.mobility_schedule_edges[edge_index][index_schedule]; + + // Select appropriate integration schedule. + auto integrator_schedule_row = schedules.local_int_schedule[current_node_idx]; + if (in_mobility_node) + integrator_schedule_row = schedules.mobility_int_schedule[current_node_idx]; + + // Locate the current schedule index. + const size_t index_current = std::distance( + integrator_schedule_row.begin(), + std::lower_bound(integrator_schedule_row.begin(), integrator_schedule_row.end(), index_schedule)); + + if (integrator_schedule_row[index_current] != index_schedule) + throw std::runtime_error("Error in schedule."); + + ScalarType dt_mobility = 0; + if (index_schedule == 99 || index_current == integrator_schedule_row.size() - 1) { + // Use minimal timestep at last iteration. + dt_mobility = (100 - index_schedule) * 0.01; + } + else { + // Compute dt based on the next integration point. + dt_mobility = round_nth_decimal((static_cast(integrator_schedule_row[index_current + 1]) - + static_cast(integrator_schedule_row[index_current])) / + 100 + + epsilon, + 2); + } + return dt_mobility; + } + + /** + * @brief Moves commuters along edges that are adressed in the current schedule. + * + * For each edge with mobility active at the current schedule index, this function either initializes mobility + * or moves commuters to the next node. + * + * @param index_schedule Current schedule index. + */ + void move_commuters(size_t index_schedule) + { + for (const auto& edge_index : schedules.edges_mobility[index_schedule]) { + auto& e = this->m_graph.edges()[edge_index]; + // Initialize mobility when reaching the first mobility timestep. + if (index_schedule == schedules.first_mobility[edge_index]) { + auto& node_from = + this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule - 1]].property.base_sim; + auto& node_to = + this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule]].property.mobility_sim; + m_mobility_functions.init_mobility(this->m_t, e.property, node_from, node_to); + } + else if (index_schedule > schedules.first_mobility[edge_index]) { + // Move commuters if there is a change in node or mobility flag. + if ((schedules.schedule_edges[edge_index][index_schedule] != + schedules.schedule_edges[edge_index][index_schedule - 1]) || + (schedules.mobility_schedule_edges[edge_index][index_schedule] != + schedules.mobility_schedule_edges[edge_index][index_schedule - 1])) { + auto& node_from = + schedules.mobility_schedule_edges[edge_index][index_schedule - 1] + ? this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule - 1]] + .property.mobility_sim + : this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule - 1]] + .property.base_sim; + + auto& node_to = schedules.mobility_schedule_edges[edge_index][index_schedule] + ? this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule]] + .property.mobility_sim + : this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule]] + .property.base_sim; + + m_mobility_functions.move_migrated(this->m_t, e.property, node_from, node_to); + } + } + } + } + + /** + * @brief Updates the status of commuters for edges active at the current schedule index. + * + * @param index_schedule Current schedule index. + * @param max_num_contacts Maximum allowed number of contacts (default 20.0). + */ + void update_status_commuters(size_t index_schedule, const double max_num_contacts = 20.0) + { + for (const auto& edge_index : schedules.edges_mobility[index_schedule]) { + auto& e = this->m_graph.edges()[edge_index]; + auto next_dt = calculate_next_dt(edge_index, index_schedule); + auto& node_to = + schedules.mobility_schedule_edges[edge_index][index_schedule] + ? this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule]].property.mobility_sim + : this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule]].property.base_sim; + + // Get and copy current contact pattern. + auto contact_pattern_curr = node_to.get_model().get_contact_pattern(); + auto contacts_copy = contact_pattern_curr; + if (schedules.mobility_schedule_edges[edge_index][index_schedule]) { + auto& contact_matrix = contact_pattern_curr.get_cont_freq_mat(); + Eigen::MatrixXd scaled_matrix = contact_matrix[0].get_baseline().eval() / e.property.travel_time; + // Cap the contact frequency at max_num_contacts. + for (int i = 0; i < scaled_matrix.rows(); ++i) { + for (int j = 0; j < scaled_matrix.cols(); ++j) { + if (scaled_matrix(i, j) > max_num_contacts) { + scaled_matrix(i, j) = max_num_contacts; + } + } + } + contact_matrix[0].get_baseline() = scaled_matrix; + node_to.get_model().set_contact_pattern(contact_pattern_curr); + } + + m_mobility_functions.update_commuters(this->m_t, next_dt, e.property, node_to, + schedules.mobility_schedule_edges[edge_index][index_schedule]); + + // Reset the contact pattern after update. + if (schedules.mobility_schedule_edges[edge_index][index_schedule]) + node_to.get_model().set_contact_pattern(contacts_copy); + } + } + + /** + * @brief Advances the local node simulations based on their integration schedules. + * + * For each node active in the current schedule, this function computes the next time step dt and advances the simulation. + * It also checks for negative or NaN values in the simulation result. + * + * @param index_schedule Current schedule index. + */ + void advance_local_nodes(size_t index_schedule) + { + for (const auto& node_index : schedules.nodes_mobility[index_schedule]) { + auto& n = this->m_graph.nodes()[node_index]; + const size_t index_current = + std::distance(schedules.local_int_schedule[node_index].begin(), + std::lower_bound(schedules.local_int_schedule[node_index].begin(), + schedules.local_int_schedule[node_index].end(), index_schedule)); + + const size_t val_next = (index_current == schedules.local_int_schedule[node_index].size() - 1) + ? 100 + : schedules.local_int_schedule[node_index][index_current + 1]; + const ScalarType next_dt = + round_nth_decimal((static_cast(val_next) - index_schedule) / 100 + epsilon, 2); + n.property.base_sim.advance(this->m_t + next_dt); + + // Check for negative or NaN values. + if (n.property.base_sim.get_result().get_last_value().minCoeff() < -1e-7 || + std::isnan(n.property.base_sim.get_result().get_last_value().sum())) { + auto last_value_as_vec = + std::vector(n.property.base_sim.get_result().get_last_value().data(), + n.property.base_sim.get_result().get_last_value().data() + + n.property.base_sim.get_result().get_last_value().size()); + log_error("Negative value in local node " + std::to_string(node_index) + " at schedule index " + + std::to_string(index_schedule)); + } + } + } + + /** + * @brief Handles the last timestep of a simulation day. + * + * When the simulation reaches the final schedule index (assumed to be 99), this function moves + * individuals back to their local simulation models and clears commuter time series. + * + * @param index_schedule Current schedule index. + */ + void handle_last_time_step(size_t index_schedule) + { + if (index_schedule == 99) { + size_t edge_index = 0; + for (auto& e : this->m_graph.edges()) { + auto& node_from = + this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule]].property.mobility_sim; + auto& node_to = + this->m_graph.nodes()[schedules.schedule_edges[edge_index][index_schedule]].property.base_sim; + + if (schedules.schedule_edges[edge_index][index_schedule] != e.start_node_idx) + log_error("Last node is not the start node in edge " + std::to_string(edge_index) + + " at schedule index " + std::to_string(index_schedule)); + // Move individuals back to the local model. + m_mobility_functions.move_migrated(this->m_t + 0.01, e.property, node_from, node_to); + m_mobility_functions.delete_migrated(e.property); + ++edge_index; + } + } + } +}; + +template +GraphSimulationExtended, ExtendedMobilityEdge>, MobilityFunctions> +make_mobility_sim(FP t0, FP dt, Graph, ExtendedMobilityEdge>&& graph) +{ + auto migration_modes = MobilityFunctions(); + return GraphSimulationExtended, ExtendedMobilityEdge>, + MobilityFunctions>(t0, dt, std::move(graph), {}, migration_modes); +} + +/** + * @brief number of migrated people when they return according to the model. + * E.g. during the time in the other node, some people who left as susceptible will return exposed. + * Implemented for general compartmentmodel simulations, overload for your custom model if necessary + * so that it can be found with argument-dependent lookup, i.e. in the same namespace as the model. + * @param migrated number of people that migrated as input, number of people that return as output + * @param sim Simulation that is used for the migration + * @param integrator Integrator that is used for the estimation. Has to be a one-step scheme. + * @param total total population in the node that the people migrated to. + * @param t time of migration + * @param dt timestep + */ + +template ::value>> +void update_status_migrated(Eigen::Ref::Vector> migrated, Sim& sim, + Eigen::Ref::Vector> total, FP t, FP dt, + Eigen::VectorXd& flows) +{ + auto y0 = migrated.eval(); + auto y1 = migrated.setZero(); + mio::EulerIntegratorCore().step( + [&](auto&& y, auto&& t_, auto&& dydt) { + sim.get_model().get_derivatives(total, y, t_, dydt); + }, + y0, t, dt, y1); + flows = sim.get_model().get_flow_values(); + flows *= dt; +} +} // namespace mio + +#endif //METAPOPULATION_MOBILITY_DETAILED_H diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index b64874d43b..20e40a1ae3 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -401,7 +401,7 @@ class MobilityEdge template void apply_mobility(FP t, FP dt, SimulationNode& node_from, SimulationNode& node_to); -private: +protected: MobilityParameters m_parameters; TimeSeries m_mobile_population; TimeSeries m_return_times; @@ -535,8 +535,8 @@ auto get_mobility_factors(const SimulationNode& node, double t, const Eigen * detect a get_mobility_factors function for the Model type. */ template -using test_commuters_expr_t = decltype(test_commuters( - std::declval(), std::declval&>(), std::declval())); +using test_commuters_expr_t = decltype( + test_commuters(std::declval(), std::declval&>(), std::declval())); /** * Test persons when moving from their source node. diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h index 76ea3992c5..93144bb438 100755 --- a/cpp/models/ode_secir/model.h +++ b/cpp/models/ode_secir/model.h @@ -209,6 +209,38 @@ class Model : public FlowModelparameters.template get>(); + } + + /** + * @brief Returns a reference to the contact pattern. + * + * @return A reference to the contact pattern. + */ + auto& get_contact_pattern() + { + return this->parameters.template get>(); + } + + /** + * @brief Sets the contact pattern. + * + * @tparam T The type of the new contact pattern. + * @param contact_pattern The new contact pattern to be set. + */ + template + void set_contact_pattern(T contact_pattern) + { + this->parameters.template get>() = std::move(contact_pattern); + } + /** * serialize this. * @see mio::serialize diff --git a/cpp/models/ode_secirts/model.h b/cpp/models/ode_secirts/model.h index ddc935f37d..bc5c1b5407 100644 --- a/cpp/models/ode_secirts/model.h +++ b/cpp/models/ode_secirts/model.h @@ -265,7 +265,7 @@ class Model // effective contact rate by contact rate between groups i and j and damping j FP season_val = (1 + params.template get>() * sin(3.141592653589793 * - (std::fmod((params.template get() + t), 365.0) / 182.5 + 0.5))); + (std::fmod((params.template get() + t), 365.0) / 182.5 + 0.5))); FP cont_freq_eff = season_val * contact_matrix.get_matrix_at(t)(static_cast((size_t)i), static_cast((size_t)j)); // without died people @@ -633,6 +633,38 @@ class Model return smoothed_vaccinations; } + /** + * @brief Returns a const reference to the contact pattern. + * + * @return A const reference to the contact pattern. + */ + const auto& get_contact_pattern() const + { + return this->parameters.template get>(); + } + + /** + * @brief Returns a reference to the contact pattern. + * + * @return A reference to the contact pattern. + */ + auto& get_contact_pattern() + { + return this->parameters.template get>(); + } + + /** + * @brief Sets the contact pattern. + * + * @tparam T The type of the new contact pattern. + * @param contact_pattern The new contact pattern to be set. + */ + template + void set_contact_pattern(T contact_pattern) + { + this->parameters.template get>() = std::move(contact_pattern); + } + /** * serialize this. * @see mio::serialize diff --git a/cpp/models/ode_secirts/parameter_space.h b/cpp/models/ode_secirts/parameter_space.h index 0d73a09ea8..e6cb02caae 100644 --- a/cpp/models/ode_secirts/parameter_space.h +++ b/cpp/models/ode_secirts/parameter_space.h @@ -21,6 +21,7 @@ #define MIO_ODE_SECIRTS_PARAMETER_SPACE_H #include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" #include "memilio/utils/memory.h" #include "memilio/utils/logging.h" #include "memilio/utils/parameter_distributions.h" @@ -263,6 +264,77 @@ Graph, MobilityParameters> draw_sample(Graph, MobilityPa return sampled_graph; } +/** + * Draws samples for each model node in a graph. + * Some parameters are shared between nodes and only sampled once. + * @tparam FP floating point type, e.g., double + * @param graph Graph to be sampled. + * @return Graph with nodes and edges from the input graph sampled. + */ +template +ExtendedGraph> draw_sample(ExtendedGraph>& graph) +{ + ExtendedGraph> sampled_graph; + + //sample global parameters + auto& shared_params_base_model = graph.nodes()[0].property.base_sim; + auto& shared_mobility_model = graph.nodes()[0].property.mobility_sim; + draw_sample_infection(shared_params_base_model); + auto& shared_contacts = shared_params_base_model.parameters.template get>(); + shared_contacts.draw_sample_dampings(); + auto& shared_dynamic_npis = shared_params_base_model.parameters.template get>(); + shared_dynamic_npis.draw_sample(); + + for (auto& params_node : graph.nodes()) { + auto& node_model_local = params_node.property.base_sim; + + //sample local parameters + draw_sample_demographics(params_node.property.base_sim); + + //copy global parameters + //save parameters so they aren't overwritten + auto local_icu_capacity = node_model_local.parameters.template get>(); + auto local_tnt_capacity = node_model_local.parameters.template get>(); + auto local_holidays = node_model_local.parameters.template get>().get_school_holidays(); + auto local_daily_v1 = node_model_local.parameters.template get>(); + auto local_daily_v2 = node_model_local.parameters.template get>(); + auto local_daily_v3 = node_model_local.parameters.template get>(); + node_model_local.parameters = shared_params_base_model.parameters; + node_model_local.parameters.template get>() = local_icu_capacity; + node_model_local.parameters.template get>() = local_tnt_capacity; + node_model_local.parameters.template get>().get_school_holidays() = local_holidays; + node_model_local.parameters.template get>() = local_daily_v1; + node_model_local.parameters.template get>() = local_daily_v2; + node_model_local.parameters.template get>() = local_daily_v3; + + node_model_local.parameters.template get>().make_matrix(); + node_model_local.apply_constraints(); + + // do the same for the mobility model + auto& node_mobility_model = params_node.property.mobility_sim; + auto node_mobility_contacts = node_mobility_model.parameters.template get>(); + node_mobility_model.parameters = shared_mobility_model.parameters; + node_mobility_model.parameters.template get>() = node_mobility_contacts; + + // set vaccination parameters to zero + node_mobility_model.parameters.template get>().array().setConstant(0); + node_mobility_model.parameters.template get>().array().setConstant(0); + node_mobility_model.parameters.template get>().array().setConstant(0); + + node_mobility_model.apply_constraints(); + sampled_graph.add_node(params_node.id, node_model_local, node_mobility_model, + params_node.property.stay_duration); + } + + for (auto& edge : graph.edges()) { + auto edge_params = edge.property.get_parameters(); + sampled_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge_params, edge.property.travel_time, + edge.property.path); + } + + return sampled_graph; +} + } // namespace osecirts } // namespace mio diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h index 1c9b6966f2..363faf6485 100644 --- a/cpp/models/ode_secirvvs/model.h +++ b/cpp/models/ode_secirvvs/model.h @@ -511,6 +511,38 @@ class Model } } + /** + * @brief Returns a const reference to the contact pattern. + * + * @return A const reference to the contact pattern. + */ + const auto& get_contact_pattern() const + { + return this->parameters.template get>(); + } + + /** + * @brief Returns a reference to the contact pattern. + * + * @return A reference to the contact pattern. + */ + auto& get_contact_pattern() + { + return this->parameters.template get>(); + } + + /** + * @brief Sets the contact pattern. + * + * @tparam T The type of the new contact pattern. + * @param contact_pattern The new contact pattern to be set. + */ + template + void set_contact_pattern(T contact_pattern) + { + this->parameters.template get>() = std::move(contact_pattern); + } + /** * serialize this. * @see mio::serialize @@ -911,6 +943,19 @@ auto test_commuters(Simulation& sim, Eigen::Ref> mo } } +template +auto get_contact_pattern(Model& model) +{ + auto& contact_pattern = model.parameters.template get>(); + return contact_pattern; +} + +template +void set_contact_pattern(Model& model, CP contact_pattern) +{ + model.parameters.template get>() = contact_pattern; +} + } // namespace osecirvvs } // namespace mio diff --git a/cpp/models/ode_secirvvs/parameter_space.h b/cpp/models/ode_secirvvs/parameter_space.h index 6e65bce23e..c73a99e195 100644 --- a/cpp/models/ode_secirvvs/parameter_space.h +++ b/cpp/models/ode_secirvvs/parameter_space.h @@ -21,6 +21,7 @@ #define MIO_ODE_SECIRVVS_PARAMETER_SPACE_H #include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" #include "memilio/utils/logging.h" #include "ode_secirvvs/model.h" #include "ode_secirvvs/infection_state.h" @@ -223,6 +224,83 @@ Graph, MobilityParameters> draw_sample(Graph, MobilityPa return sampled_graph; } + +/** + * Draws samples for each model node in a graph. + * Some parameters are shared between nodes and only sampled once. + * @tparam FP floating point type, e.g., double + * @param graph Graph to be sampled. + * @param fact_mask_transport Factor to adjust the transmission probability on contact for the mobility model. + * @return Graph with nodes and edges from the input graph sampled. + */ +template +ExtendedGraph> draw_sample(ExtendedGraph>& graph, FP fact_mask_transport) +{ + ExtendedGraph> sampled_graph; + + //sample global parameters + auto& shared_params_base_model = graph.nodes()[0].property.base_sim; + auto& shared_mobility_model = graph.nodes()[0].property.mobility_sim; + draw_sample_infection(shared_params_base_model); + auto& shared_contacts = shared_params_base_model.parameters.template get>(); + shared_contacts.draw_sample_dampings(); + auto& shared_dynamic_npis = shared_params_base_model.parameters.template get>(); + shared_dynamic_npis.draw_sample(); + + for (auto& params_node : graph.nodes()) { + auto& node_model_local = params_node.property.base_sim; + + //sample local parameters + draw_sample_demographics(params_node.property.base_sim); + + //copy global parameters + //save demographic parameters so they aren't overwritten + auto local_icu_capacity = node_model_local.parameters.template get>(); + auto local_tnt_capacity = node_model_local.parameters.template get>(); + auto local_holidays = node_model_local.parameters.template get>().get_school_holidays(); + auto local_daily_v1 = node_model_local.parameters.template get>(); + auto local_daily_v2 = node_model_local.parameters.template get>(); + node_model_local.parameters = shared_params_base_model.parameters; + node_model_local.parameters.template get>() = local_icu_capacity; + node_model_local.parameters.template get>() = local_tnt_capacity; + node_model_local.parameters.template get>().get_school_holidays() = local_holidays; + node_model_local.parameters.template get>() = local_daily_v1; + node_model_local.parameters.template get>() = local_daily_v2; + + node_model_local.parameters.template get>().make_matrix(); + node_model_local.apply_constraints(); + + // do the same for the mobility model. It should have the same parametrization as the base model but + // without vaccinations and the contact patterns are different. + auto& node_mobility_model = params_node.property.mobility_sim; + auto node_mobility_contacts = node_mobility_model.parameters.template get>(); + node_mobility_model.parameters = shared_mobility_model.parameters; + node_mobility_model.parameters.template get>() = node_mobility_contacts; + + // set vaccination parameters to zero + node_mobility_model.parameters.template get>().array().setConstant(0); + node_mobility_model.parameters.template get>().array().setConstant(0); + + // adjust the transmission factor for the mobility model based on the usage of masks + for (auto age = AgeGroup(0); age < node_mobility_model.parameters.get_num_groups(); ++age) { + node_mobility_model.parameters.template get>()[age] *= + fact_mask_transport; + } + node_mobility_model.apply_constraints(); + + sampled_graph.add_node(params_node.id, node_model_local, node_mobility_model, + params_node.property.stay_duration); + } + + for (auto& edge : graph.edges()) { + auto edge_params = edge.property.get_parameters(); + sampled_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge_params, edge.property.travel_time, + edge.property.path); + } + + return sampled_graph; +} + } // namespace osecirvvs } // namespace mio diff --git a/cpp/models/ode_seir/model.h b/cpp/models/ode_seir/model.h index d771bae94c..e1af7124a3 100644 --- a/cpp/models/ode_seir/model.h +++ b/cpp/models/ode_seir/model.h @@ -206,6 +206,37 @@ class Model return mio::success(static_cast(result)); } + /** + * @brief Returns a const reference to the contact pattern. + * + * @return A const reference to the contact pattern. + */ + const auto& get_contact_pattern() const + { + return this->parameters.template get>(); + } + + /** + * @brief Returns a reference to the contact pattern. + * + * @return A reference to the contact pattern. + */ + auto& get_contact_pattern() + { + return this->parameters.template get>(); + } + + /** + * @brief Sets the contact pattern. + * + * @tparam T The type of the new contact pattern. + * @param contact_pattern The new contact pattern to be set. + */ + template + void set_contact_pattern(T contact_pattern) + { + this->parameters.template get>() = std::move(contact_pattern); + } /** * serialize this. * @see mio::serialize diff --git a/cpp/models/ode_sir/model.h b/cpp/models/ode_sir/model.h index 5b972ac8a0..d15aabd729 100644 --- a/cpp/models/ode_sir/model.h +++ b/cpp/models/ode_sir/model.h @@ -92,6 +92,38 @@ class Model } } + /** + * @brief Returns a const reference to the contact pattern. + * + * @return A const reference to the contact pattern. + */ + const auto& get_contact_pattern() const + { + return this->parameters.template get>(); + } + + /** + * @brief Returns a reference to the contact pattern. + * + * @return A reference to the contact pattern. + */ + auto& get_contact_pattern() + { + return this->parameters.template get>(); + } + + /** + * @brief Sets the contact pattern. + * + * @tparam T The type of the new contact pattern. + * @param contact_pattern The new contact pattern to be set. + */ + template + void set_contact_pattern(T contact_pattern) + { + this->parameters.template get>() = std::move(contact_pattern); + } + /** * serialize this. * @see mio::serialize diff --git a/cpp/simulations/2022_omicron_late_phase_mobility.cpp b/cpp/simulations/2022_omicron_late_phase_mobility.cpp new file mode 100644 index 0000000000..bfb348fb73 --- /dev/null +++ b/cpp/simulations/2022_omicron_late_phase_mobility.cpp @@ -0,0 +1,953 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "memilio/compartments/parameter_studies.h" +#include "memilio/geography/regions.h" +#include "memilio/io/epi_data.h" +#include "memilio/io/result_io.h" +#include "memilio/io/mobility_io.h" +#include "memilio/utils/date.h" +#include "memilio/utils/miompi.h" +#include "memilio/utils/random_number_generator.h" +#include "ode_secirts/parameters.h" +#include "ode_secirts/parameters_io.h" +#include "ode_secirts/parameter_space.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" +#include "memilio/utils/stl_util.h" +#include "boost/filesystem.hpp" +#include +#include + +namespace fs = boost::filesystem; + +/** + * Assigns a uniform distribution to an UncertainValue with a specified range. + * The value is set to the average of min and max, and the distribution is UNIFORM(min, max). + * @param[in,out] param UncertainValue to configure. + * @param[in] min Lower bound of the uniform distribution. + * @param[in] max Upper bound of the uniform distribution. + */ +void assign_uniform_distribution(mio::UncertainValue& param, double min, double max) +{ + param = mio::UncertainValue(0.5 * (min + max)); + param.set_distribution(mio::ParameterDistributionUniform(min, max)); +} + +/** + * Assigns uniform distributions to an array of UncertainValues. + * Each element i is set to the average of min[i] and max[i] with a UNIFORM(min[i], max[i]) distribution. + * @param[in,out] array Array of UncertainValues to configure. + * @param[in] min Array of lower bounds for each element. + * @param[in] max Array of upper bounds for each element. + * @tparam N Size of the array, must match the number of elements in min and max. + */ +template +void assign_uniform_distribution_array(mio::CustomIndexArray, mio::AgeGroup>& array, + const double (&min)[N], const double (&max)[N]) +{ + assert(N == array.numel() && "Array size must match the number of elements in min and max."); + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(N); ++i) { + assign_uniform_distribution(array[i], min[static_cast(i)], max[static_cast(i)]); + } +} + +/** + * Assigns a uniform distribution to all elements of an array of UncertainValues using a single range. + * Each element is set to the average of min and max with a UNIFORM(min, max) distribution. + * @param[in,out] array Array of UncertainValues to configure. + * @param[in] min Lower bound of the uniform distribution applied to all elements. + * @param[in] max Upper bound of the uniform distribution applied to all elements. + */ +void assign_uniform_distribution_array(mio::CustomIndexArray, mio::AgeGroup>& array, + double min, double max) +{ + for (auto i = mio::AgeGroup(0); i < array.size(); ++i) { + assign_uniform_distribution(array[i], min, max); + } +} + +/** + * Configures epidemiological parameters for COVID-19 model (Omicron variant) based on literature. + * @param[in,out] params Object that the parameters will be added to. + * @return IOResult indicating success (currently no failure cases defined). + */ +mio::IOResult set_covid_parameters(mio::osecirts::Parameters& params) +{ + constexpr size_t num_age_groups = 6; + + // --- Transition Times --- + // Incubation and infectious periods sourced from literature + // Sources: doi.org/10.1016/j.lanepe.2022.100446, doi.org/10.3201/eid2806.220158 + const double time_exposed_min = 1.66; + const double time_exposed_max = 1.66; + const double time_infected_no_symptoms_min = 1.44; + const double time_infected_no_symptoms_max = 1.44; + + // Symptomatic period: doi.org/10.1016/S0140-6736(22)00327-0 + const double time_infected_symptoms_min = 6.58; + const double time_infected_symptoms_max = 7.16; + + // Severe and critical periods: doi.org/10.1186/s12879-022-07971-6 + const double time_infected_severe_min[num_age_groups] = {1.8, 1.8, 1.8, 2.5, 3.5, 4.91}; + const double time_infected_severe_max[num_age_groups] = {2.3, 2.3, 2.3, 3.67, 5, 7.01}; + const double time_infected_critical_min[num_age_groups] = {9.29, 9.29, 9.29, 10.842, 11.15, 11.07}; + const double time_infected_critical_max[num_age_groups] = {10.57, 10.57, 10.57, 12.86, 13.23, 13.25}; + + assign_uniform_distribution_array(params.get>(), time_exposed_min, + time_exposed_max); + assign_uniform_distribution_array(params.get>(), + time_infected_no_symptoms_min, time_infected_no_symptoms_max); + assign_uniform_distribution_array(params.get>(), + time_infected_symptoms_min, time_infected_symptoms_max); + assign_uniform_distribution_array(params.get>(), time_infected_severe_min, + time_infected_severe_max); + assign_uniform_distribution_array(params.get>(), + time_infected_critical_min, time_infected_critical_max); + + // --- Transmission Probabilities --- + // Adjusted for Omicron variant: doi.org/10.7554/eLife.78933 + const double variant_factor = 1.94; + const double transmission_prob_contact_min[num_age_groups] = {0.02 * variant_factor, 0.05 * variant_factor, + 0.05 * variant_factor, 0.05 * variant_factor, + 0.08 * variant_factor, 0.10 * variant_factor}; + const double transmission_prob_contact_max[num_age_groups] = {0.04 * variant_factor, 0.07 * variant_factor, + 0.07 * variant_factor, 0.07 * variant_factor, + 0.10 * variant_factor, 0.15 * variant_factor}; + + // Relative transmission from asymptomatic cases (fixed for simplicity, could use DOI: 10.1097/INF.0000000000003791) + const double rel_transmission_no_symptoms_min = 0.5; + const double rel_transmission_no_symptoms_max = 0.5; + + // Risk of infection from symptomatic cases (depends on incidence and testing capacity) + const double risk_infection_symptomatic_min = 0.0; + const double risk_infection_symptomatic_max = 0.2; + const double max_risk_infection_symptomatic_min = 0.4; + const double max_risk_infection_symptomatic_max = 0.5; + + assign_uniform_distribution_array(params.get>(), + transmission_prob_contact_min, transmission_prob_contact_max); + assign_uniform_distribution_array(params.get>(), + rel_transmission_no_symptoms_min, rel_transmission_no_symptoms_max); + assign_uniform_distribution_array(params.get>(), + risk_infection_symptomatic_min, risk_infection_symptomatic_max); + assign_uniform_distribution_array(params.get>(), + max_risk_infection_symptomatic_min, max_risk_infection_symptomatic_max); + + // --- Disease Progression Probabilities --- + // Recovery from asymptomatic infection: doi.org/10.1101/2022.05.05.22274697 + const double recovered_per_infected_no_symptoms_min[num_age_groups] = {0.2, 0.25, 0.2, 0.2, 0.175, 0.1}; + const double recovered_per_infected_no_symptoms_max[num_age_groups] = {0.4, 0.45, 0.35, 0.3, 0.25, 0.15}; + + // Severe cases from symptomatic infection (2021 data with factors): doi.org/10.1016/S0140-6736(22)00462-7 + const double severe_per_infected_symptoms_min[num_age_groups] = {1 * 0.006, 0.8 * 0.006, 0.4 * 0.015, + 0.3 * 0.049, 0.25 * 0.15, 0.35 * 0.2}; + const double severe_per_infected_symptoms_max[num_age_groups] = {1 * 0.009, 0.8 * 0.009, 0.4 * 0.023, + 0.3 * 0.074, 0.25 * 0.18, 0.35 * 0.25}; + + // Critical cases from severe (Delta-adjusted, factor 0.52): doi.org/10.1177/14034948221108548 + const double critical_factor = 0.52; + const double critical_per_severe_min[num_age_groups] = {0.05 * critical_factor, 0.05 * critical_factor, + 0.05 * critical_factor, 0.10 * critical_factor, + 0.25 * critical_factor, 0.35 * critical_factor}; + const double critical_per_severe_max[num_age_groups] = {0.10 * critical_factor, 0.10 * critical_factor, + 0.10 * critical_factor, 0.20 * critical_factor, + 0.35 * critical_factor, 0.45 * critical_factor}; + + // Deaths from critical cases (factor 0.39): doi.org/10.1136/bmjgh-2023-012328 + const double death_factor = 0.39; + const double deaths_per_critical_min[num_age_groups] = {0.00 * death_factor, 0.00 * death_factor, + 0.10 * death_factor, 0.10 * death_factor, + 0.30 * death_factor, 0.50 * death_factor}; + const double deaths_per_critical_max[num_age_groups] = {0.10 * death_factor, 0.10 * death_factor, + 0.18 * death_factor, 0.18 * death_factor, + 0.50 * death_factor, 0.70 * death_factor}; + + assign_uniform_distribution_array(params.get>(), + recovered_per_infected_no_symptoms_min, recovered_per_infected_no_symptoms_max); + assign_uniform_distribution_array(params.get>(), + severe_per_infected_symptoms_min, severe_per_infected_symptoms_max); + assign_uniform_distribution_array(params.get>(), critical_per_severe_min, + critical_per_severe_max); + assign_uniform_distribution_array(params.get>(), deaths_per_critical_min, + deaths_per_critical_max); + + // --- Immunity Parameters --- + // Exposure reduction (no reduction assumed here) + const double reduc_exposed_partial_immunity_min = 1.0; + const double reduc_exposed_partial_immunity_max = 1.0; + const double reduc_exposed_improved_immunity_min = 1.0; + const double reduc_exposed_improved_immunity_max = 1.0; + + // Symptom reduction: doi.org/10.1056/NEJMoa2119451 + const double reduc_infected_symptoms_partial_min = 0.746; + const double reduc_infected_symptoms_partial_max = 0.961; + const double reduc_infected_symptoms_improved_min = 0.295; + const double reduc_infected_symptoms_improved_max = 0.344; + + // Severe/critical/death reduction + // Partial immunity: doi.org/10.1056/NEJMoa2119451 (week 4 report) + const double reduc_severe_critical_dead_partial_min = 0.52; + const double reduc_severe_critical_dead_partial_max = 0.82; + // Improved immunity: doi.org/10.1136/bmj-2022-071502 + const double reduc_severe_critical_dead_improved_min = 0.1; + const double reduc_severe_critical_dead_improved_max = 0.19; + + // Time reduction for mild infections: doi.org/10.1101/2021.09.24.21263978 + const double reduc_time_infected_mild = 0.5; + + assign_uniform_distribution_array(params.get>(), + reduc_exposed_partial_immunity_min, reduc_exposed_partial_immunity_max); + assign_uniform_distribution_array(params.get>(), + reduc_exposed_improved_immunity_min, reduc_exposed_improved_immunity_max); + assign_uniform_distribution_array(params.get>(), + reduc_infected_symptoms_partial_min, reduc_infected_symptoms_partial_max); + assign_uniform_distribution_array(params.get>(), + reduc_infected_symptoms_improved_min, reduc_infected_symptoms_improved_max); + assign_uniform_distribution_array( + params.get>(), + reduc_severe_critical_dead_partial_min, reduc_severe_critical_dead_partial_max); + assign_uniform_distribution_array( + params.get>(), + reduc_severe_critical_dead_improved_min, reduc_severe_critical_dead_improved_max); + assign_uniform_distribution_array(params.get>(), + reduc_time_infected_mild, reduc_time_infected_mild); + + // --- Seasonality --- + // Seasonal variation + const double seasonality_min = 0.1; + const double seasonality_max = 0.3; + assign_uniform_distribution(params.get>(), seasonality_min, seasonality_max); + + // --- Variant-Specific Parameters --- + params.get() = mio::get_day_in_year(mio::Date(2022, 6, 6)); + + // --- Waning Immunity Durations --- + // Temporary immunity periods: doi.org/10.1016/S1473-3099(22)00801-5 + const double immunity_interval_partial_min = 60; + const double immunity_interval_partial_max = 60; + const double immunity_interval_improved_min = 60; + const double immunity_interval_improved_max = 60; + + assign_uniform_distribution_array(params.get>(), + immunity_interval_partial_min, immunity_interval_partial_max); + assign_uniform_distribution_array(params.get>(), + immunity_interval_improved_min, immunity_interval_improved_max); + + // Waning immunity duration: doi.org/10.1016/S1473-3099(22)00801-5 + params.get>() = 365.0; + params.get>() = 365.0; + + return mio::success(); +} + +enum class ContactLocation +{ + Home = 0, + School, + Work, + Other, + Count, +}; + +/** + * different types of NPI, used as DampingType. + */ +enum class Intervention +{ + Home, + SchoolClosure, + HomeOffice, + GatheringBanFacilitiesClosure, + PhysicalDistanceAndMasks, + SeniorAwareness, +}; + +/** + * different level of NPI, used as DampingLevel. + */ +enum class InterventionLevel +{ + Main, + PhysicalDistanceAndMasks, + SeniorAwareness, + Holidays, +}; + +static const std::map contact_locations = {{ContactLocation::Home, "home"}, + {ContactLocation::School, "school_pf_eig"}, + {ContactLocation::Work, "work"}, + {ContactLocation::Other, "other"}}; + +/** + * Configures contact matrices for the model by reading data from files and applying scaling adjustments. + * @param[in] data_dir Directory containing contact data files. + * @param[in,out] params Parameters object to store the configured contact matrices. + * @param[in] avg_transport_time Average transport time (default: 0.0), used to scale contact frequencies. + * @param[in] share_staying_local Proportion of individuals staying local (default: 1.0), affects contact scaling. + * @return IOResult indicating success or an IO error if file reading fails. + */ +mio::IOResult set_contact_matrices(const fs::path& data_dir, mio::osecirts::Parameters& params, + double avg_transport_time = 0.0, double share_staying_local = 1.0) +{ + // Read transport contact matrix + BOOST_OUTCOME_TRY(auto&& transport_matrix, + mio::read_mobility_plain((data_dir / "contacts" / "contacts_transport.txt").string())); + + // Init contact matrices for all locations + constexpr size_t num_age_groups = 6; + mio::ContactMatrixGroup contact_matrices(contact_locations.size(), num_age_groups); + + for (const auto& [location_id, location_name] : contact_locations) { + // Read baseline contact matrix for this location + BOOST_OUTCOME_TRY( + auto&& baseline, + mio::read_mobility_plain((data_dir / "contacts" / ("baseline_" + location_name + ".txt")).string())); + + // Adjust "other" location by subtracting transport contacts + if (location_name == "other") { + baseline = (baseline - transport_matrix).cwiseAbs(); // Ensure non-negative values + } + + // Scale contacts based on travel time and local staying proportion + auto scaled_baseline = + (1.0 - share_staying_local) * baseline / (1.0 - avg_transport_time) + share_staying_local * baseline; + + // Adjust contacts: only scale for commuting age groups (2-5), others remain unchanged + Eigen::MatrixXd adjusted_baseline = Eigen::MatrixXd::Zero(num_age_groups, num_age_groups); + for (size_t i = 0; i < num_age_groups; ++i) { + for (size_t j = 0; j < num_age_groups; ++j) { + if ((i >= 2 && i <= 5) || (j >= 2 && j <= 5)) { + adjusted_baseline(i, j) = scaled_baseline(i, j); + } + else { + adjusted_baseline(i, j) = baseline(i, j); + } + } + } + + // Assign adjusted baseline and zero minimum to the contact matrix + contact_matrices[static_cast(location_id)].get_baseline() = adjusted_baseline; + contact_matrices[static_cast(location_id)].get_minimum() = + Eigen::MatrixXd::Zero(num_age_groups, num_age_groups); + } + + // Store the configured contact matrices in the parameters object + params.get>() = mio::UncertainContactMatrix(contact_matrices); + + return mio::success(); +} + +/** + * Configures a transport-specific contact matrix for the model by reading data from a file. + * @param[in] data_dir Directory containing the transport contact data file ("contacts/contacts_transport.txt"). + * @param[in,out] params Parameters object to store the configured transport contact matrix. + * @return IOResult indicating success or an IO error if file reading fails. + */ +mio::IOResult set_contact_matrices_transport(const fs::path& data_dir, mio::osecirts::Parameters& params) +{ + // Define the path to the transport contacts file + fs::path transport_file = data_dir / "contacts" / "contacts_transport.txt"; + + // number of age groups + const auto num_age_groups = static_cast(params.get_num_groups()); + + // Read the transport contact matrix from the file + BOOST_OUTCOME_TRY(auto&& transport_matrix, mio::read_mobility_plain(transport_file.string())); + + // Initialize contact matrices with a single group for transport + auto contact_matrices = mio::ContactMatrixGroup(1, num_age_groups); + + // Assign the transport matrix as the baseline and set a zero minimum + contact_matrices[0].get_baseline() = transport_matrix; + contact_matrices[0].get_minimum() = Eigen::MatrixXd::Zero(num_age_groups, num_age_groups); + + // Store the configured contact matrix in the parameters object + params.get>() = mio::UncertainContactMatrix(contact_matrices); + + return mio::success(); +} + +/** + * Scales contact matrices for all nodes in a graph based on commuting data and travel time. + * Computes average travel time and commuter share, then applies these to adjust local contact patterns. + * @param[in,out] params_graph Graph containing model parameters for each node, updated with scaled contact matrices. + * @param[in] data_dir Directory containing baseline contact data files. + * @param[in] mobility_data_dir Path to the file with commuter mobility data. + * @param[in] commuting_weights Weights for each age group to compute total population (size must match age groups). + * @return IOResult indicating success or an IO error if file reading or contact matrix setup fails. + */ +mio::IOResult scale_contacts_local(mio::ExtendedGraph>& params_graph, + const fs::path& data_dir, const std::string& mobility_data_dir, + const std::vector& commuting_weights) +{ + // Read commuter mobility data + BOOST_OUTCOME_TRY(auto&& mobility_data_commuter, mio::read_mobility_plain(mobility_data_dir)); + + // Calculate average travel time across all edges + const double avg_travel_time = std::accumulate(params_graph.edges().begin(), params_graph.edges().end(), 0.0, + [](double sum, const auto& edge) { + return sum + edge.property.travel_time; + }) / + static_cast(params_graph.edges().size()); + + // Calculate total population weighted by age group commuting factors + const double total_population = + std::accumulate(params_graph.nodes().begin(), params_graph.nodes().end(), 0.0, + [&commuting_weights](double sum, const auto& node) { + const auto& populations = node.property.base_sim.populations; + return sum + populations.get_group_total(mio::AgeGroup(0)) * commuting_weights[0] + + populations.get_group_total(mio::AgeGroup(1)) * commuting_weights[1] + + populations.get_group_total(mio::AgeGroup(2)) * commuting_weights[2] + + populations.get_group_total(mio::AgeGroup(3)) * commuting_weights[3] + + populations.get_group_total(mio::AgeGroup(4)) * commuting_weights[4] + + populations.get_group_total(mio::AgeGroup(5)) * commuting_weights[5]; + }); + + // Calculate total number of commuters from mobility data + const double total_commuters = + std::accumulate(params_graph.edges().begin(), params_graph.edges().end(), 0.0, + [&mobility_data_commuter](double sum, const auto& edge) { + return sum + mobility_data_commuter(edge.start_node_idx, edge.end_node_idx); + }); + + // Compute average commuter share relative to total population + const double avg_commuter_share = total_commuters / total_population; + + // Log results on the root MPI process + if (mio::mpi::is_root()) { + std::cout << "Average commuter share: " << avg_commuter_share << "\n"; + std::cout << "Average travel time: " << avg_travel_time << "\n"; + } + + // Scale contact matrices for each node using computed averages + for (auto& node : params_graph.nodes()) { + BOOST_OUTCOME_TRY( + set_contact_matrices(data_dir, node.property.base_sim.parameters, avg_travel_time, avg_commuter_share)); + } + + return mio::success(); +} + +/** + * @brief Configures graph nodes for counties with epidemiological data. + * Reads node IDs and population data from files and initializes models for each node. + * @param[in] params Model parameters applied to all nodes. + * @param[in] start_date Start date for reading epidemiological data. + * @param[in] end_date End date for reading epidemiological data. + * @param[in] data_dir Directory containing data files (e.g., contact matrices). + * @param[in] population_data_path Path to JSON file with population data. + * @param[in] stay_times_data_path Path to TXT file with stay times for local entities. + * @param[in] is_node_for_county True for county IDs, false for district IDs. + * @param[in,out] params_graph Graph to populate with nodes. + * @param[in] read_func Function to read input data and set model compartments. + * @param[in] node_func Function to get node IDs from population data. + * @param[in] scaling_factor_inf Factors for confirmed cases to adjust for undetected infections. + * @param[in] scaling_factor_icu Factor for ICU cases to adjust for underreporting. + * @param[in] tnt_capacity_factor Factor for test-and-trace capacity relative to population. + * @param[in] immunity_population Immunity distribution per age group across. + * @param[in] num_days Number of simulation days (default: 0), used for vaccination data. + * @param[in] export_time_series If true, exports daily time series data to the input directory. + * @param[in] rki_age_groups If true, uses RKI age group definitions. + * @param[in] masks If true, applies mask-related transmission reductions in mobility nodes. + * @param[in] ffp2 If true, uses FFP2 mask factors; otherwise, uses surgical mask factors. + * @return IOResult indicating success or an IO error if file reading fails. + */ +template +mio::IOResult +set_nodes(const mio::osecirts::Parameters& params, mio::Date start_date, mio::Date end_date, + const fs::path& data_dir, const std::string& population_data_path, const std::string& stay_times_data_path, + bool is_node_for_county, mio::ExtendedGraph>& params_graph, ReadFunction&& read_func, + NodeIdFunction&& node_func, const std::vector& scaling_factor_inf, FP scaling_factor_icu, + FP tnt_capacity_factor, const std::vector>& immunity_population, int num_days = 0, + bool export_time_series = false, bool rki_age_groups = true, bool masks = false, bool ffp2 = false) +{ + // Read stay times for mobility + BOOST_OUTCOME_TRY(auto&& duration_stay, mio::read_duration_stay(stay_times_data_path)); + + // Get node IDs (counties or districts) + BOOST_OUTCOME_TRY(auto&& node_ids, node_func(population_data_path, is_node_for_county, rki_age_groups)); + + // Initialize models for each node + std::vector> nodes(node_ids.size(), + mio::osecirts::Model(static_cast(params.get_num_groups()))); + for (auto& node : nodes) { + node.parameters = params; // Assign base parameters to each node + } + + // Populate nodes with epidemiological data + BOOST_OUTCOME_TRY(read_func(nodes, start_date, node_ids, scaling_factor_inf, scaling_factor_icu, data_dir.string(), + num_days, immunity_population, export_time_series)); + + // Configure each node + for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { + auto& node = nodes[node_idx]; + + // Set test-and-trace capacity with uncertainty + const FP tnt_capacity = node.populations.get_total() * tnt_capacity_factor; + auto& tnt_value = node.parameters.template get>(); + assign_uniform_distribution(tnt_value, 0.8 * tnt_capacity, 1.2 * tnt_capacity); + + // Configure school holidays based on node's state + const int county_id = static_cast(mio::regions::CountyId(node_ids[node_idx])); + const auto holiday_periods = + mio::regions::get_holidays(mio::regions::get_state_id(county_id), start_date, end_date); + auto& contacts = node.parameters.template get>(); + contacts.get_school_holidays().resize(holiday_periods.size()); + std::transform(holiday_periods.begin(), holiday_periods.end(), contacts.get_school_holidays().begin(), + [start_date](const auto& period) { + return std::make_pair( + mio::SimulationTime(mio::get_offset_in_days(period.first, start_date)), + mio::SimulationTime(mio::get_offset_in_days(period.second, start_date))); + }); + + // Add uncertainty to population compartments + for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + for (auto j = mio::Index::Compartments>(0); + j < mio::osecirts::Model::Compartments::Count; ++j) { + auto& compartment_value = nodes[node_idx].populations[{i, j}]; + compartment_value = + mio::UncertainValue(0.5 * (1.1 * double(compartment_value) + 0.9 * double(compartment_value))); + compartment_value.set_distribution(mio::ParameterDistributionUniform(0.9 * double(compartment_value), + 1.1 * double(compartment_value))); + } + } + + // Create mobility node with zero population + auto mobility_model = node; + mobility_model.populations.set_total(0); + + // Adjust transmission probability in mobility node if masks are enabled + if (masks) { + constexpr double surgical_mask_factor = 0.1; // Reduction factor for surgical masks + constexpr double ffp2_mask_factor = 0.001; // Reduction factor for FFP2 masks + double mask_factors[] = {1.0, + surgical_mask_factor, + ffp2 ? ffp2_mask_factor : surgical_mask_factor, + ffp2 ? ffp2_mask_factor : surgical_mask_factor, + ffp2 ? ffp2_mask_factor : surgical_mask_factor, + ffp2 ? ffp2_mask_factor : surgical_mask_factor}; + + constexpr double variant_factor = 1.94; // Omicron adjustment: doi.org/10.7554/eLife.78933 + double trans_prob_min[] = {0.02 * variant_factor, 0.05 * variant_factor, 0.05 * variant_factor, + 0.05 * variant_factor, 0.08 * variant_factor, 0.10 * variant_factor}; + double trans_prob_max[] = {0.04 * variant_factor, 0.07 * variant_factor, 0.07 * variant_factor, + 0.07 * variant_factor, 0.10 * variant_factor, 0.15 * variant_factor}; + + for (int i = 0; i < 6; ++i) { + trans_prob_min[i] *= mask_factors[i]; + trans_prob_max[i] *= mask_factors[i]; + } + assign_uniform_distribution_array( + mobility_model.parameters.template get>(), + trans_prob_min, trans_prob_max); + } + + // Disable vaccinations in mobility node + for (int day = 0; day < num_days; ++day) { + auto t = mio::SimulationDay(static_cast(day)); + for (auto age = mio::AgeGroup(0); age < params.get_num_groups(); ++age) { + mobility_model.parameters.template get>()[{age, t}] = 0; + mobility_model.parameters.template get>()[{age, t}] = 0; + mobility_model.parameters.template get>()[{age, t}] = 0; + } + } + + // Add node to graph with its mobility counterpart and stay duration + params_graph.add_node(node_ids[node_idx], node, mobility_model, + duration_stay(static_cast(node_idx))); + } + + return mio::success(); +} + +/** + * @brief Configures graph edges for commuting between counties or districts. + * Reads commuting matrices, travel times, and paths from files, then creates edges between node pairs with sufficient mobility. + * @param[in] travel_times_dir Path to TXT file containing travel times between counties. + * @param[in] mobility_data_dir Path to TXT file containing commuting matrices. + * @param[in] travel_path_dir Path to TXT file containing paths between counties. + * @param[in,out] params_graph Graph to add with edges. + * @param[in] migrating_compartments Infection states that participate in commuting. + * @param[in] contact_locations_size Number of contact locations (e.g., home, work) for mobility coefficients. + * @param[in] commuting_weights Weights for each age groups commuting contribution (default: empty, filled with 1.0). + * @param[in] threshold_edges Minimum commuter coefficient to create an edge (default: 4e-5). + * @return IOResult indicating success or an IO error if file reading fails. + */ +mio::IOResult set_edges(const std::string& travel_times_dir, const std::string& mobility_data_dir, + const std::string& travel_path_dir, + mio::ExtendedGraph>& params_graph, + const std::vector& migrating_compartments, + size_t contact_locations_size, std::vector commuting_weights = {}, + const double threshold_edges = 4e-5) +{ + // Read mobility, travel times, and path data + BOOST_OUTCOME_TRY(auto&& mobility_data_commuter, mio::read_mobility_plain(mobility_data_dir)); + BOOST_OUTCOME_TRY(auto&& travel_times, mio::read_mobility_plain(travel_times_dir)); + BOOST_OUTCOME_TRY(auto&& path_mobility, mio::read_path_mobility(travel_path_dir)); + + // Iterate over all pairs of nodes to set edges + for (size_t county_idx_i = 0; county_idx_i < params_graph.nodes().size(); ++county_idx_i) { + for (size_t county_idx_j = 0; county_idx_j < params_graph.nodes().size(); ++county_idx_j) { + auto& populations = params_graph.nodes()[county_idx_i].property.base_sim.populations; + const size_t num_age_groups = + static_cast(params_graph.nodes()[county_idx_i].property.base_sim.parameters.get_num_groups()); + + // Initialize mobility coefficients aligned with contact locations + auto mobility_coeffs = mio::MobilityCoefficientGroup(contact_locations_size, populations.numel()); + + // Default commuting weights to 1.0 for all age groups if not provided + if (commuting_weights.empty()) { + commuting_weights.resize(num_age_groups, 1.0); + } + + // Calculate working population weighted by commuting factors + double working_population = 0.0; + const auto min_commuter_age = mio::AgeGroup(2); // Youngest commuting age + const auto max_commuter_age = mio::AgeGroup(4); // Partially retired group + for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { + working_population += populations.get_group_total(age) * commuting_weights[static_cast(age)]; + } + + // Compute commuter coefficient from i to j + const double commuter_coeff_ij = mobility_data_commuter(county_idx_i, county_idx_j) / working_population; + + // Assign mobility coefficients for commuting compartments at work location + for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { + for (auto compartment : migrating_compartments) { + const size_t coeff_index = populations.get_flat_index({age, compartment}); + mobility_coeffs[static_cast(ContactLocation::Work)].get_baseline()[coeff_index] = + commuter_coeff_ij * commuting_weights[static_cast(age)]; + } + } + + // Add edge if commuter coefficient exceeds threshold + if (commuter_coeff_ij > threshold_edges) { + params_graph.add_edge(county_idx_i, county_idx_j, std::move(mobility_coeffs), + travel_times(county_idx_i, county_idx_j), + path_mobility[county_idx_i][county_idx_j]); + } + } + } + + return mio::success(); +} + +/** + * @brief Creates an input graph for a parameter study by reading data from files. + * @param[in] start_date Start date of the simulation. + * @param[in] end_date End date of the simulation. + * @param[in] num_days Number of days to simulate. + * @param[in] data_dir Directory containing input data files. + * @param[in] masks If true, applies mask-related transmission reductions. + * @param[in] ffp2 If true, uses FFP2 mask factors; otherwise, uses surgical masks (requires masks=true). + * @param[in] edges If true, configures graph edges for mobility. + * @return IOResult containing the created graph or an IO error if file reading fails. + */ +mio::IOResult>> +get_graph(mio::Date start_date, const mio::Date end_date, const int num_days, const std::string& data_dir, + const bool masks, const bool ffp2, const std::vector> immunity_population, + const bool edges) +{ + // file paths + const std::string travel_times_dir = mio::path_join(data_dir, "mobility", "travel_times_pathes.txt"); + const std::string durations_dir = mio::path_join(data_dir, "mobility", "activity_duration_work.txt"); + const std::string mobility_data = mio::path_join(data_dir, "mobility", "commuter_migration_with_locals.txt"); + const std::string travel_paths = mio::path_join(data_dir, "mobility", "wegketten_ohne_komma.txt"); + const std::string population_file = mio::path_join(data_dir, "pydata", "Germany", "county_current_population.json"); + + // Initialize global parameters + constexpr int num_age_groups = 6; + mio::osecirts::Parameters params(num_age_groups); + params.get() = mio::get_day_in_year(start_date); + + // Configure parameters and contact matrices + BOOST_OUTCOME_TRY(set_covid_parameters(params)); + BOOST_OUTCOME_TRY(set_contact_matrices(data_dir, params)); + + // Create graph + mio::ExtendedGraph> params_graph; + const std::vector scaling_factor_infected(num_age_groups, 1.7); + constexpr double scaling_factor_icu = 1.0; + constexpr double tnt_capacity_factor = 1.43 / 100000.0; + + BOOST_OUTCOME_TRY(set_nodes(params, start_date, end_date, data_dir, population_file, durations_dir, true, + params_graph, mio::osecirts::read_input_data_county>, + mio::get_node_ids, scaling_factor_infected, scaling_factor_icu, tnt_capacity_factor, + immunity_population, num_days, false, true, masks, ffp2)); + + // Set transport contact matrices for mobility models in all nodes + for (auto& node : params_graph.nodes()) { + BOOST_OUTCOME_TRY(set_contact_matrices_transport(data_dir, node.property.mobility_sim.parameters)); + } + + // Add edges if enabled + if (edges) { + const auto migrating_compartments = {mio::osecirts::InfectionState::SusceptibleNaive, + mio::osecirts::InfectionState::ExposedNaive, + mio::osecirts::InfectionState::InfectedNoSymptomsNaive, + mio::osecirts::InfectionState::InfectedSymptomsNaive, + mio::osecirts::InfectionState::SusceptibleImprovedImmunity, + mio::osecirts::InfectionState::SusceptiblePartialImmunity, + mio::osecirts::InfectionState::ExposedPartialImmunity, + mio::osecirts::InfectionState::InfectedNoSymptomsPartialImmunity, + mio::osecirts::InfectionState::InfectedSymptomsPartialImmunity, + mio::osecirts::InfectionState::ExposedImprovedImmunity, + mio::osecirts::InfectionState::InfectedNoSymptomsImprovedImmunity, + mio::osecirts::InfectionState::InfectedSymptomsImprovedImmunity}; + + const std::vector commuting_weights = {0.0, 0.0, 1.0, 1.0, 0.33, 0.0}; + BOOST_OUTCOME_TRY(set_edges(travel_times_dir, mobility_data, travel_paths, params_graph, migrating_compartments, + contact_locations.size(), commuting_weights)); + BOOST_OUTCOME_TRY(scale_contacts_local(params_graph, data_dir, mobility_data, commuting_weights)); + } + + return params_graph; +} + +/** + * @brief Runs a parameter study simulation and saves results. + * Configures a graph, runs multiple simulations, and exports time series and flow data. + * @param[in] data_dir Directory containing input data files. + * @param[in] res_dir Directory to save simulation results. + * @param[in] num_runs Number of simulation runs for the parameter study. + * @param[in] num_days Number of days to simulate. + * @param[in] masks Whether masks are enabled. + * @param[in] ffp2 Whether FFP2 is enabled. + * @param[in] edges Whether edges should be used. + * @return IOResult indicating success or an error during execution. + */ +mio::IOResult run(const std::string& data_dir, std::string res_dir, int num_runs, int num_days, bool masks, + bool ffp2, bool edges) +{ + // Uncomment to suppress logs + // mio::set_log_level(mio::LogLevel::critical); + const mio::Date start_date(2022, 8, 1); + const mio::Date end_date(2022, 11, 1); + + if (!masks && ffp2) { + return mio::failure(mio::StatusCode::InvalidValue, "FFP2 only possible with masks enabled"); + } + + const auto immunity_population = + std::vector>{{0.04, 0.61, 0.35}, {0.04, 0.61, 0.35}, {0.075, 0.62, 0.305}, + {0.08, 0.62, 0.3}, {0.035, 0.58, 0.385}, {0.01, 0.41, 0.58}}; + + BOOST_OUTCOME_TRY(auto params_graph, + get_graph(start_date, end_date, num_days, data_dir, masks, ffp2, immunity_population, edges)); + + // Construct result directory name + res_dir += "/masks_" + std::to_string(masks) + (ffp2 ? "_ffp2" : "") + (!edges ? "_no_edges" : ""); + if (mio::mpi::is_root()) { + std::cout << "Result directory: " << res_dir << "\n"; + } + + // Ensure result directory exists + if (!fs::exists(res_dir)) { + fs::create_directories(res_dir); + } + + // Extract county IDs from graph nodes + std::vector county_ids; + county_ids.reserve(params_graph.nodes().size()); + std::transform(params_graph.nodes().begin(), params_graph.nodes().end(), std::back_inserter(county_ids), + [](const auto& node) { + return node.id; + }); + + // Run parameter study + auto parameter_study = mio::ParameterStudy< + mio::osecirts::Simulation>>, + mio::ExtendedGraph>, + mio::ExtendedGraph< + mio::osecirts::Simulation>>>>( + params_graph, 0.0, num_days, 0.01, num_runs); + + if (mio::mpi::is_root()) { + std::cout << "Seeds: "; + for (auto seed : parameter_study.get_rng().get_seeds()) { + std::cout << seed << ", "; + } + std::cout << "\n"; + } + + auto ensemble = parameter_study.run( + [](auto&& graph) { + return draw_sample(graph); + }, + [&](auto&& results_graph, size_t run_idx) { + // get base simulation results + std::vector> base_results; + base_results.reserve(results_graph.nodes().size()); + std::transform(results_graph.nodes().begin(), results_graph.nodes().end(), std::back_inserter(base_results), + [](auto& node) { + return mio::interpolate_simulation_result(node.property.base_sim.get_result()); + }); + + std::vector> base_params; + base_params.reserve(results_graph.nodes().size()); + std::transform(results_graph.nodes().begin(), results_graph.nodes().end(), std::back_inserter(base_params), + [](auto& node) { + return node.property.base_sim.get_model(); + }); + + std::vector> base_flows; + base_flows.reserve(results_graph.nodes().size()); + std::transform(results_graph.nodes().begin(), results_graph.nodes().end(), std::back_inserter(base_flows), + [](auto& node) { + return mio::interpolate_simulation_result(node.property.base_sim.get_flows()); + }); + + // get mobility simulation results + std::vector> mobility_params; + mobility_params.reserve(results_graph.nodes().size()); + std::transform(results_graph.nodes().begin(), results_graph.nodes().end(), + std::back_inserter(mobility_params), [](auto& node) { + return node.property.mobility_sim.get_model(); + }); + + std::vector> mobility_flows; + mobility_flows.reserve(results_graph.nodes().size()); + std::transform(results_graph.nodes().begin(), results_graph.nodes().end(), + std::back_inserter(mobility_flows), [](auto& node) { + return mio::interpolate_simulation_result(node.property.mobility_sim.get_flows()); + }); + + if (mio::mpi::is_root()) { + std::cout << "Run " << run_idx << " complete.\n"; + } + + return std::make_tuple(base_results, base_params, base_flows, mobility_params, mobility_flows); + }); + + // Save results if ensemble is non-empty + if (!ensemble.empty()) { + std::vector>> ensemble_results, ensemble_flows, ensemble_flows_mobility; + std::vector>> ensemble_params, ensemble_params_mobility; + + ensemble_results.reserve(ensemble.size()); + ensemble_params.reserve(ensemble.size()); + ensemble_flows.reserve(ensemble.size()); + ensemble_params_mobility.reserve(ensemble.size()); + ensemble_flows_mobility.reserve(ensemble.size()); + + for (auto&& run : ensemble) { + ensemble_results.emplace_back(std::move(std::get<0>(run))); + ensemble_params.emplace_back(std::move(std::get<1>(run))); + ensemble_flows.emplace_back(std::move(std::get<2>(run))); + ensemble_params_mobility.emplace_back(std::move(std::get<3>(run))); + ensemble_flows_mobility.emplace_back(std::move(std::get<4>(run))); + } + + // Save base results + BOOST_OUTCOME_TRY(mio::save_results(ensemble_results, ensemble_params, county_ids, res_dir, false)); + + // Save base flows + const std::string flows_dir = res_dir + "/flows"; + if (mio::mpi::is_root()) { + fs::create_directories(flows_dir); + std::cout << "Saving flow results to \"" << flows_dir << "\".\n"; + } + BOOST_OUTCOME_TRY(mio::save_results(ensemble_flows, ensemble_params, county_ids, flows_dir, false)); + + // Save mobility flows + const std::string mobility_dir = res_dir + "/mobility"; + const std::string mobility_flows_dir = mobility_dir + "/flows"; + if (mio::mpi::is_root()) { + fs::create_directories(mobility_flows_dir); + std::cout << "Saving mobility flow results to \"" << mobility_flows_dir << "\".\n"; + } + BOOST_OUTCOME_TRY(mio::save_results(ensemble_flows_mobility, ensemble_params_mobility, county_ids, + mobility_flows_dir, false)); + } + + return mio::success(); +} + +int main(int argc, char** argv) +{ + mio::set_log_level(mio::LogLevel::warn); + mio::mpi::init(); + + std::string data_dir = ""; + std::string results_dir = ""; + int num_days = 10; + int num_runs = 2; + bool masks = true; + bool ffp2 = true; + bool edges = true; + + if (argc == 8) { + data_dir = argv[1]; + results_dir = argv[2]; + num_runs = std::atoi(argv[3]); + num_days = std::atoi(argv[4]); + masks = std::atoi(argv[5]) != 0; + ffp2 = std::atoi(argv[6]) != 0; + edges = std::atoi(argv[7]) != 0; + + if (mio::mpi::is_root()) { + printf("Reading data from \"%s\", saving results to \"%s\".\n", data_dir.c_str(), results_dir.c_str()); + printf("Number of runs: %d, Number of days: %d.\n", num_runs, num_days); + printf("Masks: %s, FFP2: %s, Edges: %s.\n", masks ? "true" : "false", ffp2 ? "true" : "false", + edges ? "true" : "false"); + } + } + else if (argc == 5) { + data_dir = argv[1]; + results_dir = argv[2]; + num_runs = std::atoi(argv[3]); + num_days = std::atoi(argv[4]); + + if (mio::mpi::is_root()) { + printf("Reading data from \"%s\", saving results to \"%s\".\n", data_dir.c_str(), results_dir.c_str()); + printf("Using default values - Number of runs: %d, Number of days: %d, Masks: true, FFP2: true, Edges: " + "true.\n", + num_runs, num_days); + } + } + else { + if (mio::mpi::is_root()) { + printf("Usage:\n"); + printf("2022_omicron_late_phase_mobility " + "\n"); + printf("\tRun simulation with data from , saving results to .\n"); + printf("\t: Number of simulation runs.\n"); + printf("\t: Number of days to simulate.\n"); + printf("\t: 0 or 1 to disable or enable masks.\n"); + printf("\t: 0 or 1 to disable or enable FFP2 (requires masks).\n"); + printf("\t: 0 or 1 to disable or enable edges.\n"); + } + mio::mpi::finalize(); + return 0; + } + + auto result = run(data_dir, results_dir, num_runs, num_days, masks, ffp2, edges); + if (!result) { + if (mio::mpi::is_root()) { + printf("%s\n", result.error().formatted_message().c_str()); + } + mio::mpi::finalize(); + return -1; + } + + mio::mpi::finalize(); + return 0; +} diff --git a/cpp/simulations/CMakeLists.txt b/cpp/simulations/CMakeLists.txt index a784b83ba5..3d85986539 100644 --- a/cpp/simulations/CMakeLists.txt +++ b/cpp/simulations/CMakeLists.txt @@ -7,6 +7,10 @@ if(MEMILIO_HAS_JSONCPP AND MEMILIO_HAS_HDF5) target_link_libraries(2021_vaccination_delta PRIVATE memilio ode_secirvvs Boost::filesystem ${HDF5_C_LIBRARIES}) target_compile_options(2021_vaccination_delta PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(2022_omicron_late_phase_mobility 2022_omicron_late_phase_mobility.cpp) + target_link_libraries(2022_omicron_late_phase_mobility PRIVATE memilio ode_secirts Boost::filesystem ${HDF5_C_LIBRARIES}) + target_compile_options(2022_omicron_late_phase_mobility PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(abm_simulation abm.cpp) target_link_libraries(abm_simulation PRIVATE memilio abm Boost::filesystem ${HDF5_C_LIBRARIES}) target_compile_options(abm_simulation PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 44e913f442..68fac5a015 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -19,6 +19,7 @@ set(TESTSOURCES test_sde_sirs.cpp test_sde_seirvv.cpp test_mobility.cpp + test_mobility_detailed.cpp test_date.cpp test_eigen_util.cpp test_odesecir_ageres.cpp diff --git a/cpp/tests/test_math_floating_point.cpp b/cpp/tests/test_math_floating_point.cpp index 3653b919f9..3b30da0bbf 100644 --- a/cpp/tests/test_math_floating_point.cpp +++ b/cpp/tests/test_math_floating_point.cpp @@ -169,3 +169,34 @@ TYPED_TEST(TestMathFloatingPoint, floating_point_greater_equal) }}; test_fp_compare(this->a, this->b, this->params, &mio::floating_point_greater_equal, truth_table); } + +TYPED_TEST(TestMathFloatingPoint, round_nth_decimal) +{ + using FP = double; + FP value = static_cast(1.234567); + EXPECT_EQ(mio::round_nth_decimal(value, 0), static_cast(1)); + EXPECT_EQ(mio::round_nth_decimal(value, 1), static_cast(1.2)); + EXPECT_EQ(mio::round_nth_decimal(value, 2), static_cast(1.23)); + EXPECT_EQ(mio::round_nth_decimal(value, 3), static_cast(1.235)); + EXPECT_EQ(mio::round_nth_decimal(value, 4), static_cast(1.2346)); + EXPECT_EQ(mio::round_nth_decimal(value, 5), static_cast(1.23457)); + EXPECT_EQ(mio::round_nth_decimal(value, 6), static_cast(1.234567)); + + value = static_cast(-1.234567); + EXPECT_EQ(mio::round_nth_decimal(value, 0), static_cast(-1)); + EXPECT_EQ(mio::round_nth_decimal(value, 1), static_cast(-1.2)); + EXPECT_EQ(mio::round_nth_decimal(value, 2), static_cast(-1.23)); + EXPECT_EQ(mio::round_nth_decimal(value, 3), static_cast(-1.235)); + EXPECT_EQ(mio::round_nth_decimal(value, 4), static_cast(-1.2346)); + EXPECT_EQ(mio::round_nth_decimal(value, 5), static_cast(-1.23457)); + EXPECT_EQ(mio::round_nth_decimal(value, 6), static_cast(-1.234567)); + + value = static_cast(0.999999); + EXPECT_EQ(mio::round_nth_decimal(value, 0), static_cast(1)); + EXPECT_EQ(mio::round_nth_decimal(value, 1), static_cast(1.0)); + EXPECT_EQ(mio::round_nth_decimal(value, 2), static_cast(1.0)); + EXPECT_EQ(mio::round_nth_decimal(value, 3), static_cast(1.0)); + EXPECT_EQ(mio::round_nth_decimal(value, 4), static_cast(1.0)); + EXPECT_EQ(mio::round_nth_decimal(value, 5), static_cast(1.0)); + EXPECT_EQ(mio::round_nth_decimal(value, 6), static_cast(0.999999)); +} diff --git a/cpp/tests/test_mobility_detailed.cpp b/cpp/tests/test_mobility_detailed.cpp new file mode 100644 index 0000000000..05d6693311 --- /dev/null +++ b/cpp/tests/test_mobility_detailed.cpp @@ -0,0 +1,236 @@ +/* +* Copyright (C) 2020-2024 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_detailed.h" +#include "matchers.h" + +#include + +// Mock classes for Graph, Node, and Edge +struct NodeProperty { + double stay_duration; +}; + +struct EdgeProperty { + double travel_time; + std::vector path; +}; + +struct Node { + NodeProperty property; +}; + +struct Edge { + size_t start_node_idx; + size_t end_node_idx; + EdgeProperty property; +}; + +class MockGraph +{ +public: + const std::vector& nodes() const + { + return nodes_; + } + + const std::vector& edges() const + { + return edges_; + } + + void add_node(const Node& node) + { + nodes_.push_back(node); + } + + void add_edge(const Edge& edge) + { + edges_.push_back(edge); + } + +private: + std::vector nodes_; + std::vector edges_; +}; + +class ScheduleManagerTest : public ::testing::Test +{ +protected: + ScheduleManagerTest() + : manager(100) + { + } + + mio::ScheduleManager manager; +}; + +// Test case for computing the schedule with a simple graph +TEST_F(ScheduleManagerTest, ComputeSimpleSchedule) +{ + MockGraph graph; + + Node node1 = {NodeProperty{0.4}}; + Node node2 = {NodeProperty{0.4}}; + graph.add_node(node1); + graph.add_node(node2); + + Edge edge1 = {0, 1, EdgeProperty{0.1, {0, 1}}}; + Edge edge2 = {1, 0, EdgeProperty{0.1, {1, 0}}}; + graph.add_edge(edge1); + graph.add_edge(edge2); + + auto schedule = manager.compute_schedule(graph); + + // Check that the schedule edges and mobility schedule edges are correctly filled + EXPECT_EQ(schedule.schedule_edges.size(), 2); + EXPECT_EQ(schedule.mobility_schedule_edges.size(), 2); + + // Check the contents of the schedule + for (const auto& edge_schedule : schedule.schedule_edges) { + EXPECT_EQ(edge_schedule.size(), 100); // timesteps + } + + for (const auto& mobility_schedule : schedule.mobility_schedule_edges) { + EXPECT_EQ(mobility_schedule.size(), 100); // timesteps + } +} + +// Test case for computing the schedule with no edges +TEST_F(ScheduleManagerTest, ComputeEmptySchedule) +{ + MockGraph graph; + + Node node1 = {NodeProperty{1.0}}; + Node node2 = {NodeProperty{1.0}}; + graph.add_node(node1); + graph.add_node(node2); + + auto schedule = manager.compute_schedule(graph); + + // Check that the schedule edges and mobility schedule edges are empty + EXPECT_TRUE(schedule.schedule_edges.empty()); + EXPECT_TRUE(schedule.mobility_schedule_edges.empty()); +} + +// Test case for verifying first_mobility vector +TEST_F(ScheduleManagerTest, VerifyFirstMobility) +{ + MockGraph graph; + + Node node1 = {NodeProperty{0.5}}; + Node node2 = {NodeProperty{0.5}}; + graph.add_node(node1); + graph.add_node(node2); + + Edge edge1 = {0, 1, EdgeProperty{0.2, {0, 1}}}; + graph.add_edge(edge1); + + auto schedule = manager.compute_schedule(graph); + + // First mobility should be at timestep 10 + EXPECT_EQ(schedule.first_mobility.size(), 1); + EXPECT_EQ(schedule.first_mobility[0], 10); +} + +// Test case for verifying local and mobility integration schedules +TEST_F(ScheduleManagerTest, VerifyIntegrationSchedules) +{ + MockGraph graph; + + Node node1 = {NodeProperty{0.4}}; + Node node2 = {NodeProperty{0.4}}; + graph.add_node(node1); + graph.add_node(node2); + + Edge edge1 = {0, 1, EdgeProperty{0.1, {0, 1}}}; + graph.add_edge(edge1); + + auto schedule = manager.compute_schedule(graph); + + // Check that local and mobility integration schedules are filled + EXPECT_EQ(schedule.local_int_schedule.size(), 2); + EXPECT_EQ(schedule.mobility_int_schedule.size(), 2); + + // the local_int_schedule for the first node should contain only the start, end time and the start time of the mobility. + EXPECT_EQ(schedule.local_int_schedule[0].size(), 3); + EXPECT_EQ(schedule.local_int_schedule[0][0], 0); + EXPECT_EQ(schedule.local_int_schedule[0][1], 40); + EXPECT_EQ(schedule.local_int_schedule[0][2], 99); + + // the second node should have the start and end time and additionally the arrival and departure time of the + // individual from the first node. + EXPECT_EQ(schedule.local_int_schedule[1].size(), 4); + EXPECT_EQ(schedule.local_int_schedule[1][0], 0); + EXPECT_EQ(schedule.local_int_schedule[1][1], 50); + EXPECT_EQ(schedule.local_int_schedule[1][2], 90); + EXPECT_EQ(schedule.local_int_schedule[1][3], 99); + + // since the travel time is 0.1, the mobility_int_schedule for the first node should contain the start time of the mobility, + // the time where the individual arrives at the second mobility model and the time where the individual arrives again before getting back. + // the last time step isnt considered since we do this by default at the end of the day. + EXPECT_EQ(schedule.mobility_int_schedule[0].size(), 3); + EXPECT_EQ(schedule.mobility_int_schedule[0][0], 40); + EXPECT_EQ(schedule.mobility_int_schedule[0][1], 45); + EXPECT_EQ(schedule.mobility_int_schedule[0][2], 95); + + // for the second node, all times should be shifted by 5 and add the time point where the individual start their home trip + EXPECT_EQ(schedule.mobility_int_schedule[1].size(), 4); + EXPECT_EQ(schedule.mobility_int_schedule[1][0], 45); + EXPECT_EQ(schedule.mobility_int_schedule[1][1], 50); + EXPECT_EQ(schedule.mobility_int_schedule[1][2], 90); + EXPECT_EQ(schedule.mobility_int_schedule[1][3], 95); +} + +// Test case where travel time per region lower than allowed minimum of 0.01 +TEST_F(ScheduleManagerTest, MinimalTravelTimePerRegion) +{ + MockGraph graph; + + Node node1 = {NodeProperty{0.4}}; + Node node2 = {NodeProperty{0.4}}; + graph.add_node(node1); + graph.add_node(node2); + + Edge edge1 = {0, 1, EdgeProperty{0.001, {0, 1}}}; // Minimal travel time per region + graph.add_edge(edge1); + + auto schedule = manager.compute_schedule(graph); + + // Check that the schedule edges and mobility schedule edges are correctly filled + EXPECT_EQ(schedule.schedule_edges.size(), 1); + EXPECT_EQ(schedule.mobility_schedule_edges.size(), 1); + + // Check that the travel time was correctly set to the minimum value + EXPECT_EQ(schedule.mobility_schedule_edges[0][56], true); + EXPECT_EQ(schedule.mobility_schedule_edges[0][57], true); + EXPECT_EQ(schedule.mobility_schedule_edges[0][98], true); + EXPECT_EQ(schedule.mobility_schedule_edges[0][99], true); + + // check that the number of true values is 4 + size_t count = 0; + for (const auto& mobility_schedule : schedule.mobility_schedule_edges) { + for (const auto& mobility : mobility_schedule) { + if (mobility) { + count++; + } + } + } + EXPECT_EQ(count, 4); +} \ No newline at end of file