From 6b76fa0c6a315225155a2827f58e13a50024c756 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Tue, 19 Mar 2024 11:16:32 +0100 Subject: [PATCH 01/30] save mobility data. First draft --- .../metapopulation_mobility_instant.cpp | 21 ++++++ .../metapopulation_mobility_instant.h | 68 ++++++++++++++++++- cpp/models/ode_secir/model.h | 21 ++++++ cpp/models/ode_secirvvs/model.h | 36 ++++++++++ 4 files changed, 144 insertions(+), 2 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp index 700ebb2a79..03abf61850 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp @@ -21,5 +21,26 @@ namespace mio { +void MigrationEdge::condense_m_migrated(const double t, const std::vector& indices_non_symptomatic, + const std::vector& indices_symptomatic) +{ + + const auto& last_value = m_migrated.get_last_value(); + + auto num_INS = + std::accumulate(indices_non_symptomatic.begin(), indices_non_symptomatic.end(), 0, [&](auto sum, auto i) { + return sum + last_value[i]; + }); + + auto num_ISy = std::accumulate(indices_symptomatic.begin(), indices_symptomatic.end(), 0, [&](auto sum, auto i) { + return sum + last_value[i]; + }); + + double total_commuters = m_migrated.get_last_value().sum(); + + // as time point t which contains now the carriers, infectious and total over age groups + m_num_migrated.add_time_point(t, + (mio::TimeSeries::Vector(3) << num_INS, num_ISy, total_commuters).finished()); +} } // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 5178315036..e49a258432 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -241,6 +241,41 @@ class MigrationParameters DynamicNPIs m_dynamic_npis; }; +/** + * detect a get_indices_of_symptomatic_and_nonsymptomatic function for the Model type. + */ +template +using get_indices_of_symptomatic_and_nonsymptomatic_expr_t = + decltype(get_indices_of_symptomatic_and_nonsymptomatic(std::declval())); + +/** + * @brief Get the indices of symptomatic and non-symptomatic infection states. + * + * This function generates two vectors of indices, one for non-symptomatic infection states and one for symptomatic infection states. + * Each vector contains the flat indices of the corresponding infection states for each age group and compartment in the model. + * + * @tparam Base The base class for the simulation, defaults to mio::Simulation. + * @param[in] sim The simulation object from which we obtain the model. + * + * @return A tuple containing two vectors of size_t. The first vector contains the indices of non-symptomatic infection states, + * and the second vector contains the indices of symptomatic infection states. The indices are ordered first by age group, then by compartment. + */ +template ::value, + void*> = nullptr> +auto get_indices_of_symptomatic_and_nonsymptomatic(SimulationNode& /*node*/) +{ + return std::make_tuple(std::vector{}, std::vector{}); +} + +template ::value, + void*> = nullptr> +auto get_indices_of_symptomatic_and_nonsymptomatic(SimulationNode& node) +{ + return get_indices_of_symptomatic_and_nonsymptomatic(node.get_simulation()); +} + /** * represents the migration between two nodes. */ @@ -256,6 +291,7 @@ class MigrationEdge , m_migrated(params.get_coefficients().get_shape().rows()) , m_return_times(0) , m_return_migrated(false) + , m_num_migrated(3) { } @@ -268,6 +304,7 @@ class MigrationEdge , m_migrated(coeffs.rows()) , m_return_times(0) , m_return_migrated(false) + , m_num_migrated(3) { } @@ -279,6 +316,19 @@ class MigrationEdge return m_parameters; } + /** + * Retrieve the count of commuters in the infection states: InfectedNoSymptoms and InfectedSymptomsNaive, + * along with the total number of commuter. + */ + TimeSeries& get_migrated() + { + return m_num_migrated; + } + const TimeSeries& get_migrated() const + { + return m_num_migrated; + } + /** * compute migration from node_from to node_to. * migration is based on coefficients. @@ -299,6 +349,16 @@ class MigrationEdge 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)}; + TimeSeries m_num_migrated; + + /** + * Computes a condensed version of m_migrated and puts it in m_num_migrated. + * m_num_migrated then only contains commuters with infection states InfectedNoSymptoms and InfectedSymptoms. + * Additionally, the total number of commuters is stored in the last entry of m_num_migrated. + * @param[in] t current time + */ + void condense_m_migrated(const double t, const std::vector& indices_non_symptomatic, + const std::vector& indices_symptomatic); }; /** @@ -389,8 +449,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. @@ -484,6 +544,10 @@ void MigrationEdge::apply_migration(double t, double dt, SimulationNode& no node_to.get_result().get_last_value() += m_migrated.get_last_value(); node_from.get_result().get_last_value() -= m_migrated.get_last_value(); + + static auto indices_tuple = get_indices_of_symptomatic_and_nonsymptomatic(node_from); + auto& [indices_no_symptoms, indices_symptoms] = indices_tuple; + condense_m_migrated(t, indices_no_symptoms, indices_symptoms); } m_return_migrated = !m_return_migrated; } diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h index 62c8325785..0671c833d6 100755 --- a/cpp/models/ode_secir/model.h +++ b/cpp/models/ode_secir/model.h @@ -29,6 +29,7 @@ #include "memilio/math/smoother.h" #include "memilio/math/eigen_util.h" #include "memilio/math/interpolation.h" +#include namespace mio { @@ -669,6 +670,26 @@ auto test_commuters(Simulation& sim, Eigen::Ref migrated, } } +template > +auto get_indices_of_symptomatic_and_nonsymptomatic(Simulation& sim) +{ + const auto& model = sim.get_model(); + const auto num_groups = model.parameters.get_num_groups(); + std::vector indices_no_symptoms(2 * size_t(num_groups)); + std::vector indices_symptoms(2 * size_t(num_groups)); + + for (auto i = AgeGroup(0); i < num_groups; ++i) { + indices_no_symptoms.emplace_back(model.populations.get_flat_index({i, InfectionState::InfectedNoSymptoms})); + indices_no_symptoms.emplace_back( + model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsConfirmed})); + + indices_symptoms.emplace_back(model.populations.get_flat_index({i, InfectionState::InfectedSymptoms})); + indices_symptoms.emplace_back(model.populations.get_flat_index({i, InfectionState::InfectedSymptomsConfirmed})); + } + + return std::make_tuple(std::move(indices_no_symptoms), std::move(indices_symptoms)); +} + } // namespace osecir } // namespace mio diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h index 1c9c9ae2e6..5e15eb3561 100644 --- a/cpp/models/ode_secirvvs/model.h +++ b/cpp/models/ode_secirvvs/model.h @@ -871,6 +871,42 @@ auto test_commuters(Simulation& sim, Eigen::Ref migrated, } } +template > +auto get_indices_of_symptomatic_and_nonsymptomatic(Simulation& sim) +{ + const auto& model = sim.get_model(); + constexpr size_t num_compartments_per_group = 6; + const auto num_groups = model.parameters.get_num_groups(); + std::vector indices_no_symptoms(num_compartments_per_group * size_t(num_groups)); + std::vector indices_symptoms(num_compartments_per_group * size_t(num_groups)); + + static const std::array no_symptoms_states = { + InfectionState::InfectedNoSymptomsNaive, + InfectionState::InfectedNoSymptomsNaiveConfirmed, + InfectionState::InfectedNoSymptomsPartialImmunity, + InfectionState::InfectedNoSymptomsPartialImmunityConfirmed, + InfectionState::InfectedNoSymptomsImprovedImmunity, + InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}; + + static const std::array symptoms_states = { + InfectionState::InfectedSymptomsNaive, + InfectionState::InfectedSymptomsNaiveConfirmed, + InfectionState::InfectedSymptomsPartialImmunity, + InfectionState::InfectedSymptomsPartialImmunityConfirmed, + InfectionState::InfectedSymptomsImprovedImmunity, + InfectionState::InfectedSymptomsImprovedImmunityConfirmed}; + + for (auto i = AgeGroup(0); i < num_groups; ++i) { + for (size_t j = 0; j < num_compartments_per_group; ++j) { + indices_no_symptoms[size_t(i) * num_compartments_per_group + j] = + model.populations.get_flat_index({i, no_symptoms_states[j]}); + indices_symptoms[size_t(i) * num_compartments_per_group + j] = + model.populations.get_flat_index({i, symptoms_states[j]}); + } + } + return std::make_tuple(std::move(indices_no_symptoms), std::move(indices_symptoms)); +} + } // namespace osecirvvs } // namespace mio From ea67eefdc89e388ba279b82824e643b014fa7654 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Tue, 19 Mar 2024 13:15:19 +0100 Subject: [PATCH 02/30] added tests + some fixes --- .../metapopulation_mobility_instant.h | 4 +- cpp/models/ode_secir/model.h | 7 +-- cpp/models/ode_secirvvs/model.h | 12 ++--- cpp/tests/test_mobility.cpp | 46 +++++++++++++++++++ 4 files changed, 58 insertions(+), 11 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index e49a258432..641c21bb84 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -252,13 +252,13 @@ using get_indices_of_symptomatic_and_nonsymptomatic_expr_t = * @brief Get the indices of symptomatic and non-symptomatic infection states. * * This function generates two vectors of indices, one for non-symptomatic infection states and one for symptomatic infection states. - * Each vector contains the flat indices of the corresponding infection states for each age group and compartment in the model. + * Each vector contains the flat indices of the corresponding infection states for each age group in the model. * * @tparam Base The base class for the simulation, defaults to mio::Simulation. * @param[in] sim The simulation object from which we obtain the model. * * @return A tuple containing two vectors of size_t. The first vector contains the indices of non-symptomatic infection states, - * and the second vector contains the indices of symptomatic infection states. The indices are ordered first by age group, then by compartment. + * and the second vector contains the indices of symptomatic infection states. The indices are ordered first by age group. */ template ::value, diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h index 0671c833d6..d78baf1c8f 100755 --- a/cpp/models/ode_secir/model.h +++ b/cpp/models/ode_secir/model.h @@ -29,7 +29,6 @@ #include "memilio/math/smoother.h" #include "memilio/math/eigen_util.h" #include "memilio/math/interpolation.h" -#include namespace mio { @@ -675,8 +674,10 @@ auto get_indices_of_symptomatic_and_nonsymptomatic(Simulation& sim) { const auto& model = sim.get_model(); const auto num_groups = model.parameters.get_num_groups(); - std::vector indices_no_symptoms(2 * size_t(num_groups)); - std::vector indices_symptoms(2 * size_t(num_groups)); + std::vector indices_no_symptoms; + std::vector indices_symptoms; //(2 * size_t(num_groups)); + indices_no_symptoms.reserve(2 * size_t(num_groups)); + indices_symptoms.reserve(2 * size_t(num_groups)); for (auto i = AgeGroup(0); i < num_groups; ++i) { indices_no_symptoms.emplace_back(model.populations.get_flat_index({i, InfectionState::InfectedNoSymptoms})); diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h index 5e15eb3561..0e0f5d9333 100644 --- a/cpp/models/ode_secirvvs/model.h +++ b/cpp/models/ode_secirvvs/model.h @@ -877,8 +877,10 @@ auto get_indices_of_symptomatic_and_nonsymptomatic(Simulation& sim) const auto& model = sim.get_model(); constexpr size_t num_compartments_per_group = 6; const auto num_groups = model.parameters.get_num_groups(); - std::vector indices_no_symptoms(num_compartments_per_group * size_t(num_groups)); - std::vector indices_symptoms(num_compartments_per_group * size_t(num_groups)); + std::vector indices_no_symptoms; + std::vector indices_symptoms; + indices_no_symptoms.reserve(num_compartments_per_group * size_t(num_groups)); + indices_symptoms.reserve(num_compartments_per_group * size_t(num_groups)); static const std::array no_symptoms_states = { InfectionState::InfectedNoSymptomsNaive, @@ -898,10 +900,8 @@ auto get_indices_of_symptomatic_and_nonsymptomatic(Simulation& sim) for (auto i = AgeGroup(0); i < num_groups; ++i) { for (size_t j = 0; j < num_compartments_per_group; ++j) { - indices_no_symptoms[size_t(i) * num_compartments_per_group + j] = - model.populations.get_flat_index({i, no_symptoms_states[j]}); - indices_symptoms[size_t(i) * num_compartments_per_group + j] = - model.populations.get_flat_index({i, symptoms_states[j]}); + indices_no_symptoms.emplace_back(model.populations.get_flat_index({i, no_symptoms_states[j]})); + indices_symptoms.emplace_back(model.populations.get_flat_index({i, symptoms_states[j]})); } } return std::make_tuple(std::move(indices_no_symptoms), std::move(indices_symptoms)); diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp index 9a5895606c..15d7be4aa8 100644 --- a/cpp/tests/test_mobility.cpp +++ b/cpp/tests/test_mobility.cpp @@ -169,3 +169,49 @@ TEST(TestMobility, edgeApplyMigration) EXPECT_DOUBLE_EQ(node1.get_result().get_last_value().sum(), 900); EXPECT_DOUBLE_EQ(node2.get_result().get_last_value().sum(), 1100); } + +TEST(TestMobility, condense_m_migrated_with_indices) +{ + using Model = mio::osecir::Model; + + //setup nodes + Model model(1); + 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(); + + //setup different edges + double t = 0.; + mio::SimulationNode> node1(model, t); + mio::SimulationNode> node2(model, t); + mio::MigrationEdge edge1(Eigen::VectorXd::Constant(10, 0.1)); + edge1.apply_migration(t, 0.0, node1, node2); + auto migrated = edge1.get_migrated().get_last_value(); + EXPECT_NEAR(migrated[0], 1.0, 1e-12); + EXPECT_NEAR(migrated[1], 2.0, 1e-12); + EXPECT_NEAR(migrated[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::MigrationEdge edge2(Eigen::VectorXd::Constant(10, 0.1)); + edge2.apply_migration(t, 0.5, node3, node4); + migrated = edge2.get_migrated().get_last_value(); + EXPECT_NEAR(migrated[0], 11.0, 1e-12); + EXPECT_NEAR(migrated[1], 5.0, 1e-12); + EXPECT_NEAR(migrated[2], 113.0, 1e-12); +} \ No newline at end of file From 46d11fd1ecbc2e1edcf9d4527003e80e775239c9 Mon Sep 17 00:00:00 2001 From: HenrZu Date: Wed, 20 Mar 2024 08:26:44 +0100 Subject: [PATCH 03/30] add example --- cpp/examples/CMakeLists.txt | 6 +-- cpp/examples/graph.cpp | 60 ----------------------- cpp/examples/ode_secir_graph.cpp | 84 ++++++++++++++++++++++++++++++++ cpp/tests/test_mobility.cpp | 2 +- 4 files changed, 88 insertions(+), 64 deletions(-) delete mode 100644 cpp/examples/graph.cpp create mode 100644 cpp/examples/ode_secir_graph.cpp diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index ec8849ce13..28d627c118 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -36,9 +36,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) diff --git a/cpp/examples/graph.cpp b/cpp/examples/graph.cpp deleted file mode 100644 index 49d0c48463..0000000000 --- a/cpp/examples/graph.cpp +++ /dev/null @@ -1,60 +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" - -#include - -int main() -{ - const auto t0 = 0.; - const auto tmax = 10.; - const auto dt = 0.5; //time step of migration, daily migration every second step - - mio::oseir::Model model; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = 10000; - model.parameters.set(1); - model.parameters.get().get_baseline()(0, 0) = 2.7; - model.parameters.set(1); - - //two mostly identical groups - auto model_group1 = model; - auto model_group2 = model; - //some contact restrictions in group 1 - model_group1.parameters.get().add_damping(0.5, mio::SimulationTime(5)); - //infection starts in group 1 - model_group1.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = 9990; - model_group1.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 10; - - mio::Graph>, mio::MigrationEdge> 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_migration_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..540fb7b491 --- /dev/null +++ b/cpp/examples/ode_secir_graph.cpp @@ -0,0 +1,84 @@ +/* +* 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 "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 migration, daily migration every second step + + mio::osecir::Model model(1); + 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; + + mio::Graph>, mio::MigrationEdge> 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)); + g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.1)); + + auto sim = mio::make_migration_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_migrated(); + results.print_table({"Commuter INS", "Commuter ISy", "Commuter Total"}); + + return 0; +} diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp index 15d7be4aa8..5563041c28 100644 --- a/cpp/tests/test_mobility.cpp +++ b/cpp/tests/test_mobility.cpp @@ -170,7 +170,7 @@ TEST(TestMobility, edgeApplyMigration) EXPECT_DOUBLE_EQ(node2.get_result().get_last_value().sum(), 1100); } -TEST(TestMobility, condense_m_migrated_with_indices) +TEST(TestMobility, condense_m_migrated) { using Model = mio::osecir::Model; From 52423ade87f256cfbe5e0eeaddf2e93f1651414c Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Wed, 20 Mar 2024 09:19:56 +0100 Subject: [PATCH 04/30] fix msvc --- cpp/memilio/mobility/metapopulation_mobility_instant.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp index 03abf61850..d43bd38dec 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp @@ -28,11 +28,11 @@ void MigrationEdge::condense_m_migrated(const double t, const std::vector Date: Wed, 20 Mar 2024 09:44:02 +0100 Subject: [PATCH 05/30] also results at end of commuting --- .../metapopulation_mobility_instant.cpp | 4 ++-- .../metapopulation_mobility_instant.h | 22 ++++++++++--------- cpp/models/ode_secir/model.h | 2 +- cpp/models/ode_secirvvs/model.h | 2 +- 4 files changed, 16 insertions(+), 14 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp index d43bd38dec..e01bd3441f 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp @@ -39,8 +39,8 @@ void MigrationEdge::condense_m_migrated(const double t, const std::vector::Vector(3) << num_INS, num_ISy, total_commuters).finished()); + m_mobility_results.add_time_point( + t, (mio::TimeSeries::Vector(3) << num_INS, num_ISy, total_commuters).finished()); } } // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 641c21bb84..6995643511 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -291,7 +291,7 @@ class MigrationEdge , m_migrated(params.get_coefficients().get_shape().rows()) , m_return_times(0) , m_return_migrated(false) - , m_num_migrated(3) + , m_mobility_results(3) { } @@ -304,7 +304,7 @@ class MigrationEdge , m_migrated(coeffs.rows()) , m_return_times(0) , m_return_migrated(false) - , m_num_migrated(3) + , m_mobility_results(3) { } @@ -322,11 +322,11 @@ class MigrationEdge */ TimeSeries& get_migrated() { - return m_num_migrated; + return m_mobility_results; } const TimeSeries& get_migrated() const { - return m_num_migrated; + return m_mobility_results; } /** @@ -349,12 +349,12 @@ class MigrationEdge 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)}; - TimeSeries m_num_migrated; + TimeSeries m_mobility_results; /** - * Computes a condensed version of m_migrated and puts it in m_num_migrated. - * m_num_migrated then only contains commuters with infection states InfectedNoSymptoms and InfectedSymptoms. - * Additionally, the total number of commuters is stored in the last entry of m_num_migrated. + * Computes a condensed version of m_migrated and puts it in m_mobility_results. + * 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 current time */ void condense_m_migrated(const double t, const std::vector& indices_non_symptomatic, @@ -498,6 +498,9 @@ void MigrationEdge::apply_migration(double t, double dt, SimulationNode& no m_t_last_dynamic_npi_check = t; } + static auto indices_tuple = get_indices_of_symptomatic_and_nonsymptomatic(node_from); + auto& [indices_no_symptoms, indices_symptoms] = indices_tuple; + //returns for (Eigen::Index i = m_return_times.get_num_time_points() - 1; i >= 0; --i) { if (m_return_times.get_time(i) <= t) { @@ -527,6 +530,7 @@ void MigrationEdge::apply_migration(double t, double dt, SimulationNode& no } node_from.get_result().get_last_value() += m_migrated[i]; node_to.get_result().get_last_value() -= m_migrated[i]; + condense_m_migrated(t, indices_no_symptoms, indices_symptoms); m_migrated.remove_time_point(i); m_return_times.remove_time_point(i); } @@ -545,8 +549,6 @@ void MigrationEdge::apply_migration(double t, double dt, SimulationNode& no node_to.get_result().get_last_value() += m_migrated.get_last_value(); node_from.get_result().get_last_value() -= m_migrated.get_last_value(); - static auto indices_tuple = get_indices_of_symptomatic_and_nonsymptomatic(node_from); - auto& [indices_no_symptoms, indices_symptoms] = indices_tuple; condense_m_migrated(t, indices_no_symptoms, indices_symptoms); } m_return_migrated = !m_return_migrated; diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h index d78baf1c8f..3cfd48125e 100755 --- a/cpp/models/ode_secir/model.h +++ b/cpp/models/ode_secir/model.h @@ -670,7 +670,7 @@ auto test_commuters(Simulation& sim, Eigen::Ref migrated, } template > -auto get_indices_of_symptomatic_and_nonsymptomatic(Simulation& sim) +auto get_indices_of_symptomatic_and_nonsymptomatic(const Simulation& sim) { const auto& model = sim.get_model(); const auto num_groups = model.parameters.get_num_groups(); diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h index 0e0f5d9333..ceda8dd0fd 100644 --- a/cpp/models/ode_secirvvs/model.h +++ b/cpp/models/ode_secirvvs/model.h @@ -872,7 +872,7 @@ auto test_commuters(Simulation& sim, Eigen::Ref migrated, } template > -auto get_indices_of_symptomatic_and_nonsymptomatic(Simulation& sim) +auto get_indices_of_symptomatic_and_nonsymptomatic(const Simulation& sim) { const auto& model = sim.get_model(); constexpr size_t num_compartments_per_group = 6; From 3c07c5ac7b5c075eac7c06fbbd4fcc279b688ddb Mon Sep 17 00:00:00 2001 From: HenrZu Date: Thu, 21 Mar 2024 12:14:12 +0100 Subject: [PATCH 06/30] save_edges function + example graph study --- cpp/examples/CMakeLists.txt | 4 + .../ode_secir_parameter_study_graph.cpp | 297 ++++++++++++++++++ cpp/memilio/io/result_io.cpp | 114 +++++++ cpp/memilio/io/result_io.h | 25 ++ 4 files changed, 440 insertions(+) create mode 100644 cpp/examples/ode_secir_parameter_study_graph.cpp diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 28d627c118..131ecc87bf 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -83,6 +83,10 @@ 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_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/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp new file mode 100644 index 0000000000..099c1f80ab --- /dev/null +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -0,0 +1,297 @@ +/* +* 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 + +/** + * 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 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 + // TimeExposed and TimeInfectedNoSymptoms are calculated as described in + // Khailaie et al. (https://doi.org/10.1186/s12916-020-01884-4) + // given SI_min = 3.935, SI_max = 4.6, INC = 5.2 + 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); + } +} + +/** + * 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 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); + + 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)); + params_graph.add_edge(1, 0, + Eigen::VectorXd::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.1)); + params_graph.add_edge(1, 2, + Eigen::VectorXd::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.15)); + params_graph.add_edge(2, 1, + Eigen::VectorXd::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.2)); + + //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_migrated(); + }); + + 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))); + } + auto county_ids = std::vector{1001, 1002, 1003}; + auto save_results_status = + save_results(ensemble_results, ensemble_params, county_ids, "test_results", 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/result_io.cpp b/cpp/memilio/io/result_io.cpp index 39675db0db..feeedc8c23 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -89,6 +89,120 @@ IOResult save_result(const std::vector>& results, const } return success(); } +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 int num_edges = static_cast(results.size()); + mio::unused(num_edges); + int edge_indx = 0; + H5File file{H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, 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 + ")"); + const int num_timepoints = static_cast(result.get_num_time_points()); + const int num_infectionstates = 3 ; // (int)result.get_num_elements() / num_groups; + + 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, "Time DataSpace could not be created."); + 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, "Time DataSet could not be created (Time)."); + 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, "Time data could not be written."); + + int start_id = ids[edge_indx].first; + while (ids[edge_indx].first == start_id) { + const auto& result_edge = results[edge_indx]; + auto total = Eigen::Matrix::Zero( + num_timepoints, num_infectionstates) + .eval(); + + for (Eigen::Index t_idx = 0; t_idx < result.get_num_time_points(); ++t_idx) { + auto v = result_edge.get_value(t_idx).transpose().eval(); + // std::cout << "result_edge: " << v << std::endl; + // mio::slice(total, {t_idx, 1}, {0, num_infectionstates}) += + // mio::slice(v, {num_infectionstates, num_infectionstates}); + total.row(t_idx) = v; + } + + auto total_to_vector = std::vector(total.data(), total.data() + total.size()); + + auto results_edge_to_vector = std::vector(result_edge.get_times().begin(), result_edge.get_times().end()); + + + std::cout << "result_edge: " << total_to_vector.size() << std::endl; + + std::cout << "result_edge: " << results_edge_to_vector.size() << std::endl; + + hsize_t dims_values[] = {static_cast(num_timepoints), static_cast(num_infectionstates)}; + H5DataSpace dspace_values{H5Screate_simple(2, dims_values, NULL)}; + MEMILIO_H5_CHECK(dspace_values.id, StatusCode::UnknownError, "Values DataSpace could not be created."); + auto dset_name = 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, "Values DataSet could not be created."); + + MEMILIO_H5_CHECK(H5Dwrite(dset_values.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, total.data()), + StatusCode::UnknownError, "Values data could not be written."); + edge_indx++; + } + } + return success(); +} herr_t store_group_name(hid_t /*id*/, const char* name, const H5L_info_t* /*linfo*/, void* opdata) { diff --git a/cpp/memilio/io/result_io.h b/cpp/memilio/io/result_io.h index 8de5ec6461..2387665c30 100644 --- a/cpp/memilio/io/result_io.h +++ b/cpp/memilio/io/result_io.h @@ -121,6 +121,31 @@ IOResult save_result_with_params(const std::vector>& re return success(); } +/** + * @brief Save the results of the edges for a single graph simulation run. + * @param result Simulation results per edge of the graph. + * @param ids Identifiers for the start and end node of the edges. + * @param num_groups Number of groups in the results. + * @param filename Name of file + * @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); + +/** + * Save the results for the Edges obtained from the function condense_m_migrated. + * @param ensemble_edges Simulation results for each run for each edge. + * @param num_groups Number of age groups used simulation. + * @param ids Identifiers for the start and end node of the edges. + * @param result_dir Top level directory for all results of the parameter study. + * @param save_single_runs [Default: true] Defines if single run results are written. + * @param save_single_runs [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); + /** * @brief Save the results of a parameter study. * Stores different percentiles (p5, p25, p50, p75, p90) and sums of the results and parameters. From be9df03faebdb29b477ded151f241baf19d8790c Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 8 Apr 2024 15:23:26 +0200 Subject: [PATCH 07/30] rm print statements --- cpp/memilio/io/result_io.cpp | 38 ++++++++++++------------------------ 1 file changed, 13 insertions(+), 25 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index feeedc8c23..509e3bc354 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -89,7 +89,7 @@ IOResult save_result(const std::vector>& results, const } return success(); } -IOResult save_edges(const std::vector>>& ensemble_edges, +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) { @@ -123,22 +123,18 @@ IOResult save_edges(const std::vector>>& en 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())); + 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) +IOResult save_edges(const std::vector>& results, const std::vector>& ids, + const std::string& filename) { const int num_edges = static_cast(results.size()); mio::unused(num_edges); @@ -151,8 +147,8 @@ IOResult save_edges(const std::vector>& results, const 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 + ")"); - const int num_timepoints = static_cast(result.get_num_time_points()); - const int num_infectionstates = 3 ; // (int)result.get_num_elements() / num_groups; + const int num_timepoints = static_cast(result.get_num_time_points()); + constexpr int num_infectionstates = 3; // (int)result.get_num_elements() / num_groups; hsize_t dims_t[] = {static_cast(num_timepoints)}; H5DataSpace dspace_t{H5Screate_simple(1, dims_t, NULL)}; @@ -172,22 +168,14 @@ IOResult save_edges(const std::vector>& results, const .eval(); for (Eigen::Index t_idx = 0; t_idx < result.get_num_time_points(); ++t_idx) { - auto v = result_edge.get_value(t_idx).transpose().eval(); - // std::cout << "result_edge: " << v << std::endl; - // mio::slice(total, {t_idx, 1}, {0, num_infectionstates}) += - // mio::slice(v, {num_infectionstates, num_infectionstates}); + auto v = result_edge.get_value(t_idx).transpose().eval(); total.row(t_idx) = v; } auto total_to_vector = std::vector(total.data(), total.data() + total.size()); - auto results_edge_to_vector = std::vector(result_edge.get_times().begin(), result_edge.get_times().end()); - - - std::cout << "result_edge: " << total_to_vector.size() << std::endl; - - std::cout << "result_edge: " << results_edge_to_vector.size() << std::endl; - + auto results_edge_to_vector = + std::vector(result_edge.get_times().begin(), result_edge.get_times().end()); hsize_t dims_values[] = {static_cast(num_timepoints), static_cast(num_infectionstates)}; H5DataSpace dspace_values{H5Screate_simple(2, dims_values, NULL)}; MEMILIO_H5_CHECK(dspace_values.id, StatusCode::UnknownError, "Values DataSpace could not be created."); From b1243c2fe3bb6c100ac3d8212958a93ff3e9e2b6 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Tue, 9 Apr 2024 11:35:01 +0200 Subject: [PATCH 08/30] [ci skip] draft for tests --- cpp/memilio/io/result_io.cpp | 52 ++++++++++++------ cpp/tests/test_save_results.cpp | 94 +++++++++++++++++++++++++++++++++ 2 files changed, 130 insertions(+), 16 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 509e3bc354..141765d68f 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -161,31 +161,43 @@ IOResult save_edges(const std::vector>& results, const StatusCode::UnknownError, "Time data could not be written."); int start_id = ids[edge_indx].first; + auto total = Eigen::Matrix::Zero(num_timepoints, + num_infectionstates) + .eval(); while (ids[edge_indx].first == start_id) { - const auto& result_edge = results[edge_indx]; - auto total = Eigen::Matrix::Zero( - num_timepoints, num_infectionstates) - .eval(); + const auto& result_edge_indx = results[edge_indx]; + auto edge_result = Eigen::Matrix::Zero( + num_timepoints, num_infectionstates) + .eval(); for (Eigen::Index t_idx = 0; t_idx < result.get_num_time_points(); ++t_idx) { - auto v = result_edge.get_value(t_idx).transpose().eval(); - total.row(t_idx) = v; + auto v = result_edge_indx.get_value(t_idx).transpose().eval(); + edge_result.row(t_idx) = v; + total.row(t_idx) += v; } - auto total_to_vector = std::vector(total.data(), total.data() + total.size()); - - auto results_edge_to_vector = - std::vector(result_edge.get_times().begin(), result_edge.get_times().end()); hsize_t dims_values[] = {static_cast(num_timepoints), static_cast(num_infectionstates)}; H5DataSpace dspace_values{H5Screate_simple(2, dims_values, NULL)}; MEMILIO_H5_CHECK(dspace_values.id, StatusCode::UnknownError, "Values DataSpace could not be created."); - auto dset_name = std::to_string(ids[edge_indx].second); + 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, "Values DataSet could not be created."); - MEMILIO_H5_CHECK(H5Dwrite(dset_values.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, total.data()), - StatusCode::UnknownError, "Values data could not be written."); + MEMILIO_H5_CHECK( + H5Dwrite(dset_values.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, edge_result.data()), + StatusCode::UnknownError, "Values data could not be written."); + + // in the final iteration, we also save the total values + if (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, "Values DataSet could not be created."); + MEMILIO_H5_CHECK( + H5Dwrite(dset_total.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, total.data()), + StatusCode::UnknownError, "Values data could not be written."); + } edge_indx++; } } @@ -219,8 +231,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)}; @@ -271,7 +291,7 @@ IOResult> read_result(const std::string& filename) 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); + auto group_name = "/" + h5group_name + "/" + h5_key + std::to_string(group_idx + 1); 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."); diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index a98c12b0df..e81f7a3bde 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -343,3 +343,97 @@ TEST(TestSaveResult, save_percentiles_and_sums) auto results_run2_sum = mio::read_result(tmp_results_dir + "/results_run2_sum.h5"); ASSERT_TRUE(results_run2_sum); } + +TEST(TestSaveResult, save_edges) +{ + double t0 = 0; + double tmax = 5; + double dt = 0.1; + double cont_freq = 10.; + + double pop_total = 10000; + + mio::osecir::Model model(1); + auto& params = model.parameters; + mio::AgeGroup nb_groups = params.get_num_groups(); + ; + + for (auto i = mio::AgeGroup(0); i < nb_groups; i++) { + params.get()[i] = 3.2; + params.get()[i] = 2; + params.get()[i] = 5.; + params.get()[i] = 10.; + params.get()[i] = 8.; + + model.populations[{i, mio::osecir::InfectionState::Exposed}] = 100; + model.populations[{i, mio::osecir::InfectionState::InfectedNoSymptoms}] = 50; + model.populations[{i, mio::osecir::InfectionState::InfectedSymptoms}] = 50; + model.populations[{i, mio::osecir::InfectionState::InfectedSevere}] = 20; + model.populations[{i, mio::osecir::InfectionState::InfectedCritical}] = 10; + model.populations[{i, mio::osecir::InfectionState::Recovered}] = 10; + model.populations[{i, mio::osecir::InfectionState::Dead}] = 0; + model.populations.set_difference_from_total({i, mio::osecir::InfectionState::Susceptible}, pop_total); + + params.get()[i] = 0.06; + params.get()[i] = 0.67; + params.get()[i] = 0.09; + params.get()[i] = 0.25; + params.get()[i] = 0.2; + params.get()[i] = 0.25; + params.get()[i] = 0.3; + } + + mio::ContactMatrixGroup& contact_matrix = params.get(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant((size_t)nb_groups, (size_t)nb_groups, cont_freq)); + contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); + + // setup graph + mio::Graph>, mio::MigrationEdge> g; + g.add_node(0, model, t0); + g.add_node(1, model, t0); + g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.1)); + g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.1)); + + // simulation + auto sim = mio::make_migration_sim(t0, dt, std::move(g)); + sim.advance(tmax); + + // get_results and pair ids + auto& results_edge_0 = sim.get_graph().edges()[0].property.get_migrated(); + std::vector> results_edges = {results_edge_0, results_edge_0}; + + auto pairs_edges = std::vector>{}; + pairs_edges.push_back({0, 1}); + pairs_edges.push_back({1, 0}); + + 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); + + auto results_from_file = mio::read_result(results_file_path); + ASSERT_TRUE(results_from_file); + auto result_from_file = results_from_file.value()[0]; + + ASSERT_EQ(result_from_file.get_groups().get_num_time_points(), results_edge_0.get_num_time_points()); + ASSERT_EQ(result_from_file.get_totals().get_num_time_points(), results_edge_0.get_num_time_points()); + for (Eigen::Index i = 0; i < results_edge_0.get_num_time_points(); i++) { + ASSERT_EQ(result_from_file.get_groups().get_num_elements(), results_edge_0.get_num_elements()) + << "at row " << i; + ASSERT_EQ(result_from_file.get_totals().get_num_elements(), + results_edge_0.get_num_elements() / static_cast((size_t)nb_groups)) + << "at row " << i; + ASSERT_NEAR(results_edge_0.get_time(i), result_from_file.get_groups().get_time(i), 1e-10) << "at row " << i; + ASSERT_NEAR(results_edge_0.get_time(i), result_from_file.get_totals().get_time(i), 1e-10) << "at row " << i; + for (Eigen::Index l = 0; l < result_from_file.get_totals().get_num_elements(); l++) { + double total = 0.0; + for (Eigen::Index j = 0; j < Eigen::Index((size_t)nb_groups); j++) { + total += results_edge_0[i][j * (size_t)mio::osecir::InfectionState::Count + l]; + EXPECT_NEAR(result_from_file.get_groups()[i][j * (size_t)mio::osecir::InfectionState::Count + l], + results_edge_0[i][j * (size_t)mio::osecir::InfectionState::Count + l], 1e-10) + << " at row " << i << " at row " << l << " at Group " << j; + } + EXPECT_NEAR(result_from_file.get_totals()[i][l], total, 1e-10) << " at row " << i << " at row " << l; + } + } +} \ No newline at end of file From 68e0cd55b81011c3875bba2d2ea67df463644dff Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 10 Apr 2024 12:47:26 +0200 Subject: [PATCH 09/30] fix some bugs + tests --- cpp/memilio/io/result_io.cpp | 32 +++++++++++++++---------- cpp/tests/test_save_results.cpp | 41 +++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 12 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 141765d68f..191556e77c 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -18,6 +18,7 @@ * limitations under the License. */ #include "memilio/io/result_io.h" +#include #ifdef MEMILIO_HAS_HDF5 @@ -137,8 +138,7 @@ IOResult save_edges(const std::vector>& results, const const std::string& filename) { const int num_edges = static_cast(results.size()); - mio::unused(num_edges); - int edge_indx = 0; + int edge_indx = 0; H5File file{H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)}; MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, filename); while (edge_indx < num_edges) { @@ -164,14 +164,13 @@ IOResult save_edges(const std::vector>& results, const auto total = Eigen::Matrix::Zero(num_timepoints, num_infectionstates) .eval(); - while (ids[edge_indx].first == start_id) { - const auto& result_edge_indx = results[edge_indx]; - auto edge_result = Eigen::Matrix::Zero( + while (ids[edge_indx].first == start_id && edge_indx < num_edges) { + const auto& result_edge = results[edge_indx]; + auto edge_result = Eigen::Matrix::Zero( num_timepoints, num_infectionstates) .eval(); - - for (Eigen::Index t_idx = 0; t_idx < result.get_num_time_points(); ++t_idx) { - auto v = result_edge_indx.get_value(t_idx).transpose().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; } @@ -189,7 +188,7 @@ IOResult save_edges(const std::vector>& results, const StatusCode::UnknownError, "Values data could not be written."); // in the final iteration, we also save the total values - if (ids[edge_indx + 1].first != start_id) { + if (ids[edge_indx + 1].first != start_id || edge_indx == num_edges - 1) { 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)}; @@ -290,8 +289,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 + "/" + h5_key + 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."); @@ -316,7 +324,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/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index e81f7a3bde..c780c081f7 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -293,6 +293,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); @@ -306,6 +308,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_migrated(); + }); + return 0; //function needs to return something }); @@ -342,6 +352,37 @@ 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) From 1677a51b960a00fd020d70491aab667f17aa1823 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 10 Apr 2024 14:40:46 +0200 Subject: [PATCH 10/30] change order if condition --- cpp/memilio/io/result_io.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 191556e77c..5af8f51fda 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -188,7 +188,7 @@ IOResult save_edges(const std::vector>& results, const StatusCode::UnknownError, "Values data could not be written."); // in the final iteration, we also save the total values - if (ids[edge_indx + 1].first != start_id || edge_indx == num_edges - 1) { + 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)}; From 20803e0f673df21c7fd4f5975d27342972b63e03 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 10 Apr 2024 15:08:30 +0200 Subject: [PATCH 11/30] naming --- cpp/memilio/io/result_io.cpp | 2 +- cpp/memilio/io/result_io.h | 2 +- cpp/memilio/mobility/metapopulation_mobility_instant.cpp | 2 +- cpp/memilio/mobility/metapopulation_mobility_instant.h | 6 +++--- cpp/tests/test_mobility.cpp | 2 +- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 5af8f51fda..a5109856e6 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -164,7 +164,7 @@ IOResult save_edges(const std::vector>& results, const auto total = Eigen::Matrix::Zero(num_timepoints, num_infectionstates) .eval(); - while (ids[edge_indx].first == start_id && edge_indx < num_edges) { + 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_infectionstates) diff --git a/cpp/memilio/io/result_io.h b/cpp/memilio/io/result_io.h index 2387665c30..db42ad4c5a 100644 --- a/cpp/memilio/io/result_io.h +++ b/cpp/memilio/io/result_io.h @@ -133,7 +133,7 @@ IOResult save_edges(const std::vector>& results, const const std::string& filename); /** - * Save the results for the Edges obtained from the function condense_m_migrated. + * Save the results for the Edges obtained from the function condense_m_mobility. * @param ensemble_edges Simulation results for each run for each edge. * @param num_groups Number of age groups used simulation. * @param ids Identifiers for the start and end node of the edges. diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp index e01bd3441f..04218fde92 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp @@ -21,7 +21,7 @@ namespace mio { -void MigrationEdge::condense_m_migrated(const double t, const std::vector& indices_non_symptomatic, +void MigrationEdge::condense_m_mobility(const double t, const std::vector& indices_non_symptomatic, const std::vector& indices_symptomatic) { diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 6995643511..70c38ba854 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -357,7 +357,7 @@ class MigrationEdge * Additionally, the total number of commuters is stored in the last entry of m_mobility_results. * @param[in] t current time */ - void condense_m_migrated(const double t, const std::vector& indices_non_symptomatic, + void condense_m_mobility(const double t, const std::vector& indices_non_symptomatic, const std::vector& indices_symptomatic); }; @@ -530,7 +530,7 @@ void MigrationEdge::apply_migration(double t, double dt, SimulationNode& no } node_from.get_result().get_last_value() += m_migrated[i]; node_to.get_result().get_last_value() -= m_migrated[i]; - condense_m_migrated(t, indices_no_symptoms, indices_symptoms); + condense_m_mobility(t, indices_no_symptoms, indices_symptoms); m_migrated.remove_time_point(i); m_return_times.remove_time_point(i); } @@ -549,7 +549,7 @@ void MigrationEdge::apply_migration(double t, double dt, SimulationNode& no node_to.get_result().get_last_value() += m_migrated.get_last_value(); node_from.get_result().get_last_value() -= m_migrated.get_last_value(); - condense_m_migrated(t, indices_no_symptoms, indices_symptoms); + condense_m_mobility(t, indices_no_symptoms, indices_symptoms); } m_return_migrated = !m_return_migrated; } diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp index 5563041c28..0d0991d3c7 100644 --- a/cpp/tests/test_mobility.cpp +++ b/cpp/tests/test_mobility.cpp @@ -170,7 +170,7 @@ TEST(TestMobility, edgeApplyMigration) EXPECT_DOUBLE_EQ(node2.get_result().get_last_value().sum(), 1100); } -TEST(TestMobility, condense_m_migrated) +TEST(TestMobility, condense_m_mobility) { using Model = mio::osecir::Model; From e3e7c38355a797067fdd93d09ac5337332a6b6b9 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Thu, 25 Apr 2024 13:34:41 +0200 Subject: [PATCH 12/30] change structure --- cpp/examples/ode_secir_graph.cpp | 28 +++++- .../ode_secir_parameter_study_graph.cpp | 46 ++++++++-- .../metapopulation_mobility_instant.cpp | 35 ++++---- .../metapopulation_mobility_instant.h | 88 +++++++++++++++++-- cpp/models/ode_secir/model.h | 22 ----- cpp/models/ode_secirvvs/model.h | 36 -------- 6 files changed, 164 insertions(+), 91 deletions(-) diff --git a/cpp/examples/ode_secir_graph.cpp b/cpp/examples/ode_secir_graph.cpp index 540fb7b491..48ed48b0b8 100644 --- a/cpp/examples/ode_secir_graph.cpp +++ b/cpp/examples/ode_secir_graph.cpp @@ -31,7 +31,8 @@ int main() const auto tmax = 30.; const auto dt = 0.5; //time step of migration, daily migration every second step - mio::osecir::Model model(1); + 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); @@ -66,11 +67,32 @@ int main() 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::MigrationEdge> 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)); - g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.1)); + 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_migration_sim(t0, dt, std::move(g)); diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp index 099c1f80ab..895966f45c 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -33,6 +33,7 @@ #include "ode_secir/parameter_space.h" #include "boost/filesystem.hpp" #include "memilio/utils/stl_util.h" +#include #include #include @@ -187,6 +188,31 @@ void set_synthetic_population_data(mio::osecir::Model& model) } } +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. @@ -221,19 +247,24 @@ int main() // 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)); + 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)); + 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)); + 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)); + Eigen::VectorXd::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.2), + indices_save_edges); //run parameter study auto parameter_study = @@ -285,10 +316,9 @@ int main() ensemble_params.emplace_back(std::move(std::get<1>(run))); ensemble_edges.emplace_back(std::move(std::get<2>(run))); } - auto county_ids = std::vector{1001, 1002, 1003}; - auto save_results_status = - save_results(ensemble_results, ensemble_params, county_ids, "test_results", false); - auto pairs_edges = std::vector>{}; + auto county_ids = std::vector{1001, 1002, 1003}; + auto save_results_status = save_results(ensemble_results, ensemble_params, county_ids, "test_results", 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]}); } diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp index 04218fde92..f47eeccabd 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp @@ -18,29 +18,34 @@ * limitations under the License. */ #include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/utils/compiler_diagnostics.h" namespace mio { -void MigrationEdge::condense_m_mobility(const double t, const std::vector& indices_non_symptomatic, - const std::vector& indices_symptomatic) +void MigrationEdge::condense_m_mobility(const double t) { + // only call this if his->m_save_indices.size() is greater than 0. Perfect would be to define this in compile time + const size_t save_indices_size = this->m_save_indices.size(); + if (save_indices_size > 0) { - const auto& last_value = m_migrated.get_last_value(); + const auto& last_value = m_migrated.get_last_value(); + Eigen::VectorXd condensed_values = Eigen::VectorXd::Zero(save_indices_size + 1); - auto num_INS = - std::accumulate(indices_non_symptomatic.begin(), indices_non_symptomatic.end(), 0., [&](auto sum, auto i) { - return sum + last_value[i]; - }); + // sum up the values of m_save_indices for each group (e.g. Age groups) + std::transform(this->m_save_indices.begin(), this->m_save_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]; + }); + }); - auto num_ISy = std::accumulate(indices_symptomatic.begin(), indices_symptomatic.end(), 0., [&](auto sum, auto i) { - return sum + last_value[i]; - }); + // the last value is the sum of commuters + condensed_values[save_indices_size] = m_migrated.get_last_value().sum(); - double total_commuters = m_migrated.get_last_value().sum(); - - // as time point t which contains now the carriers, infectious and total over age groups - m_mobility_results.add_time_point( - t, (mio::TimeSeries::Vector(3) << num_INS, num_ISy, total_commuters).finished()); + // Move the condensed values to the m_mobility_results time series + m_mobility_results.add_time_point(t, std::move(condensed_values)); + } } } // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 70c38ba854..0723af59f1 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -127,6 +127,8 @@ class MigrationParameters */ MigrationParameters(const MigrationCoefficientGroup& coeffs) : m_coefficients(coeffs) + , m_save_indices(0) + { } @@ -136,6 +138,31 @@ class MigrationParameters */ MigrationParameters(const Eigen::VectorXd& coeffs) : m_coefficients({MigrationCoefficients(coeffs)}) + , m_save_indices(0) + { + } + + /** + * @brief Constructor for a new MigrationParameters object. + * + * @param coeffs migration coefficients + * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. + */ + MigrationParameters(const MigrationCoefficientGroup& coeffs, const std::vector>& save_indices) + : m_coefficients(coeffs) + , m_save_indices(save_indices) + { + } + + /** + * @brief Constructor for a new MigrationParameters object. + * + * @param coeffs migration coefficients + * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. + */ + MigrationParameters(const Eigen::VectorXd& coeffs, const std::vector>& save_indices) + : m_coefficients({MigrationCoefficients(coeffs)}) + , m_save_indices(save_indices) { } @@ -177,6 +204,20 @@ class MigrationParameters { m_coefficients = coeffs; } + + /** + * @return the indices we want to save for the edge + */ + const auto& get_save_indices() const + { + return m_save_indices; + } + + const auto& get_save_indices() + { + return m_save_indices; + } + /** @} */ /** @@ -239,6 +280,7 @@ class MigrationParameters private: MigrationCoefficientGroup m_coefficients; //one per group and compartment DynamicNPIs m_dynamic_npis; + std::vector> m_save_indices; // groups of indices from compartments to save }; /** @@ -291,7 +333,8 @@ class MigrationEdge , m_migrated(params.get_coefficients().get_shape().rows()) , m_return_times(0) , m_return_migrated(false) - , m_mobility_results(3) + , m_save_indices(params.get_save_indices()) + , m_mobility_results(m_save_indices.size()) { } @@ -304,7 +347,38 @@ class MigrationEdge , m_migrated(coeffs.rows()) , m_return_times(0) , m_return_migrated(false) - , m_mobility_results(3) + , m_save_indices(0) + , m_mobility_results(m_save_indices.size()) + { + } + + /** + * create edge with coefficients as MigrationParameters object and a 2d vector of indices which determine which compartments we save. + * @param coeffs % of people in each group and compartment that migrate in each time step. + * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. + */ + MigrationEdge(const MigrationParameters& params, const std::vector>& save_indices) + : m_parameters(params) + , m_migrated(params.get_coefficients().get_shape().rows()) + , m_return_times(0) + , m_return_migrated(false) + , m_save_indices(save_indices) + , m_mobility_results(m_save_indices.size() + 1) + { + } + + /** + * create edge with coefficients and a 2d vector of indices which determine which compartments we save. + * @param coeffs % of people in each group and compartment that migrate in each time step. + * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. + */ + MigrationEdge(const Eigen::VectorXd& coeffs, const std::vector>& save_indices) + : m_parameters(coeffs) + , m_migrated(coeffs.rows()) + , m_return_times(0) + , m_return_migrated(false) + , m_save_indices(save_indices) + , m_mobility_results(m_save_indices.size() + 1) { } @@ -349,7 +423,8 @@ class MigrationEdge 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)}; - TimeSeries m_mobility_results; + std::vector> m_save_indices; // groups of indices from compartments to save + TimeSeries m_mobility_results; // save results from edges + total number of commuters /** * Computes a condensed version of m_migrated and puts it in m_mobility_results. @@ -357,8 +432,7 @@ class MigrationEdge * Additionally, the total number of commuters is stored in the last entry of m_mobility_results. * @param[in] t current time */ - void condense_m_mobility(const double t, const std::vector& indices_non_symptomatic, - const std::vector& indices_symptomatic); + void condense_m_mobility(const double t); }; /** @@ -530,7 +604,7 @@ void MigrationEdge::apply_migration(double t, double dt, SimulationNode& no } node_from.get_result().get_last_value() += m_migrated[i]; node_to.get_result().get_last_value() -= m_migrated[i]; - condense_m_mobility(t, indices_no_symptoms, indices_symptoms); + condense_m_mobility(t); m_migrated.remove_time_point(i); m_return_times.remove_time_point(i); } @@ -549,7 +623,7 @@ void MigrationEdge::apply_migration(double t, double dt, SimulationNode& no node_to.get_result().get_last_value() += m_migrated.get_last_value(); node_from.get_result().get_last_value() -= m_migrated.get_last_value(); - condense_m_mobility(t, indices_no_symptoms, indices_symptoms); + condense_m_mobility(t); } m_return_migrated = !m_return_migrated; } diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h index 3cfd48125e..62c8325785 100755 --- a/cpp/models/ode_secir/model.h +++ b/cpp/models/ode_secir/model.h @@ -669,28 +669,6 @@ auto test_commuters(Simulation& sim, Eigen::Ref migrated, } } -template > -auto get_indices_of_symptomatic_and_nonsymptomatic(const Simulation& sim) -{ - const auto& model = sim.get_model(); - const auto num_groups = model.parameters.get_num_groups(); - std::vector indices_no_symptoms; - std::vector indices_symptoms; //(2 * size_t(num_groups)); - indices_no_symptoms.reserve(2 * size_t(num_groups)); - indices_symptoms.reserve(2 * size_t(num_groups)); - - for (auto i = AgeGroup(0); i < num_groups; ++i) { - indices_no_symptoms.emplace_back(model.populations.get_flat_index({i, InfectionState::InfectedNoSymptoms})); - indices_no_symptoms.emplace_back( - model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsConfirmed})); - - indices_symptoms.emplace_back(model.populations.get_flat_index({i, InfectionState::InfectedSymptoms})); - indices_symptoms.emplace_back(model.populations.get_flat_index({i, InfectionState::InfectedSymptomsConfirmed})); - } - - return std::make_tuple(std::move(indices_no_symptoms), std::move(indices_symptoms)); -} - } // namespace osecir } // namespace mio diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h index ceda8dd0fd..1c9c9ae2e6 100644 --- a/cpp/models/ode_secirvvs/model.h +++ b/cpp/models/ode_secirvvs/model.h @@ -871,42 +871,6 @@ auto test_commuters(Simulation& sim, Eigen::Ref migrated, } } -template > -auto get_indices_of_symptomatic_and_nonsymptomatic(const Simulation& sim) -{ - const auto& model = sim.get_model(); - constexpr size_t num_compartments_per_group = 6; - const auto num_groups = model.parameters.get_num_groups(); - std::vector indices_no_symptoms; - std::vector indices_symptoms; - indices_no_symptoms.reserve(num_compartments_per_group * size_t(num_groups)); - indices_symptoms.reserve(num_compartments_per_group * size_t(num_groups)); - - static const std::array no_symptoms_states = { - InfectionState::InfectedNoSymptomsNaive, - InfectionState::InfectedNoSymptomsNaiveConfirmed, - InfectionState::InfectedNoSymptomsPartialImmunity, - InfectionState::InfectedNoSymptomsPartialImmunityConfirmed, - InfectionState::InfectedNoSymptomsImprovedImmunity, - InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}; - - static const std::array symptoms_states = { - InfectionState::InfectedSymptomsNaive, - InfectionState::InfectedSymptomsNaiveConfirmed, - InfectionState::InfectedSymptomsPartialImmunity, - InfectionState::InfectedSymptomsPartialImmunityConfirmed, - InfectionState::InfectedSymptomsImprovedImmunity, - InfectionState::InfectedSymptomsImprovedImmunityConfirmed}; - - for (auto i = AgeGroup(0); i < num_groups; ++i) { - for (size_t j = 0; j < num_compartments_per_group; ++j) { - indices_no_symptoms.emplace_back(model.populations.get_flat_index({i, no_symptoms_states[j]})); - indices_symptoms.emplace_back(model.populations.get_flat_index({i, symptoms_states[j]})); - } - } - return std::make_tuple(std::move(indices_no_symptoms), std::move(indices_symptoms)); -} - } // namespace osecirvvs } // namespace mio From bf69ac0b9514e910b8ba5faab31021b065c908e3 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Thu, 25 Apr 2024 15:42:29 +0200 Subject: [PATCH 13/30] adjust test + rm unused var --- .../metapopulation_mobility_instant.h | 9 ++---- cpp/tests/test_mobility.cpp | 28 +++++++++++++++++-- cpp/tests/test_save_results.cpp | 26 +++++++++++++++-- 3 files changed, 52 insertions(+), 11 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 0723af59f1..fe66940ba3 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -334,7 +334,7 @@ class MigrationEdge , m_return_times(0) , m_return_migrated(false) , m_save_indices(params.get_save_indices()) - , m_mobility_results(m_save_indices.size()) + , m_mobility_results(m_save_indices.size() + 1) { } @@ -348,7 +348,7 @@ class MigrationEdge , m_return_times(0) , m_return_migrated(false) , m_save_indices(0) - , m_mobility_results(m_save_indices.size()) + , m_mobility_results(m_save_indices.size() + 1) { } @@ -424,7 +424,7 @@ class MigrationEdge 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_save_indices; // groups of indices from compartments to save - TimeSeries m_mobility_results; // save results from edges + total number of commuters + TimeSeries m_mobility_results; // save results from edges + entry for the total number of commuters /** * Computes a condensed version of m_migrated and puts it in m_mobility_results. @@ -572,9 +572,6 @@ void MigrationEdge::apply_migration(double t, double dt, SimulationNode& no m_t_last_dynamic_npi_check = t; } - static auto indices_tuple = get_indices_of_symptomatic_and_nonsymptomatic(node_from); - auto& [indices_no_symptoms, indices_symptoms] = indices_tuple; - //returns for (Eigen::Index i = m_return_times.get_num_time_points() - 1; i >= 0; --i) { if (m_return_times.get_time(i) <= t) { diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp index 0d0991d3c7..495249a52b 100644 --- a/cpp/tests/test_mobility.cpp +++ b/cpp/tests/test_mobility.cpp @@ -175,7 +175,8 @@ TEST(TestMobility, condense_m_mobility) using Model = mio::osecir::Model; //setup nodes - Model model(1); + 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; @@ -193,11 +194,32 @@ TEST(TestMobility, condense_m_mobility) 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::MigrationEdge edge1(Eigen::VectorXd::Constant(10, 0.1)); + mio::MigrationEdge edge1(Eigen::VectorXd::Constant(10, 0.1), indices_save_edges); edge1.apply_migration(t, 0.0, node1, node2); auto migrated = edge1.get_migrated().get_last_value(); EXPECT_NEAR(migrated[0], 1.0, 1e-12); @@ -208,7 +230,7 @@ TEST(TestMobility, condense_m_mobility) model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 30; mio::SimulationNode> node3(model, t); mio::SimulationNode> node4(model, t); - mio::MigrationEdge edge2(Eigen::VectorXd::Constant(10, 0.1)); + mio::MigrationEdge edge2(Eigen::VectorXd::Constant(10, 0.1), indices_save_edges); edge2.apply_migration(t, 0.5, node3, node4); migrated = edge2.get_migrated().get_last_value(); EXPECT_NEAR(migrated[0], 11.0, 1e-12); diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index c780c081f7..342bb4054f 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -235,7 +235,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; @@ -274,12 +274,34 @@ 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(); graph.add_node(0, model); graph.add_node(1, model); - graph.add_edge(0, 1, mio::MigrationParameters(Eigen::VectorXd::Constant(Eigen::Index(num_groups * 10), 1.0))); + graph.add_edge(0, 1, Eigen::VectorXd::Constant(num_groups * (size_t)mio::osecir::InfectionState::Count, 0.1), + indices_save_edges); auto num_runs = 3; auto parameter_study = mio::ParameterStudy>(graph, 0.0, 2.0, 0.5, num_runs); From 71cfe00608ba47aecbc645a041d2f81c53f8c4f6 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Thu, 25 Apr 2024 16:09:50 +0200 Subject: [PATCH 14/30] add feature to set_edges --- cpp/memilio/mobility/graph.h | 11 +++- .../2020_npis_sarscov2_wildtype_germany.cpp | 3 +- ...021_vaccination_sarscov2_delta_germany.cpp | 3 +- cpp/tests/test_graph.cpp | 59 ++++++++++++++++++- 4 files changed, 69 insertions(+), 7 deletions(-) diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index 7e7c093595..e9c45185d3 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -335,8 +335,8 @@ template IOResult set_edges(const fs::path& data_dir, Graph& params_graph, std::initializer_list& migrating_compartments, size_t contact_locations_size, - ReadFunction&& read_func, - std::vector commuting_weights = std::vector{}) + ReadFunction&& read_func, std::vector commuting_weights, + bool save_results_edges = false, std::vector> indices_save_edges = {}) { // mobility between nodes BOOST_OUTCOME_TRY(mobility_data_commuter, @@ -388,7 +388,12 @@ 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)); + if (save_results_edges) { + params_graph.add_edge(county_idx_i, county_idx_j, std::move(mobility_coeffs), indices_save_edges); + } + else { + params_graph.add_edge(county_idx_i, county_idx_j, std::move(mobility_coeffs)); + } } } } diff --git a/cpp/simulations/2020_npis_sarscov2_wildtype_germany.cpp b/cpp/simulations/2020_npis_sarscov2_wildtype_germany.cpp index e3a2a3bef7..64785b9aad 100644 --- a/cpp/simulations/2020_npis_sarscov2_wildtype_germany.cpp +++ b/cpp/simulations/2020_npis_sarscov2_wildtype_germany.cpp @@ -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, migrating_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.}, + false, {})); return mio::success(params_graph); } diff --git a/cpp/simulations/2021_vaccination_sarscov2_delta_germany.cpp b/cpp/simulations/2021_vaccination_sarscov2_delta_germany.cpp index 772d14214e..f1d5df910a 100644 --- a/cpp/simulations/2021_vaccination_sarscov2_delta_germany.cpp +++ b/cpp/simulations/2021_vaccination_sarscov2_delta_germany.cpp @@ -579,7 +579,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, migrating_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.}, + false, {})); return mio::success(params_graph); } diff --git a/cpp/tests/test_graph.cpp b/cpp/tests/test_graph.cpp index 4d292c6004..d695db54c5 100644 --- a/cpp/tests/test_graph.cpp +++ b/cpp/tests/test_graph.cpp @@ -34,6 +34,7 @@ #include "matchers.h" #include "memilio/utils/stl_util.h" #include "gmock/gmock-matchers.h" +#include #include #include #include @@ -269,6 +270,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 params_graph; + const fs::path& dir = " "; + auto migrating_compartments = {mio::osecir::InfectionState::Susceptible, mio::osecir::InfectionState::Exposed, + mio::osecir::InfectionState::InfectedNoSymptoms, + mio::osecir::InfectionState::InfectedSymptoms, + mio::osecir::InfectionState::Recovered}; + + bool save_results_edges = true; + // 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( + dir, params_graph, migrating_compartments, size_t(2), read_function_edges, + std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}, save_results_edges, 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; @@ -293,10 +348,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; From d30c44cd0bb82d4c2d2d27cb23e59d92eb48d6ea Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Thu, 25 Apr 2024 22:18:14 +0200 Subject: [PATCH 15/30] clean up and dir in example --- cpp/examples/ode_secir_parameter_study_graph.cpp | 9 ++++++++- cpp/memilio/mobility/metapopulation_mobility_instant.cpp | 1 - 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp index 895966f45c..4d9c16539d 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -316,8 +316,15 @@ int main() 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, "test_results", false); + 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]}); diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp index f47eeccabd..108f3f5e66 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp @@ -24,7 +24,6 @@ namespace mio { void MigrationEdge::condense_m_mobility(const double t) { - // only call this if his->m_save_indices.size() is greater than 0. Perfect would be to define this in compile time const size_t save_indices_size = this->m_save_indices.size(); if (save_indices_size > 0) { From e6b1f413a0228cd3308f7f0e222244a20dc15027 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 24 May 2024 14:54:50 +0200 Subject: [PATCH 16/30] fix build after merge --- cpp/examples/ode_secir_graph.cpp | 44 +++++++-------- .../ode_secir_parameter_study_graph.cpp | 54 +++++++++---------- .../metapopulation_mobility_instant.cpp | 25 --------- .../metapopulation_mobility_instant.h | 28 +++++++++- cpp/tests/test_graph.cpp | 4 +- cpp/tests/test_mobility.cpp | 16 +++--- cpp/tests/test_save_results.cpp | 31 ++++++----- 7 files changed, 102 insertions(+), 100 deletions(-) diff --git a/cpp/examples/ode_secir_graph.cpp b/cpp/examples/ode_secir_graph.cpp index 48ed48b0b8..51584f855d 100644 --- a/cpp/examples/ode_secir_graph.cpp +++ b/cpp/examples/ode_secir_graph.cpp @@ -17,6 +17,7 @@ * 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" @@ -35,32 +36,33 @@ int main() 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(); + 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(); + 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 @@ -88,7 +90,7 @@ int main() model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptomsConfirmed})); } - mio::Graph>, mio::MigrationEdge> g; + mio::Graph>, mio::MigrationEdge> 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); diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp index 4d9c16539d..34542f3318 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -44,9 +44,9 @@ * @param min minimum of distribution. * @param max minimum of distribution. */ -void assign_uniform_distribution(mio::UncertainValue& p, double min, double max) +void assign_uniform_distribution(mio::UncertainValue& p, double min, double max) { - p = mio::UncertainValue(0.5 * (max + min)); + p = mio::UncertainValue(0.5 * (max + min)); p.set_distribution(mio::ParameterDistributionUniform(min, max)); } @@ -59,7 +59,7 @@ void assign_uniform_distribution(mio::UncertainValue& p, double min, double max) * @param max minimum of distribution for each element of array. */ template -void array_assign_uniform_distribution(mio::CustomIndexArray& array, +void array_assign_uniform_distribution(mio::CustomIndexArray, mio::AgeGroup>& array, const double (&min)[N], const double (&max)[N]) { assert(N == array.numel()); @@ -75,8 +75,8 @@ void array_assign_uniform_distribution(mio::CustomIndexArray& array, double min, - double max) +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); @@ -89,7 +89,7 @@ void array_assign_uniform_distribution(mio::CustomIndexArray& params) { //times // TimeExposed and TimeInfectedNoSymptoms are calculated as described in @@ -107,14 +107,14 @@ void set_covid_parameters(mio::osecir::Parameters& params) 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, + 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, + array_assign_uniform_distribution(params.get>(), timeInfectedSevereMin, timeInfectedSevereMax); - array_assign_uniform_distribution(params.get(), timeInfectedCriticalMin, + array_assign_uniform_distribution(params.get>(), timeInfectedCriticalMin, timeInfectedCriticalMax); //probabilities @@ -137,28 +137,28 @@ void set_covid_parameters(mio::osecir::Parameters& params) 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(), + array_assign_uniform_distribution(params.get>(), transmissionProbabilityOnContactMin, transmissionProbabilityOnContactMax); - array_assign_uniform_distribution(params.get(), + array_assign_uniform_distribution(params.get>(), relativeTransmissionNoSymptomsMin, relativeTransmissionNoSymptomsMax); - array_assign_uniform_distribution(params.get(), + array_assign_uniform_distribution(params.get>(), riskOfInfectionFromSymptomaticMin, riskOfInfectionFromSymptomaticMax); - array_assign_uniform_distribution(params.get(), + array_assign_uniform_distribution(params.get>(), maxRiskOfInfectionFromSymptomaticMin, maxRiskOfInfectionFromSymptomaticMax); - array_assign_uniform_distribution(params.get(), + array_assign_uniform_distribution(params.get>(), recoveredPerInfectedNoSymptomsMin, recoveredPerInfectedNoSymptomsMax); - array_assign_uniform_distribution(params.get(), + array_assign_uniform_distribution(params.get>(), severePerInfectedSymptomsMin, severePerInfectedSymptomsMax); - array_assign_uniform_distribution(params.get(), criticalPerSevereMin, + array_assign_uniform_distribution(params.get>(), criticalPerSevereMin, criticalPerSevereMax); - array_assign_uniform_distribution(params.get(), deathsPerCriticalMin, + 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); + assign_uniform_distribution(params.get>(), seasonality_min, seasonality_max); params.set(0); } @@ -168,7 +168,7 @@ void set_covid_parameters(mio::osecir::Parameters& params) * 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) +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; @@ -188,7 +188,7 @@ void set_synthetic_population_data(mio::osecir::Model& model) } } -std::vector> get_indices_of_symptomatic_and_nonsymptomatic(mio::osecir::Model& model) +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()); @@ -231,7 +231,7 @@ int main() const auto num_days_sim = 30.0; const auto num_runs = 10; - mio::Graph params_graph; + mio::Graph, mio::MigrationParameters> params_graph; const int num_age_groups = 6; mio::osecir::Model model(num_age_groups); @@ -241,7 +241,7 @@ int main() // set contact matrix const auto cont_freq = 10.0; - mio::ContactMatrixGroup& contact_matrix = model.parameters.get(); + 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)); @@ -286,7 +286,7 @@ int main() [&](auto results_graph, auto&& run_id) { 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) { @@ -307,7 +307,7 @@ int main() if (ensemble.size() > 0) { auto ensemble_results = std::vector>>{}; ensemble_results.reserve(ensemble.size()); - auto ensemble_params = std::vector>{}; + auto ensemble_params = std::vector>>{}; ensemble_params.reserve(ensemble.size()); auto ensemble_edges = std::vector>>{}; ensemble_edges.reserve(ensemble.size()); diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp index 108f3f5e66..86f5ce4283 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp @@ -22,29 +22,4 @@ namespace mio { -void MigrationEdge::condense_m_mobility(const double t) -{ - const size_t save_indices_size = this->m_save_indices.size(); - if (save_indices_size > 0) { - - const auto& last_value = m_migrated.get_last_value(); - Eigen::VectorXd condensed_values = Eigen::VectorXd::Zero(save_indices_size + 1); - - // sum up the values of m_save_indices for each group (e.g. Age groups) - std::transform(this->m_save_indices.begin(), this->m_save_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_migrated.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)); - } -} - } // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 1684ec55ef..c458cdb2d7 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -360,7 +360,7 @@ class MigrationEdge * @param coeffs % of people in each group and compartment that migrate in each time step. * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. */ - MigrationEdge(const MigrationParameters& params, const std::vector>& save_indices) + MigrationEdge(const MigrationParameters& params, const std::vector>& save_indices) : m_parameters(params) , m_migrated(params.get_coefficients().get_shape().rows()) , m_return_times(0) @@ -438,6 +438,32 @@ class MigrationEdge void condense_m_mobility(const double t); }; +template +void MigrationEdge::condense_m_mobility(const double t) +{ + const size_t save_indices_size = this->m_save_indices.size(); + if (save_indices_size > 0) { + + const auto& last_value = m_migrated.get_last_value(); + Eigen::VectorXd condensed_values = Eigen::VectorXd::Zero(save_indices_size + 1); + + // sum up the values of m_save_indices for each group (e.g. Age groups) + std::transform(this->m_save_indices.begin(), this->m_save_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_migrated.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 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. diff --git a/cpp/tests/test_graph.cpp b/cpp/tests/test_graph.cpp index b0ab791b14..794d31c888 100644 --- a/cpp/tests/test_graph.cpp +++ b/cpp/tests/test_graph.cpp @@ -279,7 +279,7 @@ 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 params_graph; + mio::Graph, mio::MigrationParameters> params_graph; const fs::path& dir = " "; auto migrating_compartments = {mio::osecir::InfectionState::Susceptible, mio::osecir::InfectionState::Exposed, mio::osecir::InfectionState::InfectedNoSymptoms, @@ -314,7 +314,7 @@ TEST(TestGraph, set_edges_saving_edges) const auto& read_function_edges = mock_read_mobility; auto result = - mio::set_edges, mio::MigrationParameters, mio::MigrationCoefficientGroup, mio::osecir::InfectionState, decltype(read_function_edges)>( dir, params_graph, migrating_compartments, size_t(2), read_function_edges, std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}, save_results_edges, indices_save_edges); diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp index 3f2d37d9c4..72b5a23a17 100644 --- a/cpp/tests/test_mobility.cpp +++ b/cpp/tests/test_mobility.cpp @@ -169,13 +169,13 @@ TEST(TestMobility, edgeApplyMigration) TEST(TestMobility, condense_m_mobility) { - using Model = mio::osecir::Model; + 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()); + 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; @@ -183,12 +183,12 @@ TEST(TestMobility, condense_m_mobility) 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.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. diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index ab665f983f..0c8ef305f6 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -296,8 +296,7 @@ TEST(TestSaveResult, save_percentiles_and_sums) auto graph = mio::Graph, mio::MigrationParameters>(); graph.add_node(0, model); graph.add_node(1, model); - graph.add_edge(0, 1, - mio::MigrationParameters(Eigen::VectorXd::Constant(Eigen::Index(num_groups * 10), 1.0)), + graph.add_edge(0, 1, Eigen::VectorXd::Constant(num_groups * (size_t)mio::osecir::InfectionState::Count, 0.1), indices_save_edges); auto num_runs = 3; @@ -419,11 +418,11 @@ TEST(TestSaveResult, save_edges) ; for (auto i = mio::AgeGroup(0); i < nb_groups; i++) { - params.get()[i] = 3.2; - params.get()[i] = 2; - params.get()[i] = 5.; - params.get()[i] = 10.; - params.get()[i] = 8.; + params.get>()[i] = 3.2; + params.get>()[i] = 2; + params.get>()[i] = 5.; + params.get>()[i] = 10.; + params.get>()[i] = 8.; model.populations[{i, mio::osecir::InfectionState::Exposed}] = 100; model.populations[{i, mio::osecir::InfectionState::InfectedNoSymptoms}] = 50; @@ -434,21 +433,21 @@ TEST(TestSaveResult, save_edges) model.populations[{i, mio::osecir::InfectionState::Dead}] = 0; model.populations.set_difference_from_total({i, mio::osecir::InfectionState::Susceptible}, pop_total); - params.get()[i] = 0.06; - params.get()[i] = 0.67; - params.get()[i] = 0.09; - params.get()[i] = 0.25; - params.get()[i] = 0.2; - params.get()[i] = 0.25; - params.get()[i] = 0.3; + params.get>()[i] = 0.06; + params.get>()[i] = 0.67; + params.get>()[i] = 0.09; + params.get>()[i] = 0.25; + params.get>()[i] = 0.2; + params.get>()[i] = 0.25; + params.get>()[i] = 0.3; } - mio::ContactMatrixGroup& contact_matrix = params.get(); + mio::ContactMatrixGroup& contact_matrix = params.get>(); contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant((size_t)nb_groups, (size_t)nb_groups, cont_freq)); contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); // setup graph - mio::Graph>, mio::MigrationEdge> g; + mio::Graph>, mio::MigrationEdge> g; g.add_node(0, model, t0); g.add_node(1, model, t0); g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.1)); From ff137c59a9cd2b24c88634698efc851261eb6a4d Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 24 May 2024 15:04:35 +0200 Subject: [PATCH 17/30] remove old doc --- cpp/memilio/io/result_io.cpp | 1 - cpp/memilio/io/result_io.h | 4 +-- .../metapopulation_mobility_instant.h | 35 ------------------- 3 files changed, 1 insertion(+), 39 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index a5109856e6..c8c36153a6 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -18,7 +18,6 @@ * limitations under the License. */ #include "memilio/io/result_io.h" -#include #ifdef MEMILIO_HAS_HDF5 diff --git a/cpp/memilio/io/result_io.h b/cpp/memilio/io/result_io.h index d80aa90075..19de26f29b 100644 --- a/cpp/memilio/io/result_io.h +++ b/cpp/memilio/io/result_io.h @@ -125,7 +125,6 @@ IOResult save_result_with_params(const std::vector>& re * @brief Save the results of the edges for a single graph simulation run. * @param result Simulation results per edge of the graph. * @param ids Identifiers for the start and end node of the edges. - * @param num_groups Number of groups in the results. * @param filename Name of file * @return Any io errors that occur during writing of the files. */ @@ -135,8 +134,7 @@ IOResult save_edges(const std::vector>& results, const /** * Save the results for the Edges obtained from the function condense_m_mobility. * @param ensemble_edges Simulation results for each run for each edge. - * @param num_groups Number of age groups used simulation. - * @param ids Identifiers for the start and end node of the edges. + * @param pairs_edges Identifiers for the start and end node of the edges. * @param result_dir Top level directory for all results of the parameter study. * @param save_single_runs [Default: true] Defines if single run results are written. * @param save_single_runs [Default: true] Defines if percentiles are written. diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index c458cdb2d7..6b01e0c13c 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -285,41 +285,6 @@ class MigrationParameters std::vector> m_save_indices; // groups of indices from compartments to save }; -/** - * detect a get_indices_of_symptomatic_and_nonsymptomatic function for the Model type. - */ -template -using get_indices_of_symptomatic_and_nonsymptomatic_expr_t = - decltype(get_indices_of_symptomatic_and_nonsymptomatic(std::declval())); - -/** - * @brief Get the indices of symptomatic and non-symptomatic infection states. - * - * This function generates two vectors of indices, one for non-symptomatic infection states and one for symptomatic infection states. - * Each vector contains the flat indices of the corresponding infection states for each age group in the model. - * - * @tparam Base The base class for the simulation, defaults to mio::Simulation. - * @param[in] sim The simulation object from which we obtain the model. - * - * @return A tuple containing two vectors of size_t. The first vector contains the indices of non-symptomatic infection states, - * and the second vector contains the indices of symptomatic infection states. The indices are ordered first by age group. - */ -template ::value, - void*> = nullptr> -auto get_indices_of_symptomatic_and_nonsymptomatic(SimulationNode& /*node*/) -{ - return std::make_tuple(std::vector{}, std::vector{}); -} - -template ::value, - void*> = nullptr> -auto get_indices_of_symptomatic_and_nonsymptomatic(SimulationNode& node) -{ - return get_indices_of_symptomatic_and_nonsymptomatic(node.get_simulation()); -} - /** * represents the migration between two nodes. */ From 2732c62e29187c16617aa3f9f52247ff2e894e7c Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 10 Jun 2024 14:43:45 +0200 Subject: [PATCH 18/30] comments parameters, doc save_edges, rm int casts --- cpp/examples/ode_secir_parameter_study_graph.cpp | 3 --- cpp/memilio/io/result_io.cpp | 9 ++++++--- cpp/memilio/io/result_io.h | 2 +- cpp/memilio/mobility/metapopulation_mobility_instant.h | 2 +- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp index 34542f3318..753aa6155f 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -92,9 +92,6 @@ void array_assign_uniform_distribution(mio::CustomIndexArray& params) { //times - // TimeExposed and TimeInfectedNoSymptoms are calculated as described in - // Khailaie et al. (https://doi.org/10.1186/s12916-020-01884-4) - // given SI_min = 3.935, SI_max = 4.6, INC = 5.2 const double timeExposedMin = 2.67; const double timeExposedMax = 4.; const double timeInfectedNoSymptomsMin = 1.2; diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index c8c36153a6..39c5faa111 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -136,8 +136,11 @@ IOResult save_edges(const std::vector>>& en IOResult save_edges(const std::vector>& results, const std::vector>& ids, const std::string& filename) { - const int num_edges = static_cast(results.size()); - int edge_indx = 0; + 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, filename); while (edge_indx < num_edges) { @@ -146,7 +149,7 @@ IOResult save_edges(const std::vector>& results, const 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 + ")"); - const int num_timepoints = static_cast(result.get_num_time_points()); + const auto num_timepoints = result.get_num_time_points(); constexpr int num_infectionstates = 3; // (int)result.get_num_elements() / num_groups; hsize_t dims_t[] = {static_cast(num_timepoints)}; diff --git a/cpp/memilio/io/result_io.h b/cpp/memilio/io/result_io.h index 19de26f29b..487e7fcc1b 100644 --- a/cpp/memilio/io/result_io.h +++ b/cpp/memilio/io/result_io.h @@ -137,7 +137,7 @@ IOResult save_edges(const std::vector>& results, const * @param pairs_edges Identifiers for the start and end node of the edges. * @param result_dir Top level directory for all results of the parameter study. * @param save_single_runs [Default: true] Defines if single run results are written. - * @param save_single_runs [Default: true] Defines if percentiles are written. + * @param 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, diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 6b01e0c13c..c021b033ba 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -359,7 +359,7 @@ class MigrationEdge } /** - * Retrieve the count of commuters in the infection states: InfectedNoSymptoms and InfectedSymptomsNaive, + * Retrieve the count of commuters in selected infection states, * along with the total number of commuter. */ TimeSeries& get_migrated() From 1680ef0442ec5f5c3b4261806b8cd31403c4c6da Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 10 Jun 2024 14:54:53 +0200 Subject: [PATCH 19/30] [ci skip] dynamic num of elements --- cpp/memilio/io/result_io.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 39c5faa111..21e0bf07fa 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -149,8 +149,8 @@ IOResult save_edges(const std::vector>& results, const 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 + ")"); - const auto num_timepoints = result.get_num_time_points(); - constexpr int num_infectionstates = 3; // (int)result.get_num_elements() / num_groups; + const auto num_timepoints = result.get_num_time_points(); + 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)}; @@ -163,13 +163,13 @@ IOResult save_edges(const std::vector>& results, const StatusCode::UnknownError, "Time data could not be written."); int start_id = ids[edge_indx].first; - auto total = Eigen::Matrix::Zero(num_timepoints, - num_infectionstates) - .eval(); + 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_infectionstates) + 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(); @@ -177,7 +177,7 @@ IOResult save_edges(const std::vector>& results, const total.row(t_idx) += v; } - hsize_t dims_values[] = {static_cast(num_timepoints), static_cast(num_infectionstates)}; + 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, "Values DataSpace could not be created."); auto dset_name = "End" + std::to_string(ids[edge_indx].second); From 000f18f2a18e4cd17d15a622d158c891bbc8ff0a Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 10 Jun 2024 15:13:06 +0200 Subject: [PATCH 20/30] Error messages in save_edges --- cpp/memilio/io/result_io.cpp | 40 ++++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 21e0bf07fa..11ae623810 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -142,25 +142,33 @@ IOResult save_edges(const std::vector>& results, const // 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, filename); + 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 + ")"); + "Group could not be created (" + h5group_name + ") in the file: " + filename); + const auto num_timepoints = result.get_num_time_points(); 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, "Time DataSpace could not be created."); + 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, "Time DataSet could not be created (Time)."); + 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, "Time data could not be written."); + StatusCode::UnknownError, + "Failed to write 'Time' data in group " + h5group_name + " in the file: " + filename); int start_id = ids[edge_indx].first; auto total = @@ -179,25 +187,35 @@ IOResult save_edges(const std::vector>& results, const 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, "Values DataSpace could not be created."); + 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); + 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, "Values DataSet could not be created."); + 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, "Values data could not be written."); + 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 + // 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, "Values DataSet could not be created."); + 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, "Values data could not be written."); + StatusCode::UnknownError, + "Failed to write Total data in group " + h5group_name + " in the file: " + filename); } edge_indx++; } From ba53065c7612166d06cff06e517c963fc61e906b Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 10 Jun 2024 15:22:57 +0200 Subject: [PATCH 21/30] comment about group end, rm include, doc save_edges --- cpp/memilio/io/result_io.cpp | 2 ++ cpp/memilio/io/result_io.h | 2 +- cpp/memilio/mobility/metapopulation_mobility_instant.cpp | 1 - 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 11ae623810..4119df77d6 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -89,6 +89,7 @@ IOResult save_result(const std::vector>& results, const } return success(); } + 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) @@ -191,6 +192,7 @@ IOResult save_edges(const std::vector>& results, const "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)}; diff --git a/cpp/memilio/io/result_io.h b/cpp/memilio/io/result_io.h index 487e7fcc1b..54fe99dddf 100644 --- a/cpp/memilio/io/result_io.h +++ b/cpp/memilio/io/result_io.h @@ -132,7 +132,7 @@ IOResult save_edges(const std::vector>& results, const const std::string& filename); /** - * Save the results for the Edges obtained from the function condense_m_mobility. + * Saves the results of a simulation for each edge in the graph. * @param ensemble_edges Simulation results for each run for each edge. * @param pairs_edges Identifiers for the start and end node of the edges. * @param result_dir Top level directory for all results of the parameter study. diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp index 86f5ce4283..fc54caebec 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.cpp +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.cpp @@ -18,7 +18,6 @@ * limitations under the License. */ #include "memilio/mobility/metapopulation_mobility_instant.h" -#include "memilio/utils/compiler_diagnostics.h" namespace mio { From 1e5dbe05ab90aa64f65b7ab22a640ddca69ee3a9 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Tue, 11 Jun 2024 10:56:49 +0200 Subject: [PATCH 22/30] fix tests by setting indices and case of no indices --- cpp/memilio/io/result_io.cpp | 119 +++++++++++++++++--------------- cpp/tests/test_save_results.cpp | 11 ++- 2 files changed, 71 insertions(+), 59 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 4119df77d6..a046be3fad 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -153,72 +153,79 @@ IOResult save_edges(const std::vector>& results, const "Group could not be created (" + h5group_name + ") in the file: " + filename); const auto num_timepoints = result.get_num_time_points(); - const auto num_elements = result.get_value(0).size(); + 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); + 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); + 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; - } + 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); - 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); + 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; + } - 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 + + 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_total.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, total.data()), + H5Dwrite(dset_values.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, edge_result.data()), StatusCode::UnknownError, - "Failed to write Total data in group " + h5group_name + " in the file: " + filename); + "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 { edge_indx++; } } diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index 0c8ef305f6..84b47a0438 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -415,7 +415,6 @@ TEST(TestSaveResult, save_edges) mio::osecir::Model model(1); auto& params = model.parameters; mio::AgeGroup nb_groups = params.get_num_groups(); - ; for (auto i = mio::AgeGroup(0); i < nb_groups; i++) { params.get>()[i] = 3.2; @@ -450,8 +449,14 @@ TEST(TestSaveResult, save_edges) mio::Graph>, mio::MigrationEdge> g; g.add_node(0, model, t0); g.add_node(1, model, t0); - g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.1)); - g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.1)); + + // get indices of INS and ISy compartments. + std::vector> indices_save_edges(2); + indices_save_edges[0] = {2, 3}; + indices_save_edges[1] = {4, 5}; + + 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); // simulation auto sim = mio::make_migration_sim(t0, dt, std::move(g)); From fe1b67f26ea6b68dcb971805a665190125179b16 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Tue, 11 Jun 2024 14:38:53 +0200 Subject: [PATCH 23/30] adjust test for save_edges --- cpp/tests/test_save_results.cpp | 132 ++++++++++++-------------------- 1 file changed, 49 insertions(+), 83 deletions(-) diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index 84b47a0438..c6f7eb2bc8 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -405,99 +405,65 @@ TEST(TestSaveResult, save_percentiles_and_sums) TEST(TestSaveResult, save_edges) { - double t0 = 0; - double tmax = 5; - double dt = 0.1; - double cont_freq = 10.; + // 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)); - double pop_total = 10000; + 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)); - mio::osecir::Model model(1); - auto& params = model.parameters; - mio::AgeGroup nb_groups = params.get_num_groups(); - - for (auto i = mio::AgeGroup(0); i < nb_groups; i++) { - params.get>()[i] = 3.2; - params.get>()[i] = 2; - params.get>()[i] = 5.; - params.get>()[i] = 10.; - params.get>()[i] = 8.; - - model.populations[{i, mio::osecir::InfectionState::Exposed}] = 100; - model.populations[{i, mio::osecir::InfectionState::InfectedNoSymptoms}] = 50; - model.populations[{i, mio::osecir::InfectionState::InfectedSymptoms}] = 50; - model.populations[{i, mio::osecir::InfectionState::InfectedSevere}] = 20; - model.populations[{i, mio::osecir::InfectionState::InfectedCritical}] = 10; - model.populations[{i, mio::osecir::InfectionState::Recovered}] = 10; - model.populations[{i, mio::osecir::InfectionState::Dead}] = 0; - model.populations.set_difference_from_total({i, mio::osecir::InfectionState::Susceptible}, pop_total); - - params.get>()[i] = 0.06; - params.get>()[i] = 0.67; - params.get>()[i] = 0.09; - params.get>()[i] = 0.25; - params.get>()[i] = 0.2; - params.get>()[i] = 0.25; - params.get>()[i] = 0.3; - } - - mio::ContactMatrixGroup& contact_matrix = params.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant((size_t)nb_groups, (size_t)nb_groups, cont_freq)); - contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); - - // setup graph - mio::Graph>, mio::MigrationEdge> g; - g.add_node(0, model, t0); - g.add_node(1, model, t0); - - // get indices of INS and ISy compartments. - std::vector> indices_save_edges(2); - indices_save_edges[0] = {2, 3}; - indices_save_edges[1] = {4, 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)); - 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); - - // simulation - auto sim = mio::make_migration_sim(t0, dt, std::move(g)); - sim.advance(tmax); - - // get_results and pair ids - auto& results_edge_0 = sim.get_graph().edges()[0].property.get_migrated(); - std::vector> results_edges = {results_edge_0, results_edge_0}; - - auto pairs_edges = std::vector>{}; - pairs_edges.push_back({0, 1}); - pairs_edges.push_back({1, 0}); + 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); - auto result_from_file = results_from_file.value()[0]; - ASSERT_EQ(result_from_file.get_groups().get_num_time_points(), results_edge_0.get_num_time_points()); - ASSERT_EQ(result_from_file.get_totals().get_num_time_points(), results_edge_0.get_num_time_points()); - for (Eigen::Index i = 0; i < results_edge_0.get_num_time_points(); i++) { - ASSERT_EQ(result_from_file.get_groups().get_num_elements(), results_edge_0.get_num_elements()) - << "at row " << i; - ASSERT_EQ(result_from_file.get_totals().get_num_elements(), - results_edge_0.get_num_elements() / static_cast((size_t)nb_groups)) - << "at row " << i; - ASSERT_NEAR(results_edge_0.get_time(i), result_from_file.get_groups().get_time(i), 1e-10) << "at row " << i; - ASSERT_NEAR(results_edge_0.get_time(i), result_from_file.get_totals().get_time(i), 1e-10) << "at row " << i; - for (Eigen::Index l = 0; l < result_from_file.get_totals().get_num_elements(); l++) { - double total = 0.0; - for (Eigen::Index j = 0; j < Eigen::Index((size_t)nb_groups); j++) { - total += results_edge_0[i][j * (size_t)mio::osecir::InfectionState::Count + l]; - EXPECT_NEAR(result_from_file.get_groups()[i][j * (size_t)mio::osecir::InfectionState::Count + l], - results_edge_0[i][j * (size_t)mio::osecir::InfectionState::Count + l], 1e-10) - << " at row " << i << " at row " << l << " at Group " << j; - } - EXPECT_NEAR(result_from_file.get_totals()[i][l], total, 1e-10) << " at row " << i << " at row " << l; - } - } + // group 0 + auto result_from_file_group0 = results_from_file.value()[0]; + ASSERT_EQ(result_from_file_group0.get_groups().get_num_time_points(), 3); + ASSERT_EQ(result_from_file_group0.get_groups().get_num_elements(), 6); + ASSERT_EQ(result_from_file_group0.get_groups().get_value(0), (Eigen::VectorXd(6) << 0, 0, 0, 2, 2, 2).finished()); + ASSERT_EQ(result_from_file_group0.get_groups().get_value(1), (Eigen::VectorXd(6) << 1, 1, 1, 3, 3, 3).finished()); + ASSERT_EQ(result_from_file_group0.get_groups().get_value(2), (Eigen::VectorXd(6) << 1, 1, 1, 5, 5, 5).finished()); + ASSERT_EQ(result_from_file_group0.get_groups().get_time(0), 0.0); + ASSERT_EQ(result_from_file_group0.get_groups().get_time(1), 1.0); + ASSERT_EQ(result_from_file_group0.get_groups().get_time(2), 2.0); + + ASSERT_EQ(result_from_file_group0.get_totals().get_num_time_points(), 3); + ASSERT_EQ(result_from_file_group0.get_totals().get_num_elements(), 3); + ASSERT_EQ(result_from_file_group0.get_totals().get_value(0), Eigen::VectorXd::Constant(3, 2)); + ASSERT_EQ(result_from_file_group0.get_totals().get_value(1), Eigen::VectorXd::Constant(3, 4)); + ASSERT_EQ(result_from_file_group0.get_totals().get_value(2), Eigen::VectorXd::Constant(3, 6)); + ASSERT_EQ(result_from_file_group0.get_totals().get_time(0), 0.0); + ASSERT_EQ(result_from_file_group0.get_totals().get_time(1), 1.0); + ASSERT_EQ(result_from_file_group0.get_totals().get_time(2), 2.0); + + // group 1 + auto result_from_file_group1 = results_from_file.value()[1]; + ASSERT_EQ(result_from_file_group1.get_groups().get_num_elements(), 3); + ASSERT_EQ(result_from_file_group1.get_groups().get_value(0), Eigen::VectorXd::Constant(3, 3)); + ASSERT_EQ(result_from_file_group1.get_groups().get_value(1), Eigen::VectorXd::Constant(3, 4)); + ASSERT_EQ(result_from_file_group1.get_groups().get_value(2), Eigen::VectorXd::Constant(3, 7)); + ASSERT_EQ(result_from_file_group1.get_groups().get_time(0), 0.0); + ASSERT_EQ(result_from_file_group1.get_groups().get_time(1), 1.0); + ASSERT_EQ(result_from_file_group1.get_groups().get_time(2), 2.0); + + ASSERT_EQ(result_from_file_group1.get_totals().get_num_elements(), 3); + ASSERT_EQ(result_from_file_group1.get_totals().get_value(0), Eigen::VectorXd::Constant(3, 3)); + ASSERT_EQ(result_from_file_group1.get_totals().get_value(1), Eigen::VectorXd::Constant(3, 4)); + ASSERT_EQ(result_from_file_group1.get_totals().get_value(2), Eigen::VectorXd::Constant(3, 7)); } \ No newline at end of file From e84e8ca3cd67b454ddd15117a0bc64c35acba4fa Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 24 Jun 2024 13:51:31 +0200 Subject: [PATCH 24/30] rm non const function --- cpp/memilio/mobility/metapopulation_mobility_instant.h | 7 ------- 1 file changed, 7 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index c021b033ba..209bcedcf9 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -215,13 +215,6 @@ class MigrationParameters return m_save_indices; } - const auto& get_save_indices() - { - return m_save_indices; - } - - /** @} */ - /** * Get/Set dynamic NPIs that are implemented when relative infections exceed thresholds. * This feature is optional. The simulation model needs to overload the get_infected_relative function. From 35a12111a1bbb0bba4aea356bacfc23ed20f8cf5 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 28 Aug 2024 14:41:08 +0200 Subject: [PATCH 25/30] migration -> mobility --- cpp/examples/ode_secir_graph.cpp | 8 ++--- .../ode_secir_parameter_study_graph.cpp | 4 +-- .../metapopulation_mobility_instant.h | 36 +++++++++---------- cpp/tests/test_graph.cpp | 10 +++--- cpp/tests/test_mobility.cpp | 24 ++++++------- cpp/tests/test_save_results.cpp | 4 +-- 6 files changed, 43 insertions(+), 43 deletions(-) diff --git a/cpp/examples/ode_secir_graph.cpp b/cpp/examples/ode_secir_graph.cpp index 51584f855d..8296954621 100644 --- a/cpp/examples/ode_secir_graph.cpp +++ b/cpp/examples/ode_secir_graph.cpp @@ -30,7 +30,7 @@ int main() { const auto t0 = 0.; const auto tmax = 30.; - const auto dt = 0.5; //time step of migration, daily migration every second step + 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); @@ -90,18 +90,18 @@ int main() model.populations.get_flat_index({i, mio::osecir::InfectionState::InfectedSymptomsConfirmed})); } - mio::Graph>, mio::MigrationEdge> g; + 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_migration_sim(t0, dt, std::move(g)); + 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_migrated(); + 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 index 753aa6155f..0a04a2787b 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -228,7 +228,7 @@ int main() const auto num_days_sim = 30.0; const auto num_runs = 10; - mio::Graph, mio::MigrationParameters> params_graph; + mio::Graph, mio::MobilityParameters> params_graph; const int num_age_groups = 6; mio::osecir::Model model(num_age_groups); @@ -294,7 +294,7 @@ int main() 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_migrated(); + return edge.property.get_mobility_results(); }); std::cout << "Run " << run_id << " done\n"; diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index a64bd6e73c..4f65df87d7 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -139,7 +139,7 @@ class MobilityParameters * @param coeffs mobility coefficients */ MobilityParameters(const Eigen::VectorXd& coeffs) - : m_coefficients({MigrationCoefficients(coeffs)}) + : m_coefficients({MobilityCoefficients(coeffs)}) , m_save_indices(0) { } @@ -147,10 +147,10 @@ class MobilityParameters /** * @brief Constructor for a new MobilityParameters object. * - * @param coeffs migration coefficients + * @param coeffs Mobility coefficients * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. */ - MobilityParameters(const MigrationCoefficientGroup& coeffs, const std::vector>& save_indices) + MobilityParameters(const MobilityCoefficientGroup& coeffs, const std::vector>& save_indices) : m_coefficients(coeffs) , m_save_indices(save_indices) { @@ -159,11 +159,11 @@ class MobilityParameters /** * @brief Constructor for a new MobilityParameters object. * - * @param coeffs migration coefficients + * @param coeffs Mobility coefficients * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. */ MobilityParameters(const Eigen::VectorXd& coeffs, const std::vector>& save_indices) - : m_coefficients({MigrationCoefficients(coeffs)}) + : m_coefficients({MobilityCoefficients(coeffs)}) , m_save_indices(save_indices) { } @@ -314,13 +314,13 @@ class MobilityEdge } /** - * create edge with coefficients as MigrationParameters object and a 2d vector of indices which determine which compartments we save. + * create edge with coefficients as MobilityParameters object and a 2d vector of indices which determine which compartments we save. * @param coeffs % of people in each group and compartment that migrate in each time step. * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. */ - MigrationEdge(const MigrationParameters& params, const std::vector>& save_indices) + MobilityEdge(const MobilityParameters& params, const std::vector>& save_indices) : m_parameters(params) - , m_migrated(params.get_coefficients().get_shape().rows()) + , m_mobile_population(params.get_coefficients().get_shape().rows()) , m_return_times(0) , m_return_mobile_population(false) , m_save_indices(save_indices) @@ -333,9 +333,9 @@ class MobilityEdge * @param coeffs % of people in each group and compartment that migrate in each time step. * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. */ - MigrationEdge(const Eigen::VectorXd& coeffs, const std::vector>& save_indices) + MobilityEdge(const Eigen::VectorXd& coeffs, const std::vector>& save_indices) : m_parameters(coeffs) - , m_migrated(coeffs.rows()) + , m_mobile_population(coeffs.rows()) , m_return_times(0) , m_return_mobile_population(false) , m_save_indices(save_indices) @@ -355,11 +355,11 @@ class MobilityEdge * Retrieve the count of commuters in selected infection states, * along with the total number of commuter. */ - TimeSeries& get_migrated() + TimeSeries& get_mobility_results() { return m_mobility_results; } - const TimeSeries& get_migrated() const + const TimeSeries& get_mobility_results() const { return m_mobility_results; } @@ -388,7 +388,7 @@ class MobilityEdge TimeSeries m_mobility_results; // save results from edges + entry for the total number of commuters /** - * Computes a condensed version of m_migrated and puts it in m_mobility_results. + * Computes a condensed version of m_mobile_population and puts it in m_mobility_results. * 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 current time @@ -397,12 +397,12 @@ class MobilityEdge }; template -void MigrationEdge::condense_m_mobility(const double t) +void MobilityEdge::condense_m_mobility(const double t) { const size_t save_indices_size = this->m_save_indices.size(); if (save_indices_size > 0) { - const auto& last_value = m_migrated.get_last_value(); + 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_save_indices for each group (e.g. Age groups) @@ -415,7 +415,7 @@ void MigrationEdge::condense_m_mobility(const double t) }); // the last value is the sum of commuters - condensed_values[save_indices_size] = m_migrated.get_last_value().sum(); + 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)); @@ -588,8 +588,8 @@ void MobilityEdge::apply_mobility(FP t, FP dt, SimulationNode& node_fro m_mobile_population[i](j) += remaining_after_return(j); } } - node_from.get_result().get_last_value() += m_migrated[i]; - node_to.get_result().get_last_value() -= m_migrated[i]; + node_from.get_result().get_last_value() += m_mobile_population[i]; + node_to.get_result().get_last_value() -= m_mobile_population[i]; condense_m_mobility(t); m_mobile_population.remove_time_point(i); m_return_times.remove_time_point(i); diff --git a/cpp/tests/test_graph.cpp b/cpp/tests/test_graph.cpp index e6193823ce..c3481f0e1b 100644 --- a/cpp/tests/test_graph.cpp +++ b/cpp/tests/test_graph.cpp @@ -278,9 +278,9 @@ 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::MigrationParameters> params_graph; + mio::Graph, mio::MobilityParameters> params_graph; const fs::path& dir = " "; - auto migrating_compartments = {mio::osecir::InfectionState::Susceptible, mio::osecir::InfectionState::Exposed, + auto mobile_compartments = {mio::osecir::InfectionState::Susceptible, mio::osecir::InfectionState::Exposed, mio::osecir::InfectionState::InfectedNoSymptoms, mio::osecir::InfectionState::InfectedSymptoms, mio::osecir::InfectionState::Recovered}; @@ -313,9 +313,9 @@ TEST(TestGraph, set_edges_saving_edges) const auto& read_function_edges = mock_read_mobility; auto result = - mio::set_edges, mio::MigrationParameters, - mio::MigrationCoefficientGroup, mio::osecir::InfectionState, decltype(read_function_edges)>( - dir, params_graph, migrating_compartments, size_t(2), read_function_edges, + 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.}, save_results_edges, indices_save_edges); EXPECT_EQ(params_graph.edges().size(), 2); diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp index 70590d908d..0d7b9b8a53 100644 --- a/cpp/tests/test_mobility.cpp +++ b/cpp/tests/test_mobility.cpp @@ -216,21 +216,21 @@ TEST(TestMobility, condense_m_mobility) double t = 0.; mio::SimulationNode> node1(model, t); mio::SimulationNode> node2(model, t); - mio::MigrationEdge edge1(Eigen::VectorXd::Constant(10, 0.1), indices_save_edges); - edge1.apply_migration(t, 0.0, node1, node2); - auto migrated = edge1.get_migrated().get_last_value(); - EXPECT_NEAR(migrated[0], 1.0, 1e-12); - EXPECT_NEAR(migrated[1], 2.0, 1e-12); - EXPECT_NEAR(migrated[2], 100.0, 1e-12); + 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::MigrationEdge edge2(Eigen::VectorXd::Constant(10, 0.1), indices_save_edges); - edge2.apply_migration(t, 0.5, node3, node4); - migrated = edge2.get_migrated().get_last_value(); - EXPECT_NEAR(migrated[0], 11.0, 1e-12); - EXPECT_NEAR(migrated[1], 5.0, 1e-12); - EXPECT_NEAR(migrated[2], 113.0, 1e-12); + 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 a743e76bc3..2870618880 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -297,7 +297,7 @@ TEST(TestSaveResult, save_percentiles_and_sums) 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); @@ -331,7 +331,7 @@ TEST(TestSaveResult, save_percentiles_and_sums) 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_migrated(); + return edge.property.get_mobility_results(); }); return 0; //function needs to return something From bd524477d029beac65190f86cf8251ffc6a76405 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 28 Aug 2024 14:58:26 +0200 Subject: [PATCH 26/30] rm bool, use expect instead of assert --- cpp/memilio/mobility/graph.h | 10 +-- .../2020_npis_sarscov2_wildtype_germany.cpp | 8 +-- ...021_vaccination_sarscov2_delta_germany.cpp | 26 ++++---- cpp/tests/test_graph.cpp | 3 +- cpp/tests/test_save_results.cpp | 61 ++++++++++--------- 5 files changed, 52 insertions(+), 56 deletions(-) diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index 176be5f795..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, - bool save_results_edges = false, std::vector> indices_save_edges = {}) + std::vector> indices_of_saved_edges = {}) { // mobility between nodes BOOST_OUTCOME_TRY(auto&& mobility_data_commuter, @@ -390,12 +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) { - if (save_results_edges) { - params_graph.add_edge(county_idx_i, county_idx_j, std::move(mobility_coeffs), indices_save_edges); - } - else { - 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/simulations/2020_npis_sarscov2_wildtype_germany.cpp b/cpp/simulations/2020_npis_sarscov2_wildtype_germany.cpp index bea79e9c1a..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 @@ -499,7 +499,7 @@ get_graph(mio::Date start_date, mio::Date end_date, const fs::path& data_dir) 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.}, - false, {})); + {})); return mio::success(params_graph); } @@ -575,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 f64b7f651b..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 @@ -581,7 +581,7 @@ get_graph(mio::Date start_date, mio::Date end_date, const fs::path& data_dir, bo 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.}, - false, {})); + {})); return mio::success(params_graph); } @@ -663,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 c3481f0e1b..7ee8df12f9 100644 --- a/cpp/tests/test_graph.cpp +++ b/cpp/tests/test_graph.cpp @@ -285,7 +285,6 @@ TEST(TestGraph, set_edges_saving_edges) mio::osecir::InfectionState::InfectedSymptoms, mio::osecir::InfectionState::Recovered}; - bool save_results_edges = true; // get indices of INS and ISy compartments. std::vector> indices_save_edges(2); @@ -316,7 +315,7 @@ TEST(TestGraph, set_edges_saving_edges) 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.}, save_results_edges, indices_save_edges); + std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}, indices_save_edges); EXPECT_EQ(params_graph.edges().size(), 2); diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index 2870618880..bebb5fac2e 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -297,7 +297,8 @@ TEST(TestSaveResult, save_percentiles_and_sums) 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), indices_save_edges)); + 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); @@ -434,36 +435,36 @@ TEST(TestSaveResult, save_edges) // group 0 auto result_from_file_group0 = results_from_file.value()[0]; - ASSERT_EQ(result_from_file_group0.get_groups().get_num_time_points(), 3); - ASSERT_EQ(result_from_file_group0.get_groups().get_num_elements(), 6); - ASSERT_EQ(result_from_file_group0.get_groups().get_value(0), (Eigen::VectorXd(6) << 0, 0, 0, 2, 2, 2).finished()); - ASSERT_EQ(result_from_file_group0.get_groups().get_value(1), (Eigen::VectorXd(6) << 1, 1, 1, 3, 3, 3).finished()); - ASSERT_EQ(result_from_file_group0.get_groups().get_value(2), (Eigen::VectorXd(6) << 1, 1, 1, 5, 5, 5).finished()); - ASSERT_EQ(result_from_file_group0.get_groups().get_time(0), 0.0); - ASSERT_EQ(result_from_file_group0.get_groups().get_time(1), 1.0); - ASSERT_EQ(result_from_file_group0.get_groups().get_time(2), 2.0); - - ASSERT_EQ(result_from_file_group0.get_totals().get_num_time_points(), 3); - ASSERT_EQ(result_from_file_group0.get_totals().get_num_elements(), 3); - ASSERT_EQ(result_from_file_group0.get_totals().get_value(0), Eigen::VectorXd::Constant(3, 2)); - ASSERT_EQ(result_from_file_group0.get_totals().get_value(1), Eigen::VectorXd::Constant(3, 4)); - ASSERT_EQ(result_from_file_group0.get_totals().get_value(2), Eigen::VectorXd::Constant(3, 6)); - ASSERT_EQ(result_from_file_group0.get_totals().get_time(0), 0.0); - ASSERT_EQ(result_from_file_group0.get_totals().get_time(1), 1.0); - ASSERT_EQ(result_from_file_group0.get_totals().get_time(2), 2.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]; - ASSERT_EQ(result_from_file_group1.get_groups().get_num_elements(), 3); - ASSERT_EQ(result_from_file_group1.get_groups().get_value(0), Eigen::VectorXd::Constant(3, 3)); - ASSERT_EQ(result_from_file_group1.get_groups().get_value(1), Eigen::VectorXd::Constant(3, 4)); - ASSERT_EQ(result_from_file_group1.get_groups().get_value(2), Eigen::VectorXd::Constant(3, 7)); - ASSERT_EQ(result_from_file_group1.get_groups().get_time(0), 0.0); - ASSERT_EQ(result_from_file_group1.get_groups().get_time(1), 1.0); - ASSERT_EQ(result_from_file_group1.get_groups().get_time(2), 2.0); - - ASSERT_EQ(result_from_file_group1.get_totals().get_num_elements(), 3); - ASSERT_EQ(result_from_file_group1.get_totals().get_value(0), Eigen::VectorXd::Constant(3, 3)); - ASSERT_EQ(result_from_file_group1.get_totals().get_value(1), Eigen::VectorXd::Constant(3, 4)); - ASSERT_EQ(result_from_file_group1.get_totals().get_value(2), Eigen::VectorXd::Constant(3, 7)); + 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)); } \ No newline at end of file From e7437de6fe5e1f253e9e1694cc63d38687593776 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 28 Aug 2024 16:40:21 +0200 Subject: [PATCH 27/30] move to mobility_io, adjust namings, doc --- cpp/memilio/io/mobility_io.cpp | 143 ++++++++++++++++++ cpp/memilio/io/mobility_io.h | 26 ++++ cpp/memilio/io/result_io.cpp | 142 ----------------- cpp/memilio/io/result_io.h | 23 --- .../metapopulation_mobility_instant.h | 85 ++++++----- cpp/tests/test_mobility.cpp | 2 +- 6 files changed, 219 insertions(+), 202 deletions(-) diff --git a/cpp/memilio/io/mobility_io.cpp b/cpp/memilio/io/mobility_io.cpp index 1bd9fc99d0..5ad86a6558 100644 --- a/cpp/memilio/io/mobility_io.cpp +++ b/cpp/memilio/io/mobility_io.cpp @@ -1,6 +1,7 @@ #include "memilio/io/mobility_io.h" #include "memilio/math/eigen.h" #include "memilio/utils/logging.h" +#include "memilio/io/hdf5_cpp.h" #include #include @@ -135,4 +136,146 @@ IOResult read_mobility_plain(const std::string& filename) return success(mobility); } +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 { + edge_indx++; + } + } + return success(); +} + } // namespace mio diff --git a/cpp/memilio/io/mobility_io.h b/cpp/memilio/io/mobility_io.h index b6e840d0b9..e0e356cc59 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 result Simulation results per edge of the graph. + * @param ids Identifiers for the start and end node of the edges. + * @param filename Name of file + * @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); + +/** + * Saves the results of a simulation for each edge in the graph. + * @param ensemble_edges Simulation results for each run for each edge. + * @param pairs_edges Identifiers for the start and end node of the edges. + * @param result_dir Top level directory for all results of the parameter study. + * @param save_single_runs [Default: true] Defines if single run results are written. + * @param 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 a046be3fad..4a87217419 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -90,148 +90,6 @@ IOResult save_result(const std::vector>& results, const return success(); } -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 { - edge_indx++; - } - } - return success(); -} - herr_t store_group_name(hid_t /*id*/, const char* name, const H5L_info_t* /*linfo*/, void* opdata) { auto h5group_names = reinterpret_cast*>(opdata); diff --git a/cpp/memilio/io/result_io.h b/cpp/memilio/io/result_io.h index 9c8557eccb..0746a691f3 100644 --- a/cpp/memilio/io/result_io.h +++ b/cpp/memilio/io/result_io.h @@ -121,29 +121,6 @@ IOResult save_result_with_params(const std::vector>& re return success(); } -/** - * @brief Save the results of the edges for a single graph simulation run. - * @param result Simulation results per edge of the graph. - * @param ids Identifiers for the start and end node of the edges. - * @param filename Name of file - * @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); - -/** - * Saves the results of a simulation for each edge in the graph. - * @param ensemble_edges Simulation results for each run for each edge. - * @param pairs_edges Identifiers for the start and end node of the edges. - * @param result_dir Top level directory for all results of the parameter study. - * @param save_single_runs [Default: true] Defines if single run results are written. - * @param 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); - /** * @brief Save the results of a parameter study. * Stores different percentiles (p5, p25, p50, p75, p90) and sums of the results and parameters. diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 4f65df87d7..cf42526cb3 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -129,8 +129,7 @@ class MobilityParameters */ MobilityParameters(const MobilityCoefficientGroup& coeffs) : m_coefficients(coeffs) - , m_save_indices(0) - + , m_saved_compartment_indices(0) { } @@ -140,31 +139,38 @@ class MobilityParameters */ MobilityParameters(const Eigen::VectorXd& coeffs) : m_coefficients({MobilityCoefficients(coeffs)}) - , m_save_indices(0) + , m_saved_compartment_indices(0) { } /** - * @brief Constructor for a new MobilityParameters object. - * - * @param coeffs Mobility coefficients - * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. - */ + * @brief Constructor for initializing mobility parameters with coefficients from type MobilityCoefficientGroup and + * specific save indices. + * + * @param coeffs A group of mobility coefficients represented by a `MobilityCoefficientGroup` object, defining how + * individuals move between nodes. + * @param save_indices A 2D vector of indices. 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_save_indices(save_indices) + , m_saved_compartment_indices(save_indices) { } /** - * @brief Constructor for a new MobilityParameters object. - * - * @param coeffs Mobility coefficients - * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. - */ + * @brief Constructor for initializing mobility parameters with coefficients from an Eigen Vector + * and specific save indices. + * + * @param coeffs A Eigen::VectorXd containing mobility coefficients. + * @param save_indices A 2D vector of indices. 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_save_indices(save_indices) + , m_saved_compartment_indices(save_indices) { } @@ -208,11 +214,18 @@ class MobilityParameters } /** - * @return the indices we want to save for the edge - */ + * @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 specify 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. + * Each inner vector represents a group of compartments defined by indices. + */ const auto& get_save_indices() const { - return m_save_indices; + return m_saved_compartment_indices; } /** @@ -275,7 +288,7 @@ class MobilityParameters private: MobilityCoefficientGroup m_coefficients; //one per group and compartment DynamicNPIs m_dynamic_npis; - std::vector> m_save_indices; // groups of indices from compartments to save + std::vector> m_saved_compartment_indices; // groups of indices from compartments to save }; /** @@ -294,8 +307,8 @@ class MobilityEdge , m_mobile_population(params.get_coefficients().get_shape().rows()) , m_return_times(0) , m_return_mobile_population(false) - , m_save_indices(params.get_save_indices()) - , m_mobility_results(m_save_indices.size() + 1) + , m_saved_compartment_indices(params.get_save_indices()) + , m_mobility_results(m_saved_compartment_indices.size() + 1) { } @@ -308,8 +321,8 @@ class MobilityEdge , m_mobile_population(coeffs.rows()) , m_return_times(0) , m_return_mobile_population(false) - , m_save_indices(0) - , m_mobility_results(m_save_indices.size() + 1) + , m_saved_compartment_indices(0) + , m_mobility_results(m_saved_compartment_indices.size() + 1) { } @@ -323,8 +336,8 @@ class MobilityEdge , m_mobile_population(params.get_coefficients().get_shape().rows()) , m_return_times(0) , m_return_mobile_population(false) - , m_save_indices(save_indices) - , m_mobility_results(m_save_indices.size() + 1) + , m_saved_compartment_indices(save_indices) + , m_mobility_results(m_saved_compartment_indices.size() + 1) { } @@ -338,8 +351,8 @@ class MobilityEdge , m_mobile_population(coeffs.rows()) , m_return_times(0) , m_return_mobile_population(false) - , m_save_indices(save_indices) - , m_mobility_results(m_save_indices.size() + 1) + , m_saved_compartment_indices(save_indices) + , m_mobility_results(m_saved_compartment_indices.size() + 1) { } @@ -384,7 +397,7 @@ 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_save_indices; // groups of indices from compartments to save + 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 /** @@ -393,21 +406,21 @@ class MobilityEdge * Additionally, the total number of commuters is stored in the last entry of m_mobility_results. * @param[in] t current time */ - void condense_m_mobility(const double t); + void add_mobility_result_time_point(const double t); }; template -void MobilityEdge::condense_m_mobility(const double t) +void MobilityEdge::add_mobility_result_time_point(const double t) { - const size_t save_indices_size = this->m_save_indices.size(); + 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_save_indices for each group (e.g. Age groups) - std::transform(this->m_save_indices.begin(), this->m_save_indices.end(), condensed_values.data(), - [&last_value](const auto& indices) { + // 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]; @@ -590,7 +603,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]; - condense_m_mobility(t); + add_mobility_result_time_point(t); m_mobile_population.remove_time_point(i); m_return_times.remove_time_point(i); } @@ -609,7 +622,7 @@ 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(); - condense_m_mobility(t); + add_mobility_result_time_point(t); } m_return_mobile_population = !m_return_mobile_population; } diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp index 0d7b9b8a53..8accf2fc7e 100644 --- a/cpp/tests/test_mobility.cpp +++ b/cpp/tests/test_mobility.cpp @@ -167,7 +167,7 @@ TEST(TestMobility, edgeApplyMobility) EXPECT_DOUBLE_EQ(node2.get_result().get_last_value().sum(), 1100); } -TEST(TestMobility, condense_m_mobility) +TEST(TestMobility, add_mobility_result_time_point) { using Model = mio::osecir::Model; From eb684df5ee6afa0066744e24b2955ab03f5951f3 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Wed, 28 Aug 2024 16:48:44 +0200 Subject: [PATCH 28/30] if no hdf5 --- cpp/memilio/io/mobility_io.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cpp/memilio/io/mobility_io.cpp b/cpp/memilio/io/mobility_io.cpp index 5ad86a6558..832ac7c63e 100644 --- a/cpp/memilio/io/mobility_io.cpp +++ b/cpp/memilio/io/mobility_io.cpp @@ -1,7 +1,10 @@ #include "memilio/io/mobility_io.h" #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 @@ -136,6 +139,7 @@ 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) @@ -277,5 +281,6 @@ IOResult save_edges(const std::vector>& results, const } return success(); } +#endif // MEMILIO_HAS_HDF5 } // namespace mio From bfe52c1500b83252147574f194aa2b898660894b Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 30 Sep 2024 10:59:28 +0200 Subject: [PATCH 29/30] doc, test --- .../metapopulation_mobility_instant.h | 85 +++++++++++-------- cpp/tests/test_graph.cpp | 16 ++-- cpp/tests/test_save_results.cpp | 36 +++++++- 3 files changed, 93 insertions(+), 44 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index cf42526cb3..3fbbdbfee3 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -144,14 +144,14 @@ class MobilityParameters } /** - * @brief Constructor for initializing mobility parameters with coefficients from type MobilityCoefficientGroup and - * specific save indices. + * @brief Constructor for initializing mobility parameters with coefficients from type `MobilityCoefficientGroup` + * and specific save indices. * - * @param coeffs A group of mobility coefficients represented by a `MobilityCoefficientGroup` object, defining how - * individuals move between nodes. - * @param save_indices A 2D vector of indices. 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. + * @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) @@ -163,10 +163,10 @@ class MobilityParameters * @brief Constructor for initializing mobility parameters with coefficients from an Eigen Vector * and specific save indices. * - * @param coeffs A Eigen::VectorXd containing mobility coefficients. - * @param save_indices A 2D vector of indices. 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. + * @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)}) @@ -214,15 +214,15 @@ class MobilityParameters } /** - * @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 specify 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. - * Each inner vector represents a group of compartments defined by indices. - */ + * @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. Each inner vector + * represents a group of compartments defined by indices. + */ const auto& get_save_indices() const { return m_saved_compartment_indices; @@ -299,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) @@ -313,8 +313,10 @@ class MobilityEdge } /** - * 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) @@ -327,9 +329,12 @@ class MobilityEdge } /** - * create edge with coefficients as MobilityParameters object and a 2d vector of indices which determine which compartments we save. - * @param coeffs % of people in each group and compartment that migrate in each time step. - * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. + * @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) @@ -342,9 +347,12 @@ class MobilityEdge } /** - * create edge with coefficients and a 2d vector of indices which determine which compartments we save. - * @param coeffs % of people in each group and compartment that migrate in each time step. - * @param save_indices 2D vector of indices. Each inner vector represents a group of indices to be saved. + * @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) @@ -365,9 +373,11 @@ class MobilityEdge } /** - * Retrieve the count of commuters in selected infection states, - * along with the total number of commuter. - */ + * @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; @@ -376,6 +386,7 @@ class MobilityEdge { return m_mobility_results; } + /** @} */ /** * compute mobility from node_from to node_to. @@ -401,10 +412,12 @@ class MobilityEdge TimeSeries m_mobility_results; // save results from edges + entry for the total number of commuters /** - * Computes a condensed version of m_mobile_population and puts it in m_mobility_results. - * 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 current time + * @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); }; diff --git a/cpp/tests/test_graph.cpp b/cpp/tests/test_graph.cpp index 7ee8df12f9..2ba785ffea 100644 --- a/cpp/tests/test_graph.cpp +++ b/cpp/tests/test_graph.cpp @@ -32,6 +32,7 @@ #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 @@ -279,11 +280,12 @@ TEST(TestGraph, set_edges_saving_edges) mio::osecir::Model model(num_groups); model.populations[{mio::AgeGroup(3), mio::osecir::InfectionState::Exposed}] = 1; mio::Graph, mio::MobilityParameters> params_graph; - const fs::path& dir = " "; + 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}; + 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); @@ -315,7 +317,7 @@ TEST(TestGraph, set_edges_saving_edges) 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); + std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}, indices_save_edges); EXPECT_EQ(params_graph.edges().size(), 2); @@ -350,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_save_results.cpp b/cpp/tests/test_save_results.cpp index bebb5fac2e..c8c16309c0 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -467,4 +467,38 @@ TEST(TestSaveResult, save_edges) 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)); -} \ No newline at end of file +} + +TEST(TestSaveEdges, save_edges_empty_ts) +{ + 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}, {0, 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 + auto save_result = save_edges(results, pairs_edges, results_file_path); + EXPECT_EQ(save_result, mio::success()); + + // Read the results back in and check that they are correct for the first TimeSeries. + auto results_from_file = mio::read_result(results_file_path); + ASSERT_TRUE(results_from_file); + + auto result_from_file_group0 = results_from_file.value()[0]; + EXPECT_EQ(result_from_file_group0.get_groups().get_num_time_points(), 2); +} From 97b44fd62c6c1068ec155750f92c4b89a8abd7f8 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Fri, 4 Oct 2024 14:02:05 +0200 Subject: [PATCH 30/30] failure if ts empty, doc for some functions --- cpp/memilio/io/mobility_io.cpp | 4 +++- cpp/memilio/io/mobility_io.h | 18 +++++++++--------- .../mobility/metapopulation_mobility_instant.h | 4 ++-- cpp/tests/test_save_results.cpp | 15 ++++----------- 4 files changed, 18 insertions(+), 23 deletions(-) diff --git a/cpp/memilio/io/mobility_io.cpp b/cpp/memilio/io/mobility_io.cpp index 832ac7c63e..171f607796 100644 --- a/cpp/memilio/io/mobility_io.cpp +++ b/cpp/memilio/io/mobility_io.cpp @@ -276,7 +276,9 @@ IOResult save_edges(const std::vector>& results, const } } else { - edge_indx++; + 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(); diff --git a/cpp/memilio/io/mobility_io.h b/cpp/memilio/io/mobility_io.h index e0e356cc59..e55fe7ae76 100644 --- a/cpp/memilio/io/mobility_io.h +++ b/cpp/memilio/io/mobility_io.h @@ -189,21 +189,21 @@ IOResult>> read_graph(const std::string& dir #ifdef MEMILIO_HAS_HDF5 /** * @brief Save the results of the edges for a single graph simulation run. - * @param result Simulation results per edge of the graph. - * @param ids Identifiers for the start and end node of the edges. - * @param filename Name of file + * @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); /** - * Saves the results of a simulation for each edge in the graph. - * @param ensemble_edges Simulation results for each run for each edge. - * @param pairs_edges Identifiers for the start and end node of the edges. - * @param result_dir Top level directory for all results of the parameter study. - * @param save_single_runs [Default: true] Defines if single run results are written. - * @param save_percentiles [Default: true] Defines if percentiles are written. + * @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, diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 3fbbdbfee3..34996710d9 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -220,8 +220,8 @@ class MobilityParameters * 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. Each inner vector - * represents a group of compartments defined by indices. + * @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 { diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index c8c16309c0..5710c3c879 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -471,6 +471,7 @@ TEST(TestSaveResult, save_edges) TEST(TestSaveEdges, save_edges_empty_ts) { + mio::set_log_level(mio::LogLevel::off); std::vector> results; const auto num_elements = 2; @@ -485,20 +486,12 @@ TEST(TestSaveEdges, save_edges_empty_ts) mio::TimeSeries ts_2(num_elements); results.push_back(ts_2); - const std::vector> pairs_edges = {{0, 1}, {0, 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 - auto save_result = save_edges(results, pairs_edges, results_file_path); - EXPECT_EQ(save_result, mio::success()); - - // Read the results back in and check that they are correct for the first TimeSeries. - auto results_from_file = mio::read_result(results_file_path); - ASSERT_TRUE(results_from_file); - - auto result_from_file_group0 = results_from_file.value()[0]; - EXPECT_EQ(result_from_file_group0.get_groups().get_num_time_points(), 2); + // 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)); }