diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index f18ab284f7..b73975d79f 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -80,9 +80,9 @@ add_executable(ode_secir_ageres_example ode_secir_ageres.cpp) target_link_libraries(ode_secir_ageres_example PRIVATE memilio ode_secir) target_compile_options(ode_secir_ageres_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) -add_executable(graph_example graph.cpp) -target_link_libraries(graph_example PRIVATE memilio ode_seir) -target_compile_options(graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ode_secir_graph_example ode_secir_graph.cpp) +target_link_libraries(ode_secir_graph_example PRIVATE memilio ode_secir) +target_compile_options(ode_secir_graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) add_executable(graph_stochastic_mobility_example graph_stochastic_mobility.cpp) target_link_libraries(graph_stochastic_mobility_example PRIVATE memilio ode_secir) @@ -127,9 +127,13 @@ if(MEMILIO_HAS_JSONCPP) endif() if(MEMILIO_HAS_HDF5 AND MEMILIO_HAS_JSONCPP) - add_executable(ode_secir_parameter_study_example ode_secir_parameter_study.cpp) - target_link_libraries(ode_secir_parameter_study_example PRIVATE memilio ode_secir) - target_compile_options(ode_secir_parameter_study_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(ode_secir_parameter_study_example ode_secir_parameter_study.cpp) + target_link_libraries(ode_secir_parameter_study_example PRIVATE memilio ode_secir) + target_compile_options(ode_secir_parameter_study_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + + add_executable(ode_secir_parameter_study_graph ode_secir_parameter_study_graph.cpp) + target_link_libraries(ode_secir_parameter_study_graph PRIVATE memilio ode_secir) + target_compile_options(ode_secir_parameter_study_graph PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) endif() if(MEMILIO_HAS_JSONCPP) diff --git a/cpp/examples/graph.cpp b/cpp/examples/graph.cpp deleted file mode 100644 index ad944254af..0000000000 --- a/cpp/examples/graph.cpp +++ /dev/null @@ -1,69 +0,0 @@ -/* -* Copyright (C) 2020-2024 MEmilio -* -* Authors: Daniel Abele -* -* Contact: Martin J. Kuehn -* -* Licensed under the Apache License, Version 2.0 (the "License"); -* you may not use this file except in compliance with the License. -* You may obtain a copy of the License at -* -* http://www.apache.org/licenses/LICENSE-2.0 -* -* Unless required by applicable law or agreed to in writing, software -* distributed under the License is distributed on an "AS IS" BASIS, -* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -* See the License for the specific language governing permissions and -* limitations under the License. -*/ -#include "ode_seir/model.h" -#include "ode_seir/infection_state.h" -#include "ode_seir/parameters.h" -#include "memilio/mobility/metapopulation_mobility_instant.h" -#include "memilio/compartments/simulation.h" - -int main() -{ - const auto t0 = 0.; - const auto tmax = 10.; - const auto dt = 0.5; //time step of mobility, daily mobility every second step - - mio::oseir::Model<> model(1); - - // set population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000; - - // set transition times - model.parameters.set>(1); - model.parameters.set>(1); - - // set contact matrix - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>().get_cont_freq_mat(); - contact_matrix[0].get_baseline().setConstant(2.7); - - //two mostly identical groups - auto model_group1 = model; - auto model_group2 = model; - - //some contact restrictions in group 1 - mio::ContactMatrixGroup& contact_matrix1 = - model_group1.parameters.get>().get_cont_freq_mat(); - contact_matrix1[0].add_damping(0.5, mio::SimulationTime(5)); - - //infection starts in group 1 - model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9990; - model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 10; - - mio::Graph>>, mio::MobilityEdge<>> g; - g.add_node(1001, model_group1, t0); - g.add_node(1002, model_group2, t0); - g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01)); - g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01)); - - auto sim = mio::make_mobility_sim(t0, dt, std::move(g)); - - sim.advance(tmax); - - return 0; -} diff --git a/cpp/examples/ode_secir_graph.cpp b/cpp/examples/ode_secir_graph.cpp new file mode 100644 index 0000000000..8296954621 --- /dev/null +++ b/cpp/examples/ode_secir_graph.cpp @@ -0,0 +1,108 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Daniel Abele, 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/config.h" +#include "ode_secir/model.h" +#include "ode_secir/infection_state.h" +#include "ode_secir/parameters.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/compartments/simulation.h" + +#include + +int main() +{ + const auto t0 = 0.; + const auto tmax = 30.; + const auto dt = 0.5; //time step of Mobility, daily Mobility every second step + + const size_t num_groups = 1; + mio::osecir::Model model(num_groups); + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}] = 10000; + model.parameters.set(0); + model.parameters.set>(0.2); + + model.parameters.get>() = 3.2; + model.parameters.get>() = 2.0; + model.parameters.get>() = 5.8; + model.parameters.get>() = 9.5; + model.parameters.get>() = 7.1; + + model.parameters.get>() = 0.1; + model.parameters.get>() = 0.7; + model.parameters.get>() = 0.09; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.45; + model.parameters.get>() = 35; + model.parameters.get>() = 0.2; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.3; + + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10)); + + //two mostly identical groups + auto model_group1 = model; + auto model_group2 = model; + //some contact restrictions in model_group1 + mio::ContactMatrixGroup& contact_matrix_m1 = + model_group1.parameters.get>(); + contact_matrix_m1[0].add_damping(0.7, mio::SimulationTime(15.)); + + //infection starts in group 1 + model_group1.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}] = 9990; + model_group1.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = 100; + + // get indices of INS and ISy compartments. + std::vector> indices_save_edges(2); + + // Reserve Space. The multiplication by 2 is necessary because we have the + // base and the confirmed compartments for each age group. + for (auto& vec : indices_save_edges) { + vec.reserve(2 * num_groups); + } + + // get indices and write them to the vector + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); ++i) { + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptoms})); + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptomsConfirmed})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptoms})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptomsConfirmed})); + } + + mio::Graph>, mio::MobilityEdge> g; + g.add_node(1001, model_group1, t0); + g.add_node(1002, model_group2, t0); + g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.1), indices_save_edges); + g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.1), indices_save_edges); + + auto sim = mio::make_mobility_sim(t0, dt, std::move(g)); + + sim.advance(tmax); + + auto& edge_1_0 = sim.get_graph().edges()[1]; + auto& results = edge_1_0.property.get_mobility_results(); + results.print_table({"Commuter INS", "Commuter ISy", "Commuter Total"}); + + return 0; +} diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp new file mode 100644 index 0000000000..0a04a2787b --- /dev/null +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -0,0 +1,331 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Daniel Abele, Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "memilio/compartments/parameter_studies.h" +#include "memilio/config.h" +#include "memilio/geography/regions.h" +#include "memilio/io/epi_data.h" +#include "memilio/io/result_io.h" +#include "memilio/io/mobility_io.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/utils/logging.h" +#include "memilio/utils/miompi.h" +#include "memilio/utils/random_number_generator.h" +#include "memilio/utils/time_series.h" +#include "ode_secir/parameters_io.h" +#include "ode_secir/parameter_space.h" +#include "boost/filesystem.hpp" +#include "memilio/utils/stl_util.h" +#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, mio::AgeGroup>& 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, mio::AgeGroup>& 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 Sars-CoV-2 for a immune-naive + * population and wild type variant. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +void set_covid_parameters(mio::osecir::Parameters& params) +{ + //times + const double timeExposedMin = 2.67; + const double timeExposedMax = 4.; + const double timeInfectedNoSymptomsMin = 1.2; + const double timeInfectedNoSymptomsMax = 2.53; + + const double timeInfectedSymptomsMin[] = {5.6255, 5.6255, 5.6646, 5.5631, 5.501, 5.465}; + const double timeInfectedSymptomsMax[] = {8.427, 8.427, 8.4684, 8.3139, 8.169, 8.085}; + const double timeInfectedSevereMin[] = {3.925, 3.925, 4.85, 6.4, 7.2, 9.}; + const double timeInfectedSevereMax[] = {6.075, 6.075, 7., 8.7, 9.8, 13.}; + const double timeInfectedCriticalMin[] = {4.95, 4.95, 4.86, 14.14, 14.4, 10.}; + const double timeInfectedCriticalMax[] = {8.95, 8.95, 8.86, 20.58, 19.8, 13.2}; + + array_assign_uniform_distribution(params.get>(), timeExposedMin, timeExposedMax); + array_assign_uniform_distribution(params.get>(), + timeInfectedNoSymptomsMin, timeInfectedNoSymptomsMax); + array_assign_uniform_distribution(params.get>(), timeInfectedSymptomsMin, + timeInfectedSymptomsMax); + array_assign_uniform_distribution(params.get>(), timeInfectedSevereMin, + timeInfectedSevereMax); + array_assign_uniform_distribution(params.get>(), timeInfectedCriticalMin, + timeInfectedCriticalMax); + + //probabilities + const double transmissionProbabilityOnContactMin[] = {0.02, 0.05, 0.05, 0.05, 0.08, 0.15}; + const double transmissionProbabilityOnContactMax[] = {0.04, 0.07, 0.07, 0.07, 0.10, 0.20}; + const double relativeTransmissionNoSymptomsMin = 1; + const double relativeTransmissionNoSymptomsMax = 1; + // The precise value between Risk* (situation under control) and MaxRisk* (situation not under control) + // depends on incidence and test and trace capacity + const double riskOfInfectionFromSymptomaticMin = 0.1; + const double riskOfInfectionFromSymptomaticMax = 0.3; + const double maxRiskOfInfectionFromSymptomaticMin = 0.3; + const double maxRiskOfInfectionFromSymptomaticMax = 0.5; + const double recoveredPerInfectedNoSymptomsMin[] = {0.2, 0.2, 0.15, 0.15, 0.15, 0.15}; + const double recoveredPerInfectedNoSymptomsMax[] = {0.3, 0.3, 0.25, 0.25, 0.25, 0.25}; + const double severePerInfectedSymptomsMin[] = {0.006, 0.006, 0.015, 0.049, 0.15, 0.20}; + const double severePerInfectedSymptomsMax[] = {0.009, 0.009, 0.023, 0.074, 0.18, 0.25}; + const double criticalPerSevereMin[] = {0.05, 0.05, 0.05, 0.10, 0.25, 0.35}; + const double criticalPerSevereMax[] = {0.10, 0.10, 0.10, 0.20, 0.35, 0.45}; + const double deathsPerCriticalMin[] = {0.00, 0.00, 0.10, 0.10, 0.30, 0.5}; + const double deathsPerCriticalMax[] = {0.10, 0.10, 0.18, 0.18, 0.50, 0.7}; + + 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); + + //sasonality + const double seasonality_min = 0.1; + const double seasonality_max = 0.3; + + assign_uniform_distribution(params.get>(), seasonality_min, seasonality_max); + + params.set(0); +} + +/** + * Set synthetic population data for testing. + * Same total populaton but different spread of infection in each county. + * @param counties parameters for each county. + */ +void set_synthetic_population_data(mio::osecir::Model& model) +{ + double nb_total_t0 = 10000, nb_exp_t0 = 2, nb_inf_t0 = 0, nb_car_t0 = 0, nb_hosp_t0 = 0, nb_icu_t0 = 0, + nb_rec_t0 = 0, nb_dead_t0 = 0; + + for (mio::AgeGroup i = 0; i < model.parameters.get_num_groups(); i++) { + model.populations[{i, mio::osecir::InfectionState::Exposed}] = nb_exp_t0; + model.populations[{i, mio::osecir::InfectionState::InfectedNoSymptoms}] = nb_car_t0; + model.populations[{i, mio::osecir::InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model.populations[{i, mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0; + model.populations[{i, mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 0; + model.populations[{i, mio::osecir::InfectionState::InfectedSevere}] = nb_hosp_t0; + model.populations[{i, mio::osecir::InfectionState::InfectedCritical}] = nb_icu_t0; + model.populations[{i, mio::osecir::InfectionState::Recovered}] = nb_rec_t0; + model.populations[{i, mio::osecir::InfectionState::Dead}] = nb_dead_t0; + model.populations.set_difference_from_group_total({i, mio::osecir::InfectionState::Susceptible}, + nb_total_t0); + } +} + +std::vector> get_indices_of_symptomatic_and_nonsymptomatic(mio::osecir::Model& model) +{ + std::vector> indices_save_edges(2); + const auto num_groups = static_cast(model.parameters.get_num_groups()); + + // Reserve Space. The multiplication by 2 is necessary because we have the + // base and the confirmed compartments for each age group. + for (auto& vec : indices_save_edges) { + vec.reserve(2 * num_groups); + } + + // get indices and write them to the vector + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); ++i) { + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptoms})); + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptomsConfirmed})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptoms})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptomsConfirmed})); + } + return indices_save_edges; +} + +/** + * Run the parameter study. + * Load a previously stored graph or create a new one from data. + * The graph is the input for the parameter study. + * A newly created graph is saved and can be reused. + * @param mode Mode for running the parameter study. + * @param data_dir data directory. Not used if mode is RunMode::Load. + * @param save_dir directory where the graph is loaded from if mode is RunMOde::Load or save to if mode is RunMode::Save. + * @param result_dir directory where all results of the parameter study will be stored. + * @param save_single_runs [Default: true] Defines if single run results are written to the disk. + * @returns any io error that occurs during reading or writing of files. + */ +int main() +{ + mio::set_log_level(mio::LogLevel::warn); + const auto num_days_sim = 30.0; + const auto num_runs = 10; + + mio::Graph, mio::MobilityParameters> params_graph; + + const int num_age_groups = 6; + mio::osecir::Model model(num_age_groups); + + // set parameters + set_covid_parameters(model.parameters); + + // set contact matrix + const auto cont_freq = 10.0; + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix( + Eigen::MatrixXd::Constant((size_t)num_age_groups, (size_t)num_age_groups, (1. / num_age_groups) * cont_freq)); + + // set population data + set_synthetic_population_data(model); + auto indices_save_edges = get_indices_of_symptomatic_and_nonsymptomatic(model); + + params_graph.add_node(1001, model); + params_graph.add_node(1002, model); + params_graph.add_node(1003, model); + + params_graph.add_edge(0, 1, + Eigen::VectorXd::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.05), + indices_save_edges); + params_graph.add_edge(1, 0, + Eigen::VectorXd::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.1), + indices_save_edges); + params_graph.add_edge(1, 2, + Eigen::VectorXd::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.15), + indices_save_edges); + params_graph.add_edge(2, 1, + Eigen::VectorXd::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.2), + indices_save_edges); + + //run parameter study + auto parameter_study = + mio::ParameterStudy>{params_graph, 0.0, num_days_sim, 0.5, size_t(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); + }, + [&](auto results_graph, auto&& run_id) { + 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(); + }); + + auto edges = std::vector>{}; + edges.reserve(results_graph.edges().size()); + std::transform(results_graph.edges().begin(), results_graph.edges().end(), std::back_inserter(edges), + [](auto&& edge) { + return edge.property.get_mobility_results(); + }); + + std::cout << "Run " << run_id << " done\n"; + return std::make_tuple(std::move(interpolated_result), std::move(params), std::move(edges)); + }); + + if (ensemble.size() > 0) { + auto ensemble_results = std::vector>>{}; + ensemble_results.reserve(ensemble.size()); + auto ensemble_params = std::vector>>{}; + ensemble_params.reserve(ensemble.size()); + auto ensemble_edges = std::vector>>{}; + ensemble_edges.reserve(ensemble.size()); + for (auto&& run : ensemble) { + ensemble_results.emplace_back(std::move(std::get<0>(run))); + ensemble_params.emplace_back(std::move(std::get<1>(run))); + ensemble_edges.emplace_back(std::move(std::get<2>(run))); + } + // create directory for results. + boost::filesystem::path results_dir("test_results"); + bool created = boost::filesystem::create_directories(results_dir); + if (created) { + mio::log_info("Directory '{:s}' was created.", results_dir.string()); + } + + auto county_ids = std::vector{1001, 1002, 1003}; + auto save_results_status = save_results(ensemble_results, ensemble_params, county_ids, results_dir, false); + auto pairs_edges = std::vector>{}; + for (auto& edge : params_graph.edges()) { + pairs_edges.push_back({county_ids[edge.start_node_idx], county_ids[edge.end_node_idx]}); + } + auto save_edges_status = save_edges(ensemble_edges, pairs_edges, "test_results", false, true); + } +} diff --git a/cpp/memilio/io/mobility_io.cpp b/cpp/memilio/io/mobility_io.cpp index 1bd9fc99d0..171f607796 100644 --- a/cpp/memilio/io/mobility_io.cpp +++ b/cpp/memilio/io/mobility_io.cpp @@ -2,6 +2,10 @@ #include "memilio/math/eigen.h" #include "memilio/utils/logging.h" +#ifdef MEMILIO_HAS_HDF5 +#include "memilio/io/hdf5_cpp.h" +#endif // MEMILIO_HAS_HDF5 + #include #include #include @@ -135,4 +139,150 @@ IOResult read_mobility_plain(const std::string& filename) return success(mobility); } +#ifdef MEMILIO_HAS_HDF5 +IOResult save_edges(const std::vector>>& ensemble_edges, + const std::vector>& pairs_edges, const fs::path& result_dir, + bool save_single_runs, bool save_percentiles) +{ + //save results and sum of results over nodes + auto ensemble_edges_sum = sum_nodes(ensemble_edges); + if (save_single_runs) { + for (size_t i = 0; i < ensemble_edges_sum.size(); ++i) { + BOOST_OUTCOME_TRY(save_edges(ensemble_edges[i], pairs_edges, + (result_dir / ("Edges_run" + std::to_string(i) + ".h5")).string())); + } + } + + if (save_percentiles) { + // make directories for percentiles + auto result_dir_p05 = result_dir / "p05"; + auto result_dir_p25 = result_dir / "p25"; + auto result_dir_p50 = result_dir / "p50"; + auto result_dir_p75 = result_dir / "p75"; + auto result_dir_p95 = result_dir / "p95"; + BOOST_OUTCOME_TRY(create_directory(result_dir_p05.string())); + BOOST_OUTCOME_TRY(create_directory(result_dir_p25.string())); + BOOST_OUTCOME_TRY(create_directory(result_dir_p50.string())); + BOOST_OUTCOME_TRY(create_directory(result_dir_p75.string())); + BOOST_OUTCOME_TRY(create_directory(result_dir_p95.string())); + + // save percentiles of Edges + { + auto ensemble_edges_p05 = ensemble_percentile(ensemble_edges, 0.05); + auto ensemble_edges_p25 = ensemble_percentile(ensemble_edges, 0.25); + auto ensemble_edges_p50 = ensemble_percentile(ensemble_edges, 0.50); + auto ensemble_edges_p75 = ensemble_percentile(ensemble_edges, 0.75); + auto ensemble_edges_p95 = ensemble_percentile(ensemble_edges, 0.95); + + BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p05, pairs_edges, (result_dir_p05 / "Edges.h5").string())); + BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p25, pairs_edges, (result_dir_p25 / "Edges.h5").string())); + BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p50, pairs_edges, (result_dir_p50 / "Edges.h5").string())); + BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p75, pairs_edges, (result_dir_p75 / "Edges.h5").string())); + BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p95, pairs_edges, (result_dir_p95 / "Edges.h5").string())); + } + } + return success(); +} + +IOResult save_edges(const std::vector>& results, const std::vector>& ids, + const std::string& filename) +{ + const auto num_edges = results.size(); + size_t edge_indx = 0; + // H5Fcreate creates a new HDF5 file. + // H5F_ACC_TRUNC: If the file already exists, H5Fcreate fails. If the file does not exist, it is created and opened with read-write access. + // H5P_DEFAULT: default data transfer properties are used. + H5File file{H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, "Failed to open the HDF5 file: " + filename); + + while (edge_indx < num_edges) { + const auto& result = results[edge_indx]; + auto h5group_name = "/" + std::to_string(ids[edge_indx].first); + H5Group start_node_h5group{H5Gcreate(file.id, h5group_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(start_node_h5group.id, StatusCode::UnknownError, + "Group could not be created (" + h5group_name + ") in the file: " + filename); + + const auto num_timepoints = result.get_num_time_points(); + if (num_timepoints > 0) { + const auto num_elements = result.get_value(0).size(); + + hsize_t dims_t[] = {static_cast(num_timepoints)}; + H5DataSpace dspace_t{H5Screate_simple(1, dims_t, NULL)}; + MEMILIO_H5_CHECK(dspace_t.id, StatusCode::UnknownError, + "Failed to create the DataSpace for 'Time' in group " + h5group_name + + " in the file: " + filename); + + H5DataSet dset_t{H5Dcreate(start_node_h5group.id, "Time", H5T_NATIVE_DOUBLE, dspace_t.id, H5P_DEFAULT, + H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(dset_t.id, StatusCode::UnknownError, + "Failed to create the 'Time' DataSet in group " + h5group_name + + " in the file: " + filename); + + auto values_t = std::vector(result.get_times().begin(), result.get_times().end()); + MEMILIO_H5_CHECK(H5Dwrite(dset_t.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, values_t.data()), + StatusCode::UnknownError, + "Failed to write 'Time' data in group " + h5group_name + " in the file: " + filename); + + int start_id = ids[edge_indx].first; + auto total = Eigen::Matrix::Zero(num_timepoints, + num_elements) + .eval(); + while (edge_indx < num_edges && ids[edge_indx].first == start_id) { + const auto& result_edge = results[edge_indx]; + auto edge_result = Eigen::Matrix::Zero( + num_timepoints, num_elements) + .eval(); + for (Eigen::Index t_idx = 0; t_idx < result_edge.get_num_time_points(); ++t_idx) { + auto v = result_edge.get_value(t_idx).transpose().eval(); + edge_result.row(t_idx) = v; + total.row(t_idx) += v; + } + + hsize_t dims_values[] = {static_cast(num_timepoints), static_cast(num_elements)}; + H5DataSpace dspace_values{H5Screate_simple(2, dims_values, NULL)}; + MEMILIO_H5_CHECK(dspace_values.id, StatusCode::UnknownError, + "Failed to create the DataSpace for End" + std::to_string(ids[edge_indx].second) + + " in group " + h5group_name + " in the file: " + filename); + + // End is the target node of the edge + auto dset_name = "End" + std::to_string(ids[edge_indx].second); + H5DataSet dset_values{H5Dcreate(start_node_h5group.id, dset_name.c_str(), H5T_NATIVE_DOUBLE, + dspace_values.id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(dset_values.id, StatusCode::UnknownError, + "Failed to create the DataSet for " + dset_name + " in group " + h5group_name + + " in the file: " + filename); + + MEMILIO_H5_CHECK( + H5Dwrite(dset_values.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, edge_result.data()), + StatusCode::UnknownError, + "Failed to write data for " + dset_name + " in group " + h5group_name + + " in the file: " + filename); + + // In the final iteration, we also save the total values + if (edge_indx == num_edges - 1 || ids[edge_indx + 1].first != start_id) { + dset_name = "Total"; + H5DataSet dset_total{H5Dcreate(start_node_h5group.id, dset_name.c_str(), H5T_NATIVE_DOUBLE, + dspace_values.id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(dset_total.id, StatusCode::UnknownError, + "Failed to create the Total DataSet in group " + h5group_name + + " in the file: " + filename); + + MEMILIO_H5_CHECK( + H5Dwrite(dset_total.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, total.data()), + StatusCode::UnknownError, + "Failed to write Total data in group " + h5group_name + " in the file: " + filename); + } + edge_indx++; + } + } + else { + log_error("No time points in the TimeSeries for Edge combination " + std::to_string(ids[edge_indx].first) + + " -> " + std::to_string(ids[edge_indx].second)); + return failure(mio::StatusCode::InvalidValue); + } + } + return success(); +} +#endif // MEMILIO_HAS_HDF5 + } // namespace mio diff --git a/cpp/memilio/io/mobility_io.h b/cpp/memilio/io/mobility_io.h index b6e840d0b9..e55fe7ae76 100644 --- a/cpp/memilio/io/mobility_io.h +++ b/cpp/memilio/io/mobility_io.h @@ -22,6 +22,7 @@ #include "memilio/io/json_serializer.h" #include "memilio/mobility/graph.h" +#include "memilio/data/analyze_result.h" #include "memilio/mobility/metapopulation_mobility_instant.h" namespace mio @@ -185,6 +186,31 @@ IOResult>> read_graph(const std::string& dir } #endif //MEMILIO_HAS_JSONCPP +#ifdef MEMILIO_HAS_HDF5 +/** + * @brief Save the results of the edges for a single graph simulation run. + * @param[in] results Simulation results per edge of the graph. + * @param[in] ids Identifiers for the start and end node of the edges. + * @param[in] filename Name of the file where the results will be saved. + * @return Any io errors that occur during writing of the files. + */ +IOResult save_edges(const std::vector>& results, const std::vector>& ids, + const std::string& filename); + +/** + * @brief Saves the results of a simulation for each edge in the graph. + * @param[in] ensemble_edges Simulation results for each run for each edge. + * @param[in] pairs_edges Identifiers for the start and end node of the edges. + * @param[in] result_dir Top level directory for all results of the parameter study. + * @param[in] save_single_runs [Default: true] Defines if single run results are written. + * @param[in] save_percentiles [Default: true] Defines if percentiles are written. + * @return Any io errors that occur during writing of the files. + */ +IOResult save_edges(const std::vector>>& ensemble_edges, + const std::vector>& pairs_edges, const fs::path& result_dir, + bool save_single_runs = true, bool save_percentiles = true); + +#endif //MEMILIO_HAS_HDF5 } // namespace mio diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 39675db0db..4a87217419 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -117,8 +117,16 @@ IOResult> read_result(const std::string& filename) MEMILIO_H5_CHECK( H5Literate(region_h5group.id, H5_INDEX_NAME, H5_ITER_INC, NULL, &store_group_name, &h5dset_names), StatusCode::UnknownError, "Dataset names could not be read."); - auto num_groups = (Eigen::Index)std::count_if(h5dset_names.begin(), h5dset_names.end(), [](auto&& str) { - return str.find("Group") != std::string::npos; + + std::string h5_key = std::any_of(h5dset_names.begin(), h5dset_names.end(), + [](const std::string& str) { + return str.find("Group") == 0; + }) + ? "Group" + : "End"; + + auto num_groups = (Eigen::Index)std::count_if(h5dset_names.begin(), h5dset_names.end(), [&h5_key](auto&& str) { + return str.find(h5_key) != std::string::npos; }); H5DataSet dataset_t{H5Dopen(region_h5group.id, "Time", H5P_DEFAULT)}; @@ -168,8 +176,17 @@ IOResult> read_result(const std::string& filename) for (Eigen::Index t_idx = 0; t_idx < num_timepoints; ++t_idx) { groups.add_time_point(time[t_idx]); } - for (Eigen::Index group_idx = 0; group_idx < num_groups; ++group_idx) { - auto group_name = "/" + h5group_name + "/Group" + std::to_string(group_idx + 1); + + std::vector h5_key_indices; + // Extract group indices from h5dset_names + for (const auto& name : h5dset_names) { + if (name.find(h5_key) == 0) { + h5_key_indices.push_back(std::stoi(name.substr(h5_key.size()))); + } + } + + for (auto h5_key_indx = size_t(0); h5_key_indx < h5_key_indices.size(); h5_key_indx++) { + auto group_name = "/" + h5group_name + "/" + h5_key + std::to_string(h5_key_indices[h5_key_indx]); H5DataSet dataset_values{H5Dopen(file.id, group_name.c_str(), H5P_DEFAULT)}; MEMILIO_H5_CHECK(dataset_values.id, StatusCode::UnknownError, "Values DataSet could not be read."); @@ -194,7 +211,7 @@ IOResult> read_result(const std::string& filename) for (Eigen::Index idx_t = 0; idx_t < num_timepoints; idx_t++) { for (Eigen::Index idx_c = 0; idx_c < num_infectionstates; idx_c++) { - groups[idx_t][num_infectionstates * group_idx + idx_c] = group_values(idx_t, idx_c); + groups[idx_t][num_infectionstates * h5_key_indx + idx_c] = group_values(idx_t, idx_c); } } } diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index b6d573497a..7d9555c067 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -331,13 +331,14 @@ IOResult set_nodes(const Parameters& params, Date start_date, Date end_dat * @param[in] contact_locations_size Number of contact locations. * @param[in] read_func Function that reads commuting matrices. * @param[in] commuting_weights Vector with a commuting weight for every AgeGroup. + * @param[in] indices_of_saved_edges Vector of vectors with indices of the compartments that should be saved on the edges. */ template IOResult set_edges(const fs::path& data_dir, Graph& params_graph, std::initializer_list& mobile_compartments, size_t contact_locations_size, - ReadFunction&& read_func, - std::vector commuting_weights = std::vector{}) + ReadFunction&& read_func, std::vector commuting_weights, + std::vector> indices_of_saved_edges = {}) { // mobility between nodes BOOST_OUTCOME_TRY(auto&& mobility_data_commuter, @@ -390,7 +391,7 @@ IOResult set_edges(const fs::path& data_dir, Graph& //only add edges with mobility above thresholds for performance //thresholds are chosen empirically so that more than 99% of mobility is covered, approx. 1/3 of the edges if (commuter_coeff_ij > 4e-5 || twitter_coeff > 1e-5) { - params_graph.add_edge(county_idx_i, county_idx_j, std::move(mobility_coeffs)); + params_graph.add_edge(county_idx_i, county_idx_j, std::move(mobility_coeffs), indices_of_saved_edges); } } } diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp index 700ebb2a79..fc54caebec 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp @@ -21,5 +21,4 @@ namespace mio { - } // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index d47e655e51..34996710d9 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -129,6 +129,7 @@ class MobilityParameters */ MobilityParameters(const MobilityCoefficientGroup& coeffs) : m_coefficients(coeffs) + , m_saved_compartment_indices(0) { } @@ -138,6 +139,38 @@ class MobilityParameters */ MobilityParameters(const Eigen::VectorXd& coeffs) : m_coefficients({MobilityCoefficients(coeffs)}) + , m_saved_compartment_indices(0) + { + } + + /** + * @brief Constructor for initializing mobility parameters with coefficients from type `MobilityCoefficientGroup` + * and specific save indices. + * + * @param[in] coeffs A group of mobility coefficients represented by a `MobilityCoefficientGroup` object, defining + * how individuals move between nodes. + * @param[in] save_indices A 2D vector of indices. The outer vector represents different sets of compartments. + * Each inner vector represents a set of compartments whose data will be saved in the member `m_mobility_results` + * of the `MobilityEdge` class during the simulation using the `add_mobility_result_time_point` function. + */ + MobilityParameters(const MobilityCoefficientGroup& coeffs, const std::vector>& save_indices) + : m_coefficients(coeffs) + , m_saved_compartment_indices(save_indices) + { + } + + /** + * @brief Constructor for initializing mobility parameters with coefficients from an Eigen Vector + * and specific save indices. + * + * @param[in] coeffs An `Eigen::VectorXd` containing mobility coefficients. + * @param[in] save_indices A 2D vector of indices. The outer vector represents different sets of compartments. + * Each inner vector represents a set of compartments whose data will be saved in the member `m_mobility_results` + * of the `MobilityEdge` class during the simulation using the `add_mobility_result_time_point` function. + */ + MobilityParameters(const Eigen::VectorXd& coeffs, const std::vector>& save_indices) + : m_coefficients({MobilityCoefficients(coeffs)}) + , m_saved_compartment_indices(save_indices) { } @@ -179,7 +212,21 @@ class MobilityParameters { m_coefficients = coeffs; } - /** @} */ + + /** + * @brief Get the indices of compartments to be saved during mobility. + * + * This function returns a reference to the vector of `m_saved_compartment_indices`, which specifies the groups of + * compartments that are saved in the member `m_mobility_results` of the `MobilityEdge` class during the simulation + * using the `add_mobility_result_time_point` function. + * + * @return A reference to the 2D vector containing indices of compartments to be saved. The outer vector represents different sets of compartments. + * Each inner vector represents a group of compartments defined by indices. + */ + const auto& get_save_indices() const + { + return m_saved_compartment_indices; + } /** * Get/Set dynamic NPIs that are implemented when relative infections exceed thresholds. @@ -241,6 +288,7 @@ class MobilityParameters private: MobilityCoefficientGroup m_coefficients; //one per group and compartment DynamicNPIs m_dynamic_npis; + std::vector> m_saved_compartment_indices; // groups of indices from compartments to save }; /** @@ -251,7 +299,7 @@ class MobilityEdge { public: /** - * create edge with coefficients. + * @brief Create edge with coefficients. * @param coeffs % of people in each group and compartment that change node in each time step. */ MobilityEdge(const MobilityParameters& params) @@ -259,18 +307,60 @@ class MobilityEdge , m_mobile_population(params.get_coefficients().get_shape().rows()) , m_return_times(0) , m_return_mobile_population(false) + , m_saved_compartment_indices(params.get_save_indices()) + , m_mobility_results(m_saved_compartment_indices.size() + 1) { } /** - * create edge with coefficients. - * @param coeffs % of people in each group and compartment that change node in each time step. + * @brief Create edge with coefficients. + * + * @param[in] coeffs An `Eigen::VectorXd` representing the percentage of people in each group and compartment + * that change nodes in each time step. */ MobilityEdge(const Eigen::VectorXd& coeffs) : m_parameters(coeffs) , m_mobile_population(coeffs.rows()) , m_return_times(0) , m_return_mobile_population(false) + , m_saved_compartment_indices(0) + , m_mobility_results(m_saved_compartment_indices.size() + 1) + { + } + + /** + * @brief Create edge with coefficients as MobilityParameters object and a 2D vector of indices which determine which compartments are saved. + * + * @param[in] params A `MobilityParameters` object representing the percentage of people in each group and compartment + * that change nodes in each time step. + * @param[in] save_indices A 2D vector of indices. The outer vector represents different sets of compartments. + * Each inner vector represents a group of indices to be saved. + */ + MobilityEdge(const MobilityParameters& params, const std::vector>& save_indices) + : m_parameters(params) + , m_mobile_population(params.get_coefficients().get_shape().rows()) + , m_return_times(0) + , m_return_mobile_population(false) + , m_saved_compartment_indices(save_indices) + , m_mobility_results(m_saved_compartment_indices.size() + 1) + { + } + + /** + * @brief Create edge with coefficients and a 2D vector of indices which determine which compartments are saved. + * + * @param[in] coeffs An `Eigen::VectorXd` representing the percentage of people in each group and compartment that migrate + * in each time step. + * @param[in] save_indices A 2D vector of indices. The outer vector represents different sets of compartments, while each + * inner vector represents a group of indices for compartments to be saved. + */ + MobilityEdge(const Eigen::VectorXd& coeffs, const std::vector>& save_indices) + : m_parameters(coeffs) + , m_mobile_population(coeffs.rows()) + , m_return_times(0) + , m_return_mobile_population(false) + , m_saved_compartment_indices(save_indices) + , m_mobility_results(m_saved_compartment_indices.size() + 1) { } @@ -282,6 +372,22 @@ class MobilityEdge return m_parameters; } + /** + * @brief Get the count of commuters in selected compartments, along with the total number of commuters. + * + * @return A reference to the TimeSeries object representing the mobility results. + * @{ + */ + TimeSeries& get_mobility_results() + { + return m_mobility_results; + } + const TimeSeries& get_mobility_results() const + { + return m_mobility_results; + } + /** @} */ + /** * compute mobility from node_from to node_to. * mobility is based on coefficients. @@ -302,8 +408,46 @@ class MobilityEdge bool m_return_mobile_population; double m_t_last_dynamic_npi_check = -std::numeric_limits::infinity(); std::pair m_dynamic_npi = {-std::numeric_limits::max(), SimulationTime(0)}; + std::vector> m_saved_compartment_indices; // groups of indices from compartments to save + TimeSeries m_mobility_results; // save results from edges + entry for the total number of commuters + + /** + * @brief Computes a condensed version of `m_mobile_population` and stores it in `m_mobility_results`. + * + * The `m_mobility_results` then only contains commuters with infection states `InfectedNoSymptoms` and + * `InfectedSymptoms`. Additionally, the total number of commuters is stored in the last entry of `m_mobility_results`. + * + * @param[in] t The current time. + */ + void add_mobility_result_time_point(const double t); }; +template +void MobilityEdge::add_mobility_result_time_point(const double t) +{ + const size_t save_indices_size = this->m_saved_compartment_indices.size(); + if (save_indices_size > 0) { + + const auto& last_value = m_mobile_population.get_last_value(); + Eigen::VectorXd condensed_values = Eigen::VectorXd::Zero(save_indices_size + 1); + + // sum up the values of m_saved_compartment_indices for each group (e.g. Age groups) + std::transform(this->m_saved_compartment_indices.begin(), this->m_saved_compartment_indices.end(), + condensed_values.data(), [&last_value](const auto& indices) { + return std::accumulate(indices.begin(), indices.end(), 0.0, + [&last_value](double sum, auto i) { + return sum + last_value[i]; + }); + }); + + // the last value is the sum of commuters + condensed_values[save_indices_size] = m_mobile_population.get_last_value().sum(); + + // Move the condensed values to the m_mobility_results time series + m_mobility_results.add_time_point(t, std::move(condensed_values)); + } +} + /** * adjust number of people that changed node 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. @@ -391,8 +535,8 @@ auto get_mobility_factors(const SimulationNode& node, double t, const Eigen * detect a get_mobility_factors function for the Model type. */ template -using test_commuters_expr_t = decltype(test_commuters( - std::declval(), std::declval&>(), std::declval())); +using test_commuters_expr_t = decltype( + test_commuters(std::declval(), std::declval&>(), std::declval())); /** * Test persons when moving from their source node. @@ -472,6 +616,7 @@ void MobilityEdge::apply_mobility(FP t, FP dt, SimulationNode& node_fro } node_from.get_result().get_last_value() += m_mobile_population[i]; node_to.get_result().get_last_value() -= m_mobile_population[i]; + add_mobility_result_time_point(t); m_mobile_population.remove_time_point(i); m_return_times.remove_time_point(i); } @@ -489,6 +634,8 @@ void MobilityEdge::apply_mobility(FP t, FP dt, SimulationNode& node_fro node_to.get_result().get_last_value() += m_mobile_population.get_last_value(); node_from.get_result().get_last_value() -= m_mobile_population.get_last_value(); + + add_mobility_result_time_point(t); } m_return_mobile_population = !m_return_mobile_population; } diff --git a/cpp/simulations/2020_npis_sarscov2_wildtype_germany.cpp b/cpp/simulations/2020_npis_sarscov2_wildtype_germany.cpp index 57c905e88f..25f188cb35 100644 --- a/cpp/simulations/2020_npis_sarscov2_wildtype_germany.cpp +++ b/cpp/simulations/2020_npis_sarscov2_wildtype_germany.cpp @@ -475,8 +475,8 @@ get_graph(mio::Date start_date, mio::Date end_date, const fs::path& data_dir) auto scaling_factor_icu = 1.0; auto tnt_capacity_factor = 7.5 / 100000.; auto mobile_compartments = {mio::osecir::InfectionState::Susceptible, mio::osecir::InfectionState::Exposed, - mio::osecir::InfectionState::InfectedNoSymptoms, - mio::osecir::InfectionState::InfectedSymptoms, mio::osecir::InfectionState::Recovered}; + mio::osecir::InfectionState::InfectedNoSymptoms, + mio::osecir::InfectionState::InfectedSymptoms, mio::osecir::InfectionState::Recovered}; // graph of counties with populations and local parameters // and mobility between counties @@ -498,7 +498,8 @@ get_graph(mio::Date start_date, mio::Date end_date, const fs::path& data_dir) true, params_graph, read_function_nodes, node_id_function, scaling_factor_infected, scaling_factor_icu, tnt_capacity_factor, 0, false, true)); BOOST_OUTCOME_TRY(set_edge_function(data_dir, params_graph, mobile_compartments, contact_locations.size(), - read_function_edges, std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.})); + read_function_edges, std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}, + {})); return mio::success(params_graph); } @@ -574,7 +575,7 @@ mio::IOResult run(RunMode mode, const fs::path& data_dir, const fs::path& 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) { + [](auto&& node) { return node.property.get_simulation().get_model(); }); diff --git a/cpp/simulations/2021_vaccination_sarscov2_delta_germany.cpp b/cpp/simulations/2021_vaccination_sarscov2_delta_germany.cpp index 070d0b24d6..377890a716 100644 --- a/cpp/simulations/2021_vaccination_sarscov2_delta_germany.cpp +++ b/cpp/simulations/2021_vaccination_sarscov2_delta_germany.cpp @@ -548,17 +548,17 @@ get_graph(mio::Date start_date, mio::Date end_date, const fs::path& data_dir, bo auto scaling_factor_icu = 1.0; auto tnt_capacity_factor = 1.43 / 100000.; auto mobile_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}; + 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}; // graph of counties with populations and local parameters // and mobility between counties @@ -580,7 +580,8 @@ get_graph(mio::Date start_date, mio::Date end_date, const fs::path& data_dir, bo 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, true)); BOOST_OUTCOME_TRY(set_edge_function(data_dir, params_graph, mobile_compartments, contact_locations.size(), - read_function_edges, std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.})); + read_function_edges, std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}, + {})); return mio::success(params_graph); } @@ -662,7 +663,7 @@ mio::IOResult run(RunMode mode, const fs::path& data_dir, const fs::path& 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) { + [](auto&& node) { return node.property.get_simulation().get_model(); }); diff --git a/cpp/tests/test_graph.cpp b/cpp/tests/test_graph.cpp index 9cfb33c528..2ba785ffea 100644 --- a/cpp/tests/test_graph.cpp +++ b/cpp/tests/test_graph.cpp @@ -32,8 +32,10 @@ #include "models/ode_secirvvs/model.h" #include "memilio/io/io.h" #include "matchers.h" +#include "temp_file_register.h" #include "memilio/utils/stl_util.h" #include "gmock/gmock-matchers.h" +#include #include #include #include @@ -272,6 +274,60 @@ TEST(TestGraph, set_edges) MatrixNear(print_wrap(e_other.array().cast()), 1e-5, 1e-5)); } +TEST(TestGraph, set_edges_saving_edges) +{ + const size_t num_groups = 6; + mio::osecir::Model model(num_groups); + model.populations[{mio::AgeGroup(3), mio::osecir::InfectionState::Exposed}] = 1; + mio::Graph, mio::MobilityParameters> params_graph; + TempFileRegister file_register; + const auto dir = file_register.get_unique_path("TestGraph-%%%%-%%%%"); + + auto mobile_compartments = {mio::osecir::InfectionState::Susceptible, mio::osecir::InfectionState::Exposed, + mio::osecir::InfectionState::InfectedNoSymptoms, + mio::osecir::InfectionState::InfectedSymptoms, mio::osecir::InfectionState::Recovered}; + + // get indices of INS and ISy compartments. + std::vector> indices_save_edges(2); + + // Reserve Space. The multiplication by 2 is necessary because we have the + // base and the confirmed compartments for each age group. + for (auto& vec : indices_save_edges) { + vec.reserve(2 * num_groups); + } + + // get indices and write them to the vector + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); ++i) { + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptoms})); + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptomsConfirmed})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptoms})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptomsConfirmed})); + } + + params_graph.add_node(1001, model); + params_graph.add_node(1002, model); + + const auto& read_function_edges = mock_read_mobility; + + auto result = + mio::set_edges, mio::MobilityParameters, + mio::MobilityCoefficientGroup, mio::osecir::InfectionState, decltype(read_function_edges)>( + dir, params_graph, mobile_compartments, size_t(2), read_function_edges, + std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}, indices_save_edges); + + EXPECT_EQ(params_graph.edges().size(), 2); + + auto& indices_edge0 = params_graph.edges()[0].property.get_save_indices(); + auto& indices_edge1 = params_graph.edges()[1].property.get_save_indices(); + + EXPECT_EQ(indices_edge0, indices_save_edges); + EXPECT_EQ(indices_edge1, indices_save_edges); +} + TEST(TestGraph, ot_edges) { mio::Graph g; @@ -296,10 +352,10 @@ namespace struct MoveOnly { MoveOnly(); - MoveOnly(const MoveOnly&) = delete; + MoveOnly(const MoveOnly&) = delete; MoveOnly& operator=(const MoveOnly&) = delete; MoveOnly(MoveOnly&&) = default; - MoveOnly& operator=(MoveOnly&&) = default; + MoveOnly& operator=(MoveOnly&&) = default; }; using MoveOnlyGraph = mio::Graph; diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp index 7460c776ae..8accf2fc7e 100644 --- a/cpp/tests/test_mobility.cpp +++ b/cpp/tests/test_mobility.cpp @@ -166,3 +166,71 @@ TEST(TestMobility, edgeApplyMobility) EXPECT_DOUBLE_EQ(node1.get_result().get_last_value().sum(), 900); EXPECT_DOUBLE_EQ(node2.get_result().get_last_value().sum(), 1100); } + +TEST(TestMobility, add_mobility_result_time_point) +{ + using Model = mio::osecir::Model; + + //setup nodes + const size_t num_groups = 1; + Model model(num_groups); + auto& params = model.parameters; + auto& cm = static_cast(model.parameters.get>()); + cm[0].get_baseline()(0, 0) = 5.0; + + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = 10; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = 20; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 0; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, 1000); + params.get>()[(mio::AgeGroup)0] = 1.; + params.get>()[(mio::AgeGroup)0] = 1.; + params.get>()[(mio::AgeGroup)0] = 1.; + params.get>()[(mio::AgeGroup)0] = 0.5; + params.get>()[(mio::AgeGroup)0] = 1.; + params.get>()[(mio::AgeGroup)0] = 1.; + params.apply_constraints(); + + // get indices of INS and ISy compartments. + std::vector> indices_save_edges(2); + + // Reserve Space. The multiplication by 2 is necessary because we have the + // base and the confirmed compartments for each age group. + for (auto& vec : indices_save_edges) { + vec.reserve(2 * num_groups); + } + + // get indices and write them to the vector + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); ++i) { + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptoms})); + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptomsConfirmed})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptoms})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptomsConfirmed})); + } + + //setup different edges + double t = 0.; + mio::SimulationNode> node1(model, t); + mio::SimulationNode> node2(model, t); + mio::MobilityEdge edge1(Eigen::VectorXd::Constant(10, 0.1), indices_save_edges); + edge1.apply_mobility(t, 0.0, node1, node2); + auto mobility = edge1.get_mobility_results().get_last_value(); + EXPECT_NEAR(mobility[0], 1.0, 1e-12); + EXPECT_NEAR(mobility[1], 2.0, 1e-12); + EXPECT_NEAR(mobility[2], 100.0, 1e-12); + + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptomsConfirmed}] = 100; + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 30; + mio::SimulationNode> node3(model, t); + mio::SimulationNode> node4(model, t); + mio::MobilityEdge edge2(Eigen::VectorXd::Constant(10, 0.1), indices_save_edges); + edge2.apply_mobility(t, 0.5, node3, node4); + mobility = edge2.get_mobility_results().get_last_value(); + EXPECT_NEAR(mobility[0], 11.0, 1e-12); + EXPECT_NEAR(mobility[1], 5.0, 1e-12); + EXPECT_NEAR(mobility[2], 113.0, 1e-12); +} \ No newline at end of file diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index 3e7d715ce7..5710c3c879 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -231,7 +231,7 @@ TEST(TestSaveResult, save_percentiles_and_sums) double num_total_t0 = 10000, num_exp_t0 = 100, num_inf_t0 = 50, num_car_t0 = 50, num_hosp_t0 = 20, num_icu_t0 = 10, num_rec_t0 = 10, num_dead_t0 = 0; - size_t num_groups = 3; + const size_t num_groups = 3; mio::osecir::Model model((int)num_groups); double fact = 1.0 / (double)num_groups; @@ -270,13 +270,35 @@ TEST(TestSaveResult, save_percentiles_and_sums) mio::ContactMatrixGroup& contact_matrix = params.get>(); contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_groups, num_groups, fact * cont_freq)); + // get indices of INS and ISy compartments. + std::vector> indices_save_edges(2); + + // Reserve Space. The multiplication by 2 is necessary because we have the + // base and the confirmed compartments for each age group. + for (auto& vec : indices_save_edges) { + vec.reserve(2 * num_groups); + } + + // get indices and write them to the vector + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); ++i) { + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptoms})); + indices_save_edges[0].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedNoSymptomsConfirmed})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptoms})); + indices_save_edges[1].emplace_back( + model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptomsConfirmed})); + } + mio::osecir::set_params_distributions_normal(model, t0, tmax, 0.2); auto graph = mio::Graph, mio::MobilityParameters>(); graph.add_node(0, model); graph.add_node(1, model); graph.add_edge(0, 1, - mio::MobilityParameters(Eigen::VectorXd::Constant(Eigen::Index(num_groups * 10), 1.0))); + mio::MobilityParameters(Eigen::VectorXd::Constant(Eigen::Index(num_groups * 10), 1.0), + indices_save_edges)); auto num_runs = 3; auto parameter_study = mio::ParameterStudy>(graph, 0.0, 2.0, 0.5, num_runs); @@ -290,6 +312,8 @@ TEST(TestSaveResult, save_percentiles_and_sums) ensemble_results.reserve(size_t(num_runs)); auto ensemble_params = std::vector>>{}; ensemble_params.reserve(size_t(num_runs)); + auto ensemble_edges = std::vector>>{}; + ensemble_edges.reserve(size_t(num_runs)); parameter_study.run( [](auto&& g) { return draw_sample(g); @@ -303,6 +327,14 @@ TEST(TestSaveResult, save_percentiles_and_sums) std::back_inserter(ensemble_params.back()), [](auto&& node) { return node.property.get_simulation().get_model(); }); + + ensemble_edges.emplace_back(); + ensemble_edges.back().reserve(results_graph.edges().size()); + std::transform(results_graph.edges().begin(), results_graph.edges().end(), + std::back_inserter(ensemble_edges.back()), [](auto&& edge) { + return edge.property.get_mobility_results(); + }); + return 0; //function needs to return something }); @@ -339,4 +371,127 @@ TEST(TestSaveResult, save_percentiles_and_sums) ASSERT_TRUE(results_run2); auto results_run2_sum = mio::read_result(tmp_results_dir + "/results_run2_sum.h5"); ASSERT_TRUE(results_run2_sum); + + // test save edges (percentiles and results from single runs) + std::vector> pairs_edges = {{0, 1}}; + + auto save_edges_status = save_edges(ensemble_edges, pairs_edges, tmp_results_dir, true, true); + ASSERT_TRUE(save_edges_status); + + // percentiles + auto results_edges_from_file_p05 = mio::read_result(tmp_results_dir + "/p05/Edges.h5"); + ASSERT_TRUE(results_edges_from_file_p05); + auto results_edges_from_file_p25 = mio::read_result(tmp_results_dir + "/p25/Edges.h5"); + ASSERT_TRUE(results_edges_from_file_p25); + auto results_edges_from_file_p50 = mio::read_result(tmp_results_dir + "/p50/Edges.h5"); + ASSERT_TRUE(results_edges_from_file_p50); + auto results_edges_from_file_p75 = mio::read_result(tmp_results_dir + "/p75/Edges.h5"); + ASSERT_TRUE(results_edges_from_file_p75); + auto results_edges_from_file_p95 = mio::read_result(tmp_results_dir + "/p95/Edges.h5"); + ASSERT_TRUE(results_edges_from_file_p95); + + auto result_edges_from_file = results_edges_from_file_p25.value()[0]; + EXPECT_EQ(ensemble_edges.back().back().get_num_elements(), result_edges_from_file.get_groups().get_num_elements()); + EXPECT_EQ(ensemble_edges.back().back().get_num_time_points(), + result_edges_from_file.get_groups().get_num_time_points()); + + // single runs + auto results_edges_run0 = mio::read_result(tmp_results_dir + "/Edges_run0.h5"); + ASSERT_TRUE(results_edges_run0); + auto results_edges_run1 = mio::read_result(tmp_results_dir + "/Edges_run1.h5"); + ASSERT_TRUE(results_edges_run1); + auto results_edges_run2 = mio::read_result(tmp_results_dir + "/Edges_run2.h5"); + ASSERT_TRUE(results_edges_run2); +} + +TEST(TestSaveResult, save_edges) +{ + // create some results and pairs_edges + const auto n = Eigen::Index(3); + std::vector> results_edges(3, mio::TimeSeries(n)); + results_edges[0].add_time_point(0.0, Eigen::VectorXd::Constant(n, 0)); + results_edges[0].add_time_point(1.0, Eigen::VectorXd::Constant(n, 1)); + results_edges[0].add_time_point(2.0, Eigen::VectorXd::Constant(n, 1)); + + results_edges[1].add_time_point(0.0, Eigen::VectorXd::Constant(n, 2)); + results_edges[1].add_time_point(1.0, Eigen::VectorXd::Constant(n, 3)); + results_edges[1].add_time_point(2.0, Eigen::VectorXd::Constant(n, 5)); + + results_edges[2].add_time_point(0.0, Eigen::VectorXd::Constant(n, 3)); + results_edges[2].add_time_point(1.0, Eigen::VectorXd::Constant(n, 4)); + results_edges[2].add_time_point(2.0, Eigen::VectorXd::Constant(n, 7)); + + const std::vector> pairs_edges = {{0, 1}, {0, 2}, {1, 2}}; + + // save the results to a file + TempFileRegister file_register; + auto results_file_path = file_register.get_unique_path("test_result-%%%%-%%%%.h5"); + auto save_edges_status = mio::save_edges(results_edges, pairs_edges, results_file_path); + ASSERT_TRUE(save_edges_status); + + // read the results back in and check that they are correct. + auto results_from_file = mio::read_result(results_file_path); + ASSERT_TRUE(results_from_file); + + // group 0 + auto result_from_file_group0 = results_from_file.value()[0]; + EXPECT_EQ(result_from_file_group0.get_groups().get_num_time_points(), 3); + EXPECT_EQ(result_from_file_group0.get_groups().get_num_elements(), 6); + EXPECT_EQ(result_from_file_group0.get_groups().get_value(0), (Eigen::VectorXd(6) << 0, 0, 0, 2, 2, 2).finished()); + EXPECT_EQ(result_from_file_group0.get_groups().get_value(1), (Eigen::VectorXd(6) << 1, 1, 1, 3, 3, 3).finished()); + EXPECT_EQ(result_from_file_group0.get_groups().get_value(2), (Eigen::VectorXd(6) << 1, 1, 1, 5, 5, 5).finished()); + EXPECT_EQ(result_from_file_group0.get_groups().get_time(0), 0.0); + EXPECT_EQ(result_from_file_group0.get_groups().get_time(1), 1.0); + EXPECT_EQ(result_from_file_group0.get_groups().get_time(2), 2.0); + + EXPECT_EQ(result_from_file_group0.get_totals().get_num_time_points(), 3); + EXPECT_EQ(result_from_file_group0.get_totals().get_num_elements(), 3); + EXPECT_EQ(result_from_file_group0.get_totals().get_value(0), Eigen::VectorXd::Constant(3, 2)); + EXPECT_EQ(result_from_file_group0.get_totals().get_value(1), Eigen::VectorXd::Constant(3, 4)); + EXPECT_EQ(result_from_file_group0.get_totals().get_value(2), Eigen::VectorXd::Constant(3, 6)); + EXPECT_EQ(result_from_file_group0.get_totals().get_time(0), 0.0); + EXPECT_EQ(result_from_file_group0.get_totals().get_time(1), 1.0); + EXPECT_EQ(result_from_file_group0.get_totals().get_time(2), 2.0); + + // group 1 + auto result_from_file_group1 = results_from_file.value()[1]; + EXPECT_EQ(result_from_file_group1.get_groups().get_num_elements(), 3); + EXPECT_EQ(result_from_file_group1.get_groups().get_value(0), Eigen::VectorXd::Constant(3, 3)); + EXPECT_EQ(result_from_file_group1.get_groups().get_value(1), Eigen::VectorXd::Constant(3, 4)); + EXPECT_EQ(result_from_file_group1.get_groups().get_value(2), Eigen::VectorXd::Constant(3, 7)); + EXPECT_EQ(result_from_file_group1.get_groups().get_time(0), 0.0); + EXPECT_EQ(result_from_file_group1.get_groups().get_time(1), 1.0); + EXPECT_EQ(result_from_file_group1.get_groups().get_time(2), 2.0); + + EXPECT_EQ(result_from_file_group1.get_totals().get_num_elements(), 3); + EXPECT_EQ(result_from_file_group1.get_totals().get_value(0), Eigen::VectorXd::Constant(3, 3)); + EXPECT_EQ(result_from_file_group1.get_totals().get_value(1), Eigen::VectorXd::Constant(3, 4)); + EXPECT_EQ(result_from_file_group1.get_totals().get_value(2), Eigen::VectorXd::Constant(3, 7)); +} + +TEST(TestSaveEdges, save_edges_empty_ts) +{ + mio::set_log_level(mio::LogLevel::off); + std::vector> results; + + const auto num_elements = 2; + + // Add filled TimeSeries to the results vector + mio::TimeSeries ts_1(num_elements); + ts_1.add_time_point(0.0, Eigen::VectorXd::Constant(num_elements, 1)); + ts_1.add_time_point(1.0, Eigen::VectorXd::Constant(num_elements, 2)); + results.push_back(ts_1); + + // Add an empty TimeSeries and add it to the results vector + mio::TimeSeries ts_2(num_elements); + results.push_back(ts_2); + + const std::vector> pairs_edges = {{0, 1}, {1, 2}}; + + // Create a TempFile for HDF5 output + TempFileRegister file_register; + auto results_file_path = file_register.get_unique_path("TestEdges-%%%%-%%%%.h5"); + + // Call the save_edges function and check if it returns a failure + ASSERT_THAT(save_edges(results, pairs_edges, results_file_path), IsFailure(mio::StatusCode::InvalidValue)); }