From 655781253d2ad41aa9769754e4fa93d16529dc67 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 15 Sep 2025 12:27:38 +0000 Subject: [PATCH 1/4] Move from source to header --- cpp/memilio/io/mobility_io.cpp | 255 ---------------- cpp/memilio/io/mobility_io.h | 255 +++++++++++++++- cpp/memilio/io/parameters_io.cpp | 49 +--- cpp/memilio/io/parameters_io.h | 39 ++- cpp/memilio/io/result_io.cpp | 59 ---- cpp/memilio/io/result_io.h | 61 +++- cpp/models/ode_secir/parameters_io.cpp | 152 ---------- cpp/models/ode_secir/parameters_io.h | 164 ++++++++++- cpp/models/ode_secirvvs/parameters_io.cpp | 234 --------------- cpp/models/ode_secirvvs/parameters_io.h | 338 ++++++++++++++++++---- 10 files changed, 776 insertions(+), 830 deletions(-) diff --git a/cpp/memilio/io/mobility_io.cpp b/cpp/memilio/io/mobility_io.cpp index cdcffa762a..f73821e710 100644 --- a/cpp/memilio/io/mobility_io.cpp +++ b/cpp/memilio/io/mobility_io.cpp @@ -18,17 +18,6 @@ * limitations under the License. */ #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 -#include -#include namespace mio { @@ -60,248 +49,4 @@ IOResult count_lines(const std::string& filename) return success(count); } -IOResult read_mobility_formatted(const std::string& filename) -{ - BOOST_OUTCOME_TRY(auto&& num_lines, count_lines(filename)); - - if (num_lines == 0) { - return success(Eigen::MatrixXd(0, 0)); - } - - std::fstream file; - file.open(filename, std::ios::in); - if (!file.is_open()) { - return failure(StatusCode::FileNotFound, filename); - } - - Eigen::MatrixXd txt_matrix(num_lines - 1, 3); - std::vector ids; - - try { - std::string tp; - int linenumber = 0; - while (getline(file, tp)) { - if (linenumber > 0) { - auto line = split(tp, '\t'); - if (line.size() < 5) { - return failure(StatusCode::InvalidFileFormat, - filename + ":" + std::to_string(linenumber) + ": Not enough entries in line."); - } - ids.push_back(std::stoi(line[2])); - txt_matrix(linenumber - 1, 0) = std::stoi(line[2]); - txt_matrix(linenumber - 1, 1) = std::stoi(line[3]); - txt_matrix(linenumber - 1, 2) = std::stod(line[4]); - } - linenumber++; - } - } - catch (std::runtime_error& ex) { - return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what()); - } - - sort(ids.begin(), ids.end()); - std::vector::iterator iter = std::unique(ids.begin(), ids.end()); - ids.resize(std::distance(ids.begin(), iter)); - - Eigen::MatrixXd mobility = Eigen::MatrixXd::Zero(ids.size(), ids.size()); - - for (int k = 0; k < num_lines - 1; k++) { - int row_ind = 0; - int col_ind = 0; - while (txt_matrix(k, 0) != ids[row_ind]) { - row_ind++; - } - while (txt_matrix(k, 1) != ids[col_ind]) { - col_ind++; - } - mobility(row_ind, col_ind) = txt_matrix(k, 2); - } - - return success(mobility); -} - -IOResult read_mobility_plain(const std::string& filename) -{ - BOOST_OUTCOME_TRY(auto&& num_lines, count_lines(filename)); - - if (num_lines == 0) { - return success(Eigen::MatrixXd(0, 0)); - } - - std::fstream file; - file.open(filename, std::ios::in); - if (!file.is_open()) { - return failure(StatusCode::FileNotFound, filename); - } - - Eigen::MatrixXd mobility(num_lines, num_lines); - - try { - std::string tp; - int linenumber = 0; - while (getline(file, tp)) { - auto line = split(tp, ' '); - if (line.size() != size_t(num_lines)) { - return failure(StatusCode::InvalidFileFormat, filename + ": Not a square matrix."); - } - Eigen::Index i = static_cast(linenumber); - for (Eigen::Index j = 0; j < static_cast(line.size()); j++) { - mobility(i, j) = std::stod(line[j]); - } - linenumber++; - } - } - catch (std::runtime_error& ex) { - return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what()); - } - - return success(mobility); -} - -#ifdef MEMILIO_HAS_HDF5 -IOResult save_edges(const std::vector>>& ensemble_edges, - const std::vector>& pairs_edges, const fs::path& result_dir, - bool save_single_runs, bool save_percentiles) -{ - //save results and sum of results over nodes - auto ensemble_edges_sum = sum_nodes(ensemble_edges); - if (save_single_runs) { - for (size_t i = 0; i < ensemble_edges_sum.size(); ++i) { - BOOST_OUTCOME_TRY(save_edges(ensemble_edges[i], pairs_edges, - (result_dir / ("Edges_run" + std::to_string(i) + ".h5")).string())); - } - } - - if (save_percentiles) { - // make directories for percentiles - auto result_dir_p05 = result_dir / "p05"; - auto result_dir_p25 = result_dir / "p25"; - auto result_dir_p50 = result_dir / "p50"; - auto result_dir_p75 = result_dir / "p75"; - auto result_dir_p95 = result_dir / "p95"; - BOOST_OUTCOME_TRY(create_directory(result_dir_p05.string())); - BOOST_OUTCOME_TRY(create_directory(result_dir_p25.string())); - BOOST_OUTCOME_TRY(create_directory(result_dir_p50.string())); - BOOST_OUTCOME_TRY(create_directory(result_dir_p75.string())); - BOOST_OUTCOME_TRY(create_directory(result_dir_p95.string())); - - // save percentiles of Edges - { - auto ensemble_edges_p05 = ensemble_percentile(ensemble_edges, 0.05); - auto ensemble_edges_p25 = ensemble_percentile(ensemble_edges, 0.25); - auto ensemble_edges_p50 = ensemble_percentile(ensemble_edges, 0.50); - auto ensemble_edges_p75 = ensemble_percentile(ensemble_edges, 0.75); - auto ensemble_edges_p95 = ensemble_percentile(ensemble_edges, 0.95); - - BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p05, pairs_edges, (result_dir_p05 / "Edges.h5").string())); - BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p25, pairs_edges, (result_dir_p25 / "Edges.h5").string())); - BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p50, pairs_edges, (result_dir_p50 / "Edges.h5").string())); - BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p75, pairs_edges, (result_dir_p75 / "Edges.h5").string())); - BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p95, pairs_edges, (result_dir_p95 / "Edges.h5").string())); - } - } - return success(); -} - -IOResult save_edges(const std::vector>& results, - const std::vector>& ids, const std::string& filename) -{ - const auto num_edges = results.size(); - size_t edge_indx = 0; - // H5Fcreate creates a new HDF5 file. - // H5F_ACC_TRUNC: If the file already exists, H5Fcreate fails. If the file does not exist, it is created and opened with read-write access. - // H5P_DEFAULT: default data transfer properties are used. - H5File file{H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)}; - MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, "Failed to open the HDF5 file: " + filename); - - while (edge_indx < num_edges) { - const auto& result = results[edge_indx]; - auto h5group_name = "/" + std::to_string(ids[edge_indx].first); - H5Group start_node_h5group{H5Gcreate(file.id, h5group_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; - MEMILIO_H5_CHECK(start_node_h5group.id, StatusCode::UnknownError, - "Group could not be created (" + h5group_name + ") in the file: " + filename); - - const auto num_timepoints = result.get_num_time_points(); - if (num_timepoints > 0) { - const auto num_elements = result.get_value(0).size(); - - hsize_t dims_t[] = {static_cast(num_timepoints)}; - H5DataSpace dspace_t{H5Screate_simple(1, dims_t, NULL)}; - MEMILIO_H5_CHECK(dspace_t.id, StatusCode::UnknownError, - "Failed to create the DataSpace for 'Time' in group " + h5group_name + - " in the file: " + filename); - - H5DataSet dset_t{H5Dcreate(start_node_h5group.id, "Time", H5T_NATIVE_DOUBLE, dspace_t.id, H5P_DEFAULT, - H5P_DEFAULT, H5P_DEFAULT)}; - MEMILIO_H5_CHECK(dset_t.id, StatusCode::UnknownError, - "Failed to create the 'Time' DataSet in group " + h5group_name + - " in the file: " + filename); - - auto values_t = std::vector(result.get_times().begin(), result.get_times().end()); - MEMILIO_H5_CHECK(H5Dwrite(dset_t.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, values_t.data()), - StatusCode::UnknownError, - "Failed to write 'Time' data in group " + h5group_name + " in the file: " + filename); - - int start_id = ids[edge_indx].first; - auto total = Eigen::Matrix::Zero( - num_timepoints, num_elements) - .eval(); - while (edge_indx < num_edges && ids[edge_indx].first == start_id) { - const auto& result_edge = results[edge_indx]; - auto edge_result = Eigen::Matrix::Zero( - num_timepoints, num_elements) - .eval(); - for (Eigen::Index t_idx = 0; t_idx < result_edge.get_num_time_points(); ++t_idx) { - auto v = result_edge.get_value(t_idx).transpose().eval(); - edge_result.row(t_idx) = v; - total.row(t_idx) += v; - } - - hsize_t dims_values[] = {static_cast(num_timepoints), static_cast(num_elements)}; - H5DataSpace dspace_values{H5Screate_simple(2, dims_values, NULL)}; - MEMILIO_H5_CHECK(dspace_values.id, StatusCode::UnknownError, - "Failed to create the DataSpace for End" + std::to_string(ids[edge_indx].second) + - " in group " + h5group_name + " in the file: " + filename); - - // End is the target node of the edge - auto dset_name = "End" + std::to_string(ids[edge_indx].second); - H5DataSet dset_values{H5Dcreate(start_node_h5group.id, dset_name.c_str(), H5T_NATIVE_DOUBLE, - dspace_values.id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; - MEMILIO_H5_CHECK(dset_values.id, StatusCode::UnknownError, - "Failed to create the DataSet for " + dset_name + " in group " + h5group_name + - " in the file: " + filename); - - MEMILIO_H5_CHECK( - H5Dwrite(dset_values.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, edge_result.data()), - StatusCode::UnknownError, - "Failed to write data for " + dset_name + " in group " + h5group_name + - " in the file: " + filename); - - // In the final iteration, we also save the total values - if (edge_indx == num_edges - 1 || ids[edge_indx + 1].first != start_id) { - dset_name = "Total"; - H5DataSet dset_total{H5Dcreate(start_node_h5group.id, dset_name.c_str(), H5T_NATIVE_DOUBLE, - dspace_values.id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; - MEMILIO_H5_CHECK(dset_total.id, StatusCode::UnknownError, - "Failed to create the Total DataSet in group " + h5group_name + - " in the file: " + filename); - - MEMILIO_H5_CHECK( - H5Dwrite(dset_total.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, total.data()), - StatusCode::UnknownError, - "Failed to write Total data in group " + h5group_name + " in the file: " + filename); - } - edge_indx++; - } - } - else { - log_error("No time points in the TimeSeries for Edge combination " + std::to_string(ids[edge_indx].first) + - " -> " + std::to_string(ids[edge_indx].second)); - return failure(mio::StatusCode::InvalidValue); - } - } - return success(); -} -#endif // MEMILIO_HAS_HDF5 - } // namespace mio diff --git a/cpp/memilio/io/mobility_io.h b/cpp/memilio/io/mobility_io.h index 20ed71f3e9..9e62ce5075 100644 --- a/cpp/memilio/io/mobility_io.h +++ b/cpp/memilio/io/mobility_io.h @@ -24,6 +24,18 @@ #include "memilio/mobility/graph.h" #include "memilio/data/analyze_result.h" #include "memilio/mobility/metapopulation_mobility_instant.h" +#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 +#include +#include namespace mio { @@ -48,7 +60,66 @@ IOResult count_lines(const std::string& filename); * where N is the number of regions * @param filename name of file to be read */ -IOResult read_mobility_formatted(const std::string& filename); +template +IOResult read_mobility_formatted(const std::string& filename) +{ + BOOST_OUTCOME_TRY(auto&& num_lines, count_lines(filename)); + + if (num_lines == 0) { + return success(Eigen::MatrixXd(0, 0)); + } + + std::fstream file; + file.open(filename, std::ios::in); + if (!file.is_open()) { + return failure(StatusCode::FileNotFound, filename); + } + + Eigen::MatrixXd txt_matrix(num_lines - 1, 3); + std::vector ids; + + try { + std::string tp; + int linenumber = 0; + while (getline(file, tp)) { + if (linenumber > 0) { + auto line = split(tp, '\t'); + if (line.size() < 5) { + return failure(StatusCode::InvalidFileFormat, + filename + ":" + std::to_string(linenumber) + ": Not enough entries in line."); + } + ids.push_back(std::stoi(line[2])); + txt_matrix(linenumber - 1, 0) = std::stoi(line[2]); + txt_matrix(linenumber - 1, 1) = std::stoi(line[3]); + txt_matrix(linenumber - 1, 2) = std::stod(line[4]); + } + linenumber++; + } + } + catch (std::runtime_error& ex) { + return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what()); + } + + sort(ids.begin(), ids.end()); + std::vector::iterator iter = std::unique(ids.begin(), ids.end()); + ids.resize(std::distance(ids.begin(), iter)); + + Eigen::MatrixXd mobility = Eigen::MatrixXd::Zero(ids.size(), ids.size()); + + for (int k = 0; k < num_lines - 1; k++) { + int row_ind = 0; + int col_ind = 0; + while (txt_matrix(k, 0) != ids[row_ind]) { + row_ind++; + } + while (txt_matrix(k, 1) != ids[col_ind]) { + col_ind++; + } + mobility(row_ind, col_ind) = txt_matrix(k, 2); + } + + return success(mobility); +} /** * @brief Reads txt mobility data or contact which is given by values only @@ -56,7 +127,44 @@ IOResult read_mobility_formatted(const std::string& filename); * Matrix, where N is the number of regions * @param filename name of file to be read */ -IOResult read_mobility_plain(const std::string& filename); +template +IOResult read_mobility_plain(const std::string& filename) +{ + BOOST_OUTCOME_TRY(auto&& num_lines, count_lines(filename)); + + if (num_lines == 0) { + return success(Eigen::MatrixXd(0, 0)); + } + + std::fstream file; + file.open(filename, std::ios::in); + if (!file.is_open()) { + return failure(StatusCode::FileNotFound, filename); + } + + Eigen::MatrixXd mobility(num_lines, num_lines); + + try { + std::string tp; + int linenumber = 0; + while (getline(file, tp)) { + auto line = split(tp, ' '); + if (line.size() != size_t(num_lines)) { + return failure(StatusCode::InvalidFileFormat, filename + ": Not a square matrix."); + } + Eigen::Index i = static_cast(linenumber); + for (Eigen::Index j = 0; j < static_cast(line.size()); j++) { + mobility(i, j) = std::stod(line[j]); + } + linenumber++; + } + } + catch (std::runtime_error& ex) { + return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what()); + } + + return success(mobility); +} #ifdef MEMILIO_HAS_JSONCPP @@ -194,8 +302,106 @@ IOResult>> read_graph(const std::string& dir * @param[in] filename Name of the file where the results will be saved. * @return Any io errors that occur during writing of the files. */ +template IOResult save_edges(const std::vector>& results, - const std::vector>& ids, const std::string& filename); + const std::vector>& ids, const std::string& filename) +{ + const auto num_edges = results.size(); + size_t edge_indx = 0; + // H5Fcreate creates a new HDF5 file. + // H5F_ACC_TRUNC: If the file already exists, H5Fcreate fails. If the file does not exist, it is created and opened with read-write access. + // H5P_DEFAULT: default data transfer properties are used. + H5File file{H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, "Failed to open the HDF5 file: " + filename); + + while (edge_indx < num_edges) { + const auto& result = results[edge_indx]; + auto h5group_name = "/" + std::to_string(ids[edge_indx].first); + H5Group start_node_h5group{H5Gcreate(file.id, h5group_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(start_node_h5group.id, StatusCode::UnknownError, + "Group could not be created (" + h5group_name + ") in the file: " + filename); + + const auto num_timepoints = result.get_num_time_points(); + if (num_timepoints > 0) { + const auto num_elements = result.get_value(0).size(); + + hsize_t dims_t[] = {static_cast(num_timepoints)}; + H5DataSpace dspace_t{H5Screate_simple(1, dims_t, NULL)}; + MEMILIO_H5_CHECK(dspace_t.id, StatusCode::UnknownError, + "Failed to create the DataSpace for 'Time' in group " + h5group_name + + " in the file: " + filename); + + H5DataSet dset_t{H5Dcreate(start_node_h5group.id, "Time", H5T_NATIVE_DOUBLE, dspace_t.id, H5P_DEFAULT, + H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(dset_t.id, StatusCode::UnknownError, + "Failed to create the 'Time' DataSet in group " + h5group_name + + " in the file: " + filename); + + auto values_t = std::vector(result.get_times().begin(), result.get_times().end()); + MEMILIO_H5_CHECK(H5Dwrite(dset_t.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, values_t.data()), + StatusCode::UnknownError, + "Failed to write 'Time' data in group " + h5group_name + " in the file: " + filename); + + int start_id = ids[edge_indx].first; + auto total = Eigen::Matrix::Zero( + num_timepoints, num_elements) + .eval(); + while (edge_indx < num_edges && ids[edge_indx].first == start_id) { + const auto& result_edge = results[edge_indx]; + auto edge_result = Eigen::Matrix::Zero( + num_timepoints, num_elements) + .eval(); + for (Eigen::Index t_idx = 0; t_idx < result_edge.get_num_time_points(); ++t_idx) { + auto v = result_edge.get_value(t_idx).transpose().eval(); + edge_result.row(t_idx) = v; + total.row(t_idx) += v; + } + + hsize_t dims_values[] = {static_cast(num_timepoints), static_cast(num_elements)}; + H5DataSpace dspace_values{H5Screate_simple(2, dims_values, NULL)}; + MEMILIO_H5_CHECK(dspace_values.id, StatusCode::UnknownError, + "Failed to create the DataSpace for End" + std::to_string(ids[edge_indx].second) + + " in group " + h5group_name + " in the file: " + filename); + + // End is the target node of the edge + auto dset_name = "End" + std::to_string(ids[edge_indx].second); + H5DataSet dset_values{H5Dcreate(start_node_h5group.id, dset_name.c_str(), H5T_NATIVE_DOUBLE, + dspace_values.id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(dset_values.id, StatusCode::UnknownError, + "Failed to create the DataSet for " + dset_name + " in group " + h5group_name + + " in the file: " + filename); + + MEMILIO_H5_CHECK( + H5Dwrite(dset_values.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, edge_result.data()), + StatusCode::UnknownError, + "Failed to write data for " + dset_name + " in group " + h5group_name + + " in the file: " + filename); + + // In the final iteration, we also save the total values + if (edge_indx == num_edges - 1 || ids[edge_indx + 1].first != start_id) { + dset_name = "Total"; + H5DataSet dset_total{H5Dcreate(start_node_h5group.id, dset_name.c_str(), H5T_NATIVE_DOUBLE, + dspace_values.id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(dset_total.id, StatusCode::UnknownError, + "Failed to create the Total DataSet in group " + h5group_name + + " in the file: " + filename); + + MEMILIO_H5_CHECK( + H5Dwrite(dset_total.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, total.data()), + StatusCode::UnknownError, + "Failed to write Total data in group " + h5group_name + " in the file: " + filename); + } + edge_indx++; + } + } + else { + log_error("No time points in the TimeSeries for Edge combination " + std::to_string(ids[edge_indx].first) + + " -> " + std::to_string(ids[edge_indx].second)); + return failure(mio::StatusCode::InvalidValue); + } + } + return success(); +} /** * @brief Saves the results of a simulation for each edge in the graph. @@ -206,9 +412,50 @@ IOResult save_edges(const std::vector>& results, * @param[in] save_percentiles [Default: true] Defines if percentiles are written. * @return Any io errors that occur during writing of the files. */ +template 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); + bool save_single_runs = true, bool save_percentiles = true) +{ + //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(); +} #endif //MEMILIO_HAS_HDF5 diff --git a/cpp/memilio/io/parameters_io.cpp b/cpp/memilio/io/parameters_io.cpp index dff30087e7..019d10f3e1 100644 --- a/cpp/memilio/io/parameters_io.cpp +++ b/cpp/memilio/io/parameters_io.cpp @@ -18,55 +18,8 @@ * limitations under the License. */ -#include "memilio/config.h" - -#ifdef MEMILIO_HAS_JSONCPP - -#include "memilio/io/epi_data.h" -#include "memilio/io/result_io.h" -#include "json/value.h" -#include -#include +#include "memilio/io/parameters_io.h" namespace mio { -IOResult>> -read_population_data(const std::vector& population_data, const std::vector& vregion) -{ - std::vector> vnum_population( - vregion.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - - for (auto&& county_entry : population_data) { - //accumulate population of states or country from population of counties - if (!county_entry.county_id && !county_entry.district_id) { - return failure(StatusCode::InvalidFileFormat, "File with county population expected."); - } - //find region that this county belongs to - //all counties belong to the country (id = 0) - auto it = std::find_if(vregion.begin(), vregion.end(), [&county_entry](auto r) { - return r == 0 || - (county_entry.county_id && - regions::StateId(r) == regions::get_state_id(int(*county_entry.county_id))) || - (county_entry.county_id && regions::CountyId(r) == *county_entry.county_id) || - (county_entry.district_id && regions::DistrictId(r) == *county_entry.district_id); - }); - if (it != vregion.end()) { - auto region_idx = size_t(it - vregion.begin()); - auto& num_population = vnum_population[region_idx]; - for (size_t age = 0; age < num_population.size(); age++) { - num_population[age] += county_entry.population[AgeGroup(age)]; - } - } - } - - return success(vnum_population); -} - -IOResult>> read_population_data(const std::string& path, - const std::vector& vregion) -{ - BOOST_OUTCOME_TRY(auto&& population_data, mio::read_population_data(path)); - return read_population_data(population_data, vregion); -} } // namespace mio -#endif //MEMILIO_HAS_JSONCPP diff --git a/cpp/memilio/io/parameters_io.h b/cpp/memilio/io/parameters_io.h index 9dee5046f2..23c23b729a 100644 --- a/cpp/memilio/io/parameters_io.h +++ b/cpp/memilio/io/parameters_io.h @@ -122,8 +122,38 @@ IOResult read_divi_data(const std::string& path, const std::vector& v * @return An IOResult containing a vector of vectors, where each inner vector represents the population * distribution across age groups for a specific region, or an error if the function fails. */ +template IOResult>> -read_population_data(const std::vector& population_data, const std::vector& vregion); +read_population_data(const std::vector& population_data, const std::vector& vregion) +{ + std::vector> vnum_population( + vregion.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); + + for (auto&& county_entry : population_data) { + //accumulate population of states or country from population of counties + if (!county_entry.county_id && !county_entry.district_id) { + return failure(StatusCode::InvalidFileFormat, "File with county population expected."); + } + //find region that this county belongs to + //all counties belong to the country (id = 0) + auto it = std::find_if(vregion.begin(), vregion.end(), [&county_entry](auto r) { + return r == 0 || + (county_entry.county_id && + regions::StateId(r) == regions::get_state_id(int(*county_entry.county_id))) || + (county_entry.county_id && regions::CountyId(r) == *county_entry.county_id) || + (county_entry.district_id && regions::DistrictId(r) == *county_entry.district_id); + }); + if (it != vregion.end()) { + auto region_idx = size_t(it - vregion.begin()); + auto& num_population = vnum_population[region_idx]; + for (size_t age = 0; age < num_population.size(); age++) { + num_population[age] += county_entry.population[AgeGroup(age)]; + } + } + } + + return success(vnum_population); +} /** * @brief Reads population data from census data. @@ -133,8 +163,13 @@ read_population_data(const std::vector& population_data, co * @return An IOResult containing a vector of vectors, where each inner vector represents the population * distribution across age groups for a specific region, or an error if the function fails. */ +template IOResult>> read_population_data(const std::string& path, - const std::vector& vregion); + const std::vector& vregion) +{ + BOOST_OUTCOME_TRY(auto&& population_data, mio::read_population_data(path)); + return read_population_data(population_data, vregion); +} } // namespace mio diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 9fc3aa23b7..4c3f87971e 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -31,65 +31,6 @@ namespace mio { -IOResult save_result(const std::vector>& results, const std::vector& ids, - int num_groups, const std::string& filename) -{ - int region_idx = 0; - H5File file{H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)}; - MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, filename); - for (auto& result : results) { - auto h5group_name = "/" + std::to_string(ids[region_idx]); - H5Group region_h5group{H5Gcreate(file.id, h5group_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; - MEMILIO_H5_CHECK(region_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 = (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(region_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."); - - auto total = Eigen::Matrix::Zero( - num_timepoints, num_infectionstates) - .eval(); - - for (int group_idx = 0; group_idx <= num_groups; ++group_idx) { - auto group = Eigen::Matrix::Zero( - num_timepoints, num_infectionstates) - .eval(); - if (group_idx < num_groups) { - for (Eigen::Index t_idx = 0; t_idx < result.get_num_time_points(); ++t_idx) { - auto v = result[t_idx].transpose().eval(); - auto group_slice = mio::slice(v, {group_idx * num_infectionstates, num_infectionstates}); - mio::slice(group, {t_idx, 1}, {0, num_infectionstates}) = group_slice; - mio::slice(total, {t_idx, 1}, {0, num_infectionstates}) += group_slice; - } - } - - 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 = group_idx == num_groups ? std::string("Total") : "Group" + std::to_string(group_idx + 1); - H5DataSet dset_values{H5Dcreate(region_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, - group_idx == num_groups ? total.data() : group.data()), - StatusCode::UnknownError, "Values data could not be written."); - } - region_idx++; - } - 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 e325100e48..d8e2fd1cbc 100644 --- a/cpp/memilio/io/result_io.h +++ b/cpp/memilio/io/result_io.h @@ -43,8 +43,65 @@ namespace mio * @param filename Name of file * @return Any io errors that occur during writing of the files. */ -IOResult save_result(const std::vector>& result, const std::vector& ids, - int num_groups, const std::string& filename); +template +IOResult save_result(const std::vector>& results, const std::vector& ids, + int num_groups, const std::string& filename) +{ + int region_idx = 0; + H5File file{H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, filename); + for (auto& result : results) { + auto h5group_name = "/" + std::to_string(ids[region_idx]); + H5Group region_h5group{H5Gcreate(file.id, h5group_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(region_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 = (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(region_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."); + + auto total = Eigen::Matrix::Zero( + num_timepoints, num_infectionstates) + .eval(); + + for (int group_idx = 0; group_idx <= num_groups; ++group_idx) { + auto group = Eigen::Matrix::Zero( + num_timepoints, num_infectionstates) + .eval(); + if (group_idx < num_groups) { + for (Eigen::Index t_idx = 0; t_idx < result.get_num_time_points(); ++t_idx) { + auto v = result[t_idx].transpose().eval(); + auto group_slice = mio::slice(v, {group_idx * num_infectionstates, num_infectionstates}); + mio::slice(group, {t_idx, 1}, {0, num_infectionstates}) = group_slice; + mio::slice(total, {t_idx, 1}, {0, num_infectionstates}) += group_slice; + } + } + + 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 = group_idx == num_groups ? std::string("Total") : "Group" + std::to_string(group_idx + 1); + H5DataSet dset_values{H5Dcreate(region_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, + group_idx == num_groups ? total.data() : group.data()), + StatusCode::UnknownError, "Values data could not be written."); + } + region_idx++; + } + return success(); +} class SimulationResult { diff --git a/cpp/models/ode_secir/parameters_io.cpp b/cpp/models/ode_secir/parameters_io.cpp index 621efbf276..49ca676cf0 100644 --- a/cpp/models/ode_secir/parameters_io.cpp +++ b/cpp/models/ode_secir/parameters_io.cpp @@ -49,158 +49,6 @@ int get_region_id(int id) return id; } -IOResult read_confirmed_cases_data( - std::vector& rki_data, std::vector const& vregion, Date date, - std::vector>& vnum_Exposed, std::vector>& vnum_InfectedNoSymptoms, - std::vector>& vnum_InfectedSymptoms, std::vector>& vnum_InfectedSevere, - std::vector>& vnum_icu, std::vector>& vnum_death, - std::vector>& vnum_rec, const std::vector>& vt_Exposed, - const std::vector>& vt_InfectedNoSymptoms, - const std::vector>& vt_InfectedSymptoms, const std::vector>& vt_InfectedSevere, - const std::vector>& vt_InfectedCritical, const std::vector>& vmu_C_R, - const std::vector>& vmu_I_H, const std::vector>& vmu_H_U, - const std::vector& scaling_factor_inf) -{ - auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { - return a.date < b.date; - }); - if (max_date_entry == rki_data.end()) { - log_error("RKI data file is empty."); - return failure(StatusCode::InvalidFileFormat, "RKI file is empty."); - } - auto max_date = max_date_entry->date; - if (max_date < date) { - log_error("Specified date does not exist in RKI data"); - return failure(StatusCode::OutOfRange, "Specified date does not exist in RKI data."); - } - auto days_surplus = std::min(get_offset_in_days(max_date, date) - 6, 0); - - //this statement causes maybe-uninitialized warning for some versions of gcc. - //the error is reported in an included header, so the warning is disabled for the whole file - std::sort(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { - return std::make_tuple(get_region_id(a), a.date) < std::make_tuple(get_region_id(b), b.date); - }); - - for (auto region_idx = size_t(0); region_idx < vregion.size(); ++region_idx) { - auto region_entry_range_it = - std::equal_range(rki_data.begin(), rki_data.end(), vregion[region_idx], [](auto&& a, auto&& b) { - return get_region_id(a) < get_region_id(b); - }); - auto region_entry_range = make_range(region_entry_range_it); - if (region_entry_range.begin() == region_entry_range.end()) { - log_error("No entries found for region {}", vregion[region_idx]); - return failure(StatusCode::InvalidFileFormat, - "No entries found for region " + std::to_string(vregion[region_idx])); - } - for (auto&& region_entry : region_entry_range) { - - auto& t_Exposed = vt_Exposed[region_idx]; - auto& t_InfectedNoSymptoms = vt_InfectedNoSymptoms[region_idx]; - auto& t_InfectedSymptoms = vt_InfectedSymptoms[region_idx]; - auto& t_InfectedSevere = vt_InfectedSevere[region_idx]; - auto& t_InfectedCritical = vt_InfectedCritical[region_idx]; - - auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx]; - auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx]; - auto& num_rec = vnum_rec[region_idx]; - auto& num_Exposed = vnum_Exposed[region_idx]; - auto& num_InfectedSevere = vnum_InfectedSevere[region_idx]; - auto& num_death = vnum_death[region_idx]; - auto& num_icu = vnum_icu[region_idx]; - - auto& mu_C_R = vmu_C_R[region_idx]; - auto& mu_I_H = vmu_I_H[region_idx]; - auto& mu_H_U = vmu_H_U[region_idx]; - - auto date_df = region_entry.date; - auto age = size_t(region_entry.age_group); - - if (date_df == offset_date_by_days(date, 0)) { - num_InfectedSymptoms[age] += scaling_factor_inf[age] * region_entry.num_confirmed; - // We intentionally do NOT multiply recovered with the scaling_factor_inf here. - // If we apply the scaling factor to recovered as well, we would implicitly - // assume that this factor holds for the entire historical period up to t0, which is not valid. - num_rec[age] += region_entry.num_confirmed; - } - if (date_df == offset_date_by_days(date, days_surplus)) { - num_InfectedNoSymptoms[age] -= - 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed; - } - if (date_df == offset_date_by_days(date, t_InfectedNoSymptoms[age] + days_surplus)) { - num_InfectedNoSymptoms[age] += - 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed; - num_Exposed[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed; - } - if (date_df == offset_date_by_days(date, t_Exposed[age] + t_InfectedNoSymptoms[age] + days_surplus)) { - num_Exposed[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed; - } - if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age])) { - num_InfectedSymptoms[age] -= scaling_factor_inf[age] * region_entry.num_confirmed; - num_InfectedSevere[age] += mu_I_H[age] * scaling_factor_inf[age] * region_entry.num_confirmed; - } - if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age])) { - num_InfectedSevere[age] -= mu_I_H[age] * scaling_factor_inf[age] * region_entry.num_confirmed; - num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed; - } - if (date_df == - offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age] - t_InfectedCritical[age])) { - num_death[age] += region_entry.num_deaths; - num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed; - } - } - } - - for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) { - auto region = vregion[region_idx]; - - auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx]; - auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx]; - auto& num_rec = vnum_rec[region_idx]; - auto& num_Exposed = vnum_Exposed[region_idx]; - auto& num_InfectedSevere = vnum_InfectedSevere[region_idx]; - auto& num_death = vnum_death[region_idx]; - auto& num_icu = vnum_icu[region_idx]; - - for (size_t i = 0; i < ConfirmedCasesDataEntry::age_group_names.size(); i++) { - // R(t0) = ΣC(t0) − I(t0) − H(t0) − U(t0) − D(t0) - // subtract currently infectious/hospitalized/ICU/dead - // Note: We divide I/H/U/D by scaling_factor_inf to "unscale" these contributions back to the - // reported level before subtracting from recovered. If we also applied the scaling factor to - // recovered, we would implicitly assume that the same underreporting applies to the entire - // history up to t0, which would be wrong. The scaling factor should reflect underreporting - // around t0 only. - num_rec[i] -= - (num_InfectedSymptoms[i] / scaling_factor_inf[i] + num_InfectedSevere[i] / scaling_factor_inf[i] + - num_icu[i] / scaling_factor_inf[i] + num_death[i] / scaling_factor_inf[i]); - - auto try_fix_constraints = [region, i](double& value, double error, auto str) { - if (value < error) { - //this should probably return a failure - //but the algorithm is not robust enough to avoid large negative values and there are tests that rely on it - log_error("{:s} for age group {:s} is {:.4f} for region {:d}, exceeds expected negative value.", - str, ConfirmedCasesDataEntry::age_group_names[i], value, region); - value = 0.0; - } - else if (value < 0) { - log_info("{:s} for age group {:s} is {:.4f} for region {:d}, automatically corrected", str, - ConfirmedCasesDataEntry::age_group_names[i], value, region); - value = 0.0; - } - }; - - try_fix_constraints(num_InfectedSymptoms[i], -5, "InfectedSymptoms"); - try_fix_constraints(num_InfectedNoSymptoms[i], -5, "InfectedNoSymptoms"); - try_fix_constraints(num_Exposed[i], -5, "Exposed"); - try_fix_constraints(num_InfectedSevere[i], -5, "InfectedSevere"); - try_fix_constraints(num_death[i], -5, "Dead"); - try_fix_constraints(num_icu[i], -5, "InfectedCritical"); - try_fix_constraints(num_rec[i], -20, "Recovered"); - } - } - - return success(); -} - } // namespace details } // namespace osecir } // namespace mio diff --git a/cpp/models/ode_secir/parameters_io.h b/cpp/models/ode_secir/parameters_io.h index 0505d45563..f288abcc35 100644 --- a/cpp/models/ode_secir/parameters_io.h +++ b/cpp/models/ode_secir/parameters_io.h @@ -39,16 +39,20 @@ namespace osecir namespace details { + +// Forward declaration of get_region_id +int get_region_id(int id); /** - * @brief reads populations data from RKI. - * @param[in] rki_data Vector of ConfirmedCasesDataEntry%s. - * @param[in] vregion Vector of keys of the region of interest. - * @param[in] date Date at which the data is read. - * @param[in, out] vnum_* Output vector for number of people in the corresponding compartement. - * @param[in] vt_* vector Average time it takes to get from one compartement to another for each age group. - * @param[in] vmu_* vector Probabilities to get from one compartement to another for each age group. - * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. - */ + * @brief reads populations data from RKI. + * @param[in] rki_data Vector of ConfirmedCasesDataEntry%s. + * @param[in] vregion Vector of keys of the region of interest. + * @param[in] date Date at which the data is read. + * @param[in, out] vnum_* Output vector for number of people in the corresponding compartement. + * @param[in] vt_* vector Average time it takes to get from one compartement to another for each age group. + * @param[in] vmu_* vector Probabilities to get from one compartement to another for each age group. + * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. + */ +template IOResult read_confirmed_cases_data( std::vector& rki_data, std::vector const& vregion, Date date, std::vector>& vnum_Exposed, std::vector>& vnum_InfectedNoSymptoms, @@ -59,7 +63,147 @@ IOResult read_confirmed_cases_data( const std::vector>& vt_InfectedSymptoms, const std::vector>& vt_InfectedSevere, const std::vector>& vt_InfectedCritical, const std::vector>& vmu_C_R, const std::vector>& vmu_I_H, const std::vector>& vmu_H_U, - const std::vector& scaling_factor_inf); + const std::vector& scaling_factor_inf) +{ + auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { + return a.date < b.date; + }); + if (max_date_entry == rki_data.end()) { + log_error("RKI data file is empty."); + return failure(StatusCode::InvalidFileFormat, "RKI file is empty."); + } + auto max_date = max_date_entry->date; + if (max_date < date) { + log_error("Specified date does not exist in RKI data"); + return failure(StatusCode::OutOfRange, "Specified date does not exist in RKI data."); + } + auto days_surplus = std::min(get_offset_in_days(max_date, date) - 6, 0); + + //this statement causes maybe-uninitialized warning for some versions of gcc. + //the error is reported in an included header, so the warning is disabled for the whole file + std::sort(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { + return std::make_tuple(get_region_id(a), a.date) < std::make_tuple(get_region_id(b), b.date); + }); + + for (auto region_idx = size_t(0); region_idx < vregion.size(); ++region_idx) { + auto region_entry_range_it = + std::equal_range(rki_data.begin(), rki_data.end(), vregion[region_idx], [](auto&& a, auto&& b) { + return get_region_id(a) < get_region_id(b); + }); + auto region_entry_range = make_range(region_entry_range_it); + if (region_entry_range.begin() == region_entry_range.end()) { + log_error("No entries found for region {}", vregion[region_idx]); + return failure(StatusCode::InvalidFileFormat, + "No entries found for region " + std::to_string(vregion[region_idx])); + } + for (auto&& region_entry : region_entry_range) { + + auto& t_Exposed = vt_Exposed[region_idx]; + auto& t_InfectedNoSymptoms = vt_InfectedNoSymptoms[region_idx]; + auto& t_InfectedSymptoms = vt_InfectedSymptoms[region_idx]; + auto& t_InfectedSevere = vt_InfectedSevere[region_idx]; + auto& t_InfectedCritical = vt_InfectedCritical[region_idx]; + + auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx]; + auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx]; + auto& num_rec = vnum_rec[region_idx]; + auto& num_Exposed = vnum_Exposed[region_idx]; + auto& num_InfectedSevere = vnum_InfectedSevere[region_idx]; + auto& num_death = vnum_death[region_idx]; + auto& num_icu = vnum_icu[region_idx]; + + auto& mu_C_R = vmu_C_R[region_idx]; + auto& mu_I_H = vmu_I_H[region_idx]; + auto& mu_H_U = vmu_H_U[region_idx]; + + auto date_df = region_entry.date; + auto age = size_t(region_entry.age_group); + + if (date_df == offset_date_by_days(date, 0)) { + num_InfectedSymptoms[age] += scaling_factor_inf[age] * region_entry.num_confirmed; + // We intentionally do NOT multiply recovered with the scaling_factor_inf here. + // If we apply the scaling factor to recovered as well, we would implicitly + // assume that this factor holds for the entire historical period up to t0, which is not valid. + num_rec[age] += region_entry.num_confirmed; + } + if (date_df == offset_date_by_days(date, days_surplus)) { + num_InfectedNoSymptoms[age] -= + 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed; + } + if (date_df == offset_date_by_days(date, t_InfectedNoSymptoms[age] + days_surplus)) { + num_InfectedNoSymptoms[age] += + 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed; + num_Exposed[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed; + } + if (date_df == offset_date_by_days(date, t_Exposed[age] + t_InfectedNoSymptoms[age] + days_surplus)) { + num_Exposed[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed; + } + if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age])) { + num_InfectedSymptoms[age] -= scaling_factor_inf[age] * region_entry.num_confirmed; + num_InfectedSevere[age] += mu_I_H[age] * scaling_factor_inf[age] * region_entry.num_confirmed; + } + if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age])) { + num_InfectedSevere[age] -= mu_I_H[age] * scaling_factor_inf[age] * region_entry.num_confirmed; + num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed; + } + if (date_df == + offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age] - t_InfectedCritical[age])) { + num_death[age] += region_entry.num_deaths; + num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed; + } + } + } + + for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) { + auto region = vregion[region_idx]; + + auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx]; + auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx]; + auto& num_rec = vnum_rec[region_idx]; + auto& num_Exposed = vnum_Exposed[region_idx]; + auto& num_InfectedSevere = vnum_InfectedSevere[region_idx]; + auto& num_death = vnum_death[region_idx]; + auto& num_icu = vnum_icu[region_idx]; + + for (size_t i = 0; i < ConfirmedCasesDataEntry::age_group_names.size(); i++) { + // R(t0) = ΣC(t0) − I(t0) − H(t0) − U(t0) − D(t0) + // subtract currently infectious/hospitalized/ICU/dead + // Note: We divide I/H/U/D by scaling_factor_inf to "unscale" these contributions back to the + // reported level before subtracting from recovered. If we also applied the scaling factor to + // recovered, we would implicitly assume that the same underreporting applies to the entire + // history up to t0, which would be wrong. The scaling factor should reflect underreporting + // around t0 only. + num_rec[i] -= + (num_InfectedSymptoms[i] / scaling_factor_inf[i] + num_InfectedSevere[i] / scaling_factor_inf[i] + + num_icu[i] / scaling_factor_inf[i] + num_death[i] / scaling_factor_inf[i]); + + auto try_fix_constraints = [region, i](double& value, double error, auto str) { + if (value < error) { + //this should probably return a failure + //but the algorithm is not robust enough to avoid large negative values and there are tests that rely on it + log_error("{:s} for age group {:s} is {:.4f} for region {:d}, exceeds expected negative value.", + str, ConfirmedCasesDataEntry::age_group_names[i], value, region); + value = 0.0; + } + else if (value < 0) { + log_info("{:s} for age group {:s} is {:.4f} for region {:d}, automatically corrected", str, + ConfirmedCasesDataEntry::age_group_names[i], value, region); + value = 0.0; + } + }; + + try_fix_constraints(num_InfectedSymptoms[i], -5, "InfectedSymptoms"); + try_fix_constraints(num_InfectedNoSymptoms[i], -5, "InfectedNoSymptoms"); + try_fix_constraints(num_Exposed[i], -5, "Exposed"); + try_fix_constraints(num_InfectedSevere[i], -5, "InfectedSevere"); + try_fix_constraints(num_death[i], -5, "Dead"); + try_fix_constraints(num_icu[i], -5, "InfectedCritical"); + try_fix_constraints(num_rec[i], -20, "Recovered"); + } + } + + return success(); +} /** * @brief Sets populations data from already read case data with multiple age groups into a Model with one age group. diff --git a/cpp/models/ode_secirvvs/parameters_io.cpp b/cpp/models/ode_secirvvs/parameters_io.cpp index b1f0dcddb2..1d26785a8f 100644 --- a/cpp/models/ode_secirvvs/parameters_io.cpp +++ b/cpp/models/ode_secirvvs/parameters_io.cpp @@ -29,240 +29,6 @@ namespace osecirvvs { namespace details { -IOResult read_confirmed_cases_data( - std::string const& path, std::vector const& vregion, Date date, std::vector>& vnum_Exposed, - std::vector>& vnum_InfectedNoSymptoms, std::vector>& vnum_InfectedSymptoms, - std::vector>& vnum_InfectedSevere, std::vector>& vnum_icu, - std::vector>& vnum_death, std::vector>& vnum_rec, - const std::vector>& vt_Exposed, const std::vector>& vt_InfectedNoSymptoms, - const std::vector>& vt_InfectedSymptoms, const std::vector>& vt_InfectedSevere, - const std::vector>& vt_InfectedCritical, const std::vector>& vmu_C_R, - const std::vector>& vmu_I_H, const std::vector>& vmu_H_U, - const std::vector& scaling_factor_inf) -{ - BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path)); - return read_confirmed_cases_data(rki_data, vregion, date, vnum_Exposed, vnum_InfectedNoSymptoms, - vnum_InfectedSymptoms, vnum_InfectedSevere, vnum_icu, vnum_death, vnum_rec, - vt_Exposed, vt_InfectedNoSymptoms, vt_InfectedSymptoms, vt_InfectedSevere, - vt_InfectedCritical, vmu_C_R, vmu_I_H, vmu_H_U, scaling_factor_inf); -} - -IOResult read_confirmed_cases_data( - const std::vector& rki_data, std::vector const& vregion, Date date, - std::vector>& vnum_Exposed, std::vector>& vnum_InfectedNoSymptoms, - std::vector>& vnum_InfectedSymptoms, std::vector>& vnum_InfectedSevere, - std::vector>& vnum_icu, std::vector>& vnum_death, - std::vector>& vnum_rec, const std::vector>& vt_Exposed, - const std::vector>& vt_InfectedNoSymptoms, - const std::vector>& vt_InfectedSymptoms, const std::vector>& vt_InfectedSevere, - const std::vector>& vt_InfectedCritical, const std::vector>& vmu_C_R, - const std::vector>& vmu_I_H, const std::vector>& vmu_H_U, - const std::vector& scaling_factor_inf) -{ - auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { - return a.date < b.date; - }); - if (max_date_entry == rki_data.end()) { - log_error("RKI data file is empty."); - return failure(StatusCode::InvalidValue, "RKI data is empty."); - } - auto max_date = max_date_entry->date; - if (max_date < date) { - log_error("Specified date does not exist in RKI data"); - return failure(StatusCode::OutOfRange, "RKI data does not contain specified date."); - } - - // shifts the initilization to the recent past if simulation starts - // around current day and data of the future would be required. - // Only needed for preinfection compartments, exposed and InfectedNoSymptoms. - auto days_surplus = get_offset_in_days(max_date, date) - 6; // 6 > T_E + T_C - if (days_surplus > 0) { - days_surplus = 0; - } - - for (auto&& entry : rki_data) { - auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) { - return r == 0 || get_region_id(entry) == r; - }); - if (it != vregion.end()) { - auto region_idx = size_t(it - vregion.begin()); - - auto& t_Exposed = vt_Exposed[region_idx]; - auto& t_InfectedNoSymptoms = vt_InfectedNoSymptoms[region_idx]; - auto& t_InfectedSymptoms = vt_InfectedSymptoms[region_idx]; - auto& t_InfectedSevere = vt_InfectedSevere[region_idx]; - auto& t_InfectedCritical = vt_InfectedCritical[region_idx]; - - auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx]; - auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx]; - auto& num_rec = vnum_rec[region_idx]; - auto& num_Exposed = vnum_Exposed[region_idx]; - auto& num_InfectedSevere = vnum_InfectedSevere[region_idx]; - auto& num_death = vnum_death[region_idx]; - auto& num_icu = vnum_icu[region_idx]; - - auto& mu_C_R = vmu_C_R[region_idx]; - auto& mu_I_H = vmu_I_H[region_idx]; - auto& mu_H_U = vmu_H_U[region_idx]; - - auto age = (size_t)entry.age_group; - if (entry.date == offset_date_by_days(date, 0)) { - num_InfectedSymptoms[age] += scaling_factor_inf[age] * entry.num_confirmed; - num_rec[age] += entry.num_confirmed; - } - if (entry.date == offset_date_by_days(date, t_InfectedNoSymptoms[age] + days_surplus)) { - num_InfectedNoSymptoms[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed; - num_Exposed[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed; - } - if (entry.date == offset_date_by_days(date, days_surplus)) { - num_InfectedNoSymptoms[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed; - } - if (entry.date == offset_date_by_days(date, t_Exposed[age] + t_InfectedNoSymptoms[age] + days_surplus)) { - num_Exposed[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed; - } - if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms[age])) { - num_InfectedSymptoms[age] -= scaling_factor_inf[age] * entry.num_confirmed; - num_InfectedSevere[age] += mu_I_H[age] * scaling_factor_inf[age] * entry.num_confirmed; - } - if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age])) { - num_InfectedSevere[age] -= mu_I_H[age] * scaling_factor_inf[age] * entry.num_confirmed; - num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed; - } - if (entry.date == - offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age] - t_InfectedCritical[age])) { - num_death[age] += entry.num_deaths; - num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed; - } - } - } - - for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) { - auto region = vregion[region_idx]; - - auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx]; - auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx]; - auto& num_rec = vnum_rec[region_idx]; - auto& num_Exposed = vnum_Exposed[region_idx]; - auto& num_InfectedSevere = vnum_InfectedSevere[region_idx]; - auto& num_death = vnum_death[region_idx]; - auto& num_icu = vnum_icu[region_idx]; - - size_t num_groups = ConfirmedCasesDataEntry::age_group_names.size(); - for (size_t i = 0; i < num_groups; i++) { - // subtract infected confirmed cases which are not yet recovered - // and remove dark figure scaling factor - num_rec[i] -= num_InfectedSymptoms[i] / scaling_factor_inf[i]; - num_rec[i] -= num_InfectedSevere[i] / scaling_factor_inf[i]; - num_rec[i] -= - num_icu[i] / - scaling_factor_inf[i]; // TODO: this has to be adapted for scaling_factor_inf != 1 or != ***_icu - num_rec[i] -= num_death[i] / scaling_factor_inf[i]; - auto try_fix_constraints = [region, i](double& value, double error, auto str) { - if (value < error) { - // this should probably return a failure - // but the algorithm is not robust enough to avoid large negative - // values and there are tests that rely on it - log_error("{:s} for age group {:s} is {:.4f} for region {:d}, " - "exceeds expected negative value.", - str, ConfirmedCasesDataEntry::age_group_names[i], value, region); - value = 0.0; - } - else if (value < 0) { - log_info("{:s} for age group {:s} is {:.4f} for region {:d}, " - "automatically corrected", - str, ConfirmedCasesDataEntry::age_group_names[i], value, region); - value = 0.0; - } - }; - - try_fix_constraints(num_InfectedSymptoms[i], -5, "InfectedSymptoms"); - try_fix_constraints(num_InfectedNoSymptoms[i], -5, "InfectedNoSymptoms"); - try_fix_constraints(num_Exposed[i], -5, "Exposed"); - try_fix_constraints(num_InfectedSevere[i], -5, "InfectedSevere"); - try_fix_constraints(num_death[i], -5, "Dead"); - try_fix_constraints(num_icu[i], -5, "InfectedCritical"); - try_fix_constraints(num_rec[i], -20, "Recovered or vaccinated"); - } - } - - return success(); -} - -IOResult read_confirmed_cases_data_fix_recovered(std::string const& path, std::vector const& vregion, - Date date, std::vector>& vnum_rec, - double delay) -{ - BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path)); - return read_confirmed_cases_data_fix_recovered(rki_data, vregion, date, vnum_rec, delay); -} - -IOResult read_confirmed_cases_data_fix_recovered(const std::vector& rki_data, - std::vector const& vregion, Date date, - std::vector>& vnum_rec, double delay) -{ - auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { - return a.date < b.date; - }); - if (max_date_entry == rki_data.end()) { - log_error("RKI data is empty."); - return failure(StatusCode::InvalidValue, "RKI data is empty."); - } - auto max_date = max_date_entry->date; - if (max_date < date) { - log_error("Specified date does not exist in RKI data"); - return failure(StatusCode::OutOfRange, "RKI data does not contain specified date."); - } - - // shifts the initilization to the recent past if simulation starts - // around current day and data of the future would be required. - // Only needed for preinfection compartments, exposed and InfectedNoSymptoms. - auto days_surplus = get_offset_in_days(max_date, date) - 6; // 6 > T_E^C + T_C^I - if (days_surplus > 0) { - days_surplus = 0; - } - - for (auto&& rki_entry : rki_data) { - auto it = std::find_if(vregion.begin(), vregion.end(), [&rki_entry](auto r) { - return r == 0 || get_region_id(rki_entry) == r; - }); - if (it != vregion.end()) { - auto region_idx = size_t(it - vregion.begin()); - if (rki_entry.date == offset_date_by_days(date, int(-delay))) { - vnum_rec[region_idx][size_t(rki_entry.age_group)] = rki_entry.num_confirmed; - } - } - } - - for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) { - auto region = vregion[region_idx]; - auto& num_rec = vnum_rec[region_idx]; - - size_t num_groups = ConfirmedCasesDataEntry::age_group_names.size(); - for (size_t i = 0; i < num_groups; i++) { - auto try_fix_constraints = [region, i](double& value, double error, auto str) { - if (value < error) { - // this should probably return a failure - // but the algorithm is not robust enough to avoid large negative - // values and there are tests that rely on it - log_error("{:s} for age group {:s} is {:.4f} for region {:d}, " - "exceeds expected negative value.", - str, ConfirmedCasesDataEntry::age_group_names[i], value, region); - value = 0.0; - } - else if (value < 0) { - log_info("{:s} for age group {:s} is {:.4f} for region {:d}, " - "automatically corrected", - str, ConfirmedCasesDataEntry::age_group_names[i], value, region); - value = 0.0; - } - }; - try_fix_constraints(num_rec[i], 0, "Recovered"); - } - } - - return success(); -} - } // namespace details } // namespace osecirvvs } // namespace mio diff --git a/cpp/models/ode_secirvvs/parameters_io.h b/cpp/models/ode_secirvvs/parameters_io.h index 3cd0a412c0..1e592b8951 100644 --- a/cpp/models/ode_secirvvs/parameters_io.h +++ b/cpp/models/ode_secirvvs/parameters_io.h @@ -41,67 +41,277 @@ namespace osecirvvs namespace details { /** - * @brief Reads subpopulations of infection states from transformed RKI cases file. - * @param[in] path Path to transformed RKI cases file. - * @param[in] vregion Vector of keys of the region of interest. - * @param[in] date Date for which the arrays are initialized. - * @param[in, out] num_* Output vector for number of people in the corresponding compartement. - * @param[in] vt_* Average time it takes to get from one compartement to another (vector with one element per age group). - * @param[in] vmu_* Probabilities to get from one compartement to another (vector with one element per age group). - * @param[in] scaling_factor_inf Factor for scaling the confirmed cases to account for estimated undetected cases. - * @see mio::read_confirmed_cases_data - * @{ - */ +* @brief Reads subpopulations of infection states from transformed RKI cases file. +* @param[in] path Path to transformed RKI cases file. +* @param[in] vregion Vector of keys of the region of interest. +* @param[in] date Date for which the arrays are initialized. +* @param[in, out] num_* Output vector for number of people in the corresponding compartement. +* @param[in] vt_* Average time it takes to get from one compartement to another (vector with one element per age group). +* @param[in] vmu_* Probabilities to get from one compartement to another (vector with one element per age group). +* @param[in] scaling_factor_inf Factor for scaling the confirmed cases to account for estimated undetected cases. +* @see mio::read_confirmed_cases_data +* @{ +*/ +template IOResult read_confirmed_cases_data( - std::string const& path, std::vector const& vregion, Date date, std::vector>& num_Exposed, - std::vector>& num_InfectedNoSymptoms, std::vector>& num_InfectedSymptoms, - std::vector>& num_InfectedSevere, std::vector>& num_icu, - std::vector>& num_death, std::vector>& num_rec, - const std::vector>& t_Exposed, const std::vector>& t_InfectedNoSymptoms, - const std::vector>& t_InfectedSymptoms, const std::vector>& t_InfectedSevere, - const std::vector>& t_InfectedCritical, const std::vector>& mu_C_R, - const std::vector>& mu_I_H, const std::vector>& mu_H_U, - const std::vector& scaling_factor_inf); + const std::vector& rki_data, std::vector const& vregion, Date date, + std::vector>& vnum_Exposed, std::vector>& vnum_InfectedNoSymptoms, + std::vector>& vnum_InfectedSymptoms, std::vector>& vnum_InfectedSevere, + std::vector>& vnum_icu, std::vector>& vnum_death, + std::vector>& vnum_rec, const std::vector>& vt_Exposed, + const std::vector>& vt_InfectedNoSymptoms, + const std::vector>& vt_InfectedSymptoms, const std::vector>& vt_InfectedSevere, + const std::vector>& vt_InfectedCritical, const std::vector>& vmu_C_R, + const std::vector>& vmu_I_H, const std::vector>& vmu_H_U, + const std::vector& scaling_factor_inf) +{ + auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { + return a.date < b.date; + }); + if (max_date_entry == rki_data.end()) { + log_error("RKI data file is empty."); + return failure(StatusCode::InvalidValue, "RKI data is empty."); + } + auto max_date = max_date_entry->date; + if (max_date < date) { + log_error("Specified date does not exist in RKI data"); + return failure(StatusCode::OutOfRange, "RKI data does not contain specified date."); + } + + // shifts the initilization to the recent past if simulation starts + // around current day and data of the future would be required. + // Only needed for preinfection compartments, exposed and InfectedNoSymptoms. + auto days_surplus = get_offset_in_days(max_date, date) - 6; // 6 > T_E + T_C + if (days_surplus > 0) { + days_surplus = 0; + } + + for (auto&& entry : rki_data) { + auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) { + return r == 0 || get_region_id(entry) == r; + }); + if (it != vregion.end()) { + auto region_idx = size_t(it - vregion.begin()); + + auto& t_Exposed = vt_Exposed[region_idx]; + auto& t_InfectedNoSymptoms = vt_InfectedNoSymptoms[region_idx]; + auto& t_InfectedSymptoms = vt_InfectedSymptoms[region_idx]; + auto& t_InfectedSevere = vt_InfectedSevere[region_idx]; + auto& t_InfectedCritical = vt_InfectedCritical[region_idx]; + + auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx]; + auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx]; + auto& num_rec = vnum_rec[region_idx]; + auto& num_Exposed = vnum_Exposed[region_idx]; + auto& num_InfectedSevere = vnum_InfectedSevere[region_idx]; + auto& num_death = vnum_death[region_idx]; + auto& num_icu = vnum_icu[region_idx]; + + auto& mu_C_R = vmu_C_R[region_idx]; + auto& mu_I_H = vmu_I_H[region_idx]; + auto& mu_H_U = vmu_H_U[region_idx]; + + auto age = (size_t)entry.age_group; + if (entry.date == offset_date_by_days(date, 0)) { + num_InfectedSymptoms[age] += scaling_factor_inf[age] * entry.num_confirmed; + num_rec[age] += entry.num_confirmed; + } + if (entry.date == offset_date_by_days(date, t_InfectedNoSymptoms[age] + days_surplus)) { + num_InfectedNoSymptoms[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed; + num_Exposed[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed; + } + if (entry.date == offset_date_by_days(date, days_surplus)) { + num_InfectedNoSymptoms[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed; + } + if (entry.date == offset_date_by_days(date, t_Exposed[age] + t_InfectedNoSymptoms[age] + days_surplus)) { + num_Exposed[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed; + } + if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms[age])) { + num_InfectedSymptoms[age] -= scaling_factor_inf[age] * entry.num_confirmed; + num_InfectedSevere[age] += mu_I_H[age] * scaling_factor_inf[age] * entry.num_confirmed; + } + if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age])) { + num_InfectedSevere[age] -= mu_I_H[age] * scaling_factor_inf[age] * entry.num_confirmed; + num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed; + } + if (entry.date == + offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age] - t_InfectedCritical[age])) { + num_death[age] += entry.num_deaths; + num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed; + } + } + } + + for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) { + auto region = vregion[region_idx]; + + auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx]; + auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx]; + auto& num_rec = vnum_rec[region_idx]; + auto& num_Exposed = vnum_Exposed[region_idx]; + auto& num_InfectedSevere = vnum_InfectedSevere[region_idx]; + auto& num_death = vnum_death[region_idx]; + auto& num_icu = vnum_icu[region_idx]; + + size_t num_groups = ConfirmedCasesDataEntry::age_group_names.size(); + for (size_t i = 0; i < num_groups; i++) { + // subtract infected confirmed cases which are not yet recovered + // and remove dark figure scaling factor + num_rec[i] -= num_InfectedSymptoms[i] / scaling_factor_inf[i]; + num_rec[i] -= num_InfectedSevere[i] / scaling_factor_inf[i]; + num_rec[i] -= + num_icu[i] / + scaling_factor_inf[i]; // TODO: this has to be adapted for scaling_factor_inf != 1 or != ***_icu + num_rec[i] -= num_death[i] / scaling_factor_inf[i]; + auto try_fix_constraints = [region, i](double& value, double error, auto str) { + if (value < error) { + // this should probably return a failure + // but the algorithm is not robust enough to avoid large negative + // values and there are tests that rely on it + log_error("{:s} for age group {:s} is {:.4f} for region {:d}, " + "exceeds expected negative value.", + str, ConfirmedCasesDataEntry::age_group_names[i], value, region); + value = 0.0; + } + else if (value < 0) { + log_info("{:s} for age group {:s} is {:.4f} for region {:d}, " + "automatically corrected", + str, ConfirmedCasesDataEntry::age_group_names[i], value, region); + value = 0.0; + } + }; + + try_fix_constraints(num_InfectedSymptoms[i], -5, "InfectedSymptoms"); + try_fix_constraints(num_InfectedNoSymptoms[i], -5, "InfectedNoSymptoms"); + try_fix_constraints(num_Exposed[i], -5, "Exposed"); + try_fix_constraints(num_InfectedSevere[i], -5, "InfectedSevere"); + try_fix_constraints(num_death[i], -5, "Dead"); + try_fix_constraints(num_icu[i], -5, "InfectedCritical"); + try_fix_constraints(num_rec[i], -20, "Recovered or vaccinated"); + } + } + + return success(); +} +template IOResult read_confirmed_cases_data( - const std::vector& rki_data, std::vector const& vregion, Date date, - std::vector>& num_Exposed, std::vector>& num_InfectedNoSymptoms, - std::vector>& num_InfectedSymptoms, std::vector>& num_InfectedSevere, - std::vector>& num_icu, std::vector>& num_death, - std::vector>& num_rec, const std::vector>& t_Exposed, - const std::vector>& t_InfectedNoSymptoms, const std::vector>& t_InfectedSymptoms, - const std::vector>& t_InfectedSevere, const std::vector>& t_InfectedCritical, - const std::vector>& mu_C_R, const std::vector>& mu_I_H, - const std::vector>& mu_H_U, const std::vector& scaling_factor_inf); + std::string const& path, std::vector const& vregion, Date date, std::vector>& vnum_Exposed, + std::vector>& vnum_InfectedNoSymptoms, std::vector>& vnum_InfectedSymptoms, + std::vector>& vnum_InfectedSevere, std::vector>& vnum_icu, + std::vector>& vnum_death, std::vector>& vnum_rec, + const std::vector>& vt_Exposed, const std::vector>& vt_InfectedNoSymptoms, + const std::vector>& vt_InfectedSymptoms, const std::vector>& vt_InfectedSevere, + const std::vector>& vt_InfectedCritical, const std::vector>& vmu_C_R, + const std::vector>& vmu_I_H, const std::vector>& vmu_H_U, + const std::vector& scaling_factor_inf) +{ + BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path)); + return read_confirmed_cases_data(rki_data, vregion, date, vnum_Exposed, vnum_InfectedNoSymptoms, + vnum_InfectedSymptoms, vnum_InfectedSevere, vnum_icu, vnum_death, vnum_rec, + vt_Exposed, vt_InfectedNoSymptoms, vt_InfectedSymptoms, vt_InfectedSevere, + vt_InfectedCritical, vmu_C_R, vmu_I_H, vmu_H_U, scaling_factor_inf); +} + /**@}*/ /** - * @brief Reads confirmed cases data and translates data of day t0-delay to recovered compartment. - * @param[in] path Path to RKI confirmed cases file. - * @param[in] vregion Vector of keys of the region of interest. - * @param[in] date Date for which the arrays are initialized. - * @param[in, out] num_rec Output vector for number of people in the compartement recovered. - * @param[in] delay Number of days in the past the are used to set recovered compartment. - * @see mio::read_confirmed_cases_data - * @{ - */ +* @brief Reads confirmed cases data and translates data of day t0-delay to recovered compartment. +* @param[in] path Path to RKI confirmed cases file. +* @param[in] vregion Vector of keys of the region of interest. +* @param[in] date Date for which the arrays are initialized. +* @param[in, out] num_rec Output vector for number of people in the compartement recovered. +* @param[in] delay Number of days in the past the are used to set recovered compartment. +* @see mio::read_confirmed_cases_data +* @{ +*/ +template IOResult read_confirmed_cases_data_fix_recovered(const std::vector& rki_data, std::vector const& vregion, Date date, - std::vector>& vnum_rec, double delay = 14.); + std::vector>& vnum_rec, double delay = 14.0) +{ + auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { + return a.date < b.date; + }); + if (max_date_entry == rki_data.end()) { + log_error("RKI data is empty."); + return failure(StatusCode::InvalidValue, "RKI data is empty."); + } + auto max_date = max_date_entry->date; + if (max_date < date) { + log_error("Specified date does not exist in RKI data"); + return failure(StatusCode::OutOfRange, "RKI data does not contain specified date."); + } + + // shifts the initilization to the recent past if simulation starts + // around current day and data of the future would be required. + // Only needed for preinfection compartments, exposed and InfectedNoSymptoms. + auto days_surplus = get_offset_in_days(max_date, date) - 6; // 6 > T_E^C + T_C^I + if (days_surplus > 0) { + days_surplus = 0; + } + + for (auto&& rki_entry : rki_data) { + auto it = std::find_if(vregion.begin(), vregion.end(), [&rki_entry](auto r) { + return r == 0 || get_region_id(rki_entry) == r; + }); + if (it != vregion.end()) { + auto region_idx = size_t(it - vregion.begin()); + if (rki_entry.date == offset_date_by_days(date, int(-delay))) { + vnum_rec[region_idx][size_t(rki_entry.age_group)] = rki_entry.num_confirmed; + } + } + } + + for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) { + auto region = vregion[region_idx]; + auto& num_rec = vnum_rec[region_idx]; + + size_t num_groups = ConfirmedCasesDataEntry::age_group_names.size(); + for (size_t i = 0; i < num_groups; i++) { + auto try_fix_constraints = [region, i](double& value, double error, auto str) { + if (value < error) { + // this should probably return a failure + // but the algorithm is not robust enough to avoid large negative + // values and there are tests that rely on it + log_error("{:s} for age group {:s} is {:.4f} for region {:d}, " + "exceeds expected negative value.", + str, ConfirmedCasesDataEntry::age_group_names[i], value, region); + value = 0.0; + } + else if (value < 0) { + log_info("{:s} for age group {:s} is {:.4f} for region {:d}, " + "automatically corrected", + str, ConfirmedCasesDataEntry::age_group_names[i], value, region); + value = 0.0; + } + }; + try_fix_constraints(num_rec[i], 0, "Recovered"); + } + } + + return success(); +} + +template IOResult read_confirmed_cases_data_fix_recovered(std::string const& path, std::vector const& vregion, Date date, std::vector>& vnum_rec, - double delay = 14.); + double delay = 14.0) +{ + BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path)); + return read_confirmed_cases_data_fix_recovered(rki_data, vregion, date, vnum_rec, delay); +} /**@}*/ /** - * @brief Sets the confirmed cases data for a vector of models based on input data. - * @param[in, out] model Vector of objects in which the data is set. - * @param[in] case_data Vector of case data. Each inner vector represents a different region. - * @param[in] region Vector of keys of the region of interest. - * @param[in] date Date for which the arrays are initialized. - * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of RKI data. - * @param[in] set_death If true, set the number of deaths. - */ +* @brief Sets the confirmed cases data for a vector of models based on input data. +* @param[in, out] model Vector of objects in which the data is set. +* @param[in] case_data Vector of case data. Each inner vector represents a different region. +* @param[in] region Vector of keys of the region of interest. +* @param[in] date Date for which the arrays are initialized. +* @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of RKI data. +* @param[in] set_death If true, set the number of deaths. +*/ template IOResult set_confirmed_cases_data(std::vector& model, const std::vector& case_data, @@ -404,15 +614,15 @@ IOResult set_confirmed_cases_data(std::vector& model, } /** - * @brief sets populations data from a transformed RKI cases file into a Model. - * @param[in, out] model vector of objects in which the data is set - * @param[in] path Path to transformed RKI cases file - * @param[in] region vector of keys of the region of interest - * @param[in] date Date for which the arrays are initialized - * @param[in] scaling_factor_inf factors by which to scale the confirmed cases of - * rki data - * @param set_death[in] If true, set the number of deaths. - */ +* @brief sets populations data from a transformed RKI cases file into a Model. +* @param[in, out] model vector of objects in which the data is set +* @param[in] path Path to transformed RKI cases file +* @param[in] region vector of keys of the region of interest +* @param[in] date Date for which the arrays are initialized +* @param[in] scaling_factor_inf factors by which to scale the confirmed cases of +* rki data +* @param set_death[in] If true, set the number of deaths. +*/ template IOResult set_confirmed_cases_data(std::vector& model, const std::string& path, std::vector const& region, Date date, @@ -424,13 +634,13 @@ IOResult set_confirmed_cases_data(std::vector& model, const std::st } /** - * @brief sets populations data from DIVI register into Model - * @param[in, out] model vector of objects in which the data is set - * @param[in] path Path to transformed DIVI file - * @param[in] vregion vector of keys of the regions of interest - * @param[in] date Date for which the arrays are initialized - * @param[in] scaling_factor_icu factor by which to scale the icu cases of divi data - */ +* @brief sets populations data from DIVI register into Model +* @param[in, out] model vector of objects in which the data is set +* @param[in] path Path to transformed DIVI file +* @param[in] vregion vector of keys of the regions of interest +* @param[in] date Date for which the arrays are initialized +* @param[in] scaling_factor_icu factor by which to scale the icu cases of divi data +*/ template IOResult set_divi_data(std::vector& model, const std::string& path, const std::vector& vregion, Date date, double scaling_factor_icu) From e9c2cbced7f9bdab447644664c10a1ac990489b3 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 15 Sep 2025 12:35:08 +0000 Subject: [PATCH 2/4] I like to move it move it --- cpp/memilio/io/result_io.cpp | 125 +----------------------- cpp/memilio/io/result_io.h | 127 ++++++++++++++++++++++++- cpp/models/ode_secir/parameters_io.cpp | 2 +- 3 files changed, 128 insertions(+), 126 deletions(-) diff --git a/cpp/memilio/io/result_io.cpp b/cpp/memilio/io/result_io.cpp index 4c3f87971e..f8146ec66f 100644 --- a/cpp/memilio/io/result_io.cpp +++ b/cpp/memilio/io/result_io.cpp @@ -31,6 +31,7 @@ namespace mio { + herr_t store_group_name(hid_t /*id*/, const char* name, const H5L_info_t* /*linfo*/, void* opdata) { auto h5group_names = reinterpret_cast*>(opdata); @@ -38,130 +39,6 @@ herr_t store_group_name(hid_t /*id*/, const char* name, const H5L_info_t* /*linf return 0; } -IOResult> read_result(const std::string& filename) -{ - std::vector results; - - H5File file{H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)}; - MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, filename); - - std::vector h5group_names; - MEMILIO_H5_CHECK(H5Literate(file.id, H5_INDEX_NAME, H5_ITER_INC, NULL, &store_group_name, &h5group_names), - StatusCode::UnknownError, "Group names could not be read."); - - for (auto& h5group_name : h5group_names) { - H5Group region_h5group{H5Gopen(file.id, h5group_name.c_str(), H5P_DEFAULT)}; - MEMILIO_H5_CHECK(region_h5group.id, StatusCode::UnknownError, - "Group could not be opened (" + h5group_name + ")"); - - std::vector h5dset_names; - 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."); - - 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)}; - MEMILIO_H5_CHECK(dataset_t.id, StatusCode::UnknownError, "Time DataSet could not be read."); - - // dataset dimensions - H5DataSpace dataspace_t{H5Dget_space(dataset_t.id)}; - MEMILIO_H5_CHECK(dataspace_t.id, StatusCode::UnknownError, "Time DataSpace could not be read."); - const auto n_dims_t = 1; - hsize_t dims_t[n_dims_t]; - H5Sget_simple_extent_dims(dataspace_t.id, dims_t, NULL); - - auto num_timepoints = Eigen::Index(dims_t[0]); - auto time = std::vector(num_timepoints); - MEMILIO_H5_CHECK(H5Dread(dataset_t.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, time.data()), - StatusCode::UnknownError, "Time data could not be read."); - - auto dataset_name_total("/" + h5group_name + "/Total"); - H5DataSet dataset_total{H5Dopen(file.id, dataset_name_total.c_str(), H5P_DEFAULT)}; - MEMILIO_H5_CHECK(dataset_total.id, StatusCode::UnknownError, "Totals DataSet could not be read."); - - //read data space dimensions - H5DataSpace dataspace_total{H5Dget_space(dataset_total.id)}; - MEMILIO_H5_CHECK(dataspace_total.id, StatusCode::UnknownError, "Totals DataSpace could not be read."); - const auto n_dims_total = 2; - hsize_t dims_total[n_dims_total]; - H5Sget_simple_extent_dims(dataspace_total.id, dims_total, NULL); - if (num_timepoints != Eigen::Index(dims_total[0])) { - return failure(StatusCode::InvalidFileFormat, "Number of time points does not match."); - } - auto num_infectionstates = Eigen::Index(dims_total[1]); - - auto total_values = Eigen::Matrix( - num_timepoints, num_infectionstates); - MEMILIO_H5_CHECK( - H5Dread(dataset_total.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, total_values.data()), - StatusCode::UnknownError, "Totals data could not be read"); - - auto totals = TimeSeries(num_infectionstates); - totals.reserve(num_timepoints); - for (auto t_idx = 0; t_idx < num_timepoints; ++t_idx) { - totals.add_time_point(time[t_idx], slice(total_values, {t_idx, 1}, {0, num_infectionstates}).transpose()); - } - - auto groups = TimeSeries(num_infectionstates * num_groups); - groups.reserve(num_timepoints); - for (Eigen::Index t_idx = 0; t_idx < num_timepoints; ++t_idx) { - groups.add_time_point(time[t_idx]); - } - - 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."); - - //read data space dimensions - H5DataSpace dataspace_values{H5Dget_space(dataset_values.id)}; - MEMILIO_H5_CHECK(dataspace_values.id, StatusCode::UnknownError, "Values DataSpace could not be read."); - const auto n_dims_values = 2; - hsize_t dims_values[n_dims_values]; - H5Sget_simple_extent_dims(dataspace_values.id, dims_values, NULL); - if (num_timepoints != Eigen::Index(dims_values[0])) { - return failure(StatusCode::InvalidFileFormat, "Number of time points does not match."); - } - if (num_infectionstates != Eigen::Index(dims_values[1])) { - return failure(StatusCode::InvalidFileFormat, "Number of infection states does not match."); - } - - auto group_values = Eigen::Matrix( - num_timepoints, num_infectionstates); - MEMILIO_H5_CHECK( - H5Dread(dataset_values.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, group_values.data()), - StatusCode::UnknownError, "Values data could not be read"); - - 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 * h5_key_indx + idx_c] = group_values(idx_t, idx_c); - } - } - } - - results.push_back(SimulationResult(groups, totals)); - } - return success(results); -} - } // namespace mio #endif //MEMILIO_HAS_HDF5 diff --git a/cpp/memilio/io/result_io.h b/cpp/memilio/io/result_io.h index d8e2fd1cbc..38fb060adf 100644 --- a/cpp/memilio/io/result_io.h +++ b/cpp/memilio/io/result_io.h @@ -149,11 +149,136 @@ class SimulationResult TimeSeries m_totals; }; +// Forward declaration of store_group_name +herr_t store_group_name(hid_t /*id*/, const char* name, const H5L_info_t* /*linfo*/, void* opdata); /** * @brief Read simulation result from h5 file. * @param filename name of the file to be read. */ -IOResult> read_result(const std::string& filename); +template +IOResult> read_result(const std::string& filename) +{ + std::vector results; + + H5File file{H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)}; + MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, filename); + + std::vector h5group_names; + MEMILIO_H5_CHECK(H5Literate(file.id, H5_INDEX_NAME, H5_ITER_INC, NULL, &store_group_name, &h5group_names), + StatusCode::UnknownError, "Group names could not be read."); + + for (auto& h5group_name : h5group_names) { + H5Group region_h5group{H5Gopen(file.id, h5group_name.c_str(), H5P_DEFAULT)}; + MEMILIO_H5_CHECK(region_h5group.id, StatusCode::UnknownError, + "Group could not be opened (" + h5group_name + ")"); + + std::vector h5dset_names; + 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."); + + 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)}; + MEMILIO_H5_CHECK(dataset_t.id, StatusCode::UnknownError, "Time DataSet could not be read."); + + // dataset dimensions + H5DataSpace dataspace_t{H5Dget_space(dataset_t.id)}; + MEMILIO_H5_CHECK(dataspace_t.id, StatusCode::UnknownError, "Time DataSpace could not be read."); + const auto n_dims_t = 1; + hsize_t dims_t[n_dims_t]; + H5Sget_simple_extent_dims(dataspace_t.id, dims_t, NULL); + + auto num_timepoints = Eigen::Index(dims_t[0]); + auto time = std::vector(num_timepoints); + MEMILIO_H5_CHECK(H5Dread(dataset_t.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, time.data()), + StatusCode::UnknownError, "Time data could not be read."); + + auto dataset_name_total("/" + h5group_name + "/Total"); + H5DataSet dataset_total{H5Dopen(file.id, dataset_name_total.c_str(), H5P_DEFAULT)}; + MEMILIO_H5_CHECK(dataset_total.id, StatusCode::UnknownError, "Totals DataSet could not be read."); + + //read data space dimensions + H5DataSpace dataspace_total{H5Dget_space(dataset_total.id)}; + MEMILIO_H5_CHECK(dataspace_total.id, StatusCode::UnknownError, "Totals DataSpace could not be read."); + const auto n_dims_total = 2; + hsize_t dims_total[n_dims_total]; + H5Sget_simple_extent_dims(dataspace_total.id, dims_total, NULL); + if (num_timepoints != Eigen::Index(dims_total[0])) { + return failure(StatusCode::InvalidFileFormat, "Number of time points does not match."); + } + auto num_infectionstates = Eigen::Index(dims_total[1]); + + auto total_values = Eigen::Matrix( + num_timepoints, num_infectionstates); + MEMILIO_H5_CHECK( + H5Dread(dataset_total.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, total_values.data()), + StatusCode::UnknownError, "Totals data could not be read"); + + auto totals = TimeSeries(num_infectionstates); + totals.reserve(num_timepoints); + for (auto t_idx = 0; t_idx < num_timepoints; ++t_idx) { + totals.add_time_point(time[t_idx], slice(total_values, {t_idx, 1}, {0, num_infectionstates}).transpose()); + } + + auto groups = TimeSeries(num_infectionstates * num_groups); + groups.reserve(num_timepoints); + for (Eigen::Index t_idx = 0; t_idx < num_timepoints; ++t_idx) { + groups.add_time_point(time[t_idx]); + } + + 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."); + + //read data space dimensions + H5DataSpace dataspace_values{H5Dget_space(dataset_values.id)}; + MEMILIO_H5_CHECK(dataspace_values.id, StatusCode::UnknownError, "Values DataSpace could not be read."); + const auto n_dims_values = 2; + hsize_t dims_values[n_dims_values]; + H5Sget_simple_extent_dims(dataspace_values.id, dims_values, NULL); + if (num_timepoints != Eigen::Index(dims_values[0])) { + return failure(StatusCode::InvalidFileFormat, "Number of time points does not match."); + } + if (num_infectionstates != Eigen::Index(dims_values[1])) { + return failure(StatusCode::InvalidFileFormat, "Number of infection states does not match."); + } + + auto group_values = Eigen::Matrix( + num_timepoints, num_infectionstates); + MEMILIO_H5_CHECK( + H5Dread(dataset_values.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, group_values.data()), + StatusCode::UnknownError, "Values data could not be read"); + + 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 * h5_key_indx + idx_c] = group_values(idx_t, idx_c); + } + } + } + + results.push_back(SimulationResult(groups, totals)); + } + return success(results); +} /** * Save the results and the parameters of a single graph simulation run. diff --git a/cpp/models/ode_secir/parameters_io.cpp b/cpp/models/ode_secir/parameters_io.cpp index 49ca676cf0..0b57b297ba 100644 --- a/cpp/models/ode_secir/parameters_io.cpp +++ b/cpp/models/ode_secir/parameters_io.cpp @@ -43,7 +43,7 @@ namespace osecir namespace details { -//overload for integers, so the comparison of data entry to integers is symmetric (required by e.g. equal_range) +// Overload for integers, so the comparison of data entry to integers is symmetric (required by e.g. equal_range) int get_region_id(int id) { return id; From 24ec698ab7b5411fbb8183f581a1c59d28071240 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 15 Sep 2025 12:44:49 +0000 Subject: [PATCH 3/4] Pycode templates --- .../memilio/simulation/bindings/models/osecir.cpp | 11 ++++++----- .../memilio/simulation/bindings/models/osecirvvs.cpp | 2 +- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp index cd8d6fcbe6..ed3ad0bf97 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp @@ -289,9 +289,9 @@ PYBIND11_MODULE(_simulation_osecir, m) auto result = mio::set_edges, mio::MobilityParameters, mio::MobilityCoefficientGroup, mio::osecir::InfectionState, - decltype(mio::read_mobility_plain)>(mobility_data_file, params_graph, - mobile_comp, contact_locations_size, - mio::read_mobility_plain, weights); + decltype(mio::read_mobility_plain)>( + mobility_data_file, params_graph, mobile_comp, contact_locations_size, mio::read_mobility_plain, + weights); return pymio::check_and_throw(result); }, py::return_value_policy::move); @@ -314,8 +314,9 @@ PYBIND11_MODULE(_simulation_osecir, m) py::return_value_policy::move); #endif // MEMILIO_HAS_JSONCPP - m.def("interpolate_simulation_result", py::overload_cast( - &mio::interpolate_simulation_result>)); + m.def("interpolate_simulation_result", + py::overload_cast( + &mio::interpolate_simulation_result>)); m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results); diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp index c7919b4237..08dd9aaf31 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp @@ -342,7 +342,7 @@ PYBIND11_MODULE(_simulation_osecirvvs, m) auto result = mio::set_edges, mio::MobilityParameters, mio::MobilityCoefficientGroup, - mio::osecirvvs::InfectionState, decltype(mio::read_mobility_plain)>( + mio::osecirvvs::InfectionState, decltype(mio::read_mobility_plain)>( mobility_data_file, params_graph, mobile_comp, contact_locations_size, mio::read_mobility_plain, weights); return pymio::check_and_throw(result); From 71833e4c4ac963dca3543289d20977c28c87f9ea Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 15 Sep 2025 12:48:20 +0000 Subject: [PATCH 4/4] mio::read_mobility_plain --- .../memilio/simulation/bindings/models/osecirvvs.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp index 08dd9aaf31..e14e8b34f4 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp @@ -343,7 +343,7 @@ PYBIND11_MODULE(_simulation_osecirvvs, m) ContactLocation, mio::osecirvvs::Model, mio::MobilityParameters, mio::MobilityCoefficientGroup, mio::osecirvvs::InfectionState, decltype(mio::read_mobility_plain)>( - mobility_data_file, params_graph, mobile_comp, contact_locations_size, mio::read_mobility_plain, + mobility_data_file, params_graph, mobile_comp, contact_locations_size, mio::read_mobility_plain, weights); return pymio::check_and_throw(result); },