From 49338f1ec4ca6d570b4d54665f2ae9682822bb71 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Wed, 6 Dec 2023 13:36:24 +0100 Subject: [PATCH 01/11] init --- cpp/examples/CMakeLists.txt | 3 +- cpp/memilio/compartments/parameter_studies.h | 19 +- cpp/memilio/compartments/simulation.h | 14 + cpp/memilio/io/mobility_io.cpp | 90 ++++ cpp/memilio/io/mobility_io.h | 23 +- cpp/memilio/mobility/graph.h | 110 +++++ cpp/memilio/mobility/graph_simulation.h | 454 +++++++++++++++++- .../metapopulation_mobility_detailed.cpp | 25 + .../metapopulation_mobility_detailed.h | 176 +++++++ .../metapopulation_mobility_instant.h | 25 +- 10 files changed, 922 insertions(+), 17 deletions(-) create mode 100644 cpp/memilio/mobility/metapopulation_mobility_detailed.cpp create mode 100644 cpp/memilio/mobility/metapopulation_mobility_detailed.h diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 532796883d..b15034f6c7 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -21,7 +21,8 @@ target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNI add_executable(ode_sir_example ode_sir.cpp) target_link_libraries(ode_sir_example PRIVATE memilio ode_sir) - +# add_executable(2022_sim_waning 2022_simulation_omicron_parallel.cpp) +# target_link_libraries(2022_sim_waning PRIVATE memilio ode_secirvvs) add_executable(seir_flows_example ode_seir_flows.cpp) target_link_libraries(seir_flows_example PRIVATE memilio ode_seir) diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 73ec8fa461..9147f37860 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -140,8 +140,8 @@ class ParameterStudy #else num_procs = 1; rank = 0; -#endif - +#endif + //The ParameterDistributions used for sampling parameters use thread_local_rng() //So we set our own RNG to be used. //Assume that sampling uses the thread_local_rng() and isn't multithreaded @@ -149,7 +149,8 @@ class ParameterStudy thread_local_rng() = m_rng; auto run_distribution = distribute_runs(m_num_runs, num_procs); - auto start_run_idx = std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0)); + auto start_run_idx = + std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0)); auto end_run_idx = start_run_idx + run_distribution[size_t(rank)]; std::vector> ensemble_result; @@ -167,7 +168,8 @@ class ParameterStudy //sample auto sim = create_sampled_simulation(sample_graph); - log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", (thread_local_rng().get_counter() - run_rng_counter).get()); + log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", + (thread_local_rng().get_counter() - run_rng_counter).get()); //perform run sim.advance(m_tmax); @@ -328,7 +330,8 @@ class ParameterStudy } /** @} */ - RandomNumberGenerator& get_rng() { + RandomNumberGenerator& get_rng() + { return m_rng; } @@ -341,10 +344,10 @@ class ParameterStudy auto sampled_graph = sample_graph(m_graph); for (auto&& node : sampled_graph.nodes()) { - sim_graph.add_node(node.id, node.property, m_t0, m_dt_integration); + sim_graph.add_node(node.id, node.stay_duration, node.property, node.node_pt, m_t0, m_dt_integration); } for (auto&& edge : sampled_graph.edges()) { - sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property); + sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.path, edge.property); } return make_migration_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); @@ -355,7 +358,7 @@ class ParameterStudy //evenly distribute runs //lower processes do one more run if runs are not evenly distributable auto num_runs_local = num_runs / num_procs; //integer division! - auto remainder = num_runs % num_procs; + auto remainder = num_runs % num_procs; std::vector run_distribution(num_procs); std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1); diff --git a/cpp/memilio/compartments/simulation.h b/cpp/memilio/compartments/simulation.h index ff2b4fff70..e2179de30b 100644 --- a/cpp/memilio/compartments/simulation.h +++ b/cpp/memilio/compartments/simulation.h @@ -59,6 +59,20 @@ class Simulation { } + /** + * @brief Copy constructor. + * @param[in] other The 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) + { + m_integrator.set_integrator(m_integratorCore); + } + /** * @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 5b0169fa04..1e602927e1 100644 --- a/cpp/memilio/io/mobility_io.cpp +++ b/cpp/memilio/io/mobility_io.cpp @@ -135,4 +135,94 @@ IOResult read_mobility_plain(const std::string& filename) return success(migration); } +IOResult read_duration_stay(const std::string& filename) +{ + BOOST_OUTCOME_TRY(num_lines, count_lines(filename)); + + if (num_lines == 0) { + return success(Eigen::MatrixXd(0, 0)); + } + + std::fstream file; + file.open(filename, std::ios::in); + if (!file.is_open()) { + return failure(StatusCode::FileNotFound, filename); + } + + Eigen::VectorXd duration(num_lines); + + try { + std::string tp; + int linenumber = 0; + 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(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); + } + + // std::vector> arr(std::sqrt(num_lines), std::vector(std::sqrt(num_lines))); + std::vector>> arr(std::sqrt(num_lines), + std::vector>(std::sqrt(num_lines))); + + 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); +} + } // namespace mio diff --git a/cpp/memilio/io/mobility_io.h b/cpp/memilio/io/mobility_io.h index 7059951600..44a04e8fb9 100644 --- a/cpp/memilio/io/mobility_io.h +++ b/cpp/memilio/io/mobility_io.h @@ -63,6 +63,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 /** @@ -133,7 +149,8 @@ IOResult write_graph(const Graph& graph, const * @param read_edges boolean value that decides whether the edges of the graph should also be read in. */ template -IOResult> read_graph(const std::string& directory, int ioflags = IOF_None, bool read_edges = true) +IOResult> read_graph(const std::string& directory, int ioflags = IOF_None, + bool read_edges = true) { std::string abs_path; if (!file_exists(directory, abs_path)) { @@ -160,7 +177,7 @@ IOResult> read_graph(const std::string& direct } //read edges; nodes must already be available for that) - if(read_edges){ + if (read_edges) { for (auto inode = size_t(0); inode < graph.nodes().size(); ++inode) { //list of edges auto edge_filename = path_join(abs_path, "GraphEdges_node" + std::to_string(inode) + ".json"); @@ -177,7 +194,7 @@ IOResult> read_graph(const std::string& direct if (end_node_idx >= graph.nodes().size()) { log_error("EndNodeIndex not in range of number of graph nodes."); return failure(StatusCode::OutOfRange, - edge_filename + ", EndNodeIndex not in range of number of graph nodes."); + edge_filename + ", EndNodeIndex not in range of number of graph nodes."); } BOOST_OUTCOME_TRY(parameters, deserialize_json(e["Parameters"], Tag{}, ioflags)); graph.add_edge(start_node_idx, end_node_idx, parameters); diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index a421b0ca6e..6cd0e946c9 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -29,6 +29,8 @@ #include "memilio/epidemiology/damping.h" #include "memilio/geography/regions.h" #include +#include +#include #include "boost/filesystem.hpp" @@ -69,11 +71,42 @@ struct Node { template Node(int node_id, Args&&... args) : id{node_id} + , stay_duration(0.5) , property(std::forward(args)...) + , node_pt(std::forward(args)...) { } + + template + Node(int node_id, double duration, Model property_arg, Model node_pt_arg, double m_t0, double m_dt_integration) + : id{node_id} + , stay_duration(duration) + , property(property_arg, m_t0, m_dt_integration) + , node_pt(node_pt_arg, m_t0, m_dt_integration) + { + } + + template + Node(int node_id, double duration, Args&&... args) + : id{node_id} + , stay_duration(duration) + , property(std::forward(args)...) + , node_pt(std::forward(args)...) + { + } + + Node(int node_id, double duration, NodePropertyT property_arg, NodePropertyT node_pt_arg) + : id{node_id} + , stay_duration(duration) + , property(property_arg) + , node_pt(node_pt_arg) + { + } + int id; + double stay_duration; NodePropertyT property; + NodePropertyT node_pt; }; /** @@ -81,13 +114,36 @@ struct Node { */ template struct Edge : public EdgeBase { + template Edge(size_t start, size_t end, Args&&... args) : EdgeBase{start, end} + , traveltime(0.) + , path{static_cast(start), static_cast(end)} + , property(std::forward(args)...) + { + } + + template + Edge(size_t start, size_t end, double t_travel, Args&&... args) + : EdgeBase{start, end} + , traveltime(t_travel) + , path{static_cast(start), static_cast(end)} , property(std::forward(args)...) { } + template + Edge(size_t start, size_t end, double t_travel, std::vector path_mobility, Args&&... args) + : EdgeBase{start, end} + , traveltime(t_travel) + , path(path_mobility) + , property(std::forward(args)...) + { + } + + double traveltime; + std::vector path; EdgePropertyT property; }; @@ -162,6 +218,29 @@ class Graph return m_nodes.back(); } + template + Node& add_node(int id, double duration_stay, Args&&... args) + { + m_nodes.emplace_back(id, duration_stay, std::forward(args)...); + return m_nodes.back(); + } + + // Ändere die Funktionsdeklaration von add_node + template + Node& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2) + { + m_nodes.emplace_back(id, duration_stay, model1, model2); + return m_nodes.back(); + } + + template + Node& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2, double m_t0, + double m_dt_integration) + { + m_nodes.emplace_back(id, duration_stay, model1, model2, m_t0, m_dt_integration); + return m_nodes.back(); + } + /** * @brief add an edge to the graph. property of the edge is constructed from arguments. */ @@ -178,6 +257,37 @@ class Graph }); } + /** + * @brief add an edge to the graph. property of the edge is constructed from arguments. + */ + template + Edge& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, Args&&... args) + { + assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); + return *insert_sorted_replace( + m_edges, Edge(start_node_idx, end_node_idx, traveltime, std::forward(args)...), + [](auto&& e1, auto&& e2) { + return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx + : e1.start_node_idx < e2.start_node_idx; + }); + } + + /** + * @brief add an edge to the graph. property of the edge is constructed from arguments. + */ + template + Edge& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, std::vector path, + Args&&... args) + { + assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); + return *insert_sorted_replace( + m_edges, Edge(start_node_idx, end_node_idx, traveltime, path, std::forward(args)...), + [](auto&& e1, auto&& e2) { + return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx + : e1.start_node_idx < e2.start_node_idx; + }); + } + /** * @brief range of nodes */ diff --git a/cpp/memilio/mobility/graph_simulation.h b/cpp/memilio/mobility/graph_simulation.h index 03975b8c2f..255d9d8c21 100644 --- a/cpp/memilio/mobility/graph_simulation.h +++ b/cpp/memilio/mobility/graph_simulation.h @@ -20,6 +20,7 @@ #ifndef MIO_MOBILITY_GRAPH_SIMULATION_H #define MIO_MOBILITY_GRAPH_SIMULATION_H +#include "memilio/math/euler.h" #include "memilio/mobility/graph.h" #include "memilio/utils/random_number_generator.h" @@ -181,8 +182,8 @@ class GraphSimulationStochastic //perform transition Base::m_edge_func(event_edge.property, flat_index, - Base::m_graph.nodes()[event_edge.start_node_idx].property, - Base::m_graph.nodes()[event_edge.end_node_idx].property); + Base::m_graph.nodes()[event_edge.start_node_idx].property, + Base::m_graph.nodes()[event_edge.end_node_idx].property); //calculate new cumulative rate cumulative_rate = get_cumulative_transition_rate(); @@ -247,6 +248,455 @@ class GraphSimulationStochastic RandomNumberGenerator m_rng; }; +template > +class GraphSimulationDetailed + : public GraphSimulationBase> +{ + +public: + using edge_function = edge_f; + using Base = GraphSimulationBase; + using node_function = std::function; + + GraphSimulationDetailed(double t0, double dt, const Graph& g, const node_function& node_func, + const edge_function&& edge_func) + : Base(t0, dt, g, node_func, std::move(edge_func)) + { + } + + GraphSimulationDetailed(double t0, double dt, Graph&& g, const node_function& node_func, + const edge_function&& edge_func) + : Base(t0, dt, std::move(g), node_func, std::move(edge_func)) + { + } + + inline double round_second_decimal(double x) + { + return std::round(x * 100.0) / 100.0; + } + + void advance(double t_max = 1.0) + { + ScalarType dt_first_mobility = + std::accumulate(Base::m_graph.edges().begin(), Base::m_graph.edges().end(), + std::numeric_limits::max(), [&](ScalarType current_min, const auto& e) { + auto traveltime_per_region = round_second_decimal(e.traveltime / e.path.size()); + if (traveltime_per_region < 0.01) + traveltime_per_region = 0.01; + auto start_mobility = + round_second_decimal(1 - 2 * (traveltime_per_region * e.path.size()) - + Base::m_graph.nodes()[e.end_node_idx].stay_duration); + if (start_mobility < 0) { + start_mobility = 0.; + } + return std::min(current_min, start_mobility); + }); + + // set population to zero in mobility nodes before starting + for (auto& n : Base::m_graph.nodes()) { + n.node_pt.get_result().get_last_value().setZero(); + } + + auto min_dt = 0.01; + double t_begin = Base::m_t - 1.; + + // calc schedule for each edge + precompute_schedule(); + const double epsilon = std::numeric_limits::epsilon(); + while (Base::m_t - epsilon < t_max) { + // auto start = std::chrono::system_clock::now(); + + t_begin += 1; + if (Base::m_t + dt_first_mobility > t_max) { + dt_first_mobility = t_max - Base::m_t; + for (auto& n : Base::m_graph.nodes()) { + m_node_func(Base::m_t, dt_first_mobility, n.property); + } + break; + } + + for (auto& node : Base::m_graph.nodes()) { + node.node_pt.get_simulation().set_integrator(std::make_shared()); + } + + size_t indx_schedule = 0; + while (t_begin + 1 > Base::m_t + 1e-10) { + for (const auto& edge_indx : edges_mobility[indx_schedule]) { + auto& e = Base::m_graph.edges()[edge_indx]; + // first mobility activity + if (indx_schedule == first_mobility[edge_indx]) { + auto& node_from = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + auto& node_to = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].node_pt; + m_edge_func(Base::m_t, 0.0, e.property, node_from, node_to, 0); + } + // next mobility activity + else if (indx_schedule > first_mobility[edge_indx]) { + auto current_node_indx = schedule_edges[edge_indx][indx_schedule]; + bool in_mobility_node = mobility_schedule_edges[edge_indx][indx_schedule]; + + // determine dt, which is equal to the last integration/syncronization point in the current node + auto integrator_schedule_row = ln_int_schedule[current_node_indx]; + if (in_mobility_node) + integrator_schedule_row = mb_int_schedule[current_node_indx]; + // search the index of indx_schedule in the integrator schedule + const size_t indx_current = + std::distance(integrator_schedule_row.begin(), + std::lower_bound(integrator_schedule_row.begin(), + integrator_schedule_row.end(), indx_schedule)); + + if (integrator_schedule_row[indx_current] != indx_schedule) + std::cout << "Error in schedule." + << "\n"; + + ScalarType dt_mobility; + if (indx_current == 0) { + dt_mobility = round_second_decimal(e.traveltime / e.path.size()); + if (dt_mobility < 0.01) + dt_mobility = 0.01; + } + else { + dt_mobility = + round_second_decimal((static_cast(integrator_schedule_row[indx_current]) - + static_cast(integrator_schedule_row[indx_current - 1])) / + 100 + + epsilon); + } + + // We have two cases. Either, we want to send the individuals to the next node, or we just want + // to update their state since a syncronization step is necessary in the current node. + if ((schedule_edges[edge_indx][indx_schedule] != + schedule_edges[edge_indx][indx_schedule - 1]) || + (mobility_schedule_edges[edge_indx][indx_schedule] != + mobility_schedule_edges[edge_indx][indx_schedule - 1])) { + auto& node_from = + mobility_schedule_edges[edge_indx][indx_schedule - 1] + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].node_pt + : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + + auto& node_to = + mobility_schedule_edges[edge_indx][indx_schedule] + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].node_pt + : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; + + if (indx_schedule < mobility_schedule_edges[edge_indx].size() - 1) { + m_edge_func(Base::m_t, dt_mobility, e.property, node_from, node_to, 1); + } + else + // the last time step is handled differently since we have to close the timeseries + m_edge_func(Base::m_t, dt_mobility, e.property, node_from, + Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property, + 2); + } + else { + auto& node_from = + mobility_schedule_edges[edge_indx][indx_schedule - 1] + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].node_pt + : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + + auto& node_to = + mobility_schedule_edges[edge_indx][indx_schedule] + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].node_pt + : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; + + assert(node_from.get_result().get_last_value() == node_to.get_result().get_last_value()); + m_edge_func(Base::m_t, dt_mobility, e.property, node_from, node_to, 3); + } + } + } + + // first we integrate the nodes in time. Afterwards the update on the edges is done. + // We start with the edges since the values for t0 are given. + // iterate over all local nodes and integrate them to the syncronization point + for (const auto& n_indx : nodes_mobility[indx_schedule]) { + auto& n = Base::m_graph.nodes()[n_indx]; + const size_t indx_current = + std::distance(ln_int_schedule[n_indx].begin(), + std::lower_bound(ln_int_schedule[n_indx].begin(), ln_int_schedule[n_indx].end(), + indx_schedule)); + + const size_t val_next = (indx_current == ln_int_schedule[n_indx].size() - 1) + ? 100 + : ln_int_schedule[n_indx][indx_current + 1]; + const ScalarType next_dt = + round_second_decimal((static_cast(val_next) - indx_schedule) / 100 + epsilon); + m_node_func(Base::m_t, next_dt, n.property); + } + + // iterate over all mobility nodes and integrate them to the syncronization point + for (const size_t& n_indx : nodes_mobility_m[indx_schedule]) { + auto& n = Base::m_graph.nodes()[n_indx]; + const size_t indx_current = + std::distance(mb_int_schedule[n_indx].begin(), + std::lower_bound(mb_int_schedule[n_indx].begin(), mb_int_schedule[n_indx].end(), + indx_schedule)); + const size_t val_next = (indx_current == mb_int_schedule[n_indx].size() - 1) + ? 100 + : mb_int_schedule[n_indx][indx_current + 1]; + const ScalarType next_dt = + round_second_decimal((static_cast(val_next) - indx_schedule) / 100 + epsilon); + + // get all time points from the last integration step + auto& last_time_point = + n.node_pt.get_result().get_time(n.node_pt.get_result().get_num_time_points() - 1); + // wenn last_time_point nicht innerhalb eines intervalls von 1-e10 von t liegt, dann setzte den letzten Zeitpunkt auf m_t + if (std::fabs(last_time_point - Base::m_t) > 1e-10) { + n.node_pt.get_result().get_last_time() = Base::m_t; + } + m_node_func(Base::m_t, next_dt, n.node_pt); + Eigen::Index indx_min; + while (n.node_pt.get_result().get_last_value().minCoeff(&indx_min) < -1e-7) { + Eigen::Index indx_max; + n.node_pt.get_result().get_last_value().maxCoeff(&indx_max); + n.node_pt.get_result().get_last_value()[indx_max] -= + n.node_pt.get_result().get_last_value()[indx_min]; + n.node_pt.get_result().get_last_value()[indx_min] = 0; + } + } + indx_schedule++; + Base::m_t += min_dt; + } + // set each compartment zero for all mobility nodes since we only model daily mobility + for (auto& n : Base::m_graph.nodes()) { + n.node_pt.get_result().get_last_value().setZero(); + } + // // messe die zeit, wie lange eine iteration bis zu dieser stelle dauert + // auto end = std::chrono::system_clock::now(); + // std::chrono::duration elapsed_seconds = end - start; + // std::cout << "Time (min) needed per Iteration is " << elapsed_seconds.count() / 60 << "min\n"; + } + } + +private: + // describes the schedule for each edge, i.e. which node is visited at which time step + std::vector> schedule_edges; + + // defines if people from a edge are currently in a mobility node. This is necessary to determine the correct + // position in the schedule. Otherwise we could add a specific identifier in schedule_edges. + std::vector> mobility_schedule_edges; + + // describes the syncronization steps which are necessary for the integrator in the edges and nodes. + std::vector> mb_int_schedule; + std::vector> ln_int_schedule; + + // defines the edges on which mobility takes place for each time step + std::vector> edges_mobility; + + // same for the nodes + std::vector> nodes_mobility; + std::vector> nodes_mobility_m; + + // first mobility activites per edge + std::vector first_mobility; + + void precompute_schedule() + { + // print transi + const size_t timesteps = 100; + schedule_edges.reserve(Base::m_graph.edges().size()); + mobility_schedule_edges.reserve(Base::m_graph.edges().size()); + + const double epsilon = std::numeric_limits::epsilon(); + + for (auto& e : Base::m_graph.edges()) { + // 100 since we round to second decimal + std::vector schedule(timesteps, 0.); + std::vector is_mobility_node(timesteps, false); + + double traveltime_per_region = round_second_decimal(e.traveltime / e.path.size()); + + if (traveltime_per_region < 0.01) + traveltime_per_region = 0.01; + + double start_mobility = round_second_decimal(0 + 1 - 2 * (traveltime_per_region * e.path.size()) - + Base::m_graph.nodes()[e.end_node_idx].stay_duration); + if (start_mobility < 0.0) { + start_mobility = 0.; + } + + double arrive_at = start_mobility + traveltime_per_region * e.path.size(); + + std::fill(schedule.begin(), schedule.begin() + static_cast((start_mobility + epsilon) * 100), + e.start_node_idx); + + // all values true between start mob und arrive + std::fill(is_mobility_node.begin() + static_cast((start_mobility + epsilon) * 100), + is_mobility_node.begin() + static_cast((arrive_at + epsilon) * 100), true); + int count_true = std::count(is_mobility_node.begin(), is_mobility_node.end(), true); + size_t current_index = static_cast((start_mobility + epsilon) * 100); + for (size_t county : e.path) { + std::fill(schedule.begin() + current_index, + schedule.begin() + current_index + + static_cast((traveltime_per_region + epsilon) * 100), + county); + current_index += static_cast((traveltime_per_region + epsilon) * 100); + } + + // stay in destination county + std::fill(schedule.begin() + current_index, schedule.end() - count_true, e.path[e.path.size() - 1]); + + // revert trip for return. + std::fill(is_mobility_node.end() - count_true, is_mobility_node.end(), true); + + auto first_true = std::find(is_mobility_node.begin(), is_mobility_node.end(), true); + auto last_true = std::find(is_mobility_node.rbegin(), is_mobility_node.rend(), true); + + // If all values are false, we get an error + if (first_true != is_mobility_node.end() && last_true != is_mobility_node.rend()) { + int first_index = std::distance(is_mobility_node.begin(), first_true); + std::vector path_reversed(schedule.begin() + first_index, + schedule.begin() + first_index + count_true); + std::reverse(path_reversed.begin(), path_reversed.end()); + std::copy(path_reversed.begin(), path_reversed.end(), schedule.end() - count_true); + schedule_edges.push_back(schedule); + mobility_schedule_edges.push_back(is_mobility_node); + } + else { + log_error("Error in creating schedule."); + } + } + + // iterate over nodes + size_t indx_node = 0; + for (auto& n : Base::m_graph.nodes()) { + // local node and automatical add zero, since we want to start at t=0 and start with integrating all nodes to the next dt + std::vector order_local_node = {0}; + std::vector indx_edges; + for (size_t indx_edge = 0; indx_edge < schedule_edges.size(); ++indx_edge) { + if (schedule_edges[indx_edge][0] == indx_node && !mobility_schedule_edges[indx_edge][0]) { + indx_edges.push_back(indx_edge); + } + } + + for (size_t indx_t = 1; indx_t < timesteps; ++indx_t) { + std::vector indx_edges_new; + for (size_t indx_edge = 0; indx_edge < schedule_edges.size(); ++indx_edge) { + if (schedule_edges[indx_edge][indx_t] == indx_node && !mobility_schedule_edges[indx_edge][indx_t]) { + indx_edges_new.push_back(indx_edge); + } + } + + if (indx_edges_new.size() != indx_edges.size() || + !std::equal(indx_edges.begin(), indx_edges.end(), indx_edges_new.begin())) { + order_local_node.push_back(indx_t); + indx_edges = indx_edges_new; + } + } + if (order_local_node[order_local_node.size() - 1] != timesteps - 1) + order_local_node.push_back(timesteps - 1); + ln_int_schedule.push_back(order_local_node); + + // mobility node + std::vector order_mobility_node; + std::vector indx_edges_mobility; + for (size_t indx_edge = 0; indx_edge < schedule_edges.size(); ++indx_edge) { + if (schedule_edges[indx_edge][0] == indx_node && mobility_schedule_edges[indx_edge][0]) { + indx_edges_mobility.push_back(indx_edge); + } + } + for (size_t indx_t = 1; indx_t < timesteps; ++indx_t) { + std::vector indx_edges_mobility_new; + for (size_t indx_edge = 0; indx_edge < schedule_edges.size(); ++indx_edge) { + if (schedule_edges[indx_edge][indx_t] == indx_node && mobility_schedule_edges[indx_edge][indx_t]) { + indx_edges_mobility_new.push_back(indx_edge); + } + } + + if (indx_edges_mobility_new.size() != indx_edges_mobility.size() || + !std::equal(indx_edges_mobility.begin(), indx_edges_mobility.end(), + indx_edges_mobility_new.begin())) { + order_mobility_node.push_back(indx_t); + indx_edges_mobility = indx_edges_mobility_new; + } + } + mb_int_schedule.push_back(order_mobility_node); + indx_node++; + } + + auto indx_edge = 0; + first_mobility.reserve(Base::m_graph.edges().size()); + for (auto& e : Base::m_graph.edges()) { + this->first_mobility[indx_edge] = std::distance( + mobility_schedule_edges[indx_edge].begin(), + std::find(mobility_schedule_edges[indx_edge].begin(), mobility_schedule_edges[indx_edge].end(), true)); + indx_edge++; + } + + edges_mobility.reserve(timesteps); + nodes_mobility.reserve(timesteps); + nodes_mobility_m.reserve(timesteps); + + // we handle indx_current = 0 separately since we want have added them always in the + // intregration schedule so this would lead to wrong results + // here we iterate over first mobility activities and add the edges indx when there is a start at t = 0 + std::vector temp_edges_mobility; + indx_edge = 0; + for (auto& start_time : first_mobility) { + if (start_time == 0) { + temp_edges_mobility.push_back(indx_edge); + } + indx_edge++; + } + edges_mobility.push_back(std::move(temp_edges_mobility)); + + // same for nodes + std::vector temp_nodes_mobility(Base::m_graph.nodes().size()); + std::iota(temp_nodes_mobility.begin(), temp_nodes_mobility.end(), 0); + nodes_mobility.emplace_back(std::move(temp_nodes_mobility)); + + std::vector temp_nodes_mobility_m; + size_t node_indx = 0; + for (auto& n : Base::m_graph.nodes()) { + if (std::binary_search(mb_int_schedule[node_indx].begin(), mb_int_schedule[node_indx].end(), 0)) { + temp_nodes_mobility_m.push_back(node_indx); + node_indx++; + } + } + nodes_mobility_m.push_back(temp_nodes_mobility_m); + + for (size_t indx_current = 1; indx_current < timesteps; ++indx_current) { + std::vector temp_edge_mobility; + indx_edge = 0; + for (auto& e : Base::m_graph.edges()) { + auto current_node_indx = schedule_edges[indx_edge][indx_current]; + if (indx_current >= first_mobility[indx_edge]) { + if (mobility_schedule_edges[indx_edge][indx_current] && + std::binary_search(mb_int_schedule[current_node_indx].begin(), + mb_int_schedule[current_node_indx].end(), indx_current)) + temp_edge_mobility.push_back(indx_edge); + else if (!mobility_schedule_edges[indx_edge][indx_current] && + std::binary_search(ln_int_schedule[current_node_indx].begin(), + ln_int_schedule[current_node_indx].end(), indx_current)) + temp_edge_mobility.push_back(indx_edge); + } + indx_edge++; + } + edges_mobility.push_back(temp_edge_mobility); + + // reset temp_nodes_mobility + temp_nodes_mobility.clear(); + temp_nodes_mobility_m.clear(); + node_indx = 0; + for (auto& n : Base::m_graph.nodes()) { + if (std::binary_search(ln_int_schedule[node_indx].begin(), ln_int_schedule[node_indx].end(), + indx_current)) { + temp_nodes_mobility.push_back(node_indx); + } + + if (std::binary_search(mb_int_schedule[node_indx].begin(), mb_int_schedule[node_indx].end(), + indx_current)) { + temp_nodes_mobility_m.push_back(node_indx); + } + node_indx++; + } + nodes_mobility.push_back(temp_nodes_mobility); + nodes_mobility_m.push_back(temp_nodes_mobility_m); + } + } +}; + template auto make_graph_sim(double t0, double dt, Graph&& g, NodeF&& node_func, EdgeF&& edge_func) { diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.cpp b/cpp/memilio/mobility/metapopulation_mobility_detailed.cpp new file mode 100644 index 0000000000..b5c7d1a178 --- /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 diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h new file mode 100644 index 0000000000..3f79166968 --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -0,0 +1,176 @@ +/* +* 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/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_instant.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 "boost/filesystem.hpp" + +#include + +namespace mio +{ + +// Erbe von MigrationEdge aus metapopulation_mobility_instant.h und überschreibe die apply_migration Funktion indem diese zusätzlich das argument mode bekommt +class MigrationEdgeDetailed : public MigrationEdge +{ +public: + template + void apply_migration(double t, double dt, SimulationNode& node_from, SimulationNode& node_to, int mode); +}; + +/** + * @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, const Sim& sim, IntegratorCore& integrator, + Eigen::Ref::Vector> total, double t, double dt) +{ + auto y0 = migrated.eval(); + auto y1 = migrated; + integrator.step( + [&](auto&& y, auto&& t_, auto&& dydt) { + sim.get_model().get_derivatives(total, y, t_, dydt); + }, + y0, t, dt, y1); +} + +template +using Vector = Eigen::Matrix; + +/* + * move migrated people from one node to another. + * @param migrated number of people that migrated + * @param results_from current state of node that people migrated from + * @param results_to current state of node that people migrated to +*/ +template +void move_migrated(Eigen::Ref> migrated, Eigen::Ref> results_from, + Eigen::Ref> results_to) +{ + // The performance of the low order integrator often fails for small compartments i.e. compartments are sometimes negative. + // We handel this by checking for values below 0 in m_migrated and add them to the highest value + Eigen::Index max_index_migrated; + for (size_t i = 0; i < migrated.size(); ++i) { + migrated.maxCoeff(&max_index_migrated); + if (migrated(i) < 0) { + migrated(max_index_migrated) -= migrated(i); + migrated(i) = 0; + } + if (results_from(i) < migrated(i)) { + migrated(max_index_migrated) -= (migrated(i) - results_from(i)); + migrated(i) = results_from(i); + } + results_from(i) -= migrated(i); + results_to(i) += migrated(i); + } +} + +template +void MigrationEdgeDetailed::apply_migration(double t, double dt, SimulationNode& node_from, + SimulationNode& node_to, int mode) +{ + + if (mode == 0) { + //normal daily migration + m_migrated.add_time_point( + t, (node_from.get_last_state().array() * m_parameters.get_coefficients().get_matrix_at(t).array() * + get_migration_factors(node_from, t, node_from.get_last_state()).array()) + .matrix()); + m_return_times.add_time_point(t + dt); + move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), + node_to.get_result().get_last_value()); + } + // change county of migrated + else if (mode == 1) { + // update status of migrated before moving to next county + IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); + update_status_migrated(m_migrated.get_last_value(), node_from.get_simulation(), integrator_node, + node_from.get_result().get_last_value(), t, dt); + move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), + node_to.get_result().get_last_value()); + } + // option for last time point to remove time points + else if (mode == 2) { + Eigen::Index idx = m_return_times.get_num_time_points() - 1; + IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); + update_status_migrated(m_migrated[idx], node_from.get_simulation(), integrator_node, + node_from.get_result().get_last_value(), t, dt); + + move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), + node_to.get_result().get_last_value()); + + for (Eigen::Index i = m_return_times.get_num_time_points() - 1; i >= 0; --i) { + if (m_return_times.get_time(i) <= t) { + m_migrated.remove_time_point(i); + m_return_times.remove_time_point(i); + } + } + } + // just update status of migrated + else if (mode == 3) { + Eigen::Index idx = m_return_times.get_num_time_points() - 1; + IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); + update_status_migrated(m_migrated[idx], node_from.get_simulation(), integrator_node, + node_from.get_result().get_last_value(), t, dt); + + move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), + node_from.get_result().get_last_value()); + } + else { + std::cout << "Invalid input mode. Should be 0 or 1." + << "\n"; + } +} + +/** + * edge functor for migration simulation. + * @see MigrationEdge::apply_migration + */ +template +void apply_migration(double t, double dt, MigrationEdge& migrationEdge, SimulationNode& node_from, + SimulationNode& node_to, int mode) +{ + migrationEdge.apply_migration(t, dt, node_from, node_to, mode); +} +} // 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 5178315036..9b2fcfa953 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -39,6 +39,9 @@ namespace mio { +/** + * represents the simulation in one node of the graph. + */ /** * represents the simulation in one node of the graph. */ @@ -54,6 +57,20 @@ class SimulationNode { } + SimulationNode(const Sim& property_arg, double t0, double dt_integration) + : m_simulation(property_arg, t0, dt_integration) + , m_last_state(m_simulation.get_result().get_last_value()) + , m_t0(m_simulation.get_result().get_last_time()) + { + } + + SimulationNode(const SimulationNode& other) + : m_simulation(other.m_simulation) + , m_last_state(other.m_last_state) + , m_t0(other.m_t0) + { + } + /** * get the result of the simulation in this node. * @{ @@ -292,10 +309,12 @@ class MigrationEdge template void apply_migration(double t, double dt, SimulationNode& node_from, SimulationNode& node_to); -private: +protected: MigrationParameters m_parameters; TimeSeries m_migrated; TimeSeries m_return_times; + +private: bool m_return_migrated; double m_t_last_dynamic_npi_check = -std::numeric_limits::infinity(); std::pair m_dynamic_npi = {-std::numeric_limits::max(), SimulationTime(0)}; @@ -389,8 +408,8 @@ auto get_migration_factors(const SimulationNode& node, double t, const Eige * detect a get_migration_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 migrating from their source node. From eecaff45d00e2e4e55f39eac7812617843f2fe1b Mon Sep 17 00:00:00 2001 From: HenrZu Date: Wed, 6 Dec 2023 13:43:55 +0100 Subject: [PATCH 02/11] remove doubled doc --- cpp/memilio/mobility/metapopulation_mobility_instant.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 9b2fcfa953..6d7e6ecf5d 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -39,9 +39,6 @@ namespace mio { -/** - * represents the simulation in one node of the graph. - */ /** * represents the simulation in one node of the graph. */ From f5f78139a6d88864fe4627026c94a839f6b956f3 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Wed, 6 Dec 2023 13:46:59 +0100 Subject: [PATCH 03/11] fix msvc error --- cpp/memilio/io/mobility_io.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/cpp/memilio/io/mobility_io.cpp b/cpp/memilio/io/mobility_io.cpp index 1e602927e1..d60933447f 100644 --- a/cpp/memilio/io/mobility_io.cpp +++ b/cpp/memilio/io/mobility_io.cpp @@ -183,9 +183,8 @@ IOResult>>> read_path_mobility(const st return failure(StatusCode::FileNotFound, filename); } - // std::vector> arr(std::sqrt(num_lines), std::vector(std::sqrt(num_lines))); - std::vector>> arr(std::sqrt(num_lines), - std::vector>(std::sqrt(num_lines))); + const int num_nodes = static_cast(std::sqrt(num_lines)); + std::vector>> arr(num_nodes, std::vector>(num_nodes)); try { std::string tp; From aa14826954d63136e8e3875692da6188c8060162 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Thu, 7 Dec 2023 16:16:06 +0100 Subject: [PATCH 04/11] [wip] adjust structure. --- cpp/examples/CMakeLists.txt | 5 +- .../ode_secirvvs_mobility_detailed.cpp | 946 ++++++++++++++++++ cpp/memilio/compartments/parameter_studies.h | 2 +- cpp/memilio/mobility/graph.h | 57 -- cpp/memilio/mobility/graph_simulation.h | 32 +- .../metapopulation_mobility_detailed.h | 68 +- 6 files changed, 1033 insertions(+), 77 deletions(-) create mode 100644 cpp/examples/ode_secirvvs_mobility_detailed.cpp diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 01d2822404..e996d0a39c 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -21,8 +21,9 @@ target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNI add_executable(ode_sir_example ode_sir.cpp) target_link_libraries(ode_sir_example PRIVATE memilio ode_sir) -# add_executable(2022_sim_waning 2022_simulation_omicron_parallel.cpp) -# target_link_libraries(2022_sim_waning PRIVATE memilio ode_secirvvs) +add_executable(sim_mobility_detailed ode_secirvvs_mobility_detailed.cpp) +target_link_libraries(sim_mobility_detailed PRIVATE memilio ode_secirvvs) + add_executable(seir_flows_example ode_seir_flows.cpp) target_link_libraries(seir_flows_example PRIVATE memilio ode_seir) diff --git a/cpp/examples/ode_secirvvs_mobility_detailed.cpp b/cpp/examples/ode_secirvvs_mobility_detailed.cpp new file mode 100644 index 0000000000..ecca67783a --- /dev/null +++ b/cpp/examples/ode_secirvvs_mobility_detailed.cpp @@ -0,0 +1,946 @@ +/* +* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC) +* +* 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/flow_simulation.h" +#include "memilio/config.h" +#include "memilio/data/analyze_result.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/logging.h" +#include "ode_secirvvs/model.h" +#include "ode_secirvvs/infection_state.h" +#include "ode_secirvvs/parameters.h" +#include "ode_secirvvs/parameters_io.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" +#include "memilio/compartments/simulation.h" +#include "memilio/io/result_io.h" +#include "memilio/compartments/parameter_studies.h" +#include "memilio/utils/miompi.h" +#include "memilio/utils/random_number_generator.h" + +#include "memilio/geography/regions.h" +#include "memilio/io/epi_data.h" +#include "memilio/io/mobility_io.h" +#include "memilio/utils/date.h" +#include "ode_secirvvs/parameter_space.h" +#include "memilio/utils/stl_util.h" +#include "boost/filesystem.hpp" +#include +#include +#include + +#include +#include +#include + +/** + * Set a value and distribution of an UncertainValue. + * Assigns average of min and max as a value and UNIFORM(min, max) as a distribution. + * @param p uncertain value to set. + * @param min minimum of distribution. + * @param max minimum of distribution. + */ +void assign_uniform_distribution(mio::UncertainValue& p, double min, double max) +{ + p = mio::UncertainValue(0.5 * (max + min)); + p.set_distribution(mio::ParameterDistributionUniform(min, max)); +} + +/** + * Set a value and distribution of an array of UncertainValues. + * Assigns average of min[i] and max[i] as a value and UNIFORM(min[i], max[i]) as a distribution for + * each element i of the array. + * @param array array of UncertainValues to set. + * @param min minimum of distribution for each element of array. + * @param max minimum of distribution for each element of array. + */ +template +void array_assign_uniform_distribution(mio::CustomIndexArray& array, + const double (&min)[N], const double (&max)[N]) +{ + assert(N == array.numel()); + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(N); ++i) { + assign_uniform_distribution(array[i], min[size_t(i)], max[size_t(i)]); + } +} + +/** + * Set a value and distribution of an array of UncertainValues. + * Assigns average of min and max as a value and UNIFORM(min, max) as a distribution to every element of the array. + * @param array array of UncertainValues to set. + * @param min minimum of distribution. + * @param max minimum of distribution. + */ +void array_assign_uniform_distribution(mio::CustomIndexArray& array, double min, + double max) +{ + for (auto i = mio::AgeGroup(0); i < array.size(); ++i) { + assign_uniform_distribution(array[i], min, max); + } +} + +/** + * Set epidemiological parameters of Covid19. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +mio::IOResult set_covid_parameters(mio::osecirvvs::Parameters& params) +{ + //times + const double incubationTime = 3.1; // doi.org/10.3201/eid2806.220158 + const double serialIntervalMin = 2.38; // doi.org/10.1016/j.lanepe.2022.100446 + const double serialIntervalMax = 2.38; // doi.org/10.1016/j.lanepe.2022.100446 + const double timeInfectedSymptomsMin = 6.58; //https://doi.org/10.1016/S0140-6736(22)00327-0 + const double timeInfectedSymptomsMax = 7.16; //https://doi.org/10.1016/S0140-6736(22)00327-0 + const double timeInfectedSevereMin[] = {1.8, 1.8, 1.8, 2.5, 3.5, 4.91}; // doi.org/10.1186/s12879-022-07971-6 + const double timeInfectedSevereMax[] = {2.3, 2.3, 2.3, 3.67, 5, 7.01}; // doi.org/10.1186/s12879-022-07971-6 + + const ScalarType fact_icu = 1.; + const double timeInfectedCriticalMin[] = {fact_icu * 4.95, fact_icu * 4.95, fact_icu * 4.86, + fact_icu * 14.14, fact_icu * 14.4, fact_icu * 10.}; + const double timeInfectedCriticalMax[] = {fact_icu * 8.95, fact_icu * 8.95, fact_icu * 8.86, + fact_icu * 20.58, fact_icu * 19.8, fact_icu * 13.2}; + + array_assign_uniform_distribution(params.get(), incubationTime, incubationTime); + array_assign_uniform_distribution(params.get(), serialIntervalMin, + serialIntervalMax); + array_assign_uniform_distribution(params.get(), timeInfectedSymptomsMin, + timeInfectedSymptomsMax); + array_assign_uniform_distribution(params.get(), timeInfectedSevereMin, + timeInfectedSevereMax); + array_assign_uniform_distribution(params.get(), timeInfectedCriticalMin, + timeInfectedCriticalMax); + + double fac_variant = 1.94; //https://doi.org/10.7554/eLife.78933 + const double transmissionProbabilityOnContactMin[] = {0.02 * fac_variant, 0.05 * fac_variant, 0.05 * fac_variant, + 0.05 * fac_variant, 0.08 * fac_variant, 0.1 * fac_variant}; + + const double transmissionProbabilityOnContactMax[] = {0.04 * fac_variant, 0.07 * fac_variant, 0.07 * fac_variant, + 0.07 * fac_variant, 0.10 * fac_variant, 0.15 * fac_variant}; + + const double relativeTransmissionNoSymptomsMin = 0.5; + const double relativeTransmissionNoSymptomsMax = 0.5; + const double riskOfInfectionFromSymptomaticMin = 0.0; // beta (depends on incidence and test and trace capacity) + const double riskOfInfectionFromSymptomaticMax = 0.2; + const double maxRiskOfInfectionFromSymptomaticMin = 0.4; + const double maxRiskOfInfectionFromSymptomaticMax = 0.5; + + const double recoveredPerInfectedNoSymptomsMin[] = {0.2, 0.25, 0.2, + 0.2, 0.175, 0.1}; // doi.org/10.1101/2022.05.05.22274697 + const double recoveredPerInfectedNoSymptomsMax[] = {0.4, 0.45, 0.35, + 0.3, 0.25, 0.15}; // doi.org/10.1101/2022.05.05.22274697 + + // doi:10.1136/bmjgh-2023-0123 + + // Factors from https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(22)00462-7/fulltext + const double severePerInfectedSymptomsMin[] = {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 severePerInfectedSymptomsMax[] = {1 * 0.009, 0.8 * 0.009, 0.4 * 0.023, 0.3 * 0.074, + 0.25 * 0.18, 0.35 * 0.25}; // quelle 2021 paper + factors + + // delta paper + // risk of icu admission https://doi.org/10.1177/14034948221108548 + const double fac_icu = 0.52; + const double criticalPerSevereMin[] = {0.05 * fac_icu, 0.05 * fac_icu, 0.05 * fac_icu, + 0.10 * fac_icu, 0.25 * fac_icu, 0.35 * fac_icu}; + const double criticalPerSevereMax[] = {0.10 * fac_icu, 0.10 * fac_icu, 0.10 * fac_icu, + 0.20 * fac_icu, 0.35 * fac_icu, 0.45 * fac_icu}; + + // 61% less risk doi:10.1136/bmjgh-2023-0123 + // https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10347449/pdf/bmjgh-2023-012328.pdf + const double fac_dead = 0.39; + const double deathsPerCriticalMin[] = {fac_dead * 0.00, fac_dead * 0.00, fac_dead * 0.10, + fac_dead * 0.10, fac_dead * 0.30, fac_dead * 0.5}; // 2021 paper + const double deathsPerCriticalMax[] = {fac_dead * 0.10, fac_dead * 0.10, fac_dead * 0.18, + fac_dead * 0.18, fac_dead * 0.50, fac_dead * 0.7}; + + const double reducExposedPartialImmunityMin = 0.569; // doi.org/10.1136/bmj-2022-071502 + const double reducExposedPartialImmunityMax = 0.637; // doi.org/10.1136/bmj-2022-071502 + const double reducExposedImprovedImmunityMin = + 0.34 * reducExposedPartialImmunityMin; // https://jamanetwork.com/journals/jama/fullarticle/2788487 0.19346 + const double reducExposedImprovedImmunityMax = + 0.34 * reducExposedPartialImmunityMax; // https://jamanetwork.com/journals/jama/fullarticle/2788487 0.21658 + + const double reducInfectedSymptomsPartialImmunityMin = 0.746; // doi.org/10.1056/NEJMoa2119451 + const double reducInfectedSymptomsPartialImmunityMax = 0.961; // doi.org/10.1056/NEJMoa2119451 + const double reducInfectedSymptomsImprovedImmunityMin = 0.295; // doi.org/10.1056/NEJMoa2119451 + const double reducInfectedSymptomsImprovedImmunityMax = 0.344; // doi.org/10.1056/NEJMoa2119451 + + const double reducInfectedSevereCriticalDeadPartialImmunityMin = + 0.52; // www.assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/1050721/Vaccine-surveillance-report-week-4.pdf + const double reducInfectedSevereCriticalDeadPartialImmunityMax = + 0.82; // www.assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/1050721/Vaccine-surveillance-report-week-4.pdf + const double reducInfectedSevereCriticalDeadImprovedImmunityMin = 0.1; // doi.org/10.1136/bmj-2022-071502 + const double reducInfectedSevereCriticalDeadImprovedImmunityMax = 0.19; // doi.org/10.1136/bmj-2022-071502 + const double temp_reducTimeInfectedMild = 0.5; // doi.org/10.1101/2021.09.24.21263978 + + array_assign_uniform_distribution(params.get(), + transmissionProbabilityOnContactMin, transmissionProbabilityOnContactMax); + array_assign_uniform_distribution(params.get(), + relativeTransmissionNoSymptomsMin, relativeTransmissionNoSymptomsMax); + array_assign_uniform_distribution(params.get(), + riskOfInfectionFromSymptomaticMin, riskOfInfectionFromSymptomaticMax); + array_assign_uniform_distribution(params.get(), + maxRiskOfInfectionFromSymptomaticMin, maxRiskOfInfectionFromSymptomaticMax); + array_assign_uniform_distribution(params.get(), + recoveredPerInfectedNoSymptomsMin, recoveredPerInfectedNoSymptomsMax); + array_assign_uniform_distribution(params.get(), + severePerInfectedSymptomsMin, severePerInfectedSymptomsMax); + array_assign_uniform_distribution(params.get(), criticalPerSevereMin, + criticalPerSevereMax); + array_assign_uniform_distribution(params.get(), deathsPerCriticalMin, + deathsPerCriticalMax); + + array_assign_uniform_distribution(params.get(), + reducExposedPartialImmunityMin, reducExposedPartialImmunityMax); + array_assign_uniform_distribution(params.get(), + reducExposedImprovedImmunityMin, reducExposedImprovedImmunityMax); + array_assign_uniform_distribution(params.get(), + reducInfectedSymptomsPartialImmunityMin, reducInfectedSymptomsPartialImmunityMax); + array_assign_uniform_distribution(params.get(), + reducInfectedSymptomsImprovedImmunityMin, + reducInfectedSymptomsImprovedImmunityMax); + array_assign_uniform_distribution(params.get(), + reducInfectedSevereCriticalDeadPartialImmunityMin, + reducInfectedSevereCriticalDeadPartialImmunityMax); + array_assign_uniform_distribution(params.get(), + reducInfectedSevereCriticalDeadImprovedImmunityMin, + reducInfectedSevereCriticalDeadImprovedImmunityMax); + array_assign_uniform_distribution(params.get(), temp_reducTimeInfectedMild, + temp_reducTimeInfectedMild); + + //sasonality + const double seasonality_min = 0.1; + const double seasonality_max = 0.3; + + assign_uniform_distribution(params.get(), seasonality_min, seasonality_max); + 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"}}; + +/** + * Set contact matrices. + * Reads contact matrices from files in the data directory. + * @param data_dir data directory. + * @param params Object that the contact matrices will be added to. + * @returns any io errors that happen during reading of the files. + */ +mio::IOResult set_contact_matrices(const fs::path& data_dir, mio::osecirvvs::Parameters& params, + ScalarType scale_contacts = 1.0, ScalarType share_staying = 1.0) +{ + auto contact_transport_status = mio::read_mobility_plain(data_dir.string() + "//contacts//contacts_transport.txt"); + auto contact_matrix_transport = contact_transport_status.value(); + auto contact_matrices = mio::ContactMatrixGroup(contact_locations.size(), size_t(params.get_num_groups())); + for (auto&& contact_location : contact_locations) { + BOOST_OUTCOME_TRY(baseline, + mio::read_mobility_plain( + (data_dir / "contacts" / ("baseline_" + contact_location.second + ".txt")).string())); + + if (contact_location.second == "other") { + contact_matrices[size_t(contact_location.first)].get_baseline() = + (1 - share_staying) * abs(baseline - contact_matrix_transport) / scale_contacts + + share_staying * abs(baseline - contact_matrix_transport); + contact_matrices[size_t(contact_location.first)].get_minimum() = Eigen::MatrixXd::Zero(6, 6); + } + else { + contact_matrices[size_t(contact_location.first)].get_baseline() = + (1 - share_staying) * baseline / scale_contacts + share_staying * baseline; + contact_matrices[size_t(contact_location.first)].get_minimum() = Eigen::MatrixXd::Zero(6, 6); + } + } + params.get() = mio::UncertainContactMatrix(contact_matrices); + + return mio::success(); +} + +mio::IOResult set_contact_matrices_transport(const fs::path& data_dir, mio::osecirvvs::Parameters& params, + ScalarType scale_contacts = 1.0) +{ + auto contact_transport_status = mio::read_mobility_plain(data_dir.string() + "//contacts//contacts_transport.txt"); + auto contact_matrix_transport = contact_transport_status.value(); + auto contact_matrices = mio::ContactMatrixGroup(1, size_t(params.get_num_groups())); + // ScalarType const polymod_share_contacts_transport = 1 / 0.2770885028949545; + + contact_matrices[0].get_baseline() = contact_matrix_transport / scale_contacts; + contact_matrices[0].get_minimum() = Eigen::MatrixXd::Zero(6, 6); + + params.get() = mio::UncertainContactMatrix(contact_matrices); + + return mio::success(); +} + +// reset population in graph +void init_pop_cologne_szenario(mio::Graph& graph) +{ + std::vector> immunity = {{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}}; + + for (auto& node : graph.nodes()) { + for (auto age = mio::AgeGroup(0); age < mio::AgeGroup(6); age++) { + auto pop_age = 0.0; + for (auto inf_state = mio::Index(0); + inf_state < mio::osecirvvs::InfectionState::Count; ++inf_state) { + pop_age += node.property.populations[{age, inf_state}]; + node.property.populations[{age, inf_state}] = 0.0; + } + size_t immunity_index = static_cast(age); + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleNaive}] = + pop_age * immunity[immunity_index][0]; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptiblePartialImmunity}] = + pop_age * immunity[immunity_index][1]; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleImprovedImmunity}] = + pop_age * immunity[immunity_index][2]; + } + if (node.id == 5315) { + // infect p% of population + ScalarType p = 0.05; + for (auto age = mio::AgeGroup(0); age < graph.nodes()[0].property.parameters.get_num_groups(); age++) { + node.property.populations[{mio::AgeGroup(age), mio::osecirvvs::InfectionState::InfectedSymptomsNaive}] = + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleNaive}] * p; + node.property.populations[{age, mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity}] = + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptiblePartialImmunity}] * p; + node.property.populations[{mio::AgeGroup(age), + mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity}] = + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleImprovedImmunity}] * p; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleNaive}] *= 0.99; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptiblePartialImmunity}] *= 0.99; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleImprovedImmunity}] *= 0.99; + } + } + } +} + +/** + * Create the input graph for the parameter study. + * Reads files from the data directory. + * @param start_date start date of the simulation. + * @param end_date end date of the simulation. + * @param data_dir data directory. + * @returns created graph or any io errors that happen during reading of the files. + */ +mio::IOResult> +get_graph(const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool edges) +{ + std::string data_dir = "/localdata1/test/memilio//data"; + std::string traveltimes_path = "/localdata1/code/memilio/test/travel_times_pathes.txt"; + std::string durations_path = "/localdata1/code/memilio/test/activity_duration_work.txt"; + auto mobility_data_commuter = + mio::read_mobility_plain(("/localdata1/code/memilio/test/commuter_migration_with_locals.txt")); + auto mob_data = mobility_data_commuter.value(); + auto start_date = mio::Date(2022, 8, 1); + auto end_date = mio::Date(2022, 11, 1); + + // global parameters + const int num_age_groups = 6; + mio::Graph params_graph; + + // Nodes + mio::osecirvvs::Parameters params(num_age_groups); + + params.get() = mio::get_day_in_year(start_date); + auto params_status = set_covid_parameters(params); + auto contacts_status = set_contact_matrices(data_dir, params); + + params.get() = mio::get_day_in_year(start_date); + + // Set nodes + auto scaling_factor_infected = std::vector(size_t(params.get_num_groups()), 1.0); + auto scaling_factor_icu = 1.0; + auto tnt_capacity_factor = 1.43 / 100000.; + + auto read_duration = mio::read_duration_stay(durations_path); + if (!read_duration) { + std::cout << read_duration.error().formatted_message() << '\n'; + } + auto duration_stay = read_duration.value(); + + std::string path = "/localdata1/test/memilio/data/pydata/Germany/county_population.json"; + auto read_ids = mio::get_node_ids(path, true); + auto node_ids = read_ids.value(); + + std::vector nodes(node_ids.size(), + mio::osecirvvs::Model(int(size_t(params.get_num_groups())))); + std::vector scaling_factor_inf(6, 1.0); + + for (auto& node : nodes) { + node.parameters = params; + } + auto read_node = read_input_data_county(nodes, start_date, node_ids, scaling_factor_inf, scaling_factor_icu, + "/localdata1/test/memilio/data", num_days); + + for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { + + auto tnt_capacity = nodes[node_idx].populations.get_total() * tnt_capacity_factor; + + //local parameters + auto& tnt_value = nodes[node_idx].parameters.template get(); + tnt_value = mio::UncertainValue(0.5 * (1.2 * tnt_capacity + 0.8 * tnt_capacity)); + tnt_value.set_distribution(mio::ParameterDistributionUniform(0.8 * tnt_capacity, 1.2 * tnt_capacity)); + + //holiday periods + auto id = int(mio::regions::CountyId(node_ids[node_idx])); + auto holiday_periods = mio::regions::get_holidays(mio::regions::get_state_id(id), start_date, end_date); + auto& contacts = nodes[node_idx].parameters.template get(); + contacts.get_school_holidays() = + std::vector>(holiday_periods.size()); + std::transform( + holiday_periods.begin(), holiday_periods.end(), contacts.get_school_holidays().begin(), [=](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))); + }); + + //uncertainty in populations + for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + for (auto j = mio::Index(0); + j < mio::osecirvvs::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))); + } + } + + // Add mobility node + auto mobility = nodes[node_idx]; + mobility.populations.set_total(0); + + // reduce transmission on contact due to mask obligation in mobility node + // first age group not able to (properly) wear masks + if (masks) { + + const double fact_surgical_mask = 0.1; + const double fact_ffp2 = 0.001; + + double factor_mask[] = {1, fact_surgical_mask, fact_ffp2, fact_ffp2, fact_ffp2, fact_ffp2}; + if (!ffp2) { + std::cout << "surgical masks" + << "\n"; + // multipliziere alle außer den ersten EIntrag mit 40 + for (size_t j = 1; j < 6; j++) { + factor_mask[j] = factor_mask[j] * fact_surgical_mask / fact_ffp2; + } + } + + double fac_variant = 1.4; //1.94; //https://doi.org/10.7554/eLife.78933 + double transmissionProbabilityOnContactMin[] = {0.02 * fac_variant, 0.05 * fac_variant, 0.05 * fac_variant, + 0.05 * fac_variant, 0.08 * fac_variant, 0.1 * fac_variant}; + + double transmissionProbabilityOnContactMax[] = {0.04 * fac_variant, 0.07 * fac_variant, 0.07 * fac_variant, + 0.07 * fac_variant, 0.10 * fac_variant, 0.15 * fac_variant}; + for (int i = 0; i < 6; i++) { + transmissionProbabilityOnContactMin[i] = transmissionProbabilityOnContactMin[i] * factor_mask[i]; + transmissionProbabilityOnContactMax[i] = transmissionProbabilityOnContactMax[i] * factor_mask[i]; + } + array_assign_uniform_distribution( + mobility.parameters.get(), + transmissionProbabilityOnContactMin, transmissionProbabilityOnContactMax); + } + + for (size_t t_idx = 0; t_idx < num_days; ++t_idx) { + auto t = mio::SimulationDay((size_t)t_idx); + for (auto j = mio::AgeGroup(0); j < params.get_num_groups(); j++) { + mobility.parameters.template get()[{j, t}] = 0; + mobility.parameters.template get()[{j, t}] = 0; + } + } + + auto& params_mobility = mobility.parameters; + + params_graph.add_node(node_ids[node_idx], duration_stay((Eigen::Index)node_idx), nodes[node_idx], mobility); + } + + if (szenario_cologne) { + init_pop_cologne_szenario(params_graph); + } + + // Edges + if (edges) { + ScalarType theshold_edges = 4e-5; + auto migrating_compartments = { + mio::osecirvvs::InfectionState::SusceptibleNaive, + mio::osecirvvs::InfectionState::ExposedNaive, + mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive, + mio::osecirvvs::InfectionState::InfectedSymptomsNaive, + mio::osecirvvs::InfectionState::SusceptibleImprovedImmunity, + mio::osecirvvs::InfectionState::SusceptiblePartialImmunity, + mio::osecirvvs::InfectionState::ExposedPartialImmunity, + mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity, + mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity, + mio::osecirvvs::InfectionState::ExposedImprovedImmunity, + mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity, + mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity, + }; + + // mobility between nodes + auto read_travel_times = mio::read_mobility_plain(traveltimes_path); + auto travel_times = read_travel_times.value(); + auto read_paths_mobility = mio::read_path_mobility( + "/localdata1/code/memilio/test/wegketten_ohne_komma.txt"); // wegketten_ohne_komma.txt + auto path_mobility = read_paths_mobility.value(); + + 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) { + // if (county_idx_i == county_idx_j) + // continue; + auto& populations = nodes[county_idx_i].populations; + // mobility coefficients have the same number of components as the contact matrices. + // so that the same NPIs/dampings can be used for both (e.g. more home office => fewer commuters) + auto mobility_coeffs = mio::MigrationCoefficientGroup(contact_locations.size(), populations.numel()); + + //commuters + auto working_population = 0.0; + auto min_commuter_age = mio::AgeGroup(2); + auto max_commuter_age = mio::AgeGroup(4); //this group is partially retired, only partially commutes + for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { + working_population += populations.get_group_total(age) * (age == max_commuter_age ? 0.33 : 1.0); + } + auto commuter_coeff_ij = mob_data(county_idx_i, county_idx_j) / + working_population; //data is absolute numbers, we need relative + for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { + for (auto compartment : migrating_compartments) { + auto coeff_index = populations.get_flat_index({age, compartment}); + mobility_coeffs[size_t(ContactLocation::Work)].get_baseline()[coeff_index] = + commuter_coeff_ij * (age == max_commuter_age ? 0.33 : 1.0); + } + } + + auto path = path_mobility[county_idx_i][county_idx_j]; + if (path[0] != county_idx_i || path[path.size() - 1] != county_idx_j) + std::cout << "falsch" + << "\n"; + if (commuter_coeff_ij > theshold_edges) { + params_graph.add_edge(county_idx_i, county_idx_j, travel_times(county_idx_i, county_idx_j), + path_mobility[county_idx_i][county_idx_j], std::move(mobility_coeffs)); + } + } + } + + // average travel time + const ScalarType avg_traveltime = std::accumulate(params_graph.edges().begin(), params_graph.edges().end(), 0.0, + [](double sum, const auto& e) { + return sum + e.traveltime; + }) / + params_graph.edges().size(); + + // average share of commuters for all counties relative to the total population + const ScalarType total_population = std::accumulate(params_graph.nodes().begin(), params_graph.nodes().end(), + 0.0, [](double sum, const auto& n) { + return sum + n.property.populations.get_total(); + }); + const ScalarType num_commuters = std::accumulate(params_graph.edges().begin(), params_graph.edges().end(), 0.0, + [mob_data](double sum, const auto& e) { + auto start_node = e.start_node_idx; + auto end_node = e.end_node_idx; + return sum + mob_data(start_node, end_node); + }); + const ScalarType avg_commuter_share = num_commuters / total_population; + + if (mio::mpi::is_root()) { + std::cout << "avg commuter share = " << avg_commuter_share << "\n"; + std::cout << "avg travel time = " << avg_traveltime << "\n"; + } + + // scale all contact matrices + for (auto& node : params_graph.nodes()) { + // mobility node + auto contacts_status_mobility = + set_contact_matrices_transport(data_dir, node.mobility.parameters, avg_traveltime * avg_commuter_share); + + // local node + auto contacts_status = + set_contact_matrices(data_dir, node.property.parameters, 1 - avg_traveltime, 1 - avg_commuter_share); + } + mio::ContactMatrixGroup const& contact_matrix3 = + params_graph.nodes()[0].mobility.parameters.get(); + std::cout << "Mobility node: contact matrix at t = 0\n"; + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(6); i++) { + for (auto j = mio::AgeGroup(0); j < mio::AgeGroup(6); j++) { + std::cout << contact_matrix3.get_matrix_at(0)(static_cast((size_t)i), + static_cast((size_t)j)) + << " "; + } + std::cout << "\n"; + } + } + + if (mio::mpi::is_root()) { + std::cout << "Anzahl kanten = " << params_graph.edges().size() << "\n"; + } + return params_graph; +} +mio::IOResult run() +{ + // mio::set_log_level(mio::LogLevel::critical); + const auto num_days = 70; + const int num_runs = 150; + const bool masks = true; + const bool ffp2 = true; + const bool edges = true; + + // wenn masks false und ffp2 true, dann error ausgeben + if (!masks && ffp2) { + mio::log_error("ffp2 only possible with masks"); + } + + bool szenario_cologne = false; + + // auto params_graph = get_graph(num_days, masks); + BOOST_OUTCOME_TRY(created, get_graph(num_days, masks, ffp2, szenario_cologne, edges)); + auto params_graph = created; + + std::string res_dir = "/localdata1/code/memilio/results_paper/mask_" + std::to_string(masks); + + // if (ffp2) + // res_dir = res_dir + "_ffp2"; + // if (szenario_cologne) + // res_dir = res_dir + "_cologne"; + // if (!edges) + // res_dir = res_dir + "_no_edges"; + if (!edges && (ffp2 || masks)) + mio::log_error("no edges only possible without masks and ffp2"); + + res_dir += std::string(ffp2 ? "_ffp2" : "") + std::string(szenario_cologne ? "_cologne" : "") + + std::string(!edges ? "_no_edges" : ""); + + // check if boths dir exist, otherwise create them + if (!fs::exists(res_dir)) { + fs::create_directory(res_dir); + } + auto write_graph_status = write_graph(params_graph, "/localdata1/code/memilio/save_graph"); + + std::vector county_ids(params_graph.nodes().size()); + std::transform(params_graph.nodes().begin(), params_graph.nodes().end(), county_ids.begin(), [](auto& n) { + return n.id; + }); + + // parameter study + auto parameter_study = + mio::ParameterStudy>(params_graph, 0.0, num_days, 0.01, num_runs); + if (mio::mpi::is_root()) { + printf("Seeds: "); + for (auto s : parameter_study.get_rng().get_seeds()) { + printf("%u, ", s); + } + printf("\n"); + } + auto save_single_run_result = mio::IOResult(mio::success()); + auto ensemble = parameter_study.run( + [&](auto&& graph) { + return draw_sample(graph, false); + }, + [&](auto results_graph, auto&& run_idx) { + auto interpolated_result = mio::interpolate_simulation_result(results_graph); + + auto params = std::vector(); + // params.reserve(results_graph.nodes().size()); + // std::transform(results_graph.nodes().begin(), results_graph.nodes().end(), std::back_inserter(params), + // [](auto&& node) { + // return node.property.get_simulation().get_model(); + // }); + // const auto size_vec = interpolated_result[0].get_num_time_points(); + // std::vector num_transmissions_in_mobility_nodes(size_vec, 0.0); + // std::vector> num_transmission_per_county(results_graph.nodes().size(), + // std::vector(size_vec, 0.0)); + // std::vector> num_symptomatic_per_county(results_graph.nodes().size(), + // std::vector(size_vec, 0.0)); + + // // vaccinations + // // std::vector> vacc_county_N(results_graph.nodes().size()); + // // std::vector> vacc_county_PI(results_graph.nodes().size()); + // // std::vector> vacc_county_II(results_graph.nodes().size()); + + // std::vector waning_partial_immunity(size_vec, 0.0); + // std::vector waning_improved_immunity(size_vec, 0.0); + // auto node_idx = 0; + // for (auto& node : results_graph.nodes()) { + // auto flows = mio::interpolate_simulation_result(node.mobility.get_simulation().get_flows()); + // auto flows_local_model = mio::interpolate_simulation_result(node.property.get_simulation().get_flows()); + // for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(6); ++i) { + // auto& model = node.mobility.get_simulation().get_model(); + // auto flow_index_N = + // model.template get_flat_flow_index({i}); + // auto flow_index_N_confirmed = model.template get_flat_flow_index< + // mio::osecirvvs::InfectionState::InfectedNoSymptomsNaiveConfirmed, + // mio::osecirvvs::InfectionState::InfectedSymptomsNaiveConfirmed>({i}); + // auto flow_index_PI = model.template get_flat_flow_index< + // mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity, + // mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity>({i}); + // auto flow_index_PI_confirmed = model.template get_flat_flow_index< + // mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunityConfirmed, + // mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunityConfirmed>({i}); + // auto flow_index_II = model.template get_flat_flow_index< + // mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity, + // mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity>({i}); + // auto flow_index_II_confirmed = model.template get_flat_flow_index< + // mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed, + // mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunityConfirmed>({i}); + + // // Flows from Susceptible to Exposed + // auto flow_index_SN = + // model.template get_flat_flow_index({i}); + // auto flow_index_SII = + // model.template get_flat_flow_index( + // {i}); + // auto flow_index_SPI = + // model.template get_flat_flow_index({i}); + + // auto flow_index_waning_pi = + // model.template get_flat_flow_index({i}); + // auto flow_index_waning_ii = + // model.template get_flat_flow_index( + // {i}); + // for (Eigen::Index time_idx = 0; time_idx < flows.get_num_time_points(); ++time_idx) { + // const auto& flow_values = flows.get_value(time_idx); + // auto time = flows.get_time(time_idx); + + // auto indx_vector = static_cast(time_idx); + // // num_transmissions_in_mobility_nodes[indx_vector] += + // // flow_values(flow_index_N) + flow_values(flow_index_PI) + flow_values(flow_index_II) + + // // flow_values(flow_index_N_confirmed) + flow_values(flow_index_PI_confirmed) + + // // flow_values(flow_index_II_confirmed); + // num_transmissions_in_mobility_nodes[indx_vector] += + // flow_values(flow_index_SN) + flow_values(flow_index_SPI) + flow_values(flow_index_SII); + // waning_partial_immunity[indx_vector] += flow_values(flow_index_waning_pi); + // waning_improved_immunity[indx_vector] += flow_values(flow_index_waning_ii); + // num_symptomatic_per_county[node_idx][indx_vector] += + // flow_values(flow_index_N) + flow_values(flow_index_PI) + flow_values(flow_index_II) + + // flow_values(flow_index_N_confirmed) + flow_values(flow_index_PI_confirmed) + + // flow_values(flow_index_II_confirmed); + // num_transmission_per_county[node_idx][indx_vector] += + // flow_values(flow_index_SN) + flow_values(flow_index_SPI) + flow_values(flow_index_SII); + // } + // for (Eigen::Index time_idx = 0; time_idx < flows_local_model.get_num_time_points(); ++time_idx) { + // const auto& flow_values_local_model = flows_local_model.get_value(time_idx); + // auto time = flows_local_model.get_time(time_idx); + // auto indx_vector = static_cast(time_idx); + // waning_partial_immunity[indx_vector] += flow_values_local_model(flow_index_waning_pi); + // waning_improved_immunity[indx_vector] += flow_values_local_model(flow_index_waning_ii); + // num_transmission_per_county[node_idx][indx_vector] += flow_values_local_model(flow_index_SN) + + // flow_values_local_model(flow_index_SPI) + + // flow_values_local_model(flow_index_SII); + // num_symptomatic_per_county[node_idx][indx_vector] += + // flow_values_local_model(flow_index_N) + flow_values_local_model(flow_index_PI) + + // flow_values_local_model(flow_index_II) + flow_values_local_model(flow_index_N_confirmed) + + // flow_values_local_model(flow_index_PI_confirmed) + + // flow_values_local_model(flow_index_II_confirmed); + // } + // } + + // // // vaccinations + // // vacc_county_N[node_idx].resize(flows_local_model.get_num_time_points(), 0.0); + // // vacc_county_PI[node_idx].resize(flows_local_model.get_num_time_points(), 0.0); + // // vacc_county_II[node_idx].resize(flows_local_model.get_num_time_points(), 0.0); + // // for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(6); ++i) { + // // auto& model = node.mobility.get_simulation().get_model(); + // // auto vacc_N = + // // model.template get_flat_flow_index({i}); + // // auto flow_index_PI = + // // model.template get_flat_flow_index({i}); + // // auto flow_index_II = + // // model.template get_flat_flow_index({i}); + + // // for (Eigen::Index time_idx = 0; time_idx < flows_local_model.get_num_time_points(); ++time_idx) { + // // const auto& flow_values_local_model = flows_local_model.get_value(time_idx); + // // auto time = flows_local_model.get_time(time_idx); + // // auto indx_vector = static_cast(time_idx); + // // vacc_county_N[node_idx][indx_vector] += flow_values_local_model(vacc_N); + // // vacc_county_PI[node_idx][indx_vector] += flow_values_local_model(flow_index_PI); + // // vacc_county_II[node_idx][indx_vector] += flow_values_local_model(flow_index_II); + // // } + // // } + // node_idx++; + // } + // std::string results_path = "/localdata1/code/memilio/results_paper/mask_" + std::to_string(masks); + // if (ffp2) + // results_path = results_path + "_ffp2"; + + // if (szenario_cologne) + // results_path = results_path + "_cologne"; + + // if (!edges) + // results_path = results_path + "_no_edges"; + // const std::string results_path_mb = results_path + "/flows_mb"; + // if (!fs::exists(results_path_mb)) { + // fs::create_directory(results_path_mb); + // } + // const std::string results_filename = results_path_mb + "/results_run_" + std::to_string(run_idx) + ".txt"; + // std::ofstream outfile(results_filename); + + // outfile << "TimeIndex Num_Transmissions_in_Mobility_Nodes Waning_Partial_Immunity " + // "Waning_Improved_Immunity\n"; + // size_t vec_size = num_transmissions_in_mobility_nodes.size(); + // for (size_t i = 0; i < vec_size; ++i) { + // outfile << i << " " << num_transmissions_in_mobility_nodes[i] << " " << waning_partial_immunity[i] + // << " " << waning_improved_immunity[i] << "\n"; + // } + // outfile.close(); + + // // erzeuge ordner flows_county + // const std::string results_path_flows = results_path + "/flows_county"; + // if (!fs::exists(results_path_flows)) { + // fs::create_directory(results_path_flows); + // } + // const std::string results_filename_county = + // results_path_flows + "/transmission_county_run_" + std::to_string(run_idx) + ".txt"; + // std::ofstream outfile_county(results_filename_county); + // for (size_t i = 0; i < results_graph.nodes().size(); ++i) { + // outfile_county << i << " "; + // for (size_t j = 0; j < vec_size; ++j) { + // outfile_county << num_transmission_per_county[i][j] << " "; + // } + // outfile_county << "\n"; + // } + // outfile_county.close(); + + // const std::string results_filename_county_symptomatic = + // results_path_flows + "/symptomatic_county_run_" + std::to_string(run_idx) + ".txt"; + // std::ofstream outfile_county_symptomatic(results_filename_county_symptomatic); + // for (size_t i = 0; i < results_graph.nodes().size(); ++i) { + // outfile_county_symptomatic << i << " "; + // for (size_t j = 0; j < vec_size; ++j) { + // outfile_county_symptomatic << num_symptomatic_per_county[i][j] << " "; + // } + // outfile_county_symptomatic << "\n"; + // } + // outfile_county_symptomatic.close(); + + // // // speichere vacc_county_N, vacc_county_PI, vacc_county_II in jeweils eigene Datei + // // const std::string results_path_vacc = results_path + "/vacc_county"; + // // // check if dir exist, otherwise create them + // // if (!fs::exists(results_path_vacc)) { + // // fs::create_directory(results_path_vacc); + // // } + // // const std::string results_filename_vacc_N = + // // results_path_vacc + "/vacc_county_N_run_" + std::to_string(run_idx) + ".txt"; + // // const std::string results_filename_vacc_PI = + // // results_path_vacc + "/vacc_county_PI_run_" + std::to_string(run_idx) + ".txt"; + // // const std::string results_filename_vacc_II = + // // results_path_vacc + "/vacc_county_II_run_" + std::to_string(run_idx) + ".txt"; + + // // std::ofstream outfile_vacc_N(results_filename_vacc_N); + // // std::ofstream outfile_vacc_PI(results_filename_vacc_PI); + // // std::ofstream outfile_vacc_II(results_filename_vacc_II); + + // // for (size_t i = 0; i < results_graph.nodes().size(); ++i) { + // // outfile_vacc_N << i << " "; + // // outfile_vacc_PI << i << " "; + // // outfile_vacc_II << i << " "; + + // // // Abfrage der Größe des aktuellen inneren Vektors: + // // size_t current_vec_size_N = vacc_county_N[i].size(); + // // size_t current_vec_size_PI = vacc_county_PI[i].size(); + // // size_t current_vec_size_II = vacc_county_II[i].size(); + + // // for (size_t j = 0; j < current_vec_size_N; ++j) { + // // outfile_vacc_N << vacc_county_N[i][j] << " "; + // // } + // // for (size_t j = 0; j < current_vec_size_PI; ++j) { + // // outfile_vacc_PI << vacc_county_PI[i][j] << " "; + // // } + // // for (size_t j = 0; j < current_vec_size_II; ++j) { + // // outfile_vacc_II << vacc_county_II[i][j] << " "; + // // } + + // // outfile_vacc_N << "\n"; + // // outfile_vacc_PI << "\n"; + // // outfile_vacc_II << "\n"; + // // } + // // outfile_vacc_N.close(); + // // outfile_vacc_PI.close(); + // // outfile_vacc_II.close(); + + // std::cout << "Run " << run_idx << " complete." << std::endl; + return std::make_pair(interpolated_result, params); + }); + + if (ensemble.size() > 0) { + auto ensemble_results = std::vector>>{}; + ensemble_results.reserve(ensemble.size()); + auto ensemble_params = std::vector>{}; + ensemble_params.reserve(ensemble.size()); + for (auto&& run : ensemble) { + ensemble_results.emplace_back(std::move(run.first)); + ensemble_params.emplace_back(std::move(run.second)); + } + BOOST_OUTCOME_TRY(save_results(ensemble_results, ensemble_params, county_ids, res_dir, false)); + } + + return mio::success(); +} + +int main(int argc, char** argv) +{ + mio::set_log_level(mio::LogLevel::warn); + mio::mpi::init(); + + auto result = run(); + if (!result) { + printf("%s\n", result.error().formatted_message().c_str()); + mio::mpi::finalize(); + return -1; + } + mio::mpi::finalize(); + return 0; +} \ No newline at end of file diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 9147f37860..e8005d8e58 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -344,7 +344,7 @@ class ParameterStudy auto sampled_graph = sample_graph(m_graph); for (auto&& node : sampled_graph.nodes()) { - sim_graph.add_node(node.id, node.stay_duration, node.property, node.node_pt, m_t0, m_dt_integration); + sim_graph.add_node(node.id, node.stay_duration, node.property, node.mobility, m_t0, m_dt_integration); } for (auto&& edge : sampled_graph.edges()) { sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.path, edge.property); diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index 6cd0e946c9..92950e986e 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -71,79 +71,22 @@ struct Node { template Node(int node_id, Args&&... args) : id{node_id} - , stay_duration(0.5) , property(std::forward(args)...) - , node_pt(std::forward(args)...) { } - - template - Node(int node_id, double duration, Model property_arg, Model node_pt_arg, double m_t0, double m_dt_integration) - : id{node_id} - , stay_duration(duration) - , property(property_arg, m_t0, m_dt_integration) - , node_pt(node_pt_arg, m_t0, m_dt_integration) - { - } - - template - Node(int node_id, double duration, Args&&... args) - : id{node_id} - , stay_duration(duration) - , property(std::forward(args)...) - , node_pt(std::forward(args)...) - { - } - - Node(int node_id, double duration, NodePropertyT property_arg, NodePropertyT node_pt_arg) - : id{node_id} - , stay_duration(duration) - , property(property_arg) - , node_pt(node_pt_arg) - { - } - int id; - double stay_duration; NodePropertyT property; - NodePropertyT node_pt; }; -/** - * @brief represents an edge of the graph - */ template struct Edge : public EdgeBase { - template Edge(size_t start, size_t end, Args&&... args) : EdgeBase{start, end} - , traveltime(0.) - , path{static_cast(start), static_cast(end)} - , property(std::forward(args)...) - { - } - - template - Edge(size_t start, size_t end, double t_travel, Args&&... args) - : EdgeBase{start, end} - , traveltime(t_travel) - , path{static_cast(start), static_cast(end)} - , property(std::forward(args)...) - { - } - - template - Edge(size_t start, size_t end, double t_travel, std::vector path_mobility, Args&&... args) - : EdgeBase{start, end} - , traveltime(t_travel) - , path(path_mobility) , property(std::forward(args)...) { } - double traveltime; - std::vector path; EdgePropertyT property; }; diff --git a/cpp/memilio/mobility/graph_simulation.h b/cpp/memilio/mobility/graph_simulation.h index 255d9d8c21..59311bb2b1 100644 --- a/cpp/memilio/mobility/graph_simulation.h +++ b/cpp/memilio/mobility/graph_simulation.h @@ -297,7 +297,7 @@ class GraphSimulationDetailed // set population to zero in mobility nodes before starting for (auto& n : Base::m_graph.nodes()) { - n.node_pt.get_result().get_last_value().setZero(); + n.mobility.get_result().get_last_value().setZero(); } auto min_dt = 0.01; @@ -319,7 +319,7 @@ class GraphSimulationDetailed } for (auto& node : Base::m_graph.nodes()) { - node.node_pt.get_simulation().set_integrator(std::make_shared()); + node.mobility.get_simulation().set_integrator(std::make_shared()); } size_t indx_schedule = 0; @@ -329,7 +329,7 @@ class GraphSimulationDetailed // first mobility activity if (indx_schedule == first_mobility[edge_indx]) { auto& node_from = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; - auto& node_to = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].node_pt; + auto& node_to = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility; m_edge_func(Base::m_t, 0.0, e.property, node_from, node_to, 0); } // next mobility activity @@ -373,12 +373,12 @@ class GraphSimulationDetailed mobility_schedule_edges[edge_indx][indx_schedule - 1])) { auto& node_from = mobility_schedule_edges[edge_indx][indx_schedule - 1] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].node_pt + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; auto& node_to = mobility_schedule_edges[edge_indx][indx_schedule] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].node_pt + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; if (indx_schedule < mobility_schedule_edges[edge_indx].size() - 1) { @@ -393,12 +393,12 @@ class GraphSimulationDetailed else { auto& node_from = mobility_schedule_edges[edge_indx][indx_schedule - 1] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].node_pt + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; auto& node_to = mobility_schedule_edges[edge_indx][indx_schedule] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].node_pt + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; assert(node_from.get_result().get_last_value() == node_to.get_result().get_last_value()); @@ -440,19 +440,19 @@ class GraphSimulationDetailed // get all time points from the last integration step auto& last_time_point = - n.node_pt.get_result().get_time(n.node_pt.get_result().get_num_time_points() - 1); + n.mobility.get_result().get_time(n.mobility.get_result().get_num_time_points() - 1); // wenn last_time_point nicht innerhalb eines intervalls von 1-e10 von t liegt, dann setzte den letzten Zeitpunkt auf m_t if (std::fabs(last_time_point - Base::m_t) > 1e-10) { - n.node_pt.get_result().get_last_time() = Base::m_t; + n.mobility.get_result().get_last_time() = Base::m_t; } - m_node_func(Base::m_t, next_dt, n.node_pt); + m_node_func(Base::m_t, next_dt, n.mobility); Eigen::Index indx_min; - while (n.node_pt.get_result().get_last_value().minCoeff(&indx_min) < -1e-7) { + while (n.mobility.get_result().get_last_value().minCoeff(&indx_min) < -1e-7) { Eigen::Index indx_max; - n.node_pt.get_result().get_last_value().maxCoeff(&indx_max); - n.node_pt.get_result().get_last_value()[indx_max] -= - n.node_pt.get_result().get_last_value()[indx_min]; - n.node_pt.get_result().get_last_value()[indx_min] = 0; + n.mobility.get_result().get_last_value().maxCoeff(&indx_max); + n.mobility.get_result().get_last_value()[indx_max] -= + n.mobility.get_result().get_last_value()[indx_min]; + n.mobility.get_result().get_last_value()[indx_min] = 0; } } indx_schedule++; @@ -460,7 +460,7 @@ class GraphSimulationDetailed } // set each compartment zero for all mobility nodes since we only model daily mobility for (auto& n : Base::m_graph.nodes()) { - n.node_pt.get_result().get_last_value().setZero(); + n.mobility.get_result().get_last_value().setZero(); } // // messe die zeit, wie lange eine iteration bis zu dieser stelle dauert // auto end = std::chrono::system_clock::now(); diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h index 3f79166968..b0bfa4c129 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_detailed.h +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -32,6 +32,7 @@ #include "memilio/epidemiology/dynamic_npis.h" #include "memilio/compartments/simulation.h" #include "memilio/utils/date.h" +#include "memilio/mobility/graph.h" #include "boost/filesystem.hpp" @@ -40,7 +41,72 @@ namespace mio { -// Erbe von MigrationEdge aus metapopulation_mobility_instant.h und überschreibe die apply_migration Funktion indem diese zusätzlich das argument mode bekommt +template +struct NodeDetailed : Node { + template + NodeDetailed(int node_id, Args&&... args) + : Node(node_id, std::forward(args)...) + , stay_duration(0.5) + , mobility(std::forward(args)...) + { + } + + template + NodeDetailed(int node_id, double duration, Model property_arg, Model mobility_arg, double m_t0, + double m_dt_integration) + : Node(node_id, property_arg, m_t0, m_dt_integration) + , stay_duration(duration) + , mobility(mobility_arg, m_t0, m_dt_integration) + { + } + + template + NodeDetailed(int node_id, double duration, Args&&... args) + : Node(node_id, std::forward(args)...) + , stay_duration(duration) + , mobility(std::forward(args)...) + { + } + + NodeDetailed(int node_id, double duration, NodePropertyT property_arg, NodePropertyT mobility_pt_arg) + : Node(node_id, property_arg) + , stay_duration(duration) + , mobility(mobility_pt_arg) + { + } + + double stay_duration; + NodePropertyT mobility; +}; + +template +struct EdgeDetailed : Edge { + template + EdgeDetailed(size_t start, size_t end, Args&&... args) + : Edge(start, end, std::forward(args)...) + , traveltime(0.) + , path{static_cast(start), static_cast(end)} + { + } + + template + EdgeDetailed(size_t start, size_t end, double t_travel, Args&&... args) + : Edge(start, end, std::forward(args)...) + , traveltime(t_travel) + , path{static_cast(start), static_cast(end)} + { + } + + template + EdgeDetailed(size_t start, size_t end, double t_travel, std::vector path_mobility, Args&&... args) + : Edge(start, end, std::forward(args)...) + , traveltime(t_travel) + , path(path_mobility) + { + } + double traveltime; + std::vector path; +}; class MigrationEdgeDetailed : public MigrationEdge { public: From 441c7347eb0ff5ae9ce48db81bbbd2818d289669 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Fri, 8 Dec 2023 09:40:02 +0100 Subject: [PATCH 05/11] detailed graph --- cpp/memilio/mobility/graph.h | 54 ------------------- .../metapopulation_mobility_detailed.h | 51 ++++++++++++++++++ 2 files changed, 51 insertions(+), 54 deletions(-) diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index 92950e986e..0e59d6b389 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -161,29 +161,6 @@ class Graph return m_nodes.back(); } - template - Node& add_node(int id, double duration_stay, Args&&... args) - { - m_nodes.emplace_back(id, duration_stay, std::forward(args)...); - return m_nodes.back(); - } - - // Ändere die Funktionsdeklaration von add_node - template - Node& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2) - { - m_nodes.emplace_back(id, duration_stay, model1, model2); - return m_nodes.back(); - } - - template - Node& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2, double m_t0, - double m_dt_integration) - { - m_nodes.emplace_back(id, duration_stay, model1, model2, m_t0, m_dt_integration); - return m_nodes.back(); - } - /** * @brief add an edge to the graph. property of the edge is constructed from arguments. */ @@ -200,37 +177,6 @@ class Graph }); } - /** - * @brief add an edge to the graph. property of the edge is constructed from arguments. - */ - template - Edge& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, Args&&... args) - { - assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); - return *insert_sorted_replace( - m_edges, Edge(start_node_idx, end_node_idx, traveltime, std::forward(args)...), - [](auto&& e1, auto&& e2) { - return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx - : e1.start_node_idx < e2.start_node_idx; - }); - } - - /** - * @brief add an edge to the graph. property of the edge is constructed from arguments. - */ - template - Edge& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, std::vector path, - Args&&... args) - { - assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); - return *insert_sorted_replace( - m_edges, Edge(start_node_idx, end_node_idx, traveltime, path, std::forward(args)...), - [](auto&& e1, auto&& e2) { - return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx - : e1.start_node_idx < e2.start_node_idx; - }); - } - /** * @brief range of nodes */ diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h index b0bfa4c129..913a1ce21f 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_detailed.h +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -107,6 +107,57 @@ struct EdgeDetailed : Edge { double traveltime; std::vector path; }; + +template +class GraphDetailed : public Graph +{ +public: + using Graph::Graph; // Vererbt Konstruktoren + + // Spezifische add_node Methoden + template + NodeDetailed& add_node(int id, double duration_stay, Args&&... args) + { + this->m_nodes.emplace_back(id, duration_stay, std::forward(args)...); + return this->m_nodes.back(); + } + + template + NodeDetailed& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2) + { + this->m_nodes.emplace_back(id, duration_stay, model1, model2); + return this->m_nodes.back(); + } + + template + NodeDetailed& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2, + double m_t0, double m_dt_integration) + { + this->m_nodes.emplace_back(id, duration_stay, model1, model2, m_t0, m_dt_integration); + return this->m_nodes.back(); + } + + // Spezifische add_edge Methoden + template + EdgeDetailed& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, Args&&... args) + { + assert(this->m_nodes.size() > start_node_idx && this->m_nodes.size() > end_node_idx); + // Hier die Implementierung von insert_sorted_replace + // ... + return this->m_edges.back(); + } + + template + EdgeDetailed& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, + std::vector path, Args&&... args) + { + assert(this->m_nodes.size() > start_node_idx && this->m_nodes.size() > end_node_idx); + // Hier die Implementierung von insert_sorted_replace + // ... + return this->m_edges.back(); + } +}; + class MigrationEdgeDetailed : public MigrationEdge { public: From 84a0d6a4c1cf0ff124e628216b220ef00b91336e Mon Sep 17 00:00:00 2001 From: HenrZu Date: Fri, 8 Dec 2023 10:28:25 +0100 Subject: [PATCH 06/11] now compatible with old version --- cpp/examples/CMakeLists.txt | 5 ++-- cpp/memilio/compartments/parameter_studies.h | 30 ++++++++++++++++++-- cpp/memilio/mobility/graph.h | 3 -- 3 files changed, 30 insertions(+), 8 deletions(-) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index e996d0a39c..ca5069e1d4 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -21,9 +21,8 @@ target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNI add_executable(ode_sir_example ode_sir.cpp) target_link_libraries(ode_sir_example PRIVATE memilio ode_sir) -add_executable(sim_mobility_detailed ode_secirvvs_mobility_detailed.cpp) -target_link_libraries(sim_mobility_detailed PRIVATE memilio ode_secirvvs) - +# add_executable(sim_mobility_detailed ode_secirvvs_mobility_detailed.cpp) +# target_link_libraries(sim_mobility_detailed PRIVATE memilio ode_secirvvs) add_executable(seir_flows_example ode_seir_flows.cpp) target_link_libraries(seir_flows_example PRIVATE memilio ode_seir) diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index e8005d8e58..63e76cdd1c 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -38,6 +38,22 @@ namespace mio { +template +struct has_stay_duration : std::false_type { +}; + +template +struct has_stay_duration().stay_duration)>> : std::true_type { +}; + +template +struct has_traveltime : std::false_type { +}; + +template +struct has_traveltime().traveltime)>> : std::true_type { +}; + /** * Class that performs multiple simulation runs with randomly sampled parameters. * Can simulate migration graphs with one simulation in each node or single simulations. @@ -344,10 +360,20 @@ class ParameterStudy auto sampled_graph = sample_graph(m_graph); for (auto&& node : sampled_graph.nodes()) { - sim_graph.add_node(node.id, node.stay_duration, node.property, node.mobility, m_t0, m_dt_integration); + if constexpr (has_stay_duration::value) { + sim_graph.add_node(node.id, node.stay_duration, node.property, m_t0, m_dt_integration); + } + else { + sim_graph.add_node(node.id, node.property, m_t0, m_dt_integration); + } } for (auto&& edge : sampled_graph.edges()) { - sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.path, edge.property); + if constexpr (has_traveltime::value) { + sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.property); + } + else { + sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property); + } } return make_migration_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index 0e59d6b389..c0f248fa29 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -29,9 +29,6 @@ #include "memilio/epidemiology/damping.h" #include "memilio/geography/regions.h" #include -#include -#include - #include "boost/filesystem.hpp" //is used to provide some paths as function arguments From 8810b920419afe9ba39a014b37ce3b5254a7b189 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Fri, 8 Dec 2023 12:46:28 +0100 Subject: [PATCH 07/11] [wip] set nodes function --- .../ode_secirvvs_mobility_detailed.cpp | 2 +- cpp/memilio/mobility/graph.h | 6 +- .../metapopulation_mobility_detailed.h | 159 ++++++++++++++++-- 3 files changed, 154 insertions(+), 13 deletions(-) diff --git a/cpp/examples/ode_secirvvs_mobility_detailed.cpp b/cpp/examples/ode_secirvvs_mobility_detailed.cpp index ecca67783a..42ce21ccd7 100644 --- a/cpp/examples/ode_secirvvs_mobility_detailed.cpp +++ b/cpp/examples/ode_secirvvs_mobility_detailed.cpp @@ -930,7 +930,7 @@ mio::IOResult run() return mio::success(); } -int main(int argc, char** argv) +int main(int, char**) { mio::set_log_level(mio::LogLevel::warn); mio::mpi::init(); diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index c0f248fa29..fefc864076 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -29,6 +29,7 @@ #include "memilio/epidemiology/damping.h" #include "memilio/geography/regions.h" #include + #include "boost/filesystem.hpp" //is used to provide some paths as function arguments @@ -75,6 +76,9 @@ struct Node { NodePropertyT property; }; +/** + * @brief represents an edge of the graph + */ template struct Edge : public EdgeBase { template @@ -231,7 +235,7 @@ class Graph })); } -private: +protected: std::vector> m_nodes; std::vector> m_edges; }; // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h index 913a1ce21f..a5c2cd1b51 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_detailed.h +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -33,6 +33,7 @@ #include "memilio/compartments/simulation.h" #include "memilio/utils/date.h" #include "memilio/mobility/graph.h" +#include "memilio/io/mobility_io.h" #include "boost/filesystem.hpp" @@ -112,9 +113,8 @@ template class GraphDetailed : public Graph { public: - using Graph::Graph; // Vererbt Konstruktoren + using Graph::Graph; - // Spezifische add_node Methoden template NodeDetailed& add_node(int id, double duration_stay, Args&&... args) { @@ -137,27 +137,164 @@ class GraphDetailed : public Graph return this->m_nodes.back(); } - // Spezifische add_edge Methoden template EdgeDetailed& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, Args&&... args) { - assert(this->m_nodes.size() > start_node_idx && this->m_nodes.size() > end_node_idx); - // Hier die Implementierung von insert_sorted_replace - // ... - return this->m_edges.back(); + assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); + return *insert_sorted_replace( + this->m_edges, Edge(start_node_idx, end_node_idx, traveltime, std::forward(args)...), + [](auto&& e1, auto&& e2) { + return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx + : e1.start_node_idx < e2.start_node_idx; + }); } template EdgeDetailed& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, std::vector path, Args&&... args) { - assert(this->m_nodes.size() > start_node_idx && this->m_nodes.size() > end_node_idx); - // Hier die Implementierung von insert_sorted_replace - // ... - return this->m_edges.back(); + assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); + return *insert_sorted_replace( + this->m_edges, + Edge(start_node_idx, end_node_idx, traveltime, path, std::forward(args)...), + [](auto&& e1, auto&& e2) { + return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx + : e1.start_node_idx < e2.start_node_idx; + }); } }; +/** + * @brief Sets the graph nodes for counties or districts. + * Reads the node ids which could refer to districts or counties and the epidemiological + * data from json files and creates one node for each id. Every node contains a model. + * @param[in] params Model Parameters that are used for every node. + * @param[in] start_date Start date for which the data should be read. + * @param[in] end_data End date for which the data should be read. + * @param[in] data_dir Directory that contains the data files. + * @param[in] population_data_path Path to json file containing the population data. + * @param[in] is_node_for_county Specifies whether the node ids should be county ids (true) or district ids (false). + * @param[in, out] params_graph Graph whose nodes are set by the function. + * @param[in] read_func Function that reads input data for german counties and sets Model compartments. + * @param[in] node_func Function that returns the county ids. + * @param[in] scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county. + * @param[in] scaling_factor_icu Factor of ICU cases to account for underreporting. + * @param[in] tnt_capacity_factor Factor for test and trace capacity. + * @param[in] num_days Number of days to be simulated; required to load data for vaccinations during the simulation. + * @param[in] export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files. + */ +template +IOResult +set_nodes(const Parameters& params, Date start_date, Date end_date, const fs::path& data_dir, + const std::string& population_data_path, bool is_node_for_county, Graph& params_graph, + ReadFunction&& read_func, NodeIdFunction&& node_func, const std::vector& scaling_factor_inf, + double scaling_factor_icu, double tnt_capacity_factor, int num_days = 0, bool export_time_series = false) +{ + auto scaling_factor_infected = std::vector(size_t(params.get_num_groups()), 1.0); + + auto read_duration = mio::read_duration_stay(durations_path); + if (!read_duration) { + std::cout << read_duration.error().formatted_message() << '\n'; + } + auto duration_stay = read_duration.value(); + + std::string path = "/localdata1/test/memilio/data/pydata/Germany/county_population.json"; + auto read_ids = node_func(path, true); + auto node_ids = read_ids.value(); + + std::vector nodes(node_ids.size(), Model(int(size_t(params.get_num_groups())))); + + for (auto& node : nodes) { + node.parameters = params; + } + auto read_node = read_input_data_county(nodes, start_date, node_ids, scaling_factor_inf, scaling_factor_icu, + "/localdata1/test/memilio/data", num_days); + + for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { + + auto tnt_capacity = nodes[node_idx].populations.get_total() * tnt_capacity_factor; + + //local parameters + auto& tnt_value = nodes[node_idx].parameters.template get(); + tnt_value = mio::UncertainValue(0.5 * (1.2 * tnt_capacity + 0.8 * tnt_capacity)); + tnt_value.set_distribution(mio::ParameterDistributionUniform(0.8 * tnt_capacity, 1.2 * tnt_capacity)); + + //holiday periods + auto id = int(mio::regions::CountyId(node_ids[node_idx])); + auto holiday_periods = mio::regions::get_holidays(mio::regions::get_state_id(id), start_date, end_date); + auto& contacts = nodes[node_idx].parameters.template get(); + contacts.get_school_holidays() = + std::vector>(holiday_periods.size()); + std::transform( + holiday_periods.begin(), holiday_periods.end(), contacts.get_school_holidays().begin(), [=](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))); + }); + + //uncertainty in populations + for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + for (auto j = mio::Index(0); + j < mio::osecirvvs::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))); + } + } + + // Add mobility node + auto mobility = nodes[node_idx]; + mobility.populations.set_total(0); + + // reduce transmission on contact due to mask obligation in mobility node + // first age group not able to (properly) wear masks + if (masks) { + + const double fact_surgical_mask = 0.1; + const double fact_ffp2 = 0.001; + + double factor_mask[] = {1, fact_surgical_mask, fact_ffp2, fact_ffp2, fact_ffp2, fact_ffp2}; + if (!ffp2) { + std::cout << "surgical masks" + << "\n"; + // multipliziere alle außer den ersten EIntrag mit 40 + for (size_t j = 1; j < 6; j++) { + factor_mask[j] = factor_mask[j] * fact_surgical_mask / fact_ffp2; + } + } + + double fac_variant = 1.4; //1.94; //https://doi.org/10.7554/eLife.78933 + double transmissionProbabilityOnContactMin[] = {0.02 * fac_variant, 0.05 * fac_variant, 0.05 * fac_variant, + 0.05 * fac_variant, 0.08 * fac_variant, 0.1 * fac_variant}; + + double transmissionProbabilityOnContactMax[] = {0.04 * fac_variant, 0.07 * fac_variant, 0.07 * fac_variant, + 0.07 * fac_variant, 0.10 * fac_variant, 0.15 * fac_variant}; + for (int i = 0; i < 6; i++) { + transmissionProbabilityOnContactMin[i] = transmissionProbabilityOnContactMin[i] * factor_mask[i]; + transmissionProbabilityOnContactMax[i] = transmissionProbabilityOnContactMax[i] * factor_mask[i]; + } + array_assign_uniform_distribution( + mobility.parameters.get(), + transmissionProbabilityOnContactMin, transmissionProbabilityOnContactMax); + } + + for (size_t t_idx = 0; t_idx < num_days; ++t_idx) { + auto t = mio::SimulationDay((size_t)t_idx); + for (auto j = mio::AgeGroup(0); j < params.get_num_groups(); j++) { + mobility.parameters.template get()[{j, t}] = 0; + mobility.parameters.template get()[{j, t}] = 0; + } + } + + auto& params_mobility = mobility.parameters; + + params_graph.add_node(node_ids[node_idx], duration_stay((Eigen::Index)node_idx), nodes[node_idx], mobility); + } + return success(); +} + class MigrationEdgeDetailed : public MigrationEdge { public: From d000cbab4b1e0b782efc9d7a5d4fb0b12650b198 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Fri, 8 Dec 2023 15:00:52 +0100 Subject: [PATCH 08/11] some simple tests and some fixes --- cpp/memilio/mobility/graph.h | 2 +- .../metapopulation_mobility_detailed.h | 158 ++++++++++-------- cpp/tests/CMakeLists.txt | 17 +- cpp/tests/test_mobility_detailed.cpp | 47 ++++++ 4 files changed, 145 insertions(+), 79 deletions(-) create mode 100644 cpp/tests/test_mobility_detailed.cpp diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index fefc864076..a421b0ca6e 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -235,7 +235,7 @@ class Graph })); } -protected: +private: std::vector> m_nodes; std::vector> m_edges; }; // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h index a5c2cd1b51..15a6a93c03 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_detailed.h +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -20,6 +20,7 @@ #ifndef METAPOPULATION_MOBILITY_DETAILED_H #define METAPOPULATION_MOBILITY_DETAILED_H +#include "memilio/epidemiology/simulation_day.h" #include "memilio/mobility/graph_simulation.h" #include "memilio/mobility/metapopulation_mobility_instant.h" #include "memilio/utils/time_series.h" @@ -36,8 +37,10 @@ #include "memilio/io/mobility_io.h" #include "boost/filesystem.hpp" +#include "models/ode_secirvvs/model.h" #include +#include namespace mio { @@ -115,26 +118,19 @@ class GraphDetailed : public Graph public: using Graph::Graph; - template - NodeDetailed& add_node(int id, double duration_stay, Args&&... args) - { - this->m_nodes.emplace_back(id, duration_stay, std::forward(args)...); - return this->m_nodes.back(); - } - template NodeDetailed& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2) { - this->m_nodes.emplace_back(id, duration_stay, model1, model2); - return this->m_nodes.back(); + m_nodes.emplace_back(id, duration_stay, model1, model2); + return m_nodes.back(); } template NodeDetailed& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2, double m_t0, double m_dt_integration) { - this->m_nodes.emplace_back(id, duration_stay, model1, model2, m_t0, m_dt_integration); - return this->m_nodes.back(); + m_nodes.emplace_back(id, duration_stay, model1, model2, m_t0, m_dt_integration); + return m_nodes.back(); } template @@ -142,7 +138,7 @@ class GraphDetailed : public Graph { assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); return *insert_sorted_replace( - this->m_edges, Edge(start_node_idx, end_node_idx, traveltime, std::forward(args)...), + m_edges, EdgeDetailed(start_node_idx, end_node_idx, traveltime, std::forward(args)...), [](auto&& e1, auto&& e2) { return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx : e1.start_node_idx < e2.start_node_idx; @@ -155,13 +151,17 @@ class GraphDetailed : public Graph { assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); return *insert_sorted_replace( - this->m_edges, - Edge(start_node_idx, end_node_idx, traveltime, path, std::forward(args)...), + m_edges, + EdgeDetailed(start_node_idx, end_node_idx, traveltime, path, std::forward(args)...), [](auto&& e1, auto&& e2) { return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx : e1.start_node_idx < e2.start_node_idx; }); } + +private: + std::vector> m_nodes; + std::vector> m_edges; }; /** @@ -173,6 +173,7 @@ class GraphDetailed : public Graph * @param[in] end_data End date for which the data should be read. * @param[in] data_dir Directory that contains the data files. * @param[in] population_data_path Path to json file containing the population data. + * @param[in] stay_times_data_path Path to txt file containing the stay times for the considered local entities. * @param[in] is_node_for_county Specifies whether the node ids should be county ids (true) or district ids (false). * @param[in, out] params_graph Graph whose nodes are set by the function. * @param[in] read_func Function that reads input data for german counties and sets Model compartments. @@ -185,31 +186,24 @@ class GraphDetailed : public Graph */ template -IOResult -set_nodes(const Parameters& params, Date start_date, Date end_date, const fs::path& data_dir, - const std::string& population_data_path, bool is_node_for_county, Graph& params_graph, - ReadFunction&& read_func, NodeIdFunction&& node_func, const std::vector& scaling_factor_inf, - double scaling_factor_icu, double tnt_capacity_factor, int num_days = 0, bool export_time_series = false) +IOResult set_nodes(const Parameters& params, Date start_date, 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, Graph& params_graph, ReadFunction&& read_func, + NodeIdFunction&& node_func, const std::vector& scaling_factor_inf, + double scaling_factor_icu, double tnt_capacity_factor, int num_days = 0, + bool export_time_series = false) { - auto scaling_factor_infected = std::vector(size_t(params.get_num_groups()), 1.0); - - auto read_duration = mio::read_duration_stay(durations_path); - if (!read_duration) { - std::cout << read_duration.error().formatted_message() << '\n'; - } - auto duration_stay = read_duration.value(); - std::string path = "/localdata1/test/memilio/data/pydata/Germany/county_population.json"; - auto read_ids = node_func(path, true); - auto node_ids = read_ids.value(); + BOOST_OUTCOME_TRY(duration_stay, mio::read_duration_stay(stay_times_data_path)); + BOOST_OUTCOME_TRY(node_ids, node_func(population_data_path, is_node_for_county)); std::vector nodes(node_ids.size(), Model(int(size_t(params.get_num_groups())))); for (auto& node : nodes) { node.parameters = params; } - auto read_node = read_input_data_county(nodes, start_date, node_ids, scaling_factor_inf, scaling_factor_icu, - "/localdata1/test/memilio/data", num_days); + BOOST_OUTCOME_TRY(read_func(nodes, start_date, node_ids, scaling_factor_inf, scaling_factor_icu, data_dir.string(), + num_days, export_time_series)); for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { @@ -223,7 +217,7 @@ set_nodes(const Parameters& params, Date start_date, Date end_date, const fs::pa //holiday periods auto id = int(mio::regions::CountyId(node_ids[node_idx])); auto holiday_periods = mio::regions::get_holidays(mio::regions::get_state_id(id), start_date, end_date); - auto& contacts = nodes[node_idx].parameters.template get(); + auto& contacts = nodes[node_idx].parameters.template get(); contacts.get_school_holidays() = std::vector>(holiday_periods.size()); std::transform( @@ -234,8 +228,7 @@ set_nodes(const Parameters& params, Date start_date, Date end_date, const fs::pa //uncertainty in populations for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { - for (auto j = mio::Index(0); - j < mio::osecirvvs::Model::Compartments::Count; ++j) { + for (auto j = mio::Index(0); j < 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))); @@ -248,50 +241,75 @@ set_nodes(const Parameters& params, Date start_date, Date end_date, const fs::pa auto mobility = nodes[node_idx]; mobility.populations.set_total(0); - // reduce transmission on contact due to mask obligation in mobility node - // first age group not able to (properly) wear masks - if (masks) { - - const double fact_surgical_mask = 0.1; - const double fact_ffp2 = 0.001; + params_graph.add_node(node_ids[node_idx], duration_stay((Eigen::Index)node_idx), nodes[node_idx], mobility); + } + return success(); +} - double factor_mask[] = {1, fact_surgical_mask, fact_ffp2, fact_ffp2, fact_ffp2, fact_ffp2}; - if (!ffp2) { - std::cout << "surgical masks" - << "\n"; - // multipliziere alle außer den ersten EIntrag mit 40 - for (size_t j = 1; j < 6; j++) { - factor_mask[j] = factor_mask[j] * fact_surgical_mask / fact_ffp2; +/** + * @brief Sets the graph edges. + * Reads the commuting matrices, travel times and paths from data and creates one edge for each pair of nodes. + * @param[in] travel_times_path Path to txt file containing the travel times between counties. + * @param[in] mobility_data_path Path to txt file containing the commuting matrices. + * @param[in] travelpath_path Path to txt file containing the paths between counties. + * @param[in, out] params_graph Graph whose nodes are set by the function. + * @param[in] migrating_compartments Compartments that commute. + * @param[in] contact_locations_size Number of contact locations. + * @param[in] commuting_weights Vector with a commuting weight for every AgeGroup. + */ +template +IOResult set_edges(const std::string& travel_times_path, const std::string mobility_data_path, + const std::string& travelpath_path, Graph& params_graph, + std::initializer_list& migrating_compartments, size_t contact_locations_size, + std::vector commuting_weights = std::vector{}) +{ + ScalarType theshold_edges = 4e-5; + + BOOST_OUTCOME_TRY(mobility_data_commuter, mio::read_mobility_plain(mobility_data_path)); + BOOST_OUTCOME_TRY(travel_times, mio::read_mobility_plain(travel_times_path)); + BOOST_OUTCOME_TRY(path_mobility, mio::read_path_mobility(travelpath_path)); + + 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.populations; + + // mobility coefficients have the same number of components as the contact matrices. + // so that the same NPIs/dampings can be used for both (e.g. more home office => fewer commuters) + auto mobility_coeffs = MigrationCoefficientGroup(contact_locations_size, populations.numel()); + auto num_age_groups = (size_t)params_graph.nodes()[county_idx_i].property.parameters.get_num_groups(); + commuting_weights = + (commuting_weights.size() == 0 ? std::vector(num_age_groups, 1.0) : commuting_weights); + + //commuters + auto working_population = 0.0; + auto min_commuter_age = mio::AgeGroup(2); + auto max_commuter_age = mio::AgeGroup(4); //this group is partially retired, only partially commutes + for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { + working_population += populations.get_group_total(age) * commuting_weights[size_t(age)]; + } + auto commuter_coeff_ij = mobility_data_commuter(county_idx_i, county_idx_j) / working_population; + for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { + for (auto compartment : migrating_compartments) { + auto coeff_index = populations.get_flat_index({age, compartment}); + mobility_coeffs[size_t(ContactLocation::Work)].get_baseline()[coeff_index] = + commuter_coeff_ij * commuting_weights[size_t(age)]; } } - double fac_variant = 1.4; //1.94; //https://doi.org/10.7554/eLife.78933 - double transmissionProbabilityOnContactMin[] = {0.02 * fac_variant, 0.05 * fac_variant, 0.05 * fac_variant, - 0.05 * fac_variant, 0.08 * fac_variant, 0.1 * fac_variant}; - - double transmissionProbabilityOnContactMax[] = {0.04 * fac_variant, 0.07 * fac_variant, 0.07 * fac_variant, - 0.07 * fac_variant, 0.10 * fac_variant, 0.15 * fac_variant}; - for (int i = 0; i < 6; i++) { - transmissionProbabilityOnContactMin[i] = transmissionProbabilityOnContactMin[i] * factor_mask[i]; - transmissionProbabilityOnContactMax[i] = transmissionProbabilityOnContactMax[i] * factor_mask[i]; - } - array_assign_uniform_distribution( - mobility.parameters.get(), - transmissionProbabilityOnContactMin, transmissionProbabilityOnContactMax); - } + auto path = path_mobility[county_idx_i][county_idx_j]; + if (static_cast(path[0]) != county_idx_i || + static_cast(path[path.size() - 1]) != county_idx_j) + std::cout << "Wrong Path for edge " << county_idx_i << " " << county_idx_j << "\n"; - for (size_t t_idx = 0; t_idx < num_days; ++t_idx) { - auto t = mio::SimulationDay((size_t)t_idx); - for (auto j = mio::AgeGroup(0); j < params.get_num_groups(); j++) { - mobility.parameters.template get()[{j, t}] = 0; - mobility.parameters.template get()[{j, t}] = 0; + //only add edges with mobility above thresholds for performance + if (commuter_coeff_ij > theshold_edges) { + params_graph.add_edge(county_idx_i, county_idx_j, travel_times(county_idx_i, county_idx_j), + path_mobility[county_idx_i][county_idx_j], std::move(mobility_coeffs)); } } - - auto& params_mobility = mobility.parameters; - - params_graph.add_node(node_ids[node_idx], duration_stay((Eigen::Index)node_idx), nodes[node_idx], mobility); } + return success(); } diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index c19b9fcec2..f31153406b 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -49,6 +49,7 @@ set(TESTSOURCES test_io_framework.cpp test_binary_serializer.cpp test_compartmentsimulation.cpp + test_mobility_detailed.cpp test_mobility_io.cpp test_transform_iterator.cpp test_metaprogramming.cpp @@ -68,17 +69,17 @@ set(TESTSOURCES ) if(MEMILIO_HAS_JSONCPP) -set(TESTSOURCES ${TESTSOURCES} -test_json_serializer.cpp -test_epi_data_io.cpp -) + set(TESTSOURCES ${TESTSOURCES} + test_json_serializer.cpp + test_epi_data_io.cpp + ) endif() if(MEMILIO_HAS_JSONCPP AND MEMILIO_HAS_HDF5) -set(TESTSOURCES ${TESTSOURCES} -test_save_parameters.cpp -test_save_results.cpp -) + set(TESTSOURCES ${TESTSOURCES} + test_save_parameters.cpp + test_save_results.cpp + ) endif() add_executable(memilio-test ${TESTSOURCES}) diff --git a/cpp/tests/test_mobility_detailed.cpp b/cpp/tests/test_mobility_detailed.cpp new file mode 100644 index 0000000000..b7fa254fa9 --- /dev/null +++ b/cpp/tests/test_mobility_detailed.cpp @@ -0,0 +1,47 @@ + +#include "gtest/gtest.h" +#include "gmock/gmock.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" +#include "memilio/compartments/simulation.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "ode_secir/model.h" +#include "ode_secir/parameters.h" + +namespace mio +{ +mio::IOResult read_duration_stay(const std::string& filename) +{ + mio::unused(filename); + return mio::success(Eigen::MatrixXd(0, 0)); +} + +TEST(TestMobilityDetailed, AddNodeTest) +{ + mio::osecir::Model model(1); + GraphDetailed>, int> sim_graph; + + auto& node = sim_graph.add_node(1, 0.5, model, model, 0.0, 0.1); + EXPECT_EQ(node.stay_duration, 0.5); +} + +TEST(TestMobilityDetailed, AddNodeWithTimeIntegrationTest) +{ + mio::osecir::Model model(1); + GraphDetailed model_graph; + auto& node = model_graph.add_node(1, 0.5, model, model); + EXPECT_EQ(node.stay_duration, 0.5); +} + +TEST(TestMobilityDetailed, AddEdgeTest) +{ + mio::osecir::Model model(1); + GraphDetailed model_graph; + model_graph.add_node(0, 0.5, model, model); + model_graph.add_node(1, 0.5, model, model); + auto& edge = model_graph.add_edge(0, 1, 0.5, {1, 2}); + EXPECT_EQ(edge.traveltime, 0.5); +} + +} // namespace mio \ No newline at end of file From 6c5e4350cea62f676ee24c3bf580f66eb11938bf Mon Sep 17 00:00:00 2001 From: HenrZu Date: Mon, 11 Dec 2023 08:18:58 +0100 Subject: [PATCH 09/11] change includes --- cpp/memilio/CMakeLists.txt | 2 ++ cpp/memilio/io/mobility_io.h | 6 ++-- .../metapopulation_mobility_detailed.h | 31 ++++++++++++++++--- 3 files changed, 32 insertions(+), 7 deletions(-) diff --git a/cpp/memilio/CMakeLists.txt b/cpp/memilio/CMakeLists.txt index be65a209c4..775da2bf7f 100644 --- a/cpp/memilio/CMakeLists.txt +++ b/cpp/memilio/CMakeLists.txt @@ -59,6 +59,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/io/mobility_io.h b/cpp/memilio/io/mobility_io.h index 44a04e8fb9..90df68aaa3 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/config.h" #include "memilio/math/eigen.h" @@ -209,4 +209,4 @@ IOResult> read_graph(const std::string& direct } // namespace mio -#endif // READ_TWITTER_H +#endif // MIO_MOBILITY_IO_H diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h index 15a6a93c03..f2cf89e331 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_detailed.h +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -37,7 +37,6 @@ #include "memilio/io/mobility_io.h" #include "boost/filesystem.hpp" -#include "models/ode_secirvvs/model.h" #include #include @@ -262,10 +261,9 @@ template set_edges(const std::string& travel_times_path, const std::string mobility_data_path, const std::string& travelpath_path, Graph& params_graph, std::initializer_list& migrating_compartments, size_t contact_locations_size, - std::vector commuting_weights = std::vector{}) + std::vector commuting_weights = std::vector{}, + ScalarType theshold_edges = 4e-5) { - ScalarType theshold_edges = 4e-5; - BOOST_OUTCOME_TRY(mobility_data_commuter, mio::read_mobility_plain(mobility_data_path)); BOOST_OUTCOME_TRY(travel_times, mio::read_mobility_plain(travel_times_path)); BOOST_OUTCOME_TRY(path_mobility, mio::read_path_mobility(travelpath_path)); @@ -443,6 +441,31 @@ void apply_migration(double t, double dt, MigrationEdge& migrationEdge, Simulati { migrationEdge.apply_migration(t, dt, node_from, node_to, mode); } + +/** + * create a migration simulation. + * After every second time step, for each edge a portion of the population corresponding to the coefficients of the edge + * moves from one node to the other. In the next timestep, the migrated population return to their "home" node. + * Returns are adjusted based on the development in the target node. + * @param t0 start time of the simulation + * @param dt time step between migrations + * @param graph set up for migration simulation + * @{ + */ +template +GraphSimulationDetailed, MigrationEdge>> +make_migration_sim(double t0, double dt, const Graph, MigrationEdge>& graph) +{ + return make_graph_sim(t0, dt, graph, &evolve_model, &apply_migration); +} + +template +GraphSimulationDetailed, MigrationEdge>> +make_migration_sim(double t0, double dt, Graph, MigrationEdge>&& graph) +{ + return make_graph_sim(t0, dt, std::move(graph), &evolve_model, &apply_migration); +} + } // namespace mio #endif //METAPOPULATION_MOBILITY_DETAILED_H From 1d99402d9377d03485d3104761b35f340c7da304 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Mon, 11 Dec 2023 14:09:41 +0100 Subject: [PATCH 10/11] move to member functions --- cpp/examples/CMakeLists.txt | 5 +- .../ode_secirvvs_mobility_detailed.cpp | 395 ++---------------- cpp/memilio/mobility/graph_simulation.h | 271 ++++++------ .../metapopulation_mobility_detailed.h | 71 +++- 4 files changed, 240 insertions(+), 502 deletions(-) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index ca5069e1d4..e996d0a39c 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -21,8 +21,9 @@ target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNI add_executable(ode_sir_example ode_sir.cpp) target_link_libraries(ode_sir_example PRIVATE memilio ode_sir) -# add_executable(sim_mobility_detailed ode_secirvvs_mobility_detailed.cpp) -# target_link_libraries(sim_mobility_detailed PRIVATE memilio ode_secirvvs) +add_executable(sim_mobility_detailed ode_secirvvs_mobility_detailed.cpp) +target_link_libraries(sim_mobility_detailed PRIVATE memilio ode_secirvvs) + add_executable(seir_flows_example ode_seir_flows.cpp) target_link_libraries(seir_flows_example PRIVATE memilio ode_seir) diff --git a/cpp/examples/ode_secirvvs_mobility_detailed.cpp b/cpp/examples/ode_secirvvs_mobility_detailed.cpp index 42ce21ccd7..752a8a12e3 100644 --- a/cpp/examples/ode_secirvvs_mobility_detailed.cpp +++ b/cpp/examples/ode_secirvvs_mobility_detailed.cpp @@ -371,20 +371,17 @@ void init_pop_cologne_szenario(mio::Graph> -get_graph(const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool edges) +get_graph(const std::string data_dir, const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool edges) { - std::string data_dir = "/localdata1/test/memilio//data"; - std::string traveltimes_path = "/localdata1/code/memilio/test/travel_times_pathes.txt"; - std::string durations_path = "/localdata1/code/memilio/test/activity_duration_work.txt"; auto mobility_data_commuter = mio::read_mobility_plain(("/localdata1/code/memilio/test/commuter_migration_with_locals.txt")); auto mob_data = mobility_data_commuter.value(); - auto start_date = mio::Date(2022, 8, 1); - auto end_date = mio::Date(2022, 11, 1); + auto start_date = mio::Date(2021, 8, 1); + auto end_date = mio::Date(2021, 11, 1); // global parameters const int num_age_groups = 6; - mio::Graph params_graph; + mio::GraphDetailed params_graph; // Nodes mio::osecirvvs::Parameters params(num_age_groups); @@ -400,65 +397,21 @@ get_graph(const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool auto scaling_factor_icu = 1.0; auto tnt_capacity_factor = 1.43 / 100000.; - auto read_duration = mio::read_duration_stay(durations_path); - if (!read_duration) { - std::cout << read_duration.error().formatted_message() << '\n'; - } - auto duration_stay = read_duration.value(); - - std::string path = "/localdata1/test/memilio/data/pydata/Germany/county_population.json"; - auto read_ids = mio::get_node_ids(path, true); - auto node_ids = read_ids.value(); - - std::vector nodes(node_ids.size(), - mio::osecirvvs::Model(int(size_t(params.get_num_groups())))); - std::vector scaling_factor_inf(6, 1.0); - - for (auto& node : nodes) { - node.parameters = params; - } - auto read_node = read_input_data_county(nodes, start_date, node_ids, scaling_factor_inf, scaling_factor_icu, - "/localdata1/test/memilio/data", num_days); - - for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { - - auto tnt_capacity = nodes[node_idx].populations.get_total() * tnt_capacity_factor; - - //local parameters - auto& tnt_value = nodes[node_idx].parameters.template get(); - tnt_value = mio::UncertainValue(0.5 * (1.2 * tnt_capacity + 0.8 * tnt_capacity)); - tnt_value.set_distribution(mio::ParameterDistributionUniform(0.8 * tnt_capacity, 1.2 * tnt_capacity)); - - //holiday periods - auto id = int(mio::regions::CountyId(node_ids[node_idx])); - auto holiday_periods = mio::regions::get_holidays(mio::regions::get_state_id(id), start_date, end_date); - auto& contacts = nodes[node_idx].parameters.template get(); - contacts.get_school_holidays() = - std::vector>(holiday_periods.size()); - std::transform( - holiday_periods.begin(), holiday_periods.end(), contacts.get_school_holidays().begin(), [=](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))); - }); - - //uncertainty in populations - for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { - for (auto j = mio::Index(0); - j < mio::osecirvvs::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))); - } - } - - // Add mobility node - auto mobility = nodes[node_idx]; - mobility.populations.set_total(0); - - // reduce transmission on contact due to mask obligation in mobility node - // first age group not able to (properly) wear masks + const auto& read_function_nodes = mio::osecirvvs::read_input_data_county; + const auto& read_function_edges = mio::read_mobility_plain; + const auto& node_id_function = mio::get_node_ids; + + const auto& set_node_function = + mio::set_nodes_detailed; + BOOST_OUTCOME_TRY(set_node_function( + params, start_date, end_date, data_dir, data_dir + "//pydata//Germany//county_current_population.json", + mio::path_join(data_dir, "mobility", "activity_duration_work.txt"), true, params_graph, read_function_nodes, + node_id_function, scaling_factor_infected, scaling_factor_icu, tnt_capacity_factor, + mio::get_offset_in_days(end_date, start_date), false)); + + for (auto& node : params_graph.nodes()) { if (masks) { const double fact_surgical_mask = 0.1; @@ -474,7 +427,7 @@ get_graph(const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool } } - double fac_variant = 1.4; //1.94; //https://doi.org/10.7554/eLife.78933 + double fac_variant = 1.94; //https://doi.org/10.7554/eLife.78933 double transmissionProbabilityOnContactMin[] = {0.02 * fac_variant, 0.05 * fac_variant, 0.05 * fac_variant, 0.05 * fac_variant, 0.08 * fac_variant, 0.1 * fac_variant}; @@ -485,21 +438,17 @@ get_graph(const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool transmissionProbabilityOnContactMax[i] = transmissionProbabilityOnContactMax[i] * factor_mask[i]; } array_assign_uniform_distribution( - mobility.parameters.get(), + node.mobility.parameters.get(), transmissionProbabilityOnContactMin, transmissionProbabilityOnContactMax); } for (size_t t_idx = 0; t_idx < num_days; ++t_idx) { auto t = mio::SimulationDay((size_t)t_idx); for (auto j = mio::AgeGroup(0); j < params.get_num_groups(); j++) { - mobility.parameters.template get()[{j, t}] = 0; - mobility.parameters.template get()[{j, t}] = 0; + node.mobility.parameters.template get()[{j, t}] = 0; + node.mobility.parameters.template get()[{j, t}] = 0; } } - - auto& params_mobility = mobility.parameters; - - params_graph.add_node(node_ids[node_idx], duration_stay((Eigen::Index)node_idx), nodes[node_idx], mobility); } if (szenario_cologne) { @@ -508,7 +457,6 @@ get_graph(const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool // Edges if (edges) { - ScalarType theshold_edges = 4e-5; auto migrating_compartments = { mio::osecirvvs::InfectionState::SusceptibleNaive, mio::osecirvvs::InfectionState::ExposedNaive, @@ -524,49 +472,15 @@ get_graph(const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity, }; - // mobility between nodes - auto read_travel_times = mio::read_mobility_plain(traveltimes_path); - auto travel_times = read_travel_times.value(); - auto read_paths_mobility = mio::read_path_mobility( - "/localdata1/code/memilio/test/wegketten_ohne_komma.txt"); // wegketten_ohne_komma.txt - auto path_mobility = read_paths_mobility.value(); - - 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) { - // if (county_idx_i == county_idx_j) - // continue; - auto& populations = nodes[county_idx_i].populations; - // mobility coefficients have the same number of components as the contact matrices. - // so that the same NPIs/dampings can be used for both (e.g. more home office => fewer commuters) - auto mobility_coeffs = mio::MigrationCoefficientGroup(contact_locations.size(), populations.numel()); - - //commuters - auto working_population = 0.0; - auto min_commuter_age = mio::AgeGroup(2); - auto max_commuter_age = mio::AgeGroup(4); //this group is partially retired, only partially commutes - for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { - working_population += populations.get_group_total(age) * (age == max_commuter_age ? 0.33 : 1.0); - } - auto commuter_coeff_ij = mob_data(county_idx_i, county_idx_j) / - working_population; //data is absolute numbers, we need relative - for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { - for (auto compartment : migrating_compartments) { - auto coeff_index = populations.get_flat_index({age, compartment}); - mobility_coeffs[size_t(ContactLocation::Work)].get_baseline()[coeff_index] = - commuter_coeff_ij * (age == max_commuter_age ? 0.33 : 1.0); - } - } - - auto path = path_mobility[county_idx_i][county_idx_j]; - if (path[0] != county_idx_i || path[path.size() - 1] != county_idx_j) - std::cout << "falsch" - << "\n"; - if (commuter_coeff_ij > theshold_edges) { - params_graph.add_edge(county_idx_i, county_idx_j, travel_times(county_idx_i, county_idx_j), - path_mobility[county_idx_i][county_idx_j], std::move(mobility_coeffs)); - } - } - } + const auto& read_function_edges = mio::read_mobility_plain; + const auto& set_edge_function = + mio::set_edges_detailed; + BOOST_OUTCOME_TRY(set_edge_function(mio::path_join(data_dir, "mobility", "travel_times_pathes.txt"), + data_dir + "//mobility//commuter_migration_with_locals.txt", + mio::path_join(data_dir, "mobility", "wegketten_ohne_komma.txt"), + params_graph, migrating_compartments, contact_locations.size(), + std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}, false)); // average travel time const ScalarType avg_traveltime = std::accumulate(params_graph.edges().begin(), params_graph.edges().end(), 0.0, @@ -634,21 +548,22 @@ mio::IOResult run() if (!masks && ffp2) { mio::log_error("ffp2 only possible with masks"); } + std::string data_dir = "/localdata1/test/memilio/data"; bool szenario_cologne = false; // auto params_graph = get_graph(num_days, masks); - BOOST_OUTCOME_TRY(created, get_graph(num_days, masks, ffp2, szenario_cologne, edges)); + BOOST_OUTCOME_TRY(created, get_graph(data_dir, num_days, masks, ffp2, szenario_cologne, edges)); auto params_graph = created; std::string res_dir = "/localdata1/code/memilio/results_paper/mask_" + std::to_string(masks); - // if (ffp2) - // res_dir = res_dir + "_ffp2"; - // if (szenario_cologne) - // res_dir = res_dir + "_cologne"; - // if (!edges) - // res_dir = res_dir + "_no_edges"; + if (ffp2) + res_dir = res_dir + "_ffp2"; + if (szenario_cologne) + res_dir = res_dir + "_cologne"; + if (!edges) + res_dir = res_dir + "_no_edges"; if (!edges && (ffp2 || masks)) mio::log_error("no edges only possible without masks and ffp2"); @@ -659,7 +574,6 @@ mio::IOResult run() if (!fs::exists(res_dir)) { fs::create_directory(res_dir); } - auto write_graph_status = write_graph(params_graph, "/localdata1/code/memilio/save_graph"); std::vector county_ids(params_graph.nodes().size()); std::transform(params_graph.nodes().begin(), params_graph.nodes().end(), county_ids.begin(), [](auto& n) { @@ -683,235 +597,8 @@ mio::IOResult run() }, [&](auto results_graph, auto&& run_idx) { auto interpolated_result = mio::interpolate_simulation_result(results_graph); + auto params = std::vector(); - auto params = std::vector(); - // params.reserve(results_graph.nodes().size()); - // std::transform(results_graph.nodes().begin(), results_graph.nodes().end(), std::back_inserter(params), - // [](auto&& node) { - // return node.property.get_simulation().get_model(); - // }); - // const auto size_vec = interpolated_result[0].get_num_time_points(); - // std::vector num_transmissions_in_mobility_nodes(size_vec, 0.0); - // std::vector> num_transmission_per_county(results_graph.nodes().size(), - // std::vector(size_vec, 0.0)); - // std::vector> num_symptomatic_per_county(results_graph.nodes().size(), - // std::vector(size_vec, 0.0)); - - // // vaccinations - // // std::vector> vacc_county_N(results_graph.nodes().size()); - // // std::vector> vacc_county_PI(results_graph.nodes().size()); - // // std::vector> vacc_county_II(results_graph.nodes().size()); - - // std::vector waning_partial_immunity(size_vec, 0.0); - // std::vector waning_improved_immunity(size_vec, 0.0); - // auto node_idx = 0; - // for (auto& node : results_graph.nodes()) { - // auto flows = mio::interpolate_simulation_result(node.mobility.get_simulation().get_flows()); - // auto flows_local_model = mio::interpolate_simulation_result(node.property.get_simulation().get_flows()); - // for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(6); ++i) { - // auto& model = node.mobility.get_simulation().get_model(); - // auto flow_index_N = - // model.template get_flat_flow_index({i}); - // auto flow_index_N_confirmed = model.template get_flat_flow_index< - // mio::osecirvvs::InfectionState::InfectedNoSymptomsNaiveConfirmed, - // mio::osecirvvs::InfectionState::InfectedSymptomsNaiveConfirmed>({i}); - // auto flow_index_PI = model.template get_flat_flow_index< - // mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity, - // mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity>({i}); - // auto flow_index_PI_confirmed = model.template get_flat_flow_index< - // mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunityConfirmed, - // mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunityConfirmed>({i}); - // auto flow_index_II = model.template get_flat_flow_index< - // mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity, - // mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity>({i}); - // auto flow_index_II_confirmed = model.template get_flat_flow_index< - // mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed, - // mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunityConfirmed>({i}); - - // // Flows from Susceptible to Exposed - // auto flow_index_SN = - // model.template get_flat_flow_index({i}); - // auto flow_index_SII = - // model.template get_flat_flow_index( - // {i}); - // auto flow_index_SPI = - // model.template get_flat_flow_index({i}); - - // auto flow_index_waning_pi = - // model.template get_flat_flow_index({i}); - // auto flow_index_waning_ii = - // model.template get_flat_flow_index( - // {i}); - // for (Eigen::Index time_idx = 0; time_idx < flows.get_num_time_points(); ++time_idx) { - // const auto& flow_values = flows.get_value(time_idx); - // auto time = flows.get_time(time_idx); - - // auto indx_vector = static_cast(time_idx); - // // num_transmissions_in_mobility_nodes[indx_vector] += - // // flow_values(flow_index_N) + flow_values(flow_index_PI) + flow_values(flow_index_II) + - // // flow_values(flow_index_N_confirmed) + flow_values(flow_index_PI_confirmed) + - // // flow_values(flow_index_II_confirmed); - // num_transmissions_in_mobility_nodes[indx_vector] += - // flow_values(flow_index_SN) + flow_values(flow_index_SPI) + flow_values(flow_index_SII); - // waning_partial_immunity[indx_vector] += flow_values(flow_index_waning_pi); - // waning_improved_immunity[indx_vector] += flow_values(flow_index_waning_ii); - // num_symptomatic_per_county[node_idx][indx_vector] += - // flow_values(flow_index_N) + flow_values(flow_index_PI) + flow_values(flow_index_II) + - // flow_values(flow_index_N_confirmed) + flow_values(flow_index_PI_confirmed) + - // flow_values(flow_index_II_confirmed); - // num_transmission_per_county[node_idx][indx_vector] += - // flow_values(flow_index_SN) + flow_values(flow_index_SPI) + flow_values(flow_index_SII); - // } - // for (Eigen::Index time_idx = 0; time_idx < flows_local_model.get_num_time_points(); ++time_idx) { - // const auto& flow_values_local_model = flows_local_model.get_value(time_idx); - // auto time = flows_local_model.get_time(time_idx); - // auto indx_vector = static_cast(time_idx); - // waning_partial_immunity[indx_vector] += flow_values_local_model(flow_index_waning_pi); - // waning_improved_immunity[indx_vector] += flow_values_local_model(flow_index_waning_ii); - // num_transmission_per_county[node_idx][indx_vector] += flow_values_local_model(flow_index_SN) + - // flow_values_local_model(flow_index_SPI) + - // flow_values_local_model(flow_index_SII); - // num_symptomatic_per_county[node_idx][indx_vector] += - // flow_values_local_model(flow_index_N) + flow_values_local_model(flow_index_PI) + - // flow_values_local_model(flow_index_II) + flow_values_local_model(flow_index_N_confirmed) + - // flow_values_local_model(flow_index_PI_confirmed) + - // flow_values_local_model(flow_index_II_confirmed); - // } - // } - - // // // vaccinations - // // vacc_county_N[node_idx].resize(flows_local_model.get_num_time_points(), 0.0); - // // vacc_county_PI[node_idx].resize(flows_local_model.get_num_time_points(), 0.0); - // // vacc_county_II[node_idx].resize(flows_local_model.get_num_time_points(), 0.0); - // // for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(6); ++i) { - // // auto& model = node.mobility.get_simulation().get_model(); - // // auto vacc_N = - // // model.template get_flat_flow_index({i}); - // // auto flow_index_PI = - // // model.template get_flat_flow_index({i}); - // // auto flow_index_II = - // // model.template get_flat_flow_index({i}); - - // // for (Eigen::Index time_idx = 0; time_idx < flows_local_model.get_num_time_points(); ++time_idx) { - // // const auto& flow_values_local_model = flows_local_model.get_value(time_idx); - // // auto time = flows_local_model.get_time(time_idx); - // // auto indx_vector = static_cast(time_idx); - // // vacc_county_N[node_idx][indx_vector] += flow_values_local_model(vacc_N); - // // vacc_county_PI[node_idx][indx_vector] += flow_values_local_model(flow_index_PI); - // // vacc_county_II[node_idx][indx_vector] += flow_values_local_model(flow_index_II); - // // } - // // } - // node_idx++; - // } - // std::string results_path = "/localdata1/code/memilio/results_paper/mask_" + std::to_string(masks); - // if (ffp2) - // results_path = results_path + "_ffp2"; - - // if (szenario_cologne) - // results_path = results_path + "_cologne"; - - // if (!edges) - // results_path = results_path + "_no_edges"; - // const std::string results_path_mb = results_path + "/flows_mb"; - // if (!fs::exists(results_path_mb)) { - // fs::create_directory(results_path_mb); - // } - // const std::string results_filename = results_path_mb + "/results_run_" + std::to_string(run_idx) + ".txt"; - // std::ofstream outfile(results_filename); - - // outfile << "TimeIndex Num_Transmissions_in_Mobility_Nodes Waning_Partial_Immunity " - // "Waning_Improved_Immunity\n"; - // size_t vec_size = num_transmissions_in_mobility_nodes.size(); - // for (size_t i = 0; i < vec_size; ++i) { - // outfile << i << " " << num_transmissions_in_mobility_nodes[i] << " " << waning_partial_immunity[i] - // << " " << waning_improved_immunity[i] << "\n"; - // } - // outfile.close(); - - // // erzeuge ordner flows_county - // const std::string results_path_flows = results_path + "/flows_county"; - // if (!fs::exists(results_path_flows)) { - // fs::create_directory(results_path_flows); - // } - // const std::string results_filename_county = - // results_path_flows + "/transmission_county_run_" + std::to_string(run_idx) + ".txt"; - // std::ofstream outfile_county(results_filename_county); - // for (size_t i = 0; i < results_graph.nodes().size(); ++i) { - // outfile_county << i << " "; - // for (size_t j = 0; j < vec_size; ++j) { - // outfile_county << num_transmission_per_county[i][j] << " "; - // } - // outfile_county << "\n"; - // } - // outfile_county.close(); - - // const std::string results_filename_county_symptomatic = - // results_path_flows + "/symptomatic_county_run_" + std::to_string(run_idx) + ".txt"; - // std::ofstream outfile_county_symptomatic(results_filename_county_symptomatic); - // for (size_t i = 0; i < results_graph.nodes().size(); ++i) { - // outfile_county_symptomatic << i << " "; - // for (size_t j = 0; j < vec_size; ++j) { - // outfile_county_symptomatic << num_symptomatic_per_county[i][j] << " "; - // } - // outfile_county_symptomatic << "\n"; - // } - // outfile_county_symptomatic.close(); - - // // // speichere vacc_county_N, vacc_county_PI, vacc_county_II in jeweils eigene Datei - // // const std::string results_path_vacc = results_path + "/vacc_county"; - // // // check if dir exist, otherwise create them - // // if (!fs::exists(results_path_vacc)) { - // // fs::create_directory(results_path_vacc); - // // } - // // const std::string results_filename_vacc_N = - // // results_path_vacc + "/vacc_county_N_run_" + std::to_string(run_idx) + ".txt"; - // // const std::string results_filename_vacc_PI = - // // results_path_vacc + "/vacc_county_PI_run_" + std::to_string(run_idx) + ".txt"; - // // const std::string results_filename_vacc_II = - // // results_path_vacc + "/vacc_county_II_run_" + std::to_string(run_idx) + ".txt"; - - // // std::ofstream outfile_vacc_N(results_filename_vacc_N); - // // std::ofstream outfile_vacc_PI(results_filename_vacc_PI); - // // std::ofstream outfile_vacc_II(results_filename_vacc_II); - - // // for (size_t i = 0; i < results_graph.nodes().size(); ++i) { - // // outfile_vacc_N << i << " "; - // // outfile_vacc_PI << i << " "; - // // outfile_vacc_II << i << " "; - - // // // Abfrage der Größe des aktuellen inneren Vektors: - // // size_t current_vec_size_N = vacc_county_N[i].size(); - // // size_t current_vec_size_PI = vacc_county_PI[i].size(); - // // size_t current_vec_size_II = vacc_county_II[i].size(); - - // // for (size_t j = 0; j < current_vec_size_N; ++j) { - // // outfile_vacc_N << vacc_county_N[i][j] << " "; - // // } - // // for (size_t j = 0; j < current_vec_size_PI; ++j) { - // // outfile_vacc_PI << vacc_county_PI[i][j] << " "; - // // } - // // for (size_t j = 0; j < current_vec_size_II; ++j) { - // // outfile_vacc_II << vacc_county_II[i][j] << " "; - // // } - - // // outfile_vacc_N << "\n"; - // // outfile_vacc_PI << "\n"; - // // outfile_vacc_II << "\n"; - // // } - // // outfile_vacc_N.close(); - // // outfile_vacc_PI.close(); - // // outfile_vacc_II.close(); - - // std::cout << "Run " << run_idx << " complete." << std::endl; return std::make_pair(interpolated_result, params); }); diff --git a/cpp/memilio/mobility/graph_simulation.h b/cpp/memilio/mobility/graph_simulation.h index 59311bb2b1..9f59d00737 100644 --- a/cpp/memilio/mobility/graph_simulation.h +++ b/cpp/memilio/mobility/graph_simulation.h @@ -23,6 +23,7 @@ #include "memilio/math/euler.h" #include "memilio/mobility/graph.h" #include "memilio/utils/random_number_generator.h" +#include namespace mio { @@ -305,7 +306,6 @@ class GraphSimulationDetailed // calc schedule for each edge precompute_schedule(); - const double epsilon = std::numeric_limits::epsilon(); while (Base::m_t - epsilon < t_max) { // auto start = std::chrono::system_clock::now(); @@ -324,148 +324,22 @@ class GraphSimulationDetailed size_t indx_schedule = 0; while (t_begin + 1 > Base::m_t + 1e-10) { - for (const auto& edge_indx : edges_mobility[indx_schedule]) { - auto& e = Base::m_graph.edges()[edge_indx]; - // first mobility activity - if (indx_schedule == first_mobility[edge_indx]) { - auto& node_from = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; - auto& node_to = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility; - m_edge_func(Base::m_t, 0.0, e.property, node_from, node_to, 0); - } - // next mobility activity - else if (indx_schedule > first_mobility[edge_indx]) { - auto current_node_indx = schedule_edges[edge_indx][indx_schedule]; - bool in_mobility_node = mobility_schedule_edges[edge_indx][indx_schedule]; - - // determine dt, which is equal to the last integration/syncronization point in the current node - auto integrator_schedule_row = ln_int_schedule[current_node_indx]; - if (in_mobility_node) - integrator_schedule_row = mb_int_schedule[current_node_indx]; - // search the index of indx_schedule in the integrator schedule - const size_t indx_current = - std::distance(integrator_schedule_row.begin(), - std::lower_bound(integrator_schedule_row.begin(), - integrator_schedule_row.end(), indx_schedule)); - - if (integrator_schedule_row[indx_current] != indx_schedule) - std::cout << "Error in schedule." - << "\n"; - - ScalarType dt_mobility; - if (indx_current == 0) { - dt_mobility = round_second_decimal(e.traveltime / e.path.size()); - if (dt_mobility < 0.01) - dt_mobility = 0.01; - } - else { - dt_mobility = - round_second_decimal((static_cast(integrator_schedule_row[indx_current]) - - static_cast(integrator_schedule_row[indx_current - 1])) / - 100 + - epsilon); - } - - // We have two cases. Either, we want to send the individuals to the next node, or we just want - // to update their state since a syncronization step is necessary in the current node. - if ((schedule_edges[edge_indx][indx_schedule] != - schedule_edges[edge_indx][indx_schedule - 1]) || - (mobility_schedule_edges[edge_indx][indx_schedule] != - mobility_schedule_edges[edge_indx][indx_schedule - 1])) { - auto& node_from = - mobility_schedule_edges[edge_indx][indx_schedule - 1] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility - : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; - - auto& node_to = - mobility_schedule_edges[edge_indx][indx_schedule] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility - : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; - - if (indx_schedule < mobility_schedule_edges[edge_indx].size() - 1) { - m_edge_func(Base::m_t, dt_mobility, e.property, node_from, node_to, 1); - } - else - // the last time step is handled differently since we have to close the timeseries - m_edge_func(Base::m_t, dt_mobility, e.property, node_from, - Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property, - 2); - } - else { - auto& node_from = - mobility_schedule_edges[edge_indx][indx_schedule - 1] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility - : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; - - auto& node_to = - mobility_schedule_edges[edge_indx][indx_schedule] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility - : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; - - assert(node_from.get_result().get_last_value() == node_to.get_result().get_last_value()); - m_edge_func(Base::m_t, dt_mobility, e.property, node_from, node_to, 3); - } - } - } + advance_edges(indx_schedule); // first we integrate the nodes in time. Afterwards the update on the edges is done. // We start with the edges since the values for t0 are given. - // iterate over all local nodes and integrate them to the syncronization point - for (const auto& n_indx : nodes_mobility[indx_schedule]) { - auto& n = Base::m_graph.nodes()[n_indx]; - const size_t indx_current = - std::distance(ln_int_schedule[n_indx].begin(), - std::lower_bound(ln_int_schedule[n_indx].begin(), ln_int_schedule[n_indx].end(), - indx_schedule)); - - const size_t val_next = (indx_current == ln_int_schedule[n_indx].size() - 1) - ? 100 - : ln_int_schedule[n_indx][indx_current + 1]; - const ScalarType next_dt = - round_second_decimal((static_cast(val_next) - indx_schedule) / 100 + epsilon); - m_node_func(Base::m_t, next_dt, n.property); - } + advance_local_nodes(indx_schedule); + advance_mobility_nodes(indx_schedule); - // iterate over all mobility nodes and integrate them to the syncronization point - for (const size_t& n_indx : nodes_mobility_m[indx_schedule]) { - auto& n = Base::m_graph.nodes()[n_indx]; - const size_t indx_current = - std::distance(mb_int_schedule[n_indx].begin(), - std::lower_bound(mb_int_schedule[n_indx].begin(), mb_int_schedule[n_indx].end(), - indx_schedule)); - const size_t val_next = (indx_current == mb_int_schedule[n_indx].size() - 1) - ? 100 - : mb_int_schedule[n_indx][indx_current + 1]; - const ScalarType next_dt = - round_second_decimal((static_cast(val_next) - indx_schedule) / 100 + epsilon); - - // get all time points from the last integration step - auto& last_time_point = - n.mobility.get_result().get_time(n.mobility.get_result().get_num_time_points() - 1); - // wenn last_time_point nicht innerhalb eines intervalls von 1-e10 von t liegt, dann setzte den letzten Zeitpunkt auf m_t - if (std::fabs(last_time_point - Base::m_t) > 1e-10) { - n.mobility.get_result().get_last_time() = Base::m_t; - } - m_node_func(Base::m_t, next_dt, n.mobility); - Eigen::Index indx_min; - while (n.mobility.get_result().get_last_value().minCoeff(&indx_min) < -1e-7) { - Eigen::Index indx_max; - n.mobility.get_result().get_last_value().maxCoeff(&indx_max); - n.mobility.get_result().get_last_value()[indx_max] -= - n.mobility.get_result().get_last_value()[indx_min]; - n.mobility.get_result().get_last_value()[indx_min] = 0; - } - } indx_schedule++; Base::m_t += min_dt; } - // set each compartment zero for all mobility nodes since we only model daily mobility + // At the end of the day. we set each compartment zero for all mobility nodes since we have to estimate + // the state of the indivuals moving between different counties. + // Therefore there can be differences with the states given by the integrator used for the mobility node. for (auto& n : Base::m_graph.nodes()) { n.mobility.get_result().get_last_value().setZero(); } - // // messe die zeit, wie lange eine iteration bis zu dieser stelle dauert - // auto end = std::chrono::system_clock::now(); - // std::chrono::duration elapsed_seconds = end - start; - // std::cout << "Time (min) needed per Iteration is " << elapsed_seconds.count() / 60 << "min\n"; } } @@ -491,6 +365,8 @@ class GraphSimulationDetailed // first mobility activites per edge std::vector first_mobility; + const double epsilon = std::numeric_limits::epsilon(); + void precompute_schedule() { // print transi @@ -498,8 +374,6 @@ class GraphSimulationDetailed schedule_edges.reserve(Base::m_graph.edges().size()); mobility_schedule_edges.reserve(Base::m_graph.edges().size()); - const double epsilon = std::numeric_limits::epsilon(); - for (auto& e : Base::m_graph.edges()) { // 100 since we round to second decimal std::vector schedule(timesteps, 0.); @@ -695,6 +569,133 @@ class GraphSimulationDetailed nodes_mobility_m.push_back(temp_nodes_mobility_m); } } + + void advance_edges(size_t indx_schedule) + { + for (const auto& edge_indx : edges_mobility[indx_schedule]) { + auto& e = Base::m_graph.edges()[edge_indx]; + // first mobility activity + if (indx_schedule == first_mobility[edge_indx]) { + auto& node_from = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + auto& node_to = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility; + m_edge_func(Base::m_t, 0.0, e.property, node_from, node_to, 0); + } + // next mobility activity + else if (indx_schedule > first_mobility[edge_indx]) { + auto current_node_indx = schedule_edges[edge_indx][indx_schedule]; + bool in_mobility_node = mobility_schedule_edges[edge_indx][indx_schedule]; + + // determine dt, which is equal to the last integration/syncronization point in the current node + auto integrator_schedule_row = ln_int_schedule[current_node_indx]; + if (in_mobility_node) + integrator_schedule_row = mb_int_schedule[current_node_indx]; + // search the index of indx_schedule in the integrator schedule + const size_t indx_current = std::distance( + integrator_schedule_row.begin(), + std::lower_bound(integrator_schedule_row.begin(), integrator_schedule_row.end(), indx_schedule)); + + if (integrator_schedule_row[indx_current] != indx_schedule) + std::cout << "Error in schedule." + << "\n"; + + ScalarType dt_mobility; + if (indx_current == 0) { + dt_mobility = round_second_decimal(e.traveltime / e.path.size()); + if (dt_mobility < 0.01) + dt_mobility = 0.01; + } + else { + dt_mobility = + round_second_decimal((static_cast(integrator_schedule_row[indx_current]) - + static_cast(integrator_schedule_row[indx_current - 1])) / + 100 + + epsilon); + } + + // We have two cases. Either, we want to send the individuals to the next node, or we just want + // to update their state since a syncronization step is necessary in the current node. + if ((schedule_edges[edge_indx][indx_schedule] != schedule_edges[edge_indx][indx_schedule - 1]) || + (mobility_schedule_edges[edge_indx][indx_schedule] != + mobility_schedule_edges[edge_indx][indx_schedule - 1])) { + auto& node_from = + mobility_schedule_edges[edge_indx][indx_schedule - 1] + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility + : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + + auto& node_to = mobility_schedule_edges[edge_indx][indx_schedule] + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility + : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; + + if (indx_schedule < mobility_schedule_edges[edge_indx].size() - 1) { + m_edge_func(Base::m_t, dt_mobility, e.property, node_from, node_to, 1); + } + else + // the last time step is handled differently since we have to close the timeseries + m_edge_func(Base::m_t, dt_mobility, e.property, node_from, + Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property, 2); + } + else { + auto& node_from = + mobility_schedule_edges[edge_indx][indx_schedule - 1] + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility + : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + + auto& node_to = mobility_schedule_edges[edge_indx][indx_schedule] + ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility + : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; + + assert(node_from.get_result().get_last_value() == node_to.get_result().get_last_value()); + m_edge_func(Base::m_t, dt_mobility, e.property, node_from, node_to, 3); + } + } + } + } + + void advance_local_nodes(size_t indx_schedule) + { + for (const auto& n_indx : nodes_mobility[indx_schedule]) { + auto& n = Base::m_graph.nodes()[n_indx]; + const size_t indx_current = std::distance( + ln_int_schedule[n_indx].begin(), + std::lower_bound(ln_int_schedule[n_indx].begin(), ln_int_schedule[n_indx].end(), indx_schedule)); + + const size_t val_next = + (indx_current == ln_int_schedule[n_indx].size() - 1) ? 100 : ln_int_schedule[n_indx][indx_current + 1]; + const ScalarType next_dt = + round_second_decimal((static_cast(val_next) - indx_schedule) / 100 + epsilon); + m_node_func(Base::m_t, next_dt, n.property); + } + } + + void advance_mobility_nodes(size_t indx_schedule) + { + for (const size_t& n_indx : nodes_mobility_m[indx_schedule]) { + auto& n = Base::m_graph.nodes()[n_indx]; + const size_t indx_current = std::distance( + mb_int_schedule[n_indx].begin(), + std::lower_bound(mb_int_schedule[n_indx].begin(), mb_int_schedule[n_indx].end(), indx_schedule)); + const size_t val_next = + (indx_current == mb_int_schedule[n_indx].size() - 1) ? 100 : mb_int_schedule[n_indx][indx_current + 1]; + const ScalarType next_dt = + round_second_decimal((static_cast(val_next) - indx_schedule) / 100 + epsilon); + + // get all time points from the last integration step + auto& last_time_point = n.mobility.get_result().get_time(n.mobility.get_result().get_num_time_points() - 1); + // wenn last_time_point nicht innerhalb eines intervalls von 1-e10 von t liegt, dann setzte den letzten Zeitpunkt auf m_t + if (std::fabs(last_time_point - Base::m_t) > 1e-10) { + n.mobility.get_result().get_last_time() = Base::m_t; + } + m_node_func(Base::m_t, next_dt, n.mobility); + Eigen::Index indx_min; + while (n.mobility.get_result().get_last_value().minCoeff(&indx_min) < -1e-7) { + Eigen::Index indx_max; + n.mobility.get_result().get_last_value().maxCoeff(&indx_max); + n.mobility.get_result().get_last_value()[indx_max] -= + n.mobility.get_result().get_last_value()[indx_min]; + n.mobility.get_result().get_last_value()[indx_min] = 0; + } + } + } }; template diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h index f2cf89e331..1134c223c0 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_detailed.h +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -117,6 +117,54 @@ class GraphDetailed : public Graph public: using Graph::Graph; + /** + * @brief range of nodes + */ + auto nodes() + { + return make_range(begin(m_nodes), end(m_nodes)); + } + + /** + * @brief range of nodes + */ + auto nodes() const + { + return make_range(begin(m_nodes), end(m_nodes)); + }; + + /** + * @brief range of edges + */ + auto edges() + { + return make_range(begin(m_edges), end(m_edges)); + } + + /** + * @brief range of edges + */ + auto edges() const + { + return make_range(begin(m_edges), end(m_edges)); + } + + /** + * @brief range of edges going out from a specific node + */ + auto out_edges(size_t node_idx) + { + return out_edges(begin(m_edges), end(m_edges), node_idx); + } + + /** + * @brief range of edges going out from a specific node + */ + auto out_edges(size_t node_idx) const + { + return out_edges(begin(m_edges), end(m_edges), node_idx); + } + template NodeDetailed& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2) { @@ -185,12 +233,12 @@ class GraphDetailed : public Graph */ template -IOResult set_nodes(const Parameters& params, Date start_date, 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, Graph& params_graph, ReadFunction&& read_func, - NodeIdFunction&& node_func, const std::vector& scaling_factor_inf, - double scaling_factor_icu, double tnt_capacity_factor, int num_days = 0, - bool export_time_series = false) +IOResult set_nodes_detailed(const Parameters& params, Date start_date, 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, GraphDetailed& params_graph, + ReadFunction&& read_func, NodeIdFunction&& node_func, + const std::vector& scaling_factor_inf, double scaling_factor_icu, + double tnt_capacity_factor, int num_days = 0, bool export_time_series = false) { BOOST_OUTCOME_TRY(duration_stay, mio::read_duration_stay(stay_times_data_path)); @@ -258,11 +306,12 @@ IOResult set_nodes(const Parameters& params, Date start_date, Date end_dat */ template -IOResult set_edges(const std::string& travel_times_path, const std::string mobility_data_path, - const std::string& travelpath_path, Graph& params_graph, - std::initializer_list& migrating_compartments, size_t contact_locations_size, - std::vector commuting_weights = std::vector{}, - ScalarType theshold_edges = 4e-5) +IOResult +set_edges_detailed(const std::string& travel_times_path, const std::string mobility_data_path, + const std::string& travelpath_path, GraphDetailed& params_graph, + std::initializer_list& migrating_compartments, size_t contact_locations_size, + std::vector commuting_weights = std::vector{}, + ScalarType theshold_edges = 4e-5) { BOOST_OUTCOME_TRY(mobility_data_commuter, mio::read_mobility_plain(mobility_data_path)); BOOST_OUTCOME_TRY(travel_times, mio::read_mobility_plain(travel_times_path)); From 3eda22ad6b4946e8d231e85bceca26a3aaa342c9 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Tue, 12 Dec 2023 14:00:01 +0100 Subject: [PATCH 11/11] compiling now, but apply_migration still wrong --- .../ode_secirvvs_mobility_detailed.cpp | 9 +- cpp/memilio/compartments/parameter_studies.h | 38 +- cpp/memilio/data/analyze_result.h | 13 + cpp/memilio/mobility/graph_simulation.h | 164 +++++---- .../metapopulation_mobility_detailed.h | 336 +++++++++++++----- cpp/models/ode_secirvvs/parameter_space.cpp | 50 +++ cpp/models/ode_secirvvs/parameter_space.h | 4 + 7 files changed, 443 insertions(+), 171 deletions(-) diff --git a/cpp/examples/ode_secirvvs_mobility_detailed.cpp b/cpp/examples/ode_secirvvs_mobility_detailed.cpp index 752a8a12e3..ad63b1c309 100644 --- a/cpp/examples/ode_secirvvs_mobility_detailed.cpp +++ b/cpp/examples/ode_secirvvs_mobility_detailed.cpp @@ -370,7 +370,7 @@ void init_pop_cologne_szenario(mio::Graph> +mio::IOResult> get_graph(const std::string data_dir, const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool edges) { auto mobility_data_commuter = @@ -552,7 +552,6 @@ mio::IOResult run() bool szenario_cologne = false; - // auto params_graph = get_graph(num_days, masks); BOOST_OUTCOME_TRY(created, get_graph(data_dir, num_days, masks, ffp2, szenario_cologne, edges)); auto params_graph = created; @@ -581,8 +580,8 @@ mio::IOResult run() }); // parameter study - auto parameter_study = - mio::ParameterStudy>(params_graph, 0.0, num_days, 0.01, num_runs); + auto parameter_study = mio::ParameterStudyDetailed>( + params_graph, 0.0, num_days, 0.01, num_runs); if (mio::mpi::is_root()) { printf("Seeds: "); for (auto s : parameter_study.get_rng().get_seeds()) { @@ -593,7 +592,7 @@ mio::IOResult run() auto save_single_run_result = mio::IOResult(mio::success()); auto ensemble = parameter_study.run( [&](auto&& graph) { - return draw_sample(graph, false); + return draw_sample(graph); }, [&](auto results_graph, auto&& run_idx) { auto interpolated_result = mio::interpolate_simulation_result(results_graph); diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 63e76cdd1c..893dadd14f 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -351,6 +351,26 @@ class ParameterStudy return m_rng; } +protected: + // adaptive time step of the integrator (will be corrected if too large/small) + double m_dt_integration = 0.1; + // + RandomNumberGenerator m_rng; + + std::vector distribute_runs(size_t num_runs, int num_procs) + { + //evenly distribute runs + //lower processes do one more run if runs are not evenly distributable + auto num_runs_local = num_runs / num_procs; //integer division! + auto remainder = num_runs % num_procs; + + std::vector run_distribution(num_procs); + std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1); + std::fill(run_distribution.begin() + remainder, run_distribution.end(), num_runs_local); + + return run_distribution; + } + private: //sample parameters and create simulation template @@ -379,20 +399,6 @@ class ParameterStudy return make_migration_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); } - std::vector distribute_runs(size_t num_runs, int num_procs) - { - //evenly distribute runs - //lower processes do one more run if runs are not evenly distributable - auto num_runs_local = num_runs / num_procs; //integer division! - auto remainder = num_runs % num_procs; - - std::vector run_distribution(num_procs); - std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1); - std::fill(run_distribution.begin() + remainder, run_distribution.end(), num_runs_local); - - return run_distribution; - } - private: // Stores Graph with the names and ranges of all parameters ParametersGraph m_graph; @@ -405,10 +411,6 @@ class ParameterStudy double m_tmax; // time step of the graph double m_dt_graph_sim; - // adaptive time step of the integrator (will be corrected if too large/small) - double m_dt_integration = 0.1; - // - RandomNumberGenerator m_rng; }; } // namespace mio diff --git a/cpp/memilio/data/analyze_result.h b/cpp/memilio/data/analyze_result.h index ecb501247a..ab9577cbdc 100644 --- a/cpp/memilio/data/analyze_result.h +++ b/cpp/memilio/data/analyze_result.h @@ -22,6 +22,7 @@ #include "memilio/utils/time_series.h" #include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" #include #include @@ -120,6 +121,18 @@ interpolate_simulation_result(const Graph, MigrationE return interpolated; } +template +std::vector> +interpolate_simulation_result(const GraphDetailed, MigrationEdge>& graph_result) +{ + std::vector> interpolated; + interpolated.reserve(graph_result.nodes().size()); + std::transform(graph_result.nodes().begin(), graph_result.nodes().end(), std::back_inserter(interpolated), + [](auto& n) { + return interpolate_simulation_result(n.property.get_result()); + }); + return interpolated; +} /** * Compute the distance between two SECIR simulation results. * The distance is the 2-norm of the element-wise difference of the two results. diff --git a/cpp/memilio/mobility/graph_simulation.h b/cpp/memilio/mobility/graph_simulation.h index 9f59d00737..7a8328fb67 100644 --- a/cpp/memilio/mobility/graph_simulation.h +++ b/cpp/memilio/mobility/graph_simulation.h @@ -249,28 +249,33 @@ class GraphSimulationStochastic RandomNumberGenerator m_rng; }; -template > +template > class GraphSimulationDetailed - : public GraphSimulationBase> { public: using edge_function = edge_f; - using Base = GraphSimulationBase; - using node_function = std::function; + using node_function = std::function; - GraphSimulationDetailed(double t0, double dt, const Graph& g, const node_function& node_func, + GraphSimulationDetailed(double t0, double dt, const GraphDetailed& g, const node_function& node_func, const edge_function&& edge_func) - : Base(t0, dt, g, node_func, std::move(edge_func)) + : m_dt(dt) + , m_t(double(t0)) + , m_graph(g) + , m_node_func(node_func) + , m_edge_func(edge_func) { } - GraphSimulationDetailed(double t0, double dt, Graph&& g, const node_function& node_func, + GraphSimulationDetailed(double t0, double dt, GraphDetailed&& g, const node_function& node_func, const edge_function&& edge_func) - : Base(t0, dt, std::move(g), node_func, std::move(edge_func)) + : m_dt(dt) + , m_t(double(t0)) + , m_graph(std::move(g)) + , m_node_func(node_func) + , m_edge_func(edge_func) { } @@ -279,17 +284,37 @@ class GraphSimulationDetailed return std::round(x * 100.0) / 100.0; } + double get_t() const + { + return m_t; + } + + GraphDetailed& get_graph() & + { + return m_graph; + } + + const GraphDetailed& get_graph() const& + { + return m_graph; + } + + GraphDetailed&& get_graph() && + { + return std::move(m_graph); + } + void advance(double t_max = 1.0) { ScalarType dt_first_mobility = - std::accumulate(Base::m_graph.edges().begin(), Base::m_graph.edges().end(), - std::numeric_limits::max(), [&](ScalarType current_min, const auto& e) { + std::accumulate(m_graph.edges().begin(), m_graph.edges().end(), std::numeric_limits::max(), + [&](ScalarType current_min, const auto& e) { auto traveltime_per_region = round_second_decimal(e.traveltime / e.path.size()); if (traveltime_per_region < 0.01) traveltime_per_region = 0.01; auto start_mobility = round_second_decimal(1 - 2 * (traveltime_per_region * e.path.size()) - - Base::m_graph.nodes()[e.end_node_idx].stay_duration); + m_graph.nodes()[e.end_node_idx].stay_duration); if (start_mobility < 0) { start_mobility = 0.; } @@ -297,33 +322,33 @@ class GraphSimulationDetailed }); // set population to zero in mobility nodes before starting - for (auto& n : Base::m_graph.nodes()) { + for (auto& n : m_graph.nodes()) { n.mobility.get_result().get_last_value().setZero(); } auto min_dt = 0.01; - double t_begin = Base::m_t - 1.; + double t_begin = m_t - 1.; // calc schedule for each edge precompute_schedule(); - while (Base::m_t - epsilon < t_max) { + while (m_t - epsilon < t_max) { // auto start = std::chrono::system_clock::now(); t_begin += 1; - if (Base::m_t + dt_first_mobility > t_max) { - dt_first_mobility = t_max - Base::m_t; - for (auto& n : Base::m_graph.nodes()) { - m_node_func(Base::m_t, dt_first_mobility, n.property); + if (m_t + dt_first_mobility > t_max) { + dt_first_mobility = t_max - m_t; + for (auto& n : m_graph.nodes()) { + m_node_func(m_t, dt_first_mobility, n.property); } break; } - for (auto& node : Base::m_graph.nodes()) { + for (auto& node : m_graph.nodes()) { node.mobility.get_simulation().set_integrator(std::make_shared()); } size_t indx_schedule = 0; - while (t_begin + 1 > Base::m_t + 1e-10) { + while (t_begin + 1 > m_t + 1e-10) { advance_edges(indx_schedule); // first we integrate the nodes in time. Afterwards the update on the edges is done. @@ -332,18 +357,24 @@ class GraphSimulationDetailed advance_mobility_nodes(indx_schedule); indx_schedule++; - Base::m_t += min_dt; + m_t += min_dt; } // At the end of the day. we set each compartment zero for all mobility nodes since we have to estimate // the state of the indivuals moving between different counties. // Therefore there can be differences with the states given by the integrator used for the mobility node. - for (auto& n : Base::m_graph.nodes()) { + for (auto& n : m_graph.nodes()) { n.mobility.get_result().get_last_value().setZero(); } } } private: + GraphDetailed m_graph; + double m_t; + double m_dt; + node_function m_node_func; + edge_function m_edge_func; + // describes the schedule for each edge, i.e. which node is visited at which time step std::vector> schedule_edges; @@ -371,10 +402,10 @@ class GraphSimulationDetailed { // print transi const size_t timesteps = 100; - schedule_edges.reserve(Base::m_graph.edges().size()); - mobility_schedule_edges.reserve(Base::m_graph.edges().size()); + schedule_edges.reserve(m_graph.edges().size()); + mobility_schedule_edges.reserve(m_graph.edges().size()); - for (auto& e : Base::m_graph.edges()) { + for (auto& e : m_graph.edges()) { // 100 since we round to second decimal std::vector schedule(timesteps, 0.); std::vector is_mobility_node(timesteps, false); @@ -385,7 +416,7 @@ class GraphSimulationDetailed traveltime_per_region = 0.01; double start_mobility = round_second_decimal(0 + 1 - 2 * (traveltime_per_region * e.path.size()) - - Base::m_graph.nodes()[e.end_node_idx].stay_duration); + m_graph.nodes()[e.end_node_idx].stay_duration); if (start_mobility < 0.0) { start_mobility = 0.; } @@ -434,7 +465,7 @@ class GraphSimulationDetailed // iterate over nodes size_t indx_node = 0; - for (auto& n : Base::m_graph.nodes()) { + for (auto& n : m_graph.nodes()) { // local node and automatical add zero, since we want to start at t=0 and start with integrating all nodes to the next dt std::vector order_local_node = {0}; std::vector indx_edges; @@ -490,8 +521,8 @@ class GraphSimulationDetailed } auto indx_edge = 0; - first_mobility.reserve(Base::m_graph.edges().size()); - for (auto& e : Base::m_graph.edges()) { + first_mobility.reserve(m_graph.edges().size()); + for (auto& e : m_graph.edges()) { this->first_mobility[indx_edge] = std::distance( mobility_schedule_edges[indx_edge].begin(), std::find(mobility_schedule_edges[indx_edge].begin(), mobility_schedule_edges[indx_edge].end(), true)); @@ -516,13 +547,13 @@ class GraphSimulationDetailed edges_mobility.push_back(std::move(temp_edges_mobility)); // same for nodes - std::vector temp_nodes_mobility(Base::m_graph.nodes().size()); + std::vector temp_nodes_mobility(m_graph.nodes().size()); std::iota(temp_nodes_mobility.begin(), temp_nodes_mobility.end(), 0); nodes_mobility.emplace_back(std::move(temp_nodes_mobility)); std::vector temp_nodes_mobility_m; size_t node_indx = 0; - for (auto& n : Base::m_graph.nodes()) { + for (auto& n : m_graph.nodes()) { if (std::binary_search(mb_int_schedule[node_indx].begin(), mb_int_schedule[node_indx].end(), 0)) { temp_nodes_mobility_m.push_back(node_indx); node_indx++; @@ -533,7 +564,7 @@ class GraphSimulationDetailed for (size_t indx_current = 1; indx_current < timesteps; ++indx_current) { std::vector temp_edge_mobility; indx_edge = 0; - for (auto& e : Base::m_graph.edges()) { + for (auto& e : m_graph.edges()) { auto current_node_indx = schedule_edges[indx_edge][indx_current]; if (indx_current >= first_mobility[indx_edge]) { if (mobility_schedule_edges[indx_edge][indx_current] && @@ -553,7 +584,7 @@ class GraphSimulationDetailed temp_nodes_mobility.clear(); temp_nodes_mobility_m.clear(); node_indx = 0; - for (auto& n : Base::m_graph.nodes()) { + for (auto& n : m_graph.nodes()) { if (std::binary_search(ln_int_schedule[node_indx].begin(), ln_int_schedule[node_indx].end(), indx_current)) { temp_nodes_mobility.push_back(node_indx); @@ -573,12 +604,12 @@ class GraphSimulationDetailed void advance_edges(size_t indx_schedule) { for (const auto& edge_indx : edges_mobility[indx_schedule]) { - auto& e = Base::m_graph.edges()[edge_indx]; + auto& e = m_graph.edges()[edge_indx]; // first mobility activity if (indx_schedule == first_mobility[edge_indx]) { - auto& node_from = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; - auto& node_to = Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility; - m_edge_func(Base::m_t, 0.0, e.property, node_from, node_to, 0); + auto& node_from = m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + auto& node_to = m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility; + // m_edge_func(m_t, 0.0, e.property, node_from, node_to, 0); } // next mobility activity else if (indx_schedule > first_mobility[edge_indx]) { @@ -617,35 +648,33 @@ class GraphSimulationDetailed if ((schedule_edges[edge_indx][indx_schedule] != schedule_edges[edge_indx][indx_schedule - 1]) || (mobility_schedule_edges[edge_indx][indx_schedule] != mobility_schedule_edges[edge_indx][indx_schedule - 1])) { - auto& node_from = - mobility_schedule_edges[edge_indx][indx_schedule - 1] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility - : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + auto& node_from = mobility_schedule_edges[edge_indx][indx_schedule - 1] + ? m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility + : m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; auto& node_to = mobility_schedule_edges[edge_indx][indx_schedule] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility - : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; + ? m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility + : m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; if (indx_schedule < mobility_schedule_edges[edge_indx].size() - 1) { - m_edge_func(Base::m_t, dt_mobility, e.property, node_from, node_to, 1); + // m_edge_func(m_t, dt_mobility, e.property, node_from, node_to, 1); } - else - // the last time step is handled differently since we have to close the timeseries - m_edge_func(Base::m_t, dt_mobility, e.property, node_from, - Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property, 2); + // else + // the last time step is handled differently since we have to close the timeseries + // m_edge_func(m_t, dt_mobility, e.property, node_from, + // m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property, 2); } else { - auto& node_from = - mobility_schedule_edges[edge_indx][indx_schedule - 1] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility - : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + auto& node_from = mobility_schedule_edges[edge_indx][indx_schedule - 1] + ? m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility + : m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; auto& node_to = mobility_schedule_edges[edge_indx][indx_schedule] - ? Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility - : Base::m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; + ? m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility + : m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; assert(node_from.get_result().get_last_value() == node_to.get_result().get_last_value()); - m_edge_func(Base::m_t, dt_mobility, e.property, node_from, node_to, 3); + // m_edge_func(m_t, dt_mobility, e.property, node_from, node_to, 3); } } } @@ -654,7 +683,7 @@ class GraphSimulationDetailed void advance_local_nodes(size_t indx_schedule) { for (const auto& n_indx : nodes_mobility[indx_schedule]) { - auto& n = Base::m_graph.nodes()[n_indx]; + auto& n = m_graph.nodes()[n_indx]; const size_t indx_current = std::distance( ln_int_schedule[n_indx].begin(), std::lower_bound(ln_int_schedule[n_indx].begin(), ln_int_schedule[n_indx].end(), indx_schedule)); @@ -663,14 +692,14 @@ class GraphSimulationDetailed (indx_current == ln_int_schedule[n_indx].size() - 1) ? 100 : ln_int_schedule[n_indx][indx_current + 1]; const ScalarType next_dt = round_second_decimal((static_cast(val_next) - indx_schedule) / 100 + epsilon); - m_node_func(Base::m_t, next_dt, n.property); + m_node_func(m_t, next_dt, n.property); } } void advance_mobility_nodes(size_t indx_schedule) { for (const size_t& n_indx : nodes_mobility_m[indx_schedule]) { - auto& n = Base::m_graph.nodes()[n_indx]; + auto& n = m_graph.nodes()[n_indx]; const size_t indx_current = std::distance( mb_int_schedule[n_indx].begin(), std::lower_bound(mb_int_schedule[n_indx].begin(), mb_int_schedule[n_indx].end(), indx_schedule)); @@ -682,10 +711,10 @@ class GraphSimulationDetailed // get all time points from the last integration step auto& last_time_point = n.mobility.get_result().get_time(n.mobility.get_result().get_num_time_points() - 1); // wenn last_time_point nicht innerhalb eines intervalls von 1-e10 von t liegt, dann setzte den letzten Zeitpunkt auf m_t - if (std::fabs(last_time_point - Base::m_t) > 1e-10) { - n.mobility.get_result().get_last_time() = Base::m_t; + if (std::fabs(last_time_point - m_t) > 1e-10) { + n.mobility.get_result().get_last_time() = m_t; } - m_node_func(Base::m_t, next_dt, n.mobility); + m_node_func(m_t, next_dt, n.mobility); Eigen::Index indx_min; while (n.mobility.get_result().get_last_value().minCoeff(&indx_min) < -1e-7) { Eigen::Index indx_max; @@ -712,5 +741,12 @@ auto make_graph_sim_stochastic(double t0, double dt, Graph&& g, NodeF&& node_fun t0, dt, std::forward(g), std::forward(node_func), std::forward(edge_func)); } +template +auto make_graph_sim_detailed(double t0, double dt, GraphDetailed&& g, NodeF&& node_func, EdgeF&& edge_func) +{ + return GraphSimulationDetailed>( + t0, dt, std::forward(g), std::forward(node_func), std::forward(edge_func)); +} + } // namespace mio #endif //MIO_MOBILITY_GRAPH_SIMULATION_H diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h index 1134c223c0..af21ebcb9a 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_detailed.h +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -20,6 +20,7 @@ #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" @@ -180,16 +181,17 @@ class GraphDetailed : public Graph return m_nodes.back(); } - template - EdgeDetailed& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, Args&&... args) + EdgeDetailed& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, + mio::MigrationParameters& args) { assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); - return *insert_sorted_replace( - m_edges, EdgeDetailed(start_node_idx, end_node_idx, traveltime, std::forward(args)...), - [](auto&& e1, auto&& e2) { - return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx - : e1.start_node_idx < e2.start_node_idx; - }); + return *insert_sorted_replace(m_edges, + EdgeDetailed(start_node_idx, end_node_idx, traveltime, args), + [](auto&& e1, auto&& e2) { + return e1.start_node_idx == e2.start_node_idx + ? e1.end_node_idx < e2.end_node_idx + : e1.start_node_idx < e2.start_node_idx; + }); } template @@ -360,11 +362,188 @@ set_edges_detailed(const std::string& travel_times_path, const std::string mobil return success(); } -class MigrationEdgeDetailed : public MigrationEdge +// class MigrationEdgeDetailed : public MigrationEdge +// { +// public: +// template +// void apply_migration(double t, double dt, SimulationNode& node_from, SimulationNode& node_to, int mode); +// }; + +// /** +// * edge functor for migration simulation. +// * @see MigrationEdge::apply_migration +// */ +// template +// void apply_migration(double t, double dt, MigrationEdgeDetailed& migrationEdgeDetailed, SimulationNode& node_from, +// SimulationNode& node_to, int mode) +// { +// migrationEdgeDetailed.apply_migration(t, dt, node_from, node_to, mode); +// } + +template +class ParameterStudyDetailed : public ParameterStudy { public: - template - void apply_migration(double t, double dt, SimulationNode& node_from, SimulationNode& node_to, int mode); + /** + * 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 ParametersGraphDetailed = mio::GraphDetailed; + /** + * The Graph type that stores simulations and their results of each run. + * This is the output of ParameterStudies for each run. + */ + using SimulationGraphDetailed = mio::GraphDetailed, mio::MigrationEdge>; + + /** + * create study for graph of compartment models. + * @param graph graph of parameters + * @param t0 start time of simulations + * @param tmax end time of simulations + * @param graph_sim_dt time step of graph simulation + * @param num_runs number of runs + */ + ParameterStudyDetailed(const ParametersGraphDetailed& graph, double t0, double tmax, double graph_sim_dt, + size_t num_runs) + : ParameterStudy(graph, t0, tmax, graph_sim_dt, num_runs) + , m_graph(graph) + , m_num_runs(num_runs) + , m_t0{t0} + , m_tmax{tmax} + , m_dt_graph_sim(graph_sim_dt) + { + } + + /* + * @brief Carry out all simulations in the parameter study. + * Save memory and enable more runs by immediately processing and/or discarding the result. + * The result processing function is called when a run is finished. It receives the result of the run + * (a SimulationGraphDetailed object) and an ordered index. The values returned by the result processing function + * are gathered and returned as a list. + * This function is parallelized if memilio is configured with MEMILIO_ENABLE_MPI. + * The MPI processes each contribute a share of the runs. The sample function and result processing function + * are called in the same process that performs the run. The results returned by the result processing function are + * gathered at the root process and returned as a list by the root in the same order as if the programm + * were running sequentially. Processes other than the root return an empty list. + * @param sample_graph Function that receives the ParametersGraph and returns a sampled copy. + * @param result_processing_function Processing function for simulation results, e.g., output function. + * @returns At the root process, a list of values per run that have been returned from the result processing function. + * At all other processes, an empty list. + * @tparam SampleGraphFunction Callable type, accepts instance of ParametersGraph. + * @tparam HandleSimulationResultFunction Callable type, accepts instance of SimulationGraphDetailed and an index of type size_t. + */ + template + std::vector> + run(SampleGraphFunction sample_graph, HandleSimulationResultFunction result_processing_function) + { + int num_procs, rank; +#ifdef MEMILIO_ENABLE_MPI + MPI_Comm_size(mpi::get_world(), &num_procs); + MPI_Comm_rank(mpi::get_world(), &rank); +#else + num_procs = 1; + rank = 0; +#endif + + //The ParameterDistributions used for sampling parameters use thread_local_rng() + //So we set our own RNG to be used. + //Assume that sampling uses the thread_local_rng() and isn't multithreaded + this->m_rng.synchronize(); + thread_local_rng() = this->m_rng; + + auto run_distribution = this->distribute_runs(m_num_runs, num_procs); + auto start_run_idx = + std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0)); + auto end_run_idx = start_run_idx + run_distribution[size_t(rank)]; + + std::vector> + ensemble_result; + ensemble_result.reserve(m_num_runs); + + for (size_t run_idx = start_run_idx; run_idx < end_run_idx; run_idx++) { + log(LogLevel::info, "ParameterStudies: run {}", run_idx); + + //prepare rng for this run by setting the counter to the right offset + //Add the old counter so that this call of run() produces different results + //from the previous call + auto run_rng_counter = + this->m_rng.get_counter() + + rng_totalsequence_counter(static_cast(run_idx), Counter(0)); + thread_local_rng().set_counter(run_rng_counter); + + //sample + auto sim = create_sampled_sim(sample_graph); + log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", + (thread_local_rng().get_counter() - run_rng_counter).get()); + + //perform run + sim.advance(m_tmax); + + //handle result and store + ensemble_result.emplace_back(result_processing_function(std::move(sim).get_graph(), run_idx)); + } + + //Set the counter of our RNG so that future calls of run() produce different parameters. + this->m_rng.set_counter(this->m_rng.get_counter() + + rng_totalsequence_counter(m_num_runs, Counter(0))); + +#ifdef MEMILIO_ENABLE_MPI + //gather results + if (rank == 0) { + for (int src_rank = 1; src_rank < num_procs; ++src_rank) { + int bytes_size; + MPI_Recv(&bytes_size, 1, MPI_INT, src_rank, 0, mpi::get_world(), MPI_STATUS_IGNORE); + ByteStream bytes(bytes_size); + MPI_Recv(bytes.data(), bytes.data_size(), MPI_BYTE, src_rank, 0, mpi::get_world(), MPI_STATUS_IGNORE); + + auto src_ensemble_results = deserialize_binary(bytes, Tag{}); + if (!src_ensemble_results) { + log_error("Error receiving ensemble results from rank {}.", src_rank); + } + std::copy(src_ensemble_results.value().begin(), src_ensemble_results.value().end(), + std::back_inserter(ensemble_result)); + } + } + else { + auto bytes = serialize_binary(ensemble_result); + auto bytes_size = int(bytes.data_size()); + MPI_Send(&bytes_size, 1, MPI_INT, 0, 0, mpi::get_world()); + MPI_Send(bytes.data(), bytes.data_size(), MPI_BYTE, 0, 0, mpi::get_world()); + ensemble_result.clear(); //only return root process + } +#endif + + return ensemble_result; + } + +private: + //sample parameters and create simulation + template + mio::GraphSimulationDetailed create_sampled_sim(SampleGraphFunction sample_graph) + { + SimulationGraphDetailed sim_graph; + + auto sampled_graph = sample_graph(m_graph); + for (auto&& node : sampled_graph.nodes()) { + sim_graph.add_node(node.id, node.stay_duration, node.property, node.mobility, m_t0, this->m_dt_integration); + } + for (auto&& edge : sampled_graph.edges()) { + sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.property); + } + return make_migration_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); + } + +private: + // Stores Graph with the names and ranges of all parameters + ParametersGraphDetailed m_graph; + size_t m_num_runs; + double m_t0; + double m_tmax; + double m_dt_graph_sim; }; /** @@ -423,73 +602,62 @@ void move_migrated(Eigen::Ref> migrated, Eigen::Ref> resul } } -template -void MigrationEdgeDetailed::apply_migration(double t, double dt, SimulationNode& node_from, - SimulationNode& node_to, int mode) -{ - - if (mode == 0) { - //normal daily migration - m_migrated.add_time_point( - t, (node_from.get_last_state().array() * m_parameters.get_coefficients().get_matrix_at(t).array() * - get_migration_factors(node_from, t, node_from.get_last_state()).array()) - .matrix()); - m_return_times.add_time_point(t + dt); - move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), - node_to.get_result().get_last_value()); - } - // change county of migrated - else if (mode == 1) { - // update status of migrated before moving to next county - IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); - update_status_migrated(m_migrated.get_last_value(), node_from.get_simulation(), integrator_node, - node_from.get_result().get_last_value(), t, dt); - move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), - node_to.get_result().get_last_value()); - } - // option for last time point to remove time points - else if (mode == 2) { - Eigen::Index idx = m_return_times.get_num_time_points() - 1; - IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); - update_status_migrated(m_migrated[idx], node_from.get_simulation(), integrator_node, - node_from.get_result().get_last_value(), t, dt); - - move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), - node_to.get_result().get_last_value()); - - for (Eigen::Index i = m_return_times.get_num_time_points() - 1; i >= 0; --i) { - if (m_return_times.get_time(i) <= t) { - m_migrated.remove_time_point(i); - m_return_times.remove_time_point(i); - } - } - } - // just update status of migrated - else if (mode == 3) { - Eigen::Index idx = m_return_times.get_num_time_points() - 1; - IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); - update_status_migrated(m_migrated[idx], node_from.get_simulation(), integrator_node, - node_from.get_result().get_last_value(), t, dt); - - move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), - node_from.get_result().get_last_value()); - } - else { - std::cout << "Invalid input mode. Should be 0 or 1." - << "\n"; - } -} - -/** - * edge functor for migration simulation. - * @see MigrationEdge::apply_migration - */ -template -void apply_migration(double t, double dt, MigrationEdge& migrationEdge, SimulationNode& node_from, - SimulationNode& node_to, int mode) -{ - migrationEdge.apply_migration(t, dt, node_from, node_to, mode); -} +// template +// void MigrationEdge::apply_migration(double t, double dt, SimulationNode& node_from, SimulationNode& node_to, +// int mode) // +// { + +// if (mode == 0) { +// //normal daily migration +// m_migrated.add_time_point( +// t, (node_from.get_last_state().array() * m_parameters.get_coefficients().get_matrix_at(t).array() * +// get_migration_factors(node_from, t, node_from.get_last_state()).array()) +// .matrix()); +// m_return_times.add_time_point(t + dt); +// move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), +// node_to.get_result().get_last_value()); +// } +// // change county of migrated +// else if (mode == 1) { +// // update status of migrated before moving to next county +// IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); +// update_status_migrated(m_migrated.get_last_value(), node_from.get_simulation(), integrator_node, +// node_from.get_result().get_last_value(), t, dt); +// move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), +// node_to.get_result().get_last_value()); +// } +// // option for last time point to remove time points +// else if (mode == 2) { +// Eigen::Index idx = m_return_times.get_num_time_points() - 1; +// IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); +// update_status_migrated(m_migrated[idx], node_from.get_simulation(), integrator_node, +// node_from.get_result().get_last_value(), t, dt); + +// move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), +// node_to.get_result().get_last_value()); + +// for (Eigen::Index i = m_return_times.get_num_time_points() - 1; i >= 0; --i) { +// if (m_return_times.get_time(i) <= t) { +// m_migrated.remove_time_point(i); +// m_return_times.remove_time_point(i); +// } +// } +// } +// // just update status of migrated +// else if (mode == 3) { +// Eigen::Index idx = m_return_times.get_num_time_points() - 1; +// IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); +// update_status_migrated(m_migrated[idx], node_from.get_simulation(), integrator_node, +// node_from.get_result().get_last_value(), t, dt); + +// move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), +// node_from.get_result().get_last_value()); +// } +// else { +// std::cout << "Invalid input mode. Should be 0 or 1." +// << "\n"; +// } +// } /** * create a migration simulation. @@ -502,17 +670,17 @@ void apply_migration(double t, double dt, MigrationEdge& migrationEdge, Simulati * @{ */ template -GraphSimulationDetailed, MigrationEdge>> -make_migration_sim(double t0, double dt, const Graph, MigrationEdge>& graph) +GraphSimulationDetailed, MigrationEdge>> +make_migration_sim(double t0, double dt, const GraphDetailed, MigrationEdge>& graph) { - return make_graph_sim(t0, dt, graph, &evolve_model, &apply_migration); + return make_graph_sim_detailed(t0, dt, graph, &evolve_model, &apply_migration); } template -GraphSimulationDetailed, MigrationEdge>> -make_migration_sim(double t0, double dt, Graph, MigrationEdge>&& graph) +GraphSimulationDetailed, MigrationEdge>> +make_migration_sim(double t0, double dt, GraphDetailed, MigrationEdge>&& graph) { - return make_graph_sim(t0, dt, std::move(graph), &evolve_model, &apply_migration); + return make_graph_sim_detailed(t0, dt, std::move(graph), &evolve_model, &apply_migration); } } // namespace mio diff --git a/cpp/models/ode_secirvvs/parameter_space.cpp b/cpp/models/ode_secirvvs/parameter_space.cpp index 1748c307ae..9e26a39715 100644 --- a/cpp/models/ode_secirvvs/parameter_space.cpp +++ b/cpp/models/ode_secirvvs/parameter_space.cpp @@ -23,6 +23,8 @@ #include "ode_secirvvs/infection_state.h" #include "ode_secirvvs/model.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" + namespace mio { namespace osecirvvs @@ -191,5 +193,53 @@ Graph draw_sample(Graph& return sampled_graph; } +GraphDetailed draw_sample(GraphDetailed& graph) +{ + GraphDetailed sampled_graph; + + //sample global parameters + auto& shared_params_model = graph.nodes()[0].property; + draw_sample_infection(shared_params_model); + auto& shared_contacts = shared_params_model.parameters.template get(); + shared_contacts.draw_sample_dampings(); + auto& shared_dynamic_npis = shared_params_model.parameters.template get(); + shared_dynamic_npis.draw_sample(); + + for (auto& params_node : graph.nodes()) { + auto& node_model = params_node.property; + + //sample local parameters + draw_sample_demographics(params_node.property); + + //copy global parameters + //save demographic parameters so they aren't overwritten + auto local_icu_capacity = node_model.parameters.template get(); + auto local_tnt_capacity = node_model.parameters.template get(); + auto local_holidays = node_model.parameters.template get().get_school_holidays(); + auto local_daily_v1 = node_model.parameters.template get(); + auto local_daily_v2 = node_model.parameters.template get(); + node_model.parameters = shared_params_model.parameters; + node_model.parameters.template get() = local_icu_capacity; + node_model.parameters.template get() = local_tnt_capacity; + node_model.parameters.template get().get_school_holidays() = local_holidays; + node_model.parameters.template get() = local_daily_v1; + node_model.parameters.template get() = local_daily_v2; + + node_model.parameters.template get().make_matrix(); + node_model.apply_constraints(); + + sampled_graph.add_node(params_node.id, params_node.stay_duration, node_model, params_node.mobility); + } + + for (auto& edge : graph.edges()) { + auto edge_params = edge.property; + //no dynamic NPIs + //TODO: add switch to optionally enable dynamic NPIs to edges + sampled_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.path, edge_params); + } + + return sampled_graph; +} + } // namespace osecirvvs } // namespace mio diff --git a/cpp/models/ode_secirvvs/parameter_space.h b/cpp/models/ode_secirvvs/parameter_space.h index a883cc482a..1943806f45 100644 --- a/cpp/models/ode_secirvvs/parameter_space.h +++ b/cpp/models/ode_secirvvs/parameter_space.h @@ -21,6 +21,8 @@ #define ODESECIRVVS_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" @@ -63,6 +65,8 @@ void draw_sample(Model& model); */ Graph draw_sample(Graph& graph, bool variant_high); +GraphDetailed draw_sample(GraphDetailed& graph); + } // namespace osecirvvs } // namespace mio