From 8e9f7006888039ae7a5e960d19a172644d97469c Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 15 Sep 2025 16:15:55 +0000 Subject: [PATCH 1/7] Add support for reading templated SECIR models from JSON --- cpp/examples/ode_secir_parameter_study.cpp | 4 +- .../ode_secir_parameter_study_graph.cpp | 9 +- cpp/examples/ode_secir_read_graph.cpp | 4 +- cpp/examples/ode_secir_save_results.cpp | 3 +- cpp/memilio/ad/ad.h | 68 ++ cpp/memilio/data/analyze_result.h | 4 +- cpp/memilio/io/mobility_io.h | 72 +- cpp/memilio/io/parameters_io.h | 23 +- cpp/memilio/io/result_io.h | 162 ++--- cpp/memilio/mobility/graph.h | 11 +- cpp/models/ode_secir/parameters_io.h | 192 +++--- cpp/models/ode_secirts/analyze_result.h | 6 +- cpp/models/ode_secirts/parameters_io.h | 161 ++--- cpp/models/ode_secirvvs/analyze_result.h | 84 +-- cpp/models/ode_secirvvs/parameters_io.h | 637 +++++++++--------- cpp/tests/test_analyze_result.cpp | 8 +- cpp/tests/test_mobility_io.cpp | 4 +- cpp/tests/test_odesecir.cpp | 88 +-- cpp/tests/test_odesecirts.cpp | 49 +- cpp/tests/test_odesecirvvs.cpp | 70 +- cpp/tests/test_save_parameters.cpp | 21 +- cpp/tests/test_save_results.cpp | 62 +- .../simulation/bindings/io/result_io.h | 4 +- .../simulation/bindings/models/osecir.cpp | 20 +- .../simulation/bindings/models/osecirvvs.cpp | 6 +- .../simulation/bindings/simulation.cpp | 18 +- 26 files changed, 948 insertions(+), 842 deletions(-) diff --git a/cpp/examples/ode_secir_parameter_study.cpp b/cpp/examples/ode_secir_parameter_study.cpp index 8ad8c1cd2b..aaf1ffecff 100644 --- a/cpp/examples/ode_secir_parameter_study.cpp +++ b/cpp/examples/ode_secir_parameter_study.cpp @@ -71,8 +71,8 @@ write_single_run_result(const size_t run, return node.id; }); auto num_groups = (int)(size_t)graph.nodes()[0].property.get_simulation().get_model().parameters.get_num_groups(); - BOOST_OUTCOME_TRY(mio::save_result(all_results, ids, num_groups, - mio::path_join(abs_path, ("Results_run" + std::to_string(run) + ".h5")))); + BOOST_OUTCOME_TRY(mio::save_result( + all_results, ids, num_groups, mio::path_join(abs_path, ("Results_run" + std::to_string(run) + ".h5")))); return mio::success(); } diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp index b2e8260cc3..c9df8d1781 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -322,12 +322,13 @@ int main() mio::log_info("Directory '{:s}' was created.", results_dir.string()); } - auto county_ids = std::vector{1001, 1002, 1003}; - auto save_results_status = save_results(ensemble_results, ensemble_params, county_ids, results_dir, false); - auto pairs_edges = std::vector>{}; + auto county_ids = std::vector{1001, 1002, 1003}; + auto save_results_status = + save_results(ensemble_results, ensemble_params, county_ids, results_dir, false); + auto pairs_edges = std::vector>{}; for (auto& edge : params_graph.edges()) { pairs_edges.push_back({county_ids[edge.start_node_idx], county_ids[edge.end_node_idx]}); } - auto save_edges_status = save_edges(ensemble_edges, pairs_edges, "test_results", false, true); + auto save_edges_status = save_edges(ensemble_edges, pairs_edges, "test_results", false, true); } } diff --git a/cpp/examples/ode_secir_read_graph.cpp b/cpp/examples/ode_secir_read_graph.cpp index 8832aec976..5a56fcc973 100644 --- a/cpp/examples/ode_secir_read_graph.cpp +++ b/cpp/examples/ode_secir_read_graph.cpp @@ -110,7 +110,7 @@ int main(int argc, char** argv) mio::osecir::set_params_distributions_normal(model, t0, tmax, 0.2); std::cout << "Reading Mobility File..." << std::flush; - auto read_mobility_result = mio::read_mobility_plain(filename); + auto read_mobility_result = mio::read_mobility_plain(filename); if (!read_mobility_result) { std::cout << read_mobility_result.error().formatted_message() << '\n'; std::cout << "Create the mobility file with MEmilio Epidata's getCommuterMobility.py file." << '\n'; @@ -135,7 +135,7 @@ int main(int argc, char** argv) std::cout << "Done" << std::endl; std::cout << "Writing Json Files..." << std::flush; - auto write_status = mio::write_graph(graph, "graph_parameters"); + auto write_status = mio::write_graph(graph, "graph_parameters"); if (!write_status) { std::cout << "Error: " << write_status.error().formatted_message(); } diff --git a/cpp/examples/ode_secir_save_results.cpp b/cpp/examples/ode_secir_save_results.cpp index 5c631fd949..8351c6bd35 100644 --- a/cpp/examples/ode_secir_save_results.cpp +++ b/cpp/examples/ode_secir_save_results.cpp @@ -83,5 +83,6 @@ int main() std::vector> results_from_sim = {result_from_sim, result_from_sim}; std::vector ids = {1, 2}; - auto save_result_status = mio::save_result(results_from_sim, ids, (int)(size_t)nb_groups, "test_result.h5"); + auto save_result_status = + mio::save_result(results_from_sim, ids, (int)(size_t)nb_groups, "test_result.h5"); } diff --git a/cpp/memilio/ad/ad.h b/cpp/memilio/ad/ad.h index 96c93bcd29..92cb90e98c 100644 --- a/cpp/memilio/ad/ad.h +++ b/cpp/memilio/ad/ad.h @@ -27,6 +27,74 @@ #include #include +// Extend automatic differentiation (AD) library to support std::round. +namespace ad +{ +namespace internal +{ +using std::round; +template +static inline double round(const ad::internal::active_type& x) +{ + return round(x._value()); +} +template +static inline double round(const ad::internal::binary_intermediate_aa& x) +{ + return round(x._value()); +} +template +static inline double round(const ad::internal::binary_intermediate_ap& x) +{ + return round(x._value()); +} +template +static inline double round(const ad::internal::binary_intermediate_pa& x) +{ + return round(x._value()); +} +template +static inline double round(const ad::internal::unary_intermediate& x) +{ + return round(x._value()); +} +} // namespace internal +} // namespace ad + +// Extend automatic differentiation (AD) library to support std::trunc. +namespace ad +{ +namespace internal +{ +using std::trunc; +template +static inline double trunc(const ad::internal::active_type& x) +{ + return trunc(x._value()); +} +template +static inline double trunc(const ad::internal::binary_intermediate_aa& x) +{ + return trunc(x._value()); +} +template +static inline double trunc(const ad::internal::binary_intermediate_ap& x) +{ + return trunc(x._value()); +} +template +static inline double trunc(const ad::internal::binary_intermediate_pa& x) +{ + return trunc(x._value()); +} +template +static inline double trunc(const ad::internal::unary_intermediate& x) +{ + return trunc(x._value()); +} +} // namespace internal +} // namespace ad + // Allow std::numeric_limits to work with AD types. template struct std::numeric_limits> : public numeric_limits { diff --git a/cpp/memilio/data/analyze_result.h b/cpp/memilio/data/analyze_result.h index 5e80a05773..406f097d1e 100644 --- a/cpp/memilio/data/analyze_result.h +++ b/cpp/memilio/data/analyze_result.h @@ -309,6 +309,8 @@ std::vector> ensemble_mean(const std::vector std::vector> ensemble_percentile(const std::vector>>& ensemble_result, FP p) { + using std::floor; + assert(p > 0.0 && p < 1.0 && "Invalid percentile value."); auto num_runs = ensemble_result.size(); @@ -328,7 +330,7 @@ std::vector> ensemble_percentile(const std::vector(num_runs * p)]; + percentile[node][time][elem] = single_element_ensemble[static_cast(floor(num_runs * p))]; } } } diff --git a/cpp/memilio/io/mobility_io.h b/cpp/memilio/io/mobility_io.h index 9e62ce5075..c886820a7e 100644 --- a/cpp/memilio/io/mobility_io.h +++ b/cpp/memilio/io/mobility_io.h @@ -60,13 +60,13 @@ IOResult count_lines(const std::string& filename); * where N is the number of regions * @param filename name of file to be read */ -template -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)); + return success(Eigen::MatrixX(0, 0)); } std::fstream file; @@ -75,7 +75,7 @@ IOResult read_mobility_formatted(const std::string& filename) return failure(StatusCode::FileNotFound, filename); } - Eigen::MatrixXd txt_matrix(num_lines - 1, 3); + Eigen::MatrixX txt_matrix(num_lines - 1, 3); std::vector ids; try { @@ -89,9 +89,9 @@ IOResult read_mobility_formatted(const std::string& filename) 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]); + txt_matrix(linenumber - 1, 0) = static_cast(std::stoi(line[2])); + txt_matrix(linenumber - 1, 1) = static_cast(std::stoi(line[3])); + txt_matrix(linenumber - 1, 2) = static_cast(std::stod(line[4])); } linenumber++; } @@ -104,7 +104,7 @@ IOResult read_mobility_formatted(const std::string& filename) 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()); + Eigen::MatrixX mobility = Eigen::MatrixX::Zero(ids.size(), ids.size()); for (int k = 0; k < num_lines - 1; k++) { int row_ind = 0; @@ -127,13 +127,13 @@ IOResult read_mobility_formatted(const std::string& filename) * Matrix, where N is the number of regions * @param filename name of file to be read */ -template -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)); + return success(Eigen::MatrixX(0, 0)); } std::fstream file; @@ -142,7 +142,7 @@ IOResult read_mobility_plain(const std::string& filename) return failure(StatusCode::FileNotFound, filename); } - Eigen::MatrixXd mobility(num_lines, num_lines); + Eigen::MatrixX mobility(num_lines, num_lines); try { std::string tp; @@ -154,7 +154,7 @@ IOResult read_mobility_plain(const std::string& filename) } Eigen::Index i = static_cast(linenumber); for (Eigen::Index j = 0; j < static_cast(line.size()); j++) { - mobility(i, j) = std::stod(line[j]); + mobility(i, j) = static_cast(std::stod(line[j])); } linenumber++; } @@ -302,9 +302,9 @@ 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) +template +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; @@ -337,18 +337,18 @@ IOResult save_edges(const std::vector>& results, "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()); + 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(); + 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( + 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) { @@ -412,8 +412,8 @@ 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, +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) { @@ -421,8 +421,8 @@ IOResult save_edges(const std::vector>> 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())); + BOOST_OUTCOME_TRY(save_edges(ensemble_edges[i], pairs_edges, + (result_dir / ("Edges_run" + std::to_string(i) + ".h5")).string())); } } @@ -441,17 +441,17 @@ IOResult save_edges(const std::vector>> // 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())); + 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(); diff --git a/cpp/memilio/io/parameters_io.h b/cpp/memilio/io/parameters_io.h index 23c23b729a..393ed7c845 100644 --- a/cpp/memilio/io/parameters_io.h +++ b/cpp/memilio/io/parameters_io.h @@ -63,7 +63,7 @@ int get_region_id(const EpiDataEntry& data_entry) * * @return An IOResult indicating success or failure. */ -template +template IOResult compute_divi_data(const std::vector& divi_data, const std::vector& vregion, Date date, std::vector& vnum_icu) { @@ -106,12 +106,12 @@ IOResult compute_divi_data(const std::vector& divi_data, const * * @return An IOResult indicating success or failure. */ -template +template IOResult read_divi_data(const std::string& path, const std::vector& vregion, Date date, std::vector& vnum_icu) { BOOST_OUTCOME_TRY(auto&& divi_data, mio::read_divi_data(path)); - return compute_divi_data(divi_data, vregion, date, vnum_icu); + return compute_divi_data(divi_data, vregion, date, vnum_icu); } /** @@ -122,12 +122,12 @@ 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) +template +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)); + 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 @@ -163,12 +163,11 @@ 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) +template +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); + return read_population_data(population_data, vregion); } } // namespace mio diff --git a/cpp/memilio/io/result_io.h b/cpp/memilio/io/result_io.h index 38fb060adf..07c9274fc7 100644 --- a/cpp/memilio/io/result_io.h +++ b/cpp/memilio/io/result_io.h @@ -43,9 +43,9 @@ namespace mio * @param filename Name of file * @return Any io errors that occur during writing of the files. */ -template -IOResult save_result(const std::vector>& results, 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)}; @@ -65,17 +65,17 @@ IOResult save_result(const std::vector>& results, c 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()); + 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) + 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) + 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) { @@ -103,6 +103,7 @@ IOResult save_result(const std::vector>& results, c return success(); } +template class SimulationResult { public: @@ -122,7 +123,7 @@ class SimulationResult * @param groups Simulation results of individual groups. * @param total Simulation results as the sum over all groups. */ - SimulationResult(const TimeSeries& groups, const TimeSeries& totals) + SimulationResult(const TimeSeries& groups, const TimeSeries& totals) : m_groups(groups) , m_totals(totals) { @@ -131,7 +132,7 @@ class SimulationResult /** * @brief Simulation results of individual groups. */ - const TimeSeries& get_groups() const + const TimeSeries& get_groups() const { return m_groups; } @@ -139,14 +140,14 @@ class SimulationResult /** * @brief Simulation results of the sum over all groups. */ - const TimeSeries& get_totals() const + const TimeSeries& get_totals() const { return m_totals; } private: - TimeSeries m_groups; - TimeSeries m_totals; + TimeSeries m_groups; + TimeSeries m_totals; }; // Forward declaration of store_group_name @@ -155,10 +156,10 @@ herr_t store_group_name(hid_t /*id*/, const char* name, const H5L_info_t* /*linf * @brief Read simulation result from h5 file. * @param filename name of the file to be read. */ -template -IOResult> read_result(const std::string& filename) +template +IOResult>> read_result(const std::string& filename) { - std::vector results; + std::vector> results; H5File file{H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)}; MEMILIO_H5_CHECK(file.id, StatusCode::FileNotFound, filename); @@ -199,7 +200,7 @@ IOResult> read_result(const std::string& filename) H5Sget_simple_extent_dims(dataspace_t.id, dims_t, NULL); auto num_timepoints = Eigen::Index(dims_t[0]); - auto time = std::vector(num_timepoints); + 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."); @@ -218,19 +219,19 @@ IOResult> read_result(const std::string& filename) } auto num_infectionstates = Eigen::Index(dims_total[1]); - auto total_values = Eigen::Matrix( - num_timepoints, num_infectionstates); + 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); + 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); + 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]); @@ -262,8 +263,8 @@ IOResult> read_result(const std::string& filename) return failure(StatusCode::InvalidFileFormat, "Number of infection states does not match."); } - auto group_values = Eigen::Matrix( - num_timepoints, num_infectionstates); + 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"); @@ -275,7 +276,7 @@ IOResult> read_result(const std::string& filename) } } - results.push_back(SimulationResult(groups, totals)); + results.push_back(SimulationResult(groups, totals)); } return success(results); } @@ -290,17 +291,16 @@ IOResult> read_result(const std::string& filename) * @param run_idx Index of the run; used in directory name. * @return Any io errors that occur during writing of the files. */ -template -IOResult save_result_with_params(const std::vector>& result, - const std::vector& params, const std::vector& county_ids, - const fs::path& result_dir, size_t run_idx) +template +IOResult save_result_with_params(const std::vector>& result, const std::vector& params, + const std::vector& county_ids, const fs::path& result_dir, size_t run_idx) { auto result_dir_run = result_dir / ("run" + std::to_string(run_idx)); BOOST_OUTCOME_TRY(create_directory(result_dir_run.string())); - BOOST_OUTCOME_TRY(save_result(result, county_ids, (int)(size_t)params[0].parameters.get_num_groups(), - (result_dir_run / "Result.h5").string())); - BOOST_OUTCOME_TRY(write_graph(create_graph_without_edges>(params, county_ids), - result_dir_run.string(), IOF_OmitDistributions)); + BOOST_OUTCOME_TRY(save_result(result, county_ids, (int)(size_t)params[0].parameters.get_num_groups(), + (result_dir_run / "Result.h5").string())); + BOOST_OUTCOME_TRY(write_graph(create_graph_without_edges>(params, county_ids), + result_dir_run.string(), IOF_OmitDistributions)); return success(); } @@ -315,8 +315,8 @@ IOResult save_result_with_params(const std::vector> * @param save_single_runs [Default: true] Defines if percentiles are written to the disk. * @return Any io errors that occur during writing of the files. */ -template -IOResult save_results(const std::vector>>& ensemble_results, +template +IOResult save_results(const std::vector>>& ensemble_results, const std::vector>& ensemble_params, const std::vector& county_ids, const fs::path& result_dir, bool save_single_runs = true, bool save_percentiles = true) { @@ -325,10 +325,10 @@ IOResult save_results(const std::vector auto num_groups = (int)(size_t)ensemble_params[0][0].parameters.get_num_groups(); if (save_single_runs) { for (size_t i = 0; i < ensemble_result_sum.size(); ++i) { - BOOST_OUTCOME_TRY(save_result(ensemble_result_sum[i], {0}, num_groups, - (result_dir / ("results_run" + std::to_string(i) + "_sum.h5")).string())); - BOOST_OUTCOME_TRY(save_result(ensemble_results[i], county_ids, num_groups, - (result_dir / ("results_run" + std::to_string(i) + ".h5")).string())); + BOOST_OUTCOME_TRY(save_result(ensemble_result_sum[i], {0}, num_groups, + (result_dir / ("results_run" + std::to_string(i) + "_sum.h5")).string())); + BOOST_OUTCOME_TRY(save_result(ensemble_results[i], county_ids, num_groups, + (result_dir / ("results_run" + std::to_string(i) + ".h5")).string())); } } @@ -347,65 +347,65 @@ IOResult save_results(const std::vector // save percentiles of results, summed over nodes { - auto ensemble_results_sum_p05 = ensemble_percentile(ensemble_result_sum, 0.05); - auto ensemble_results_sum_p25 = ensemble_percentile(ensemble_result_sum, 0.25); - auto ensemble_results_sum_p50 = ensemble_percentile(ensemble_result_sum, 0.50); - auto ensemble_results_sum_p75 = ensemble_percentile(ensemble_result_sum, 0.75); - auto ensemble_results_sum_p95 = ensemble_percentile(ensemble_result_sum, 0.95); - - BOOST_OUTCOME_TRY( - save_result(ensemble_results_sum_p05, {0}, num_groups, (result_dir_p05 / "Results_sum.h5").string())); - BOOST_OUTCOME_TRY( - save_result(ensemble_results_sum_p25, {0}, num_groups, (result_dir_p25 / "Results_sum.h5").string())); - BOOST_OUTCOME_TRY( - save_result(ensemble_results_sum_p50, {0}, num_groups, (result_dir_p50 / "Results_sum.h5").string())); - BOOST_OUTCOME_TRY( - save_result(ensemble_results_sum_p75, {0}, num_groups, (result_dir_p75 / "Results_sum.h5").string())); - BOOST_OUTCOME_TRY( - save_result(ensemble_results_sum_p95, {0}, num_groups, (result_dir_p95 / "Results_sum.h5").string())); + auto ensemble_results_sum_p05 = ensemble_percentile(ensemble_result_sum, 0.05); + auto ensemble_results_sum_p25 = ensemble_percentile(ensemble_result_sum, 0.25); + auto ensemble_results_sum_p50 = ensemble_percentile(ensemble_result_sum, 0.50); + auto ensemble_results_sum_p75 = ensemble_percentile(ensemble_result_sum, 0.75); + auto ensemble_results_sum_p95 = ensemble_percentile(ensemble_result_sum, 0.95); + + BOOST_OUTCOME_TRY(save_result(ensemble_results_sum_p05, {0}, num_groups, + (result_dir_p05 / "Results_sum.h5").string())); + BOOST_OUTCOME_TRY(save_result(ensemble_results_sum_p25, {0}, num_groups, + (result_dir_p25 / "Results_sum.h5").string())); + BOOST_OUTCOME_TRY(save_result(ensemble_results_sum_p50, {0}, num_groups, + (result_dir_p50 / "Results_sum.h5").string())); + BOOST_OUTCOME_TRY(save_result(ensemble_results_sum_p75, {0}, num_groups, + (result_dir_p75 / "Results_sum.h5").string())); + BOOST_OUTCOME_TRY(save_result(ensemble_results_sum_p95, {0}, num_groups, + (result_dir_p95 / "Results_sum.h5").string())); } // save percentiles of results { - auto ensemble_results_p05 = ensemble_percentile(ensemble_results, 0.05); - auto ensemble_results_p25 = ensemble_percentile(ensemble_results, 0.25); - auto ensemble_results_p50 = ensemble_percentile(ensemble_results, 0.50); - auto ensemble_results_p75 = ensemble_percentile(ensemble_results, 0.75); - auto ensemble_results_p95 = ensemble_percentile(ensemble_results, 0.95); - - BOOST_OUTCOME_TRY( - save_result(ensemble_results_p05, county_ids, num_groups, (result_dir_p05 / "Results.h5").string())); - BOOST_OUTCOME_TRY( - save_result(ensemble_results_p25, county_ids, num_groups, (result_dir_p25 / "Results.h5").string())); - BOOST_OUTCOME_TRY( - save_result(ensemble_results_p50, county_ids, num_groups, (result_dir_p50 / "Results.h5").string())); - BOOST_OUTCOME_TRY( - save_result(ensemble_results_p75, county_ids, num_groups, (result_dir_p75 / "Results.h5").string())); - BOOST_OUTCOME_TRY( - save_result(ensemble_results_p95, county_ids, num_groups, (result_dir_p95 / "Results.h5").string())); + auto ensemble_results_p05 = ensemble_percentile(ensemble_results, 0.05); + auto ensemble_results_p25 = ensemble_percentile(ensemble_results, 0.25); + auto ensemble_results_p50 = ensemble_percentile(ensemble_results, 0.50); + auto ensemble_results_p75 = ensemble_percentile(ensemble_results, 0.75); + auto ensemble_results_p95 = ensemble_percentile(ensemble_results, 0.95); + + BOOST_OUTCOME_TRY(save_result(ensemble_results_p05, county_ids, num_groups, + (result_dir_p05 / "Results.h5").string())); + BOOST_OUTCOME_TRY(save_result(ensemble_results_p25, county_ids, num_groups, + (result_dir_p25 / "Results.h5").string())); + BOOST_OUTCOME_TRY(save_result(ensemble_results_p50, county_ids, num_groups, + (result_dir_p50 / "Results.h5").string())); + BOOST_OUTCOME_TRY(save_result(ensemble_results_p75, county_ids, num_groups, + (result_dir_p75 / "Results.h5").string())); + BOOST_OUTCOME_TRY(save_result(ensemble_results_p95, county_ids, num_groups, + (result_dir_p95 / "Results.h5").string())); } // save percentiles of parameters { - auto ensemble_params_p05 = ensemble_params_percentile(ensemble_params, 0.05); - auto ensemble_params_p25 = ensemble_params_percentile(ensemble_params, 0.25); - auto ensemble_params_p50 = ensemble_params_percentile(ensemble_params, 0.50); - auto ensemble_params_p75 = ensemble_params_percentile(ensemble_params, 0.75); - auto ensemble_params_p95 = ensemble_params_percentile(ensemble_params, 0.95); + auto ensemble_params_p05 = ensemble_params_percentile(ensemble_params, 0.05); + auto ensemble_params_p25 = ensemble_params_percentile(ensemble_params, 0.25); + auto ensemble_params_p50 = ensemble_params_percentile(ensemble_params, 0.50); + auto ensemble_params_p75 = ensemble_params_percentile(ensemble_params, 0.75); + auto ensemble_params_p95 = ensemble_params_percentile(ensemble_params, 0.95); auto make_graph = [&county_ids](auto&& params) { - return create_graph_without_edges>(params, county_ids); + return create_graph_without_edges>(params, county_ids); }; BOOST_OUTCOME_TRY( - write_graph(make_graph(ensemble_params_p05), result_dir_p05.string(), IOF_OmitDistributions)); + write_graph(make_graph(ensemble_params_p05), result_dir_p05.string(), IOF_OmitDistributions)); BOOST_OUTCOME_TRY( - write_graph(make_graph(ensemble_params_p25), result_dir_p25.string(), IOF_OmitDistributions)); + write_graph(make_graph(ensemble_params_p25), result_dir_p25.string(), IOF_OmitDistributions)); BOOST_OUTCOME_TRY( - write_graph(make_graph(ensemble_params_p50), result_dir_p50.string(), IOF_OmitDistributions)); + write_graph(make_graph(ensemble_params_p50), result_dir_p50.string(), IOF_OmitDistributions)); BOOST_OUTCOME_TRY( - write_graph(make_graph(ensemble_params_p75), result_dir_p75.string(), IOF_OmitDistributions)); + write_graph(make_graph(ensemble_params_p75), result_dir_p75.string(), IOF_OmitDistributions)); BOOST_OUTCOME_TRY( - write_graph(make_graph(ensemble_params_p95), result_dir_p95.string(), IOF_OmitDistributions)); + write_graph(make_graph(ensemble_params_p95), result_dir_p95.string(), IOF_OmitDistributions)); } } return success(); diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index 26a8bb5a6a..d736af0f25 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -281,12 +281,13 @@ IOResult set_nodes(const Parameters& params, Date start_date, Date end_dat for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { - auto tnt_capacity = nodes[node_idx].populations.get_total() * tnt_capacity_factor; + FP tnt_capacity = nodes[node_idx].populations.get_total() * tnt_capacity_factor; //local parameters auto& tnt_value = nodes[node_idx].parameters.template get(); tnt_value = UncertainValue(0.5 * (1.2 * tnt_capacity + 0.8 * tnt_capacity)); - tnt_value.set_distribution(mio::ParameterDistributionUniform(0.8 * tnt_capacity, 1.2 * tnt_capacity)); + tnt_value.set_distribution( + mio::ParameterDistributionUniform(0.8 * ad::value(tnt_capacity), 1.2 * ad::value(tnt_capacity))); auto id = 0; if (is_node_for_county) { @@ -312,8 +313,8 @@ IOResult set_nodes(const Parameters& params, Date start_date, Date end_dat auto& compartment_value = nodes[node_idx].populations[{i, j}]; compartment_value = UncertainValue(0.5 * (1.1 * compartment_value.value() + 0.9 * compartment_value.value())); - compartment_value.set_distribution(mio::ParameterDistributionUniform(0.9 * compartment_value.value(), - 1.1 * compartment_value.value())); + compartment_value.set_distribution(mio::ParameterDistributionUniform( + 0.9 * ad::value(compartment_value.value()), 1.1 * ad::value(compartment_value.value()))); } } @@ -359,7 +360,7 @@ IOResult set_edges(const fs::path& mobility_data_file, Graph(num_age_groups, 1.0) : commuting_weights); //commuters - auto working_population = 0.0; + FP working_population = 0.0; for (auto age = AgeGroup(0); age < populations.template size(); ++age) { working_population += populations.get_group_total(age) * commuting_weights[size_t(age)]; } diff --git a/cpp/models/ode_secir/parameters_io.h b/cpp/models/ode_secir/parameters_io.h index f288abcc35..fc71756de9 100644 --- a/cpp/models/ode_secir/parameters_io.h +++ b/cpp/models/ode_secir/parameters_io.h @@ -52,18 +52,18 @@ int get_region_id(int id); * @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 +template 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, + 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) + 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; @@ -77,7 +77,7 @@ IOResult read_confirmed_cases_data( 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); + int 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 @@ -177,7 +177,7 @@ IOResult read_confirmed_cases_data( (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) { + auto try_fix_constraints = [region, i](FP& value, FP 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 @@ -217,14 +217,14 @@ IOResult read_confirmed_cases_data( template IOResult set_confirmed_cases_data(std::vector>& model, std::vector& case_data, const std::vector& region, Date date, - const std::vector& scaling_factor_inf) + const std::vector& scaling_factor_inf) { const size_t num_age_groups = ConfirmedCasesDataEntry::age_group_names.size(); // allow single scalar scaling that is broadcast to all age groups assert(scaling_factor_inf.size() == 1 || scaling_factor_inf.size() == num_age_groups); // Set scaling factors to match num age groups - std::vector scaling_factor_inf_full; + std::vector scaling_factor_inf_full; if (scaling_factor_inf.size() == 1) { scaling_factor_inf_full.assign(num_age_groups, scaling_factor_inf[0]); } @@ -238,10 +238,10 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect std::vector> t_InfectedSevere{model.size()}; std::vector> t_InfectedCritical{model.size()}; - std::vector> mu_C_R{model.size()}; - std::vector> mu_I_H{model.size()}; - std::vector> mu_H_U{model.size()}; - std::vector> mu_U_D{model.size()}; + std::vector> mu_C_R{model.size()}; + std::vector> mu_I_H{model.size()}; + std::vector> mu_H_U{model.size()}; + std::vector> mu_U_D{model.size()}; for (size_t node = 0; node < model.size(); ++node) { const size_t model_groups = (size_t)model[node].parameters.get_num_groups(); @@ -270,18 +270,18 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect mu_U_D[node].push_back(model[node].parameters.template get>()[(AgeGroup)group]); } } - std::vector> num_InfectedSymptoms(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_death(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_rec(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_Exposed(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedNoSymptoms(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedSevere(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_icu(model.size(), std::vector(num_age_groups, 0.0)); - - BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms, - num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec, - t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, - t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf_full)); + std::vector> num_InfectedSymptoms(model.size(), std::vector(num_age_groups, 0.0)); + std::vector> num_death(model.size(), std::vector(num_age_groups, 0.0)); + std::vector> num_rec(model.size(), std::vector(num_age_groups, 0.0)); + std::vector> num_Exposed(model.size(), std::vector(num_age_groups, 0.0)); + std::vector> num_InfectedNoSymptoms(model.size(), std::vector(num_age_groups, 0.0)); + std::vector> num_InfectedSevere(model.size(), std::vector(num_age_groups, 0.0)); + std::vector> num_icu(model.size(), std::vector(num_age_groups, 0.0)); + + BOOST_OUTCOME_TRY(read_confirmed_cases_data( + case_data, region, date, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, + num_death, num_rec, t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, t_InfectedCritical, + mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf_full)); for (size_t node = 0; node < model.size(); node++) { @@ -353,25 +353,25 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect template IOResult set_confirmed_cases_data(std::vector>& model, const std::string& path, std::vector const& region, Date date, - const std::vector& scaling_factor_inf) + const std::vector& scaling_factor_inf) { BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path)); - BOOST_OUTCOME_TRY(set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf)); + BOOST_OUTCOME_TRY(set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf)); return success(); } /** - * @brief Sets populations data from DIVI register into Model. - * @tparam FP floating point data type, e.g., double. - * @param[in, out] model Vector of models in which the data is set. - * @param[in] path Path to 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. + * @tparam FP floating point data type, e.g., double. + * @param[in, out] model Vector of models in which the data is set. + * @param[in] path Path to 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) + Date date, FP scaling_factor_icu) { // DIVI dataset will no longer be updated from CW29 2024 on. if (!is_divi_data_available(date)) { @@ -380,8 +380,8 @@ IOResult set_divi_data(std::vector>& model, const std::string& p date); return success(); } - std::vector sum_mu_I_U(vregion.size(), 0); - std::vector> mu_I_U{model.size()}; + std::vector sum_mu_I_U(vregion.size(), 0); + std::vector> mu_I_U{model.size()}; for (size_t region = 0; region < vregion.size(); region++) { auto num_groups = model[region].parameters.get_num_groups(); for (auto i = AgeGroup(0); i < num_groups; i++) { @@ -391,8 +391,8 @@ IOResult set_divi_data(std::vector>& model, const std::string& p model[region].parameters.template get>()[i]); } } - std::vector num_icu(model.size(), 0.0); - BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu)); + std::vector num_icu(model.size(), 0.0); + BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu)); for (size_t region = 0; region < vregion.size(); region++) { auto num_groups = model[region].parameters.get_num_groups(); @@ -452,8 +452,8 @@ template IOResult set_population_data(std::vector>& model, const std::string& path, const std::vector& vregion) { - BOOST_OUTCOME_TRY(const auto&& num_population, read_population_data(path, vregion)); - BOOST_OUTCOME_TRY(set_population_data(model, num_population, vregion)); + BOOST_OUTCOME_TRY(const auto&& num_population, read_population_data(path, vregion)); + BOOST_OUTCOME_TRY(set_population_data(model, num_population, vregion)); return success(); } @@ -477,48 +477,53 @@ IOResult set_population_data(std::vector>& model, const std::str * @param[in] confirmed_cases_path Path to confirmed cases file. * @param[in] population_data_path Path to population data file. */ -template -IOResult export_input_data_county_timeseries( - std::vector models, const std::string& results_dir, std::vector const& region, Date date, - const std::vector& scaling_factor_inf, double scaling_factor_icu, int num_days, - const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path) +template +IOResult export_input_data_county_timeseries(std::vector models, const std::string& results_dir, + std::vector const& region, Date date, + const std::vector& scaling_factor_inf, FP scaling_factor_icu, + int num_days, const std::string& divi_data_path, + const std::string& confirmed_cases_path, + const std::string& population_data_path) { const auto num_age_groups = (size_t)models[0].parameters.get_num_groups(); // allow scalar scaling factor as convenience for 1-group models assert(scaling_factor_inf.size() == 1 || scaling_factor_inf.size() == num_age_groups); assert(models.size() == region.size()); - std::vector> extrapolated_data( - region.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_age_groups)); + std::vector> extrapolated_data( + region.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_age_groups)); - BOOST_OUTCOME_TRY(auto&& num_population, mio::read_population_data(population_data_path, region)); + BOOST_OUTCOME_TRY(auto&& num_population, mio::read_population_data(population_data_path, region)); BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(confirmed_cases_path)); for (int t = 0; t <= num_days; ++t) { auto offset_day = offset_date_by_days(date, t); - BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, region, offset_day, scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(models, case_data, region, offset_day, scaling_factor_inf)); - BOOST_OUTCOME_TRY(details::set_population_data(models, num_population, region)); + BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, region, offset_day, scaling_factor_icu)); + BOOST_OUTCOME_TRY( + details::set_confirmed_cases_data(models, case_data, region, offset_day, scaling_factor_inf)); + BOOST_OUTCOME_TRY(details::set_population_data(models, num_population, region)); for (size_t r = 0; r < region.size(); r++) { extrapolated_data[r][t] = models[r].get_initial_values(); } } BOOST_OUTCOME_TRY( - save_result(extrapolated_data, region, (int)num_age_groups, path_join(results_dir, "Results_rki.h5"))); + save_result(extrapolated_data, region, (int)num_age_groups, path_join(results_dir, "Results_rki.h5"))); - auto rki_data_sum = sum_nodes(std::vector>>{extrapolated_data}); + auto rki_data_sum = sum_nodes(std::vector>>{extrapolated_data}); BOOST_OUTCOME_TRY( - save_result({rki_data_sum[0][0]}, {0}, (int)num_age_groups, path_join(results_dir, "Results_rki_sum.h5"))); + save_result({rki_data_sum[0][0]}, {0}, (int)num_age_groups, path_join(results_dir, "Results_rki_sum.h5"))); return success(); } #else -template -IOResult export_input_data_county_timeseries( - std::vector models, const std::string& results_dir, std::vector const& region, Date date, - const std::vector& scaling_factor_inf, double scaling_factor_icu, int num_days, - const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path) +template +IOResult export_input_data_county_timeseries(std::vector models, const std::string& results_dir, + std::vector const& region, Date date, + const std::vector& scaling_factor_inf, FP scaling_factor_icu, + int num_days, const std::string& divi_data_path, + const std::string& confirmed_cases_path, + const std::string& population_data_path) { mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data."); return success(); @@ -533,17 +538,16 @@ IOResult export_input_data_county_timeseries( * @param[in] scaling_factor_icu Factor by which to scale the icu cases of divi data. * @param[in] pydata_dir Directory of files. */ -template -IOResult read_input_data_germany(std::vector& model, Date date, - const std::vector& scaling_factor_inf, double scaling_factor_icu, - const std::string& pydata_dir) +template +IOResult read_input_data_germany(std::vector& model, Date date, const std::vector& scaling_factor_inf, + FP scaling_factor_icu, const std::string& pydata_dir) { BOOST_OUTCOME_TRY( - details::set_divi_data(model, path_join(pydata_dir, "germany_divi.json"), {0}, date, scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_age_ma7.json"), {0}, - date, scaling_factor_inf)); + details::set_divi_data(model, path_join(pydata_dir, "germany_divi.json"), {0}, date, scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_age_ma7.json"), {0}, + date, scaling_factor_inf)); BOOST_OUTCOME_TRY( - details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), {0})); + details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), {0})); return success(); } @@ -556,18 +560,18 @@ IOResult read_input_data_germany(std::vector& model, Date date, * @param[in] scaling_factor_icu Factor by which to scale the icu cases of divi data. * @param[in] pydata_dir Directory of files. */ -template +template IOResult read_input_data_state(std::vector& model, Date date, std::vector& state, - const std::vector& scaling_factor_inf, double scaling_factor_icu, + const std::vector& scaling_factor_inf, FP scaling_factor_icu, const std::string& pydata_dir) { BOOST_OUTCOME_TRY( - details::set_divi_data(model, path_join(pydata_dir, "state_divi.json"), state, date, scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_state_age_ma7.json"), - state, date, scaling_factor_inf)); + details::set_divi_data(model, path_join(pydata_dir, "state_divi.json"), state, date, scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data( + model, path_join(pydata_dir, "cases_all_state_age_ma7.json"), state, date, scaling_factor_inf)); BOOST_OUTCOME_TRY( - details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), state)); + details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), state)); return success(); } @@ -582,17 +586,17 @@ IOResult read_input_data_state(std::vector& model, Date date, std:: * @param[in] num_days [Default: 0] Number of days to be simulated; required to extrapolate real data. * @param[in] export_time_series [Default: false] If true, reads data for each day of simulation and writes it in the same directory as the input files. */ -template +template IOResult read_input_data_county(std::vector& model, Date date, const std::vector& county, - const std::vector& scaling_factor_inf, double scaling_factor_icu, + const std::vector& scaling_factor_inf, FP scaling_factor_icu, const std::string& pydata_dir, int num_days = 0, bool export_time_series = false) { + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "county_divi_ma7.json"), county, date, + scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data( + model, path_join(pydata_dir, "cases_all_county_age_ma7.json"), county, date, scaling_factor_inf)); BOOST_OUTCOME_TRY( - details::set_divi_data(model, path_join(pydata_dir, "county_divi_ma7.json"), county, date, scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_county_age_ma7.json"), - county, date, scaling_factor_inf)); - BOOST_OUTCOME_TRY( - details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), county)); + details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), county)); if (export_time_series) { // Use only if extrapolated real data is needed for comparison. EXPENSIVE ! @@ -600,7 +604,7 @@ IOResult read_input_data_county(std::vector& model, Date date, cons // (This only represents the vectorization of the previous function over all simulation days...) log_warning("Exporting time series of extrapolated real data. This may take some minutes. " "For simulation runs over the same time period, deactivate it."); - BOOST_OUTCOME_TRY(export_input_data_county_timeseries( + BOOST_OUTCOME_TRY(export_input_data_county_timeseries( model, pydata_dir, county, date, scaling_factor_inf, scaling_factor_icu, num_days, path_join(pydata_dir, "county_divi_ma7.json"), path_join(pydata_dir, "cases_all_county_age_ma7.json"), path_join(pydata_dir, "county_current_population.json"))); @@ -618,16 +622,16 @@ IOResult read_input_data_county(std::vector& model, Date date, cons * @param[in] pydata_dir directory of files * @param[in] age_group_names strings specifying age group names */ -template +template IOResult read_input_data(std::vector& model, Date date, const std::vector& node_ids, - const std::vector& scaling_factor_inf, double scaling_factor_icu, + const std::vector& scaling_factor_inf, FP scaling_factor_icu, const std::string& pydata_dir, int num_days = 0, bool export_time_series = false) { - BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "critical_cases.json"), node_ids, date, - scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "confirmed_cases.json"), node_ids, - date, scaling_factor_inf)); - BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "population_data.json"), node_ids)); + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "critical_cases.json"), node_ids, date, + scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "confirmed_cases.json"), + node_ids, date, scaling_factor_inf)); + BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "population_data.json"), node_ids)); if (export_time_series) { // Use only if extrapolated real data is needed for comparison. EXPENSIVE ! @@ -635,7 +639,7 @@ IOResult read_input_data(std::vector& model, Date date, const std:: // (This only represents the vectorization of the previous function over all simulation days...) log_warning("Exporting time series of extrapolated real data. This may take some minutes. " "For simulation runs over the same time period, deactivate it."); - BOOST_OUTCOME_TRY(export_input_data_county_timeseries( + BOOST_OUTCOME_TRY(export_input_data_county_timeseries( model, pydata_dir, node_ids, date, scaling_factor_inf, scaling_factor_icu, num_days, path_join(pydata_dir, "critical_cases.json"), path_join(pydata_dir, "confirmed_cases.json"), path_join(pydata_dir, "population_data.json"))); diff --git a/cpp/models/ode_secirts/analyze_result.h b/cpp/models/ode_secirts/analyze_result.h index 857658ffa3..2f885b80e9 100644 --- a/cpp/models/ode_secirts/analyze_result.h +++ b/cpp/models/ode_secirts/analyze_result.h @@ -40,6 +40,8 @@ namespace osecirts template std::vector ensemble_params_percentile(const std::vector>& ensemble_params, FP p) { + using std::floor; + assert(p > 0.0 && p < 1.0 && "Invalid percentile value."); auto num_runs = ensemble_params.size(); @@ -61,7 +63,7 @@ std::vector ensemble_params_percentile(const std::vector(num_runs * p)]; + new_params = single_element[static_cast(floor(num_runs * p))]; }; for (size_t node = 0; node < num_nodes; node++) { @@ -184,7 +186,7 @@ std::vector ensemble_params_percentile(const std::vector>( - single_element_ensemble[static_cast(num_runs * p)]); + single_element_ensemble[static_cast(floor(num_runs * p))]); } return percentile; } diff --git a/cpp/models/ode_secirts/parameters_io.h b/cpp/models/ode_secirts/parameters_io.h index c86263e717..5af585c396 100644 --- a/cpp/models/ode_secirts/parameters_io.h +++ b/cpp/models/ode_secirts/parameters_io.h @@ -79,6 +79,8 @@ IOResult compute_confirmed_cases_data( std::vector const& vregion, Date date, const std::vector& model, const std::vector& scaling_factor_inf, const size_t layer) { + using std::round; + auto max_date_entry = std::max_element(case_data.begin(), case_data.end(), [](auto&& a, auto&& b) { return a.date < b.date; }); @@ -119,17 +121,17 @@ IOResult compute_confirmed_cases_data( auto age = (size_t)entry.age_group; // (rounded) transition times const int t_exposed = - static_cast(std::round(params_region.template get>()[entry.age_group])); + static_cast(round(params_region.template get>()[entry.age_group])); int t_InfectedNoSymptoms = - static_cast(std::round(params_region.template get>()[entry.age_group])); + static_cast(round(params_region.template get>()[entry.age_group])); int t_InfectedSymptoms = - static_cast(std::round(params_region.template get>()[entry.age_group])); + static_cast(round(params_region.template get>()[entry.age_group])); const int t_InfectedSevere = - static_cast(std::round(params_region.template get>()[entry.age_group])); + static_cast(round(params_region.template get>()[entry.age_group])); const int t_InfectedCritical = - static_cast(std::round(params_region.template get>()[entry.age_group])); - const int t_imm_interval_i = static_cast( - std::round(params_region.template get>()[entry.age_group])); + static_cast(round(params_region.template get>()[entry.age_group])); + const int t_imm_interval_i = + static_cast(round(params_region.template get>()[entry.age_group])); // transition probabilities FP recoveredPerInfectedNoSymptoms = @@ -139,9 +141,9 @@ IOResult compute_confirmed_cases_data( // if we select a layer with better immunity (layer > 0), we need to adjust the times and transition rates if (layer > 0) { - t_InfectedNoSymptoms = static_cast(std::round( + t_InfectedNoSymptoms = static_cast(round( t_InfectedNoSymptoms * params_region.template get>()[entry.age_group])); - t_InfectedSymptoms = static_cast(std::round( + t_InfectedSymptoms = static_cast(round( t_InfectedSymptoms * params_region.template get>()[entry.age_group])); const FP red_fact_exp = @@ -281,9 +283,9 @@ IOResult read_confirmed_cases_data( const std::vector& model, const std::vector& scaling_factor_inf, const size_t layer) { BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path)); - return compute_confirmed_cases_data(case_data, vnum_Exposed, vnum_InfectedNoSymptoms, vnum_InfectedSymptoms, - vnum_InfectedSevere, vnum_icu, vnum_death, vnum_timm_i, vregion, date, model, - scaling_factor_inf, layer); + return compute_confirmed_cases_data(case_data, vnum_Exposed, vnum_InfectedNoSymptoms, vnum_InfectedSymptoms, + vnum_InfectedSevere, vnum_icu, vnum_death, vnum_timm_i, vregion, date, + model, scaling_factor_inf, layer); } /** @@ -374,9 +376,9 @@ set_confirmed_cases_data(std::vector& model, const std::vector(case_data, num_Exposed, num_InfectedNoSymptoms, + num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, + num_timm1, region, date, model, scaling_factor_inf, 0)); for (size_t county = 0; county < model.size(); county++) { size_t num_groups = (size_t)model[county].parameters.get_num_groups(); @@ -417,9 +419,9 @@ set_confirmed_cases_data(std::vector& model, const std::vector(num_age_groups, 0.0); } - BOOST_OUTCOME_TRY(compute_confirmed_cases_data(case_data, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, - num_InfectedSevere, num_icu, num_death, num_timm1, region, date, - model, scaling_factor_inf, 1)); + BOOST_OUTCOME_TRY(compute_confirmed_cases_data(case_data, num_Exposed, num_InfectedNoSymptoms, + num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, + num_timm1, region, date, model, scaling_factor_inf, 1)); for (size_t county = 0; county < model.size(); county++) { size_t num_groups = (size_t)model[county].parameters.get_num_groups(); @@ -478,9 +480,9 @@ set_confirmed_cases_data(std::vector& model, const std::vector(num_age_groups, 0.0); } - BOOST_OUTCOME_TRY(compute_confirmed_cases_data(case_data, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, - num_InfectedSevere, num_icu, num_death, num_timm2, region, date, - model, scaling_factor_inf, 2)); + BOOST_OUTCOME_TRY(compute_confirmed_cases_data(case_data, num_Exposed, num_InfectedNoSymptoms, + num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, + num_timm2, region, date, model, scaling_factor_inf, 2)); for (size_t county = 0; county < model.size(); county++) { size_t num_groups = (size_t)model[county].parameters.get_num_groups(); @@ -558,7 +560,7 @@ IOResult set_confirmed_cases_data(std::vector& model, const std::st { BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path)); BOOST_OUTCOME_TRY( - set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf, immunity_population)); + set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf, immunity_population)); return success(); } @@ -603,7 +605,7 @@ IOResult set_divi_data(std::vector& model, const std::string& path, } } std::vector num_icu(model.size(), 0.0); - BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu)); + BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu)); for (size_t region = 0; region < vregion.size(); region++) { auto num_groups = model[region].parameters.get_num_groups(); @@ -635,6 +637,9 @@ IOResult set_population_data(std::vector& model, const std::vector< const std::vector& vregion, const std::vector> immunity_population) { + using std::max; + using std::min; + for (size_t region = 0; region < vregion.size(); region++) { if (std::accumulate(num_population[region].begin(), num_population[region].end(), FP(0.0), [](const FP& a, const FP& b) { @@ -647,7 +652,7 @@ IOResult set_population_data(std::vector& model, const std::vector< FP SPI = num_population[region][size_t(i)] * immunity_population[1][size_t(i)]; FP SII = num_population[region][size_t(i)] - SN - SPI; - model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] = std::max( + model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] = max( 0.0, FP(SII - (model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] + @@ -660,7 +665,7 @@ IOResult set_population_data(std::vector& model, const std::vector< model[region].populations[{i, InfectionState::DeadImprovedImmunity}] + model[region].populations[{i, InfectionState::TemporaryImmuneImprovedImmunity}]))); - model[region].populations[{i, InfectionState::SusceptiblePartialImmunity}] = std::max( + model[region].populations[{i, InfectionState::SusceptiblePartialImmunity}] = max( 0.0, SPI - model[region].populations[{i, InfectionState::ExposedPartialImmunity}] - model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunity}] - @@ -707,12 +712,12 @@ IOResult set_population_data(std::vector& model, const std::vector< * * @return An IOResult indicating success or failure. */ -template +template IOResult set_population_data(std::vector& model, const std::string& path, const std::vector& vregion, - const std::vector> immunity_population) + const std::vector> immunity_population) { - BOOST_OUTCOME_TRY(auto&& num_population, mio::read_population_data(path, vregion)); - BOOST_OUTCOME_TRY(set_population_data(model, num_population, vregion, immunity_population)); + BOOST_OUTCOME_TRY(auto&& num_population, mio::read_population_data(path, vregion)); + BOOST_OUTCOME_TRY(set_population_data(model, num_population, vregion, immunity_population)); return success(); } @@ -734,14 +739,16 @@ template IOResult set_vaccination_data(std::vector>& model, const std::vector& vacc_data, Date date, const std::vector& vregion, int num_days) { + using std::floor; + auto num_groups = model[0].parameters.get_num_groups(); auto days_until_effective_n = - (int)(double)model[0].parameters.template get>()[AgeGroup(0)]; + (int)(floor(model[0].parameters.template get>()[AgeGroup(0)])); auto days_until_effective_pi = - (int)(double)model[0].parameters.template get>()[AgeGroup(0)]; + (int)(floor(model[0].parameters.template get>()[AgeGroup(0)])); auto days_until_effective_ii = - (int)(double)model[0].parameters.template get>()[AgeGroup(0)]; + (int)(floor(model[0].parameters.template get>()[AgeGroup(0)])); // iterate over regions (e.g., counties) for (size_t i = 0; i < model.size(); ++i) { // iterate over age groups in region @@ -881,7 +888,7 @@ IOResult set_vaccination_data(std::vector>& model, const std::st return success(); } BOOST_OUTCOME_TRY(auto&& vacc_data, read_vaccination_data(path)); - BOOST_OUTCOME_TRY(set_vaccination_data(model, vacc_data, date, vregion, num_days)); + BOOST_OUTCOME_TRY(set_vaccination_data(model, vacc_data, date, vregion, num_days)); return success(); } @@ -913,22 +920,22 @@ IOResult set_vaccination_data(std::vector>& model, const std::st * * @return An IOResult indicating success or failure. */ -template +template IOResult export_input_data_county_timeseries( std::vector models, const std::string& results_dir, const std::vector& counties, Date date, - const std::vector& scaling_factor_inf, const double scaling_factor_icu, const int num_days, + const std::vector& scaling_factor_inf, const FP scaling_factor_icu, const int num_days, const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path, - const std::vector> immunity_population, const std::string& vaccination_data_path = "") + const std::vector> immunity_population, const std::string& vaccination_data_path = "") { const auto num_groups = (size_t)models[0].parameters.get_num_groups(); assert(scaling_factor_inf.size() == num_groups); assert(num_groups == ConfirmedCasesDataEntry::age_group_names.size()); assert(models.size() == counties.size()); - std::vector> extrapolated_data( - models.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_groups)); + std::vector> extrapolated_data( + models.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_groups)); BOOST_OUTCOME_TRY(auto&& case_data, read_confirmed_cases_data(confirmed_cases_path)); - BOOST_OUTCOME_TRY(auto&& population_data, read_population_data(population_data_path, counties)); + BOOST_OUTCOME_TRY(auto&& population_data, read_population_data(population_data_path, counties)); // empty vector if set_vaccination_data is not set std::vector vacc_data; @@ -940,17 +947,17 @@ IOResult export_input_data_county_timeseries( auto offset_day = offset_date_by_days(date, t); if (!vaccination_data_path.empty()) { - BOOST_OUTCOME_TRY(details::set_vaccination_data(models, vacc_data, offset_day, counties, num_days)); + BOOST_OUTCOME_TRY(details::set_vaccination_data(models, vacc_data, offset_day, counties, num_days)); } // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType. // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available. - BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(models, case_data, counties, offset_day, scaling_factor_inf, - immunity_population)); + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(models, case_data, counties, offset_day, + scaling_factor_inf, immunity_population)); - BOOST_OUTCOME_TRY(details::set_population_data(models, population_data, counties, immunity_population)); + BOOST_OUTCOME_TRY(details::set_population_data(models, population_data, counties, immunity_population)); for (size_t r = 0; r < counties.size(); r++) { extrapolated_data[r][t] = models[r].get_initial_values(); @@ -963,22 +970,22 @@ IOResult export_input_data_county_timeseries( } } } - BOOST_OUTCOME_TRY(save_result(extrapolated_data, counties, static_cast(num_groups), - path_join(results_dir, "Results_rki.h5"))); + BOOST_OUTCOME_TRY(save_result(extrapolated_data, counties, static_cast(num_groups), + path_join(results_dir, "Results_rki.h5"))); - auto extrapolated_rki_data_sum = sum_nodes(std::vector>>{extrapolated_data}); - BOOST_OUTCOME_TRY(save_result({extrapolated_rki_data_sum[0][0]}, {0}, static_cast(num_groups), - path_join(results_dir, "Results_rki_sum.h5"))); + auto extrapolated_rki_data_sum = sum_nodes(std::vector>>{extrapolated_data}); + BOOST_OUTCOME_TRY(save_result({extrapolated_rki_data_sum[0][0]}, {0}, static_cast(num_groups), + path_join(results_dir, "Results_rki_sum.h5"))); return success(); } #else -template +template IOResult export_input_data_county_timeseries(std::vector, const std::string&, const std::vector&, - Date, const std::vector&, const double, const int, + Date, const std::vector&, const FP, const int, const std::string&, const std::string&, const std::string&, - const std::vector>, const std::string&) + const std::vector>, const std::string&) { mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data."); return success(); @@ -1007,25 +1014,26 @@ IOResult export_input_data_county_timeseries(std::vector, const std * * @return An IOResult indicating success or failure. */ -template +template IOResult read_input_data_county(std::vector& model, Date date, const std::vector& county, - const std::vector& scaling_factor_inf, double scaling_factor_icu, + const std::vector& scaling_factor_inf, FP scaling_factor_icu, const std::string& pydata_dir, int num_days, - const std::vector> immunity_population, + const std::vector> immunity_population, bool export_time_series = false) { - BOOST_OUTCOME_TRY(details::set_vaccination_data(model, path_join(pydata_dir, "vacc_county_ageinf_ma7.json"), date, - county, num_days)); + BOOST_OUTCOME_TRY(details::set_vaccination_data(model, path_join(pydata_dir, "vacc_county_ageinf_ma7.json"), + date, county, num_days)); // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType. // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available. - BOOST_OUTCOME_TRY( - details::set_divi_data(model, path_join(pydata_dir, "county_divi_ma7.json"), county, date, scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "county_divi_ma7.json"), county, date, + scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_county_age_ma7.json"), - county, date, scaling_factor_inf, immunity_population)); - BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), - county, immunity_population)); + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, + path_join(pydata_dir, "cases_all_county_age_ma7.json"), + county, date, scaling_factor_inf, immunity_population)); + BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), + county, immunity_population)); if (export_time_series) { // Use only if extrapolated real data is needed for comparison. EXPENSIVE ! @@ -1033,7 +1041,7 @@ IOResult read_input_data_county(std::vector& model, Date date, cons // (This only represents the vectorization of the previous function over all simulation days...) log_info("Exporting time series of extrapolated real data. This may take some minutes. " "For simulation runs over the same time period, deactivate it."); - BOOST_OUTCOME_TRY(export_input_data_county_timeseries( + BOOST_OUTCOME_TRY(export_input_data_county_timeseries( model, pydata_dir, county, date, scaling_factor_inf, scaling_factor_icu, num_days, path_join(pydata_dir, "county_divi_ma7.json"), path_join(pydata_dir, "cases_all_county_age_ma7.json"), path_join(pydata_dir, "county_current_population.json"), immunity_population, @@ -1062,26 +1070,25 @@ IOResult read_input_data_county(std::vector& model, Date date, cons * * @return An IOResult indicating success or failure. */ -template +template IOResult read_input_data(std::vector& model, Date date, const std::vector& node_ids, - const std::vector& scaling_factor_inf, double scaling_factor_icu, + const std::vector& scaling_factor_inf, FP scaling_factor_icu, const std::string& pydata_dir, int num_days, - const std::vector> immunity_population, - bool export_time_series = false) + const std::vector> immunity_population, bool export_time_series = false) { - BOOST_OUTCOME_TRY( - details::set_vaccination_data(model, path_join(pydata_dir, "vaccination_data.json"), date, node_ids, num_days)); + BOOST_OUTCOME_TRY(details::set_vaccination_data(model, path_join(pydata_dir, "vaccination_data.json"), date, + node_ids, num_days)); // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType. // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available. - BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "critical_cases.json"), node_ids, date, - scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "critical_cases.json"), node_ids, date, + scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "confirmed_cases.json"), node_ids, - date, scaling_factor_inf, immunity_population)); - BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "population_data.json"), node_ids, - immunity_population)); + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "confirmed_cases.json"), + node_ids, date, scaling_factor_inf, immunity_population)); + BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "population_data.json"), node_ids, + immunity_population)); if (export_time_series) { // Use only if extrapolated real data is needed for comparison. EXPENSIVE ! @@ -1089,7 +1096,7 @@ IOResult read_input_data(std::vector& model, Date date, const std:: // (This only represents the vectorization of the previous function over all simulation days...) log_info("Exporting time series of extrapolated real data. This may take some minutes. " "For simulation runs over the same time period, deactivate it."); - BOOST_OUTCOME_TRY(export_input_data_county_timeseries( + BOOST_OUTCOME_TRY(export_input_data_county_timeseries( model, pydata_dir, node_ids, date, scaling_factor_inf, scaling_factor_icu, num_days, path_join(pydata_dir, "critical_cases.json"), path_join(pydata_dir, "confirmed_cases.json"), path_join(pydata_dir, "population_data.json"), immunity_population, diff --git a/cpp/models/ode_secirvvs/analyze_result.h b/cpp/models/ode_secirvvs/analyze_result.h index 9d603338b1..097bc614e2 100644 --- a/cpp/models/ode_secirvvs/analyze_result.h +++ b/cpp/models/ode_secirvvs/analyze_result.h @@ -33,36 +33,38 @@ namespace osecirvvs * @param p percentile value in open interval (0, 1) * @return p percentile of the parameters over all runs */ -template -std::vector ensemble_params_percentile(const std::vector>& ensemble_params, double p) +template +std::vector ensemble_params_percentile(const std::vector>& ensemble_params, FP p) { + using std::floor; + assert(p > 0.0 && p < 1.0 && "Invalid percentile value."); auto num_runs = ensemble_params.size(); auto num_nodes = ensemble_params[0].size(); auto num_groups = (size_t)ensemble_params[0][0].parameters.get_num_groups(); auto num_days = ensemble_params[0][0] - .parameters.template get>() + .parameters.template get>() .template size(); - std::vector single_element_ensemble(num_runs); + std::vector single_element_ensemble(num_runs); // lambda function that calculates the percentile of a single parameter std::vector percentile(num_nodes, Model((int)num_groups)); auto param_percentil = [&ensemble_params, p, num_runs, &percentile](auto n, auto get_param) mutable { - std::vector single_element(num_runs); + std::vector single_element(num_runs); for (size_t run = 0; run < num_runs; run++) { auto const& params = ensemble_params[run][n]; single_element[run] = get_param(params); } std::sort(single_element.begin(), single_element.end()); auto& new_params = get_param(percentile[n]); - new_params = single_element[static_cast(num_runs * p)]; + new_params = single_element[static_cast(floor(num_runs * p))]; }; for (size_t node = 0; node < num_nodes; node++) { - percentile[node].parameters.template get>().resize(num_days); - percentile[node].parameters.template get>().resize(num_days); + percentile[node].parameters.template get>().resize(num_days); + percentile[node].parameters.template get>().resize(num_days); for (auto i = AgeGroup(0); i < AgeGroup(num_groups); i++) { //Population @@ -73,112 +75,112 @@ std::vector ensemble_params_percentile(const std::vector auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); //probs param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); //vaccinations param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); for (auto day = SimulationDay(0); day < num_days; ++day) { param_percentil(node, [i, day](auto&& model) -> auto& { - return model.parameters.template get>()[{i, day}]; + return model.parameters.template get>()[{i, day}]; }); param_percentil(node, [i, day](auto&& model) -> auto& { - return model.parameters.template get>()[{i, day}]; + return model.parameters.template get>()[{i, day}]; }); } //virus variants param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; + return model.parameters.template get>()[i]; }); } // group independent params param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); + return model.parameters.template get>(); }); param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); + return model.parameters.template get>(); }); param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); + return model.parameters.template get>(); }); param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); + return model.parameters.template get>(); }); for (size_t run = 0; run < num_runs; run++) { auto const& params = ensemble_params[run][node]; single_element_ensemble[run] = - params.parameters.template get>() * params.populations.get_total(); + params.parameters.template get>() * params.populations.get_total(); } std::sort(single_element_ensemble.begin(), single_element_ensemble.end()); - percentile[node].parameters.template set>( - single_element_ensemble[static_cast(num_runs * p)]); + percentile[node].parameters.template set>( + single_element_ensemble[static_cast(floor(num_runs * p))]); } return percentile; } diff --git a/cpp/models/ode_secirvvs/parameters_io.h b/cpp/models/ode_secirvvs/parameters_io.h index 1e592b8951..83023529e6 100644 --- a/cpp/models/ode_secirvvs/parameters_io.h +++ b/cpp/models/ode_secirvvs/parameters_io.h @@ -52,18 +52,18 @@ namespace details * @see mio::read_confirmed_cases_data * @{ */ -template +template 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, + 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) + 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; @@ -163,7 +163,7 @@ IOResult read_confirmed_cases_data( 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) { + auto try_fix_constraints = [region, i](FP& value, FP error, auto str) { if (value < error) { // this should probably return a failure // but the algorithm is not robust enough to avoid large negative @@ -194,25 +194,24 @@ IOResult read_confirmed_cases_data( return success(); } -template +template 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, + 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) + 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); + 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); } - /**@}*/ /** @@ -225,11 +224,13 @@ IOResult read_confirmed_cases_data( * @see mio::read_confirmed_cases_data * @{ */ -template +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.0) + std::vector>& vnum_rec, FP delay = 14.0) { + using std::trunc; + auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { return a.date < b.date; }); @@ -257,7 +258,7 @@ IOResult read_confirmed_cases_data_fix_recovered(const std::vector read_confirmed_cases_data_fix_recovered(const std::vector read_confirmed_cases_data_fix_recovered(const std::vector +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.0) + Date date, std::vector>& vnum_rec, + FP 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); + 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. -*/ -template + * @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, std::vector const& region, Date date, - const std::vector& scaling_factor_inf, bool set_death = false) + const std::vector& scaling_factor_inf, bool set_death = false) { + using std::round; + auto num_age_groups = (size_t)model[0].parameters.get_num_groups(); assert(scaling_factor_inf.size() == num_age_groups); //TODO: allow vector or scalar valued scaling factors assert(ConfirmedCasesDataEntry::age_group_names.size() == num_age_groups); @@ -328,60 +331,59 @@ IOResult set_confirmed_cases_data(std::vector& model, std::vector> t_InfectedSevere{model.size()}; std::vector> t_InfectedCritical{model.size()}; - std::vector> mu_C_R{model.size()}; - std::vector> mu_I_H{model.size()}; - std::vector> mu_H_U{model.size()}; + std::vector> mu_C_R{model.size()}; + std::vector> mu_I_H{model.size()}; + std::vector> mu_H_U{model.size()}; - std::vector> num_InfectedSymptoms(model.size()); - std::vector> num_death(model.size()); - std::vector> num_rec(model.size()); - std::vector> num_Exposed(model.size()); - std::vector> num_InfectedNoSymptoms(model.size()); - std::vector> num_InfectedSevere(model.size()); - std::vector> num_icu(model.size()); + std::vector> num_InfectedSymptoms(model.size()); + std::vector> num_death(model.size()); + std::vector> num_rec(model.size()); + std::vector> num_Exposed(model.size()); + std::vector> num_InfectedNoSymptoms(model.size()); + std::vector> num_InfectedSevere(model.size()); + std::vector> num_icu(model.size()); /*----------- UNVACCINATED -----------*/ for (size_t county = 0; county < model.size(); county++) { - num_InfectedSymptoms[county] = std::vector(num_age_groups, 0.0); - num_death[county] = std::vector(num_age_groups, 0.0); - num_rec[county] = std::vector(num_age_groups, 0.0); - num_Exposed[county] = std::vector(num_age_groups, 0.0); - num_InfectedNoSymptoms[county] = std::vector(num_age_groups, 0.0); - num_InfectedSevere[county] = std::vector(num_age_groups, 0.0); - num_icu[county] = std::vector(num_age_groups, 0.0); + num_InfectedSymptoms[county] = std::vector(num_age_groups, 0.0); + num_death[county] = std::vector(num_age_groups, 0.0); + num_rec[county] = std::vector(num_age_groups, 0.0); + num_Exposed[county] = std::vector(num_age_groups, 0.0); + num_InfectedNoSymptoms[county] = std::vector(num_age_groups, 0.0); + num_InfectedSevere[county] = std::vector(num_age_groups, 0.0); + num_icu[county] = std::vector(num_age_groups, 0.0); for (size_t group = 0; group < num_age_groups; group++) { t_Exposed[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedNoSymptoms[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedSymptoms[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedSevere[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedCritical[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); + round(static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); + t_InfectedNoSymptoms[county].push_back(static_cast(round(static_cast( + model[county].parameters.template get>()[(AgeGroup)group])))); + t_InfectedSymptoms[county].push_back(static_cast(round( + static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); + t_InfectedSevere[county].push_back(static_cast(round( + static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); + t_InfectedCritical[county].push_back(static_cast(round( + static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); mu_C_R[county].push_back( - model[county].parameters.template get>()[(AgeGroup)group]); + model[county].parameters.template get>()[(AgeGroup)group]); mu_I_H[county].push_back( - model[county].parameters.template get>()[(AgeGroup)group]); - mu_H_U[county].push_back( - model[county].parameters.template get>()[(AgeGroup)group]); + model[county].parameters.template get>()[(AgeGroup)group]); + mu_H_U[county].push_back(model[county].parameters.template get>()[(AgeGroup)group]); } } - BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms, - num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec, - t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, - t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf)); + BOOST_OUTCOME_TRY(read_confirmed_cases_data( + case_data, region, date, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, + num_death, num_rec, t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, t_InfectedCritical, + mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf)); for (size_t county = 0; county < model.size(); county++) { // if (std::accumulate( // num_InfectedSymptoms[county].begin(), // num_InfectedSymptoms[county].end(), - // double(0.0), - // [](const double& a, const double& b) { return evaluate_intermediate(a + b); } + // FP(0.0), + // [](const FP& a, const FP& b) { return evaluate_intermediate(a + b); } // ) > 0) // { size_t num_groups = (size_t)model[county].parameters.get_num_groups(); @@ -413,9 +415,9 @@ IOResult set_confirmed_cases_data(std::vector& model, } // } - if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), double(0.0), - [](const double& a, const double& b) { - return evaluate_intermediate(a + b); + if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), FP(0.0), + [](const FP& a, const FP& b) { + return evaluate_intermediate(a + b); }) == 0.0) { log_warning( "No infections for unvaccinated reported on date {} for region {}. Population data has not been set.", @@ -435,61 +437,61 @@ IOResult set_confirmed_cases_data(std::vector& model, mu_I_H[county].clear(); mu_H_U[county].clear(); - num_InfectedSymptoms[county] = std::vector(num_age_groups, 0.0); - num_death[county] = std::vector(num_age_groups, 0.0); - num_rec[county] = std::vector(num_age_groups, 0.0); - num_Exposed[county] = std::vector(num_age_groups, 0.0); - num_InfectedNoSymptoms[county] = std::vector(num_age_groups, 0.0); - num_InfectedSevere[county] = std::vector(num_age_groups, 0.0); - num_icu[county] = std::vector(num_age_groups, 0.0); + num_InfectedSymptoms[county] = std::vector(num_age_groups, 0.0); + num_death[county] = std::vector(num_age_groups, 0.0); + num_rec[county] = std::vector(num_age_groups, 0.0); + num_Exposed[county] = std::vector(num_age_groups, 0.0); + num_InfectedNoSymptoms[county] = std::vector(num_age_groups, 0.0); + num_InfectedSevere[county] = std::vector(num_age_groups, 0.0); + num_icu[county] = std::vector(num_age_groups, 0.0); for (size_t group = 0; group < num_age_groups; group++) { - double reduc_t = model[0].parameters.template get>()[(AgeGroup)group]; + FP reduc_t = model[0].parameters.template get>()[(AgeGroup)group]; t_Exposed[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedNoSymptoms[county].push_back(static_cast(std::round( - model[county].parameters.template get>()[(AgeGroup)group] * reduc_t))); - t_InfectedSymptoms[county].push_back(static_cast(std::round( - model[county].parameters.template get>()[(AgeGroup)group] * reduc_t))); - t_InfectedSevere[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedCritical[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - - double exp_fac_part_immune = - model[county].parameters.template get>()[(AgeGroup)group]; - double inf_fac_part_immune = - model[county].parameters.template get>()[(AgeGroup)group]; - double hosp_fac_part_immune = + round(static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); + t_InfectedNoSymptoms[county].push_back(static_cast(round(static_cast( + model[county].parameters.template get>()[(AgeGroup)group] * reduc_t)))); + t_InfectedSymptoms[county].push_back(static_cast(round(static_cast( + model[county].parameters.template get>()[(AgeGroup)group] * reduc_t)))); + t_InfectedSevere[county].push_back(static_cast(round( + static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); + t_InfectedCritical[county].push_back(static_cast(round( + static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); + + FP exp_fac_part_immune = + model[county].parameters.template get>()[(AgeGroup)group]; + FP inf_fac_part_immune = + model[county].parameters.template get>()[(AgeGroup)group]; + FP hosp_fac_part_immune = model[county] - .parameters.template get>()[(AgeGroup)group]; - double icu_fac_part_immune = + .parameters.template get>()[(AgeGroup)group]; + FP icu_fac_part_immune = model[county] - .parameters.template get>()[(AgeGroup)group]; - mu_C_R[county].push_back(( - 1 - inf_fac_part_immune / exp_fac_part_immune * - (1 - model[county] - .parameters.template get>()[(AgeGroup)group]))); + .parameters.template get>()[(AgeGroup)group]; + mu_C_R[county].push_back( + (1 - + inf_fac_part_immune / exp_fac_part_immune * + (1 - + model[county].parameters.template get>()[(AgeGroup)group]))); mu_I_H[county].push_back( hosp_fac_part_immune / inf_fac_part_immune * - model[county].parameters.template get>()[(AgeGroup)group]); + model[county].parameters.template get>()[(AgeGroup)group]); // transfer from H to U, D unchanged. - mu_H_U[county].push_back( - icu_fac_part_immune / hosp_fac_part_immune * - model[county].parameters.template get>()[(AgeGroup)group]); + mu_H_U[county].push_back(icu_fac_part_immune / hosp_fac_part_immune * + model[county].parameters.template get>()[(AgeGroup)group]); } } - BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms, - num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec, - t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, - t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf)); + BOOST_OUTCOME_TRY(read_confirmed_cases_data( + case_data, region, date, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, + num_death, num_rec, t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, t_InfectedCritical, + mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf)); for (size_t county = 0; county < model.size(); county++) { // if (std::accumulate( // num_InfectedSymptoms[county].begin(), // num_InfectedSymptoms[county].end(), - // double(0.0), - // [](const double& a, const double& b) { return evaluate_intermediate(a + b); } + // FP(0.0), + // [](const FP& a, const FP& b) { return evaluate_intermediate(a + b); } // ) > 0) // { size_t num_groups = (size_t)model[county].parameters.get_num_groups(); @@ -510,9 +512,9 @@ IOResult set_confirmed_cases_data(std::vector& model, } } // } - if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), double(0.0), - [](const double& a, const double& b) { - return evaluate_intermediate(a + b); + if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), FP(0.0), + [](const FP& a, const FP& b) { + return evaluate_intermediate(a + b); }) == 0.0) { log_warning("No infections for partially vaccinated reported on date {} for region {}. " "Population data has not been set.", @@ -532,54 +534,54 @@ IOResult set_confirmed_cases_data(std::vector& model, mu_I_H[county].clear(); mu_H_U[county].clear(); - num_InfectedSymptoms[county] = std::vector(num_age_groups, 0.0); - num_death[county] = std::vector(num_age_groups, 0.0); - num_rec[county] = std::vector(num_age_groups, 0.0); - num_Exposed[county] = std::vector(num_age_groups, 0.0); - num_InfectedNoSymptoms[county] = std::vector(num_age_groups, 0.0); - num_InfectedSevere[county] = std::vector(num_age_groups, 0.0); - num_icu[county] = std::vector(num_age_groups, 0.0); + num_InfectedSymptoms[county] = std::vector(num_age_groups, 0.0); + num_death[county] = std::vector(num_age_groups, 0.0); + num_rec[county] = std::vector(num_age_groups, 0.0); + num_Exposed[county] = std::vector(num_age_groups, 0.0); + num_InfectedNoSymptoms[county] = std::vector(num_age_groups, 0.0); + num_InfectedSevere[county] = std::vector(num_age_groups, 0.0); + num_icu[county] = std::vector(num_age_groups, 0.0); for (size_t group = 0; group < num_age_groups; group++) { - double reduc_t = model[0].parameters.template get>()[(AgeGroup)group]; + FP reduc_t = model[0].parameters.template get>()[(AgeGroup)group]; t_Exposed[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedNoSymptoms[county].push_back(static_cast(std::round( - model[county].parameters.template get>()[(AgeGroup)group] * reduc_t))); - t_InfectedSymptoms[county].push_back(static_cast(std::round( - model[county].parameters.template get>()[(AgeGroup)group] * reduc_t))); - t_InfectedSevere[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedCritical[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - - double reduc_immune_exp = - model[county].parameters.template get>()[(AgeGroup)group]; - double reduc_immune_inf = - model[county].parameters.template get>()[(AgeGroup)group]; - double reduc_immune_hosp = - model[county].parameters.template get>()[( - AgeGroup)group]; - double reduc_immune_icu = - model[county].parameters.template get>()[( - AgeGroup)group]; - mu_C_R[county].push_back(( - 1 - reduc_immune_inf / reduc_immune_exp * - (1 - model[county] - .parameters.template get>()[(AgeGroup)group]))); + round(static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); + t_InfectedNoSymptoms[county].push_back(static_cast(round(static_cast( + model[county].parameters.template get>()[(AgeGroup)group] * reduc_t)))); + t_InfectedSymptoms[county].push_back(static_cast(round(static_cast( + model[county].parameters.template get>()[(AgeGroup)group] * reduc_t)))); + t_InfectedSevere[county].push_back(static_cast(round( + static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); + t_InfectedCritical[county].push_back(static_cast(round( + static_cast(model[county].parameters.template get>()[(AgeGroup)group])))); + + FP reduc_immune_exp = + model[county].parameters.template get>()[(AgeGroup)group]; + FP reduc_immune_inf = + model[county].parameters.template get>()[(AgeGroup)group]; + FP reduc_immune_hosp = + model[county] + .parameters.template get>()[(AgeGroup)group]; + FP reduc_immune_icu = + model[county] + .parameters.template get>()[(AgeGroup)group]; + mu_C_R[county].push_back( + (1 - + reduc_immune_inf / reduc_immune_exp * + (1 - + model[county].parameters.template get>()[(AgeGroup)group]))); mu_I_H[county].push_back( reduc_immune_hosp / reduc_immune_inf * - model[county].parameters.template get>()[(AgeGroup)group]); + model[county].parameters.template get>()[(AgeGroup)group]); // transfer from H to U, D unchanged. - mu_H_U[county].push_back( - reduc_immune_icu / reduc_immune_hosp * - model[county].parameters.template get>()[(AgeGroup)group]); + mu_H_U[county].push_back(reduc_immune_icu / reduc_immune_hosp * + model[county].parameters.template get>()[(AgeGroup)group]); } } - BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms, - num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec, - t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, - t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf)); + BOOST_OUTCOME_TRY(read_confirmed_cases_data( + case_data, region, date, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, + num_death, num_rec, t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, t_InfectedCritical, + mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf)); for (size_t county = 0; county < model.size(); county++) { size_t num_groups = (size_t)model[county].parameters.get_num_groups(); @@ -600,9 +602,9 @@ IOResult set_confirmed_cases_data(std::vector& model, } } - if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), double(0.0), - [](const double& a, const double& b) { - return evaluate_intermediate(a + b); + if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), FP(0.0), + [](const FP& a, const FP& b) { + return evaluate_intermediate(a + b); }) == 0.0) { log_warning("No infections for vaccinated reported on date {} for region {}. " "Population data has not been set.", @@ -614,36 +616,36 @@ 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. -*/ -template + * @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, - const std::vector& scaling_factor_inf, bool set_death = false) + const std::vector& scaling_factor_inf, bool set_death = false) { BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path)); - BOOST_OUTCOME_TRY(set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf, set_death)); + BOOST_OUTCOME_TRY(set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf, set_death)); return success(); } /** -* @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 + * @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) + Date date, FP scaling_factor_icu) { // DIVI dataset will no longer be updated from CW29 2024 on. if (!is_divi_data_available(date)) { @@ -652,19 +654,19 @@ IOResult set_divi_data(std::vector& model, const std::string& path, date); return success(); } - std::vector sum_mu_I_U(vregion.size(), 0); - std::vector> mu_I_U{model.size()}; + std::vector sum_mu_I_U(vregion.size(), 0); + std::vector> mu_I_U{model.size()}; for (size_t region = 0; region < vregion.size(); region++) { auto num_groups = model[region].parameters.get_num_groups(); for (auto i = AgeGroup(0); i < num_groups; i++) { - sum_mu_I_U[region] += model[region].parameters.template get>()[i] * - model[region].parameters.template get>()[i]; - mu_I_U[region].push_back(model[region].parameters.template get>()[i] * - model[region].parameters.template get>()[i]); + sum_mu_I_U[region] += model[region].parameters.template get>()[i] * + model[region].parameters.template get>()[i]; + mu_I_U[region].push_back(model[region].parameters.template get>()[i] * + model[region].parameters.template get>()[i]); } } - std::vector num_icu(model.size(), 0.0); - BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu)); + std::vector num_icu(model.size(), 0.0); + BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu)); for (size_t region = 0; region < vregion.size(); region++) { auto num_groups = model[region].parameters.get_num_groups(); @@ -685,34 +687,37 @@ IOResult set_divi_data(std::vector& model, const std::string& path, * @param[in] vregion vector of keys of the regions of interest * @param[in] date Date for which the arrays are initialized */ -template -IOResult set_population_data(std::vector& model, const std::vector>& num_population, +template +IOResult set_population_data(std::vector& model, const std::vector>& num_population, const std::vector& case_data, const std::vector& vregion, Date date) { + using std::max; + using std::min; + auto num_age_groups = ConfirmedCasesDataEntry::age_group_names.size(); - std::vector> vnum_rec(model.size(), std::vector(num_age_groups, 0.0)); + std::vector> vnum_rec(model.size(), std::vector(num_age_groups, 0.0)); - BOOST_OUTCOME_TRY(read_confirmed_cases_data_fix_recovered(case_data, vregion, date, vnum_rec, 14.)); + BOOST_OUTCOME_TRY(read_confirmed_cases_data_fix_recovered(case_data, vregion, date, vnum_rec, 14.)); for (size_t region = 0; region < vregion.size(); region++) { - if (std::accumulate(num_population[region].begin(), num_population[region].end(), double(0.0), - [](const double& a, const double& b) { - return evaluate_intermediate(a + b); + if (std::accumulate(num_population[region].begin(), num_population[region].end(), FP(0.0), + [](const FP& a, const FP& b) { + return evaluate_intermediate(a + b); }) > 0) { auto num_groups = model[region].parameters.get_num_groups(); for (auto i = AgeGroup(0); i < num_groups; i++) { - double S_v = std::min( - model[region].parameters.template get>()[{i, SimulationDay(0)}] + - vnum_rec[region][size_t(i)], - num_population[region][size_t(i)]); - double S_pv = std::max( - model[region].parameters.template get>()[{i, SimulationDay(0)}] - - model[region].parameters.template get>()[{i, SimulationDay(0)}], + FP S_v = + min(model[region].parameters.template get>()[{i, SimulationDay(0)}] + + vnum_rec[region][size_t(i)], + num_population[region][size_t(i)]); + FP S_pv = max( + model[region].parameters.template get>()[{i, SimulationDay(0)}] - + model[region].parameters.template get>()[{i, SimulationDay(0)}], 0.0); // use std::max with 0 - double S; + FP S; if (num_population[region][size_t(i)] - S_pv - S_v < 0.0) { log_warning("Number of vaccinated persons greater than population in county {}, age group {}.", region, size_t(i)); @@ -723,32 +728,30 @@ IOResult set_population_data(std::vector& model, const std::vector< S = num_population[region][size_t(i)] - S_pv - S_v; } - double denom_E = - 1 / (S + S_pv * model[region].parameters.template get>()[i] + - S_v * model[region].parameters.template get>()[i]); - double denom_C = - 1 / (S + S_pv * model[region].parameters.template get>()[i] + - S_v * model[region].parameters.template get>()[i]); - double denom_I = - 1 / - (S + - S_pv * model[region].parameters.template get>()[i] + - S_v * model[region].parameters.template get>()[i]); - double denom_HU = + FP denom_E = + 1 / (S + S_pv * model[region].parameters.template get>()[i] + + S_v * model[region].parameters.template get>()[i]); + FP denom_C = + 1 / (S + S_pv * model[region].parameters.template get>()[i] + + S_v * model[region].parameters.template get>()[i]); + FP denom_I = 1 / - (S + - S_pv * model[region] - .parameters.template get>()[i] + - S_v * model[region] - .parameters.template get>()[i]); + (S + S_pv * model[region].parameters.template get>()[i] + + S_v * model[region].parameters.template get>()[i]); + FP denom_HU = + 1 / (S + + S_pv * model[region] + .parameters.template get>()[i] + + S_v * model[region] + .parameters.template get>()[i]); model[region].populations[{i, InfectionState::ExposedNaive}] = S * model[region].populations[{i, InfectionState::ExposedNaive}] * denom_E; model[region].populations[{i, InfectionState::ExposedPartialImmunity}] = - S_pv * model[region].parameters.template get>()[i] * + S_pv * model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::ExposedPartialImmunity}] * denom_E; model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] = - S_v * model[region].parameters.template get>()[i] * + S_v * model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] * denom_E; model[region].populations[{i, InfectionState::InfectedNoSymptomsNaive}] = @@ -770,47 +773,45 @@ IOResult set_population_data(std::vector& model, const std::vector< model[region].populations[{i, InfectionState::InfectedSymptomsNaive}] = S * model[region].populations[{i, InfectionState::InfectedSymptomsNaive}] * denom_I; model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] = - S_pv * model[region].parameters.template get>()[i] * + S_pv * model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] * denom_I; model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] = - S_v * model[region].parameters.template get>()[i] * + S_v * model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] * denom_I; model[region].populations[{i, InfectionState::InfectedSymptomsNaiveConfirmed}] = S * model[region].populations[{i, InfectionState::InfectedSymptomsNaiveConfirmed}] * denom_I; model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] = - S_pv * model[region].parameters.template get>()[i] * + S_pv * model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] * denom_I; model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] = - S_v * model[region].parameters.template get>()[i] * + S_v * model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] * denom_I; model[region].populations[{i, InfectionState::InfectedSevereNaive}] = S * model[region].populations[{i, InfectionState::InfectedSevereNaive}] * denom_HU; model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] = S_pv * - model[region].parameters.template get>()[i] * + model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] * denom_HU; model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] = S_v * - model[region] - .parameters.template get>()[i] * + model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] * denom_HU; model[region].populations[{i, InfectionState::InfectedCriticalPartialImmunity}] = S_pv * - model[region].parameters.template get>()[i] * + model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::InfectedCriticalNaive}] * denom_HU; model[region].populations[{i, InfectionState::InfectedCriticalImprovedImmunity}] = S_v * - model[region] - .parameters.template get>()[i] * + model[region].parameters.template get>()[i] * model[region].populations[{i, InfectionState::InfectedCriticalNaive}] * denom_HU; model[region].populations[{i, InfectionState::InfectedCriticalNaive}] = S * model[region].populations[{i, InfectionState::InfectedCriticalNaive}] * denom_HU; model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] = - model[region].parameters.template get>()[{i, SimulationDay(0)}] + + model[region].parameters.template get>()[{i, SimulationDay(0)}] + model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] - (model[region].populations[{i, InfectionState::InfectedSymptomsNaive}] + model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] + @@ -828,7 +829,7 @@ IOResult set_population_data(std::vector& model, const std::vector< model[region].populations[{i, InfectionState::DeadPartialImmunity}] + model[region].populations[{i, InfectionState::DeadImprovedImmunity}]); - model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] = std::min( + model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] = min( S_v - model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] - model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunity}] - model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] - @@ -837,9 +838,9 @@ IOResult set_population_data(std::vector& model, const std::vector< model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] - model[region].populations[{i, InfectionState::InfectedCriticalImprovedImmunity}] - model[region].populations[{i, InfectionState::DeadImprovedImmunity}], - std::max(0.0, double(model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}]))); + max(0.0, FP(model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}]))); - model[region].populations[{i, InfectionState::SusceptiblePartialImmunity}] = std::max( + model[region].populations[{i, InfectionState::SusceptiblePartialImmunity}] = max( 0.0, S_pv - model[region].populations[{i, InfectionState::ExposedPartialImmunity}] - model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunity}] - @@ -880,14 +881,14 @@ IOResult set_population_data(std::vector& model, const std::vector< * @param[in] vregion Vector of keys of the regions of interest. * @param[in] date Date for which the arrays are initialized. */ -template +template IOResult set_population_data(std::vector& model, const std::string& path, const std::string& path_rki, const std::vector& vregion, Date date) { - BOOST_OUTCOME_TRY(auto&& num_population, read_population_data(path, vregion)); + BOOST_OUTCOME_TRY(auto&& num_population, read_population_data(path, vregion)); BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path_rki)); - BOOST_OUTCOME_TRY(set_population_data(model, num_population, rki_data, vregion, date)); + BOOST_OUTCOME_TRY(set_population_data(model, num_population, rki_data, vregion, date)); return success(); } @@ -901,19 +902,21 @@ IOResult set_population_data(std::vector& model, const std::string& * @param[in] vregion Vector of region identifiers. * @param[in] num_days Number of days for which the simulation is run. */ -template +template IOResult set_vaccination_data(std::vector>& model, const std::vector& vacc_data, Date date, const std::vector& vregion, int num_days) { + using std::floor; + auto num_groups = model[0].parameters.get_num_groups(); // type conversion from UncertainValue -> FP -> int auto days_until_effective1 = static_cast( - static_cast(model[0].parameters.template get>()[AgeGroup(0)])); - auto days_until_effective2 = static_cast( - static_cast(model[0].parameters.template get>()[AgeGroup(0)])); + floor(static_cast(model[0].parameters.template get>()[AgeGroup(0)]))); + auto days_until_effective2 = static_cast(floor( + static_cast(model[0].parameters.template get>()[AgeGroup(0)]))); auto vaccination_distance = - static_cast(static_cast(model[0].parameters.template get>()[AgeGroup(0)])); + static_cast(floor(static_cast(model[0].parameters.template get>()[AgeGroup(0)]))); // iterate over regions (e.g., counties) for (size_t i = 0; i < model.size(); ++i) { @@ -1050,7 +1053,7 @@ IOResult set_vaccination_data(std::vector>& model, const std::ve * @param[in] vregion Vector of region identifiers. * @param[in] num_days Number of days for which the simulation is run. */ -template +template IOResult set_vaccination_data(std::vector>& model, const std::string& path, Date date, const std::vector& vregion, int num_days) { @@ -1075,7 +1078,7 @@ IOResult set_vaccination_data(std::vector>& model, const std::st return success(); } BOOST_OUTCOME_TRY(auto&& vacc_data, read_vaccination_data(path)); - BOOST_OUTCOME_TRY(set_vaccination_data(model, vacc_data, date, vregion, num_days)); + BOOST_OUTCOME_TRY(set_vaccination_data(model, vacc_data, date, vregion, num_days)); return success(); } @@ -1100,10 +1103,10 @@ IOResult set_vaccination_data(std::vector>& model, const std::st * @param[in] population_data_path Path to population data file. * @param[in] vaccination_data_path Path to vaccination data file. */ -template +template IOResult export_input_data_county_timeseries( std::vector models, const std::string& results_dir, const std::vector& counties, Date date, - const std::vector& scaling_factor_inf, const double scaling_factor_icu, const int num_days, + const std::vector& scaling_factor_inf, const FP scaling_factor_icu, const int num_days, const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path, const std::string& vaccination_data_path = "") { @@ -1111,11 +1114,11 @@ IOResult export_input_data_county_timeseries( assert(scaling_factor_inf.size() == num_groups); assert(num_groups == ConfirmedCasesDataEntry::age_group_names.size()); assert(models.size() == counties.size()); - std::vector> extrapolated_data( - models.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_groups)); + std::vector> extrapolated_data( + models.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_groups)); BOOST_OUTCOME_TRY(auto&& case_data, read_confirmed_cases_data(confirmed_cases_path)); - BOOST_OUTCOME_TRY(auto&& population_data, read_population_data(population_data_path, counties)); + BOOST_OUTCOME_TRY(auto&& population_data, read_population_data(population_data_path, counties)); // empty vector if set_vaccination_data is not set std::vector vacc_data; @@ -1127,16 +1130,16 @@ IOResult export_input_data_county_timeseries( auto offset_day = offset_date_by_days(date, t); if (!vaccination_data_path.empty()) { - BOOST_OUTCOME_TRY(details::set_vaccination_data(models, vacc_data, offset_day, counties, num_days)); + BOOST_OUTCOME_TRY(details::set_vaccination_data(models, vacc_data, offset_day, counties, num_days)); } // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType. // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available. - BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu)); BOOST_OUTCOME_TRY( - details::set_confirmed_cases_data(models, case_data, counties, offset_day, scaling_factor_inf, true)); - BOOST_OUTCOME_TRY(details::set_population_data(models, population_data, case_data, counties, offset_day)); + details::set_confirmed_cases_data(models, case_data, counties, offset_day, scaling_factor_inf, true)); + BOOST_OUTCOME_TRY(details::set_population_data(models, population_data, case_data, counties, offset_day)); for (size_t r = 0; r < counties.size(); r++) { extrapolated_data[r][t] = models[r].get_initial_values(); @@ -1149,19 +1152,19 @@ IOResult export_input_data_county_timeseries( } } } - BOOST_OUTCOME_TRY(save_result(extrapolated_data, counties, static_cast(num_groups), - path_join(results_dir, "Results_rki.h5"))); + BOOST_OUTCOME_TRY(save_result(extrapolated_data, counties, static_cast(num_groups), + path_join(results_dir, "Results_rki.h5"))); - auto extrapolated_rki_data_sum = sum_nodes(std::vector>>{extrapolated_data}); - BOOST_OUTCOME_TRY(save_result({extrapolated_rki_data_sum[0][0]}, {0}, static_cast(num_groups), - path_join(results_dir, "Results_rki_sum.h5"))); + auto extrapolated_rki_data_sum = sum_nodes(std::vector>>{extrapolated_data}); + BOOST_OUTCOME_TRY(save_result({extrapolated_rki_data_sum[0][0]}, {0}, static_cast(num_groups), + path_join(results_dir, "Results_rki_sum.h5"))); return success(); } #else -template +template IOResult export_input_data_county_timeseries(std::vector, const std::string&, const std::vector&, - Date, const std::vector&, const double, const int, + Date, const std::vector&, const FP, const int, const std::string&, const std::string&, const std::string&, const std::string&) { @@ -1172,37 +1175,37 @@ IOResult export_input_data_county_timeseries(std::vector, const std #endif //MEMILIO_HAS_HDF5 /** - * Reads compartments for German counties at a specified date from data files. - * Estimates all compartments from available data using the model parameters, so the - * model parameters must be set before calling this function. - * Uses data files that contain centered 7-day moving average. - * @param[in, out] model Vector of SECIRVVS models, one per county. - * @param[in] date Date for which the data should be read. - * @param[in] county Ids of the counties. - * @param[in] scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county. - * @param[in] scaling_factor_icu Factor of ICU cases to account for underreporting. - * @param[in] pydata_dir Directory that contains the data files. - * @param[in] num_days Number of days to be simulated; required to load data for vaccinations during the simulation. - * @param[in] export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files. - */ -template + * Reads compartments for German counties at a specified date from data files. + * Estimates all compartments from available data using the model parameters, so the + * model parameters must be set before calling this function. + * Uses data files that contain centered 7-day moving average. + * @param[in, out] model Vector of SECIRVVS models, one per county. + * @param[in] date Date for which the data should be read. + * @param[in] county Ids of the counties. + * @param[in] scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county. + * @param[in] scaling_factor_icu Factor of ICU cases to account for underreporting. + * @param[in] pydata_dir Directory that contains the data files. + * @param[in] num_days Number of days to be simulated; required to load data for vaccinations during the simulation. + * @param[in] export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files. + */ +template IOResult read_input_data_county(std::vector& model, Date date, const std::vector& county, - const std::vector& scaling_factor_inf, double scaling_factor_icu, + const std::vector& scaling_factor_inf, FP scaling_factor_icu, const std::string& pydata_dir, int num_days, bool export_time_series = false) { - BOOST_OUTCOME_TRY(details::set_vaccination_data(model, path_join(pydata_dir, "vacc_county_ageinf_ma7.json"), date, - county, num_days)); + BOOST_OUTCOME_TRY(details::set_vaccination_data(model, path_join(pydata_dir, "vacc_county_ageinf_ma7.json"), + date, county, num_days)); // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType. // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available. - BOOST_OUTCOME_TRY( - details::set_divi_data(model, path_join(pydata_dir, "county_divi_ma7.json"), county, date, scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "county_divi_ma7.json"), county, date, + scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_county_age_ma7.json"), - county, date, scaling_factor_inf)); - BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), - path_join(pydata_dir, "cases_all_county_age_ma7.json"), county, - date)); + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data( + model, path_join(pydata_dir, "cases_all_county_age_ma7.json"), county, date, scaling_factor_inf)); + BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), + path_join(pydata_dir, "cases_all_county_age_ma7.json"), county, + date)); if (export_time_series) { // Use only if extrapolated real data is needed for comparison. EXPENSIVE ! @@ -1210,7 +1213,7 @@ IOResult read_input_data_county(std::vector& model, Date date, cons // (This only represents the vectorization of the previous function over all simulation days...) log_warning("Exporting time series of extrapolated real data. This may take some minutes. " "For simulation runs over the same time period, deactivate it."); - BOOST_OUTCOME_TRY(export_input_data_county_timeseries( + BOOST_OUTCOME_TRY(export_input_data_county_timeseries( model, pydata_dir, county, date, scaling_factor_inf, scaling_factor_icu, num_days, path_join(pydata_dir, "county_divi_ma7.json"), path_join(pydata_dir, "cases_all_county_age_ma7.json"), path_join(pydata_dir, "county_current_population.json"), @@ -1221,37 +1224,37 @@ IOResult read_input_data_county(std::vector& model, Date date, cons } /** - * Reads compartments for German counties at a specified date from data files. - * Estimates all compartments from available data using the model parameters, so the - * model parameters must be set before calling this function. - * Uses data files that contain centered 7-day moving average. - * @param[in, out] model Vector of SECIRVVS models, one per county. - * @param[in] date Date for which the data should be read. - * @param[in] node_ids Ids of the nodes. - * @param[in] scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county. - * @param[in] scaling_factor_icu Factor of ICU cases to account for underreporting. - * @param[in] pydata_dir Directory that contains the data files. - * @param[in] num_days Number of days to be simulated; required to load data for vaccinations during the simulation. - * @param[in] export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files. - */ -template + * Reads compartments for German counties at a specified date from data files. + * Estimates all compartments from available data using the model parameters, so the + * model parameters must be set before calling this function. + * Uses data files that contain centered 7-day moving average. + * @param[in, out] model Vector of SECIRVVS models, one per county. + * @param[in] date Date for which the data should be read. + * @param[in] node_ids Ids of the nodes. + * @param[in] scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county. + * @param[in] scaling_factor_icu Factor of ICU cases to account for underreporting. + * @param[in] pydata_dir Directory that contains the data files. + * @param[in] num_days Number of days to be simulated; required to load data for vaccinations during the simulation. + * @param[in] export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files. + */ +template IOResult read_input_data(std::vector& model, Date date, const std::vector& node_ids, - const std::vector& scaling_factor_inf, double scaling_factor_icu, + const std::vector& scaling_factor_inf, FP scaling_factor_icu, const std::string& pydata_dir, int num_days, bool export_time_series = false) { - BOOST_OUTCOME_TRY( - details::set_vaccination_data(model, path_join(pydata_dir, "vaccination_data.json"), date, node_ids, num_days)); + BOOST_OUTCOME_TRY(details::set_vaccination_data(model, path_join(pydata_dir, "vaccination_data.json"), date, + node_ids, num_days)); // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType. // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available. - BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "critical_cases.json"), node_ids, date, - scaling_factor_icu)); + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "critical_cases.json"), node_ids, date, + scaling_factor_icu)); - BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "confirmed_cases.json"), node_ids, - date, scaling_factor_inf)); - BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "population_data.json"), - path_join(pydata_dir, "confirmed_cases.json"), node_ids, date)); + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "confirmed_cases.json"), + node_ids, date, scaling_factor_inf)); + BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "population_data.json"), + path_join(pydata_dir, "confirmed_cases.json"), node_ids, date)); if (export_time_series) { // Use only if extrapolated real data is needed for comparison. EXPENSIVE ! @@ -1259,7 +1262,7 @@ IOResult read_input_data(std::vector& model, Date date, const std:: // (This only represents the vectorization of the previous function over all simulation days...) log_warning("Exporting time series of extrapolated real data. This may take some minutes. " "For simulation runs over the same time period, deactivate it."); - BOOST_OUTCOME_TRY(export_input_data_county_timeseries( + BOOST_OUTCOME_TRY(export_input_data_county_timeseries( model, pydata_dir, node_ids, date, scaling_factor_inf, scaling_factor_icu, num_days, path_join(pydata_dir, "divi_data.json"), path_join(pydata_dir, "confirmed_cases.json"), path_join(pydata_dir, "population_data.json"), path_join(pydata_dir, "vaccination_data.json"))); diff --git a/cpp/tests/test_analyze_result.cpp b/cpp/tests/test_analyze_result.cpp index 5360c7c3de..2cc35b92b9 100644 --- a/cpp/tests/test_analyze_result.cpp +++ b/cpp/tests/test_analyze_result.cpp @@ -351,10 +351,10 @@ TEST(TestEnsemblePercentile, basic) ensemble.back()[1].add_time_point(2.0, (Vec(2) << 0.0, 0.0).finished()); ensemble.back()[1].add_time_point(3.0, (Vec(2) << 0.0, 0.0).finished()); - auto q1 = mio::ensemble_percentile(ensemble, 0.2); - auto q2 = mio::ensemble_percentile(ensemble, 0.4); - auto q3 = mio::ensemble_percentile(ensemble, 0.7); - auto q4 = mio::ensemble_percentile(ensemble, 0.9); + auto q1 = mio::ensemble_percentile(ensemble, 0.2); + auto q2 = mio::ensemble_percentile(ensemble, 0.4); + auto q3 = mio::ensemble_percentile(ensemble, 0.7); + auto q4 = mio::ensemble_percentile(ensemble, 0.9); //checking only a few elements ASSERT_EQ(q1.size(), 2); diff --git a/cpp/tests/test_mobility_io.cpp b/cpp/tests/test_mobility_io.cpp index b1df3625da..252c29d8b1 100644 --- a/cpp/tests/test_mobility_io.cpp +++ b/cpp/tests/test_mobility_io.cpp @@ -61,7 +61,7 @@ TEST(TestReadMobility, readFormatted) file.close(); } - auto matrix_read = mio::read_mobility_formatted("test_mobility.txt"); + auto matrix_read = mio::read_mobility_formatted("test_mobility.txt"); ASSERT_TRUE(matrix_read); ASSERT_EQ(test_matrix.rows(), matrix_read.value().rows()); ASSERT_EQ(test_matrix.cols(), matrix_read.value().cols()); @@ -109,7 +109,7 @@ TEST(TestReadMobility, readPlain) test_matrix(5, 4) = 0.3497; test_matrix(5, 5) = 0.1544; - auto matrix_read = mio::read_mobility_plain(get_test_data_file_path("contacts.txt")); + auto matrix_read = mio::read_mobility_plain(get_test_data_file_path("contacts.txt")); ASSERT_TRUE(matrix_read); ASSERT_EQ(test_matrix.rows(), matrix_read.value().rows()); diff --git a/cpp/tests/test_odesecir.cpp b/cpp/tests/test_odesecir.cpp index 24533cc68d..da46a87ba2 100755 --- a/cpp/tests/test_odesecir.cpp +++ b/cpp/tests/test_odesecir.cpp @@ -1278,17 +1278,17 @@ TEST_F(ModelTestOdeSecir, read_input_data) auto model2 = std::vector>{model}; auto model3 = std::vector>{model}; - auto read_result1 = mio::osecir::read_input_data_county(model1, {2020, 12, 01}, {1002}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - TEST_GERMANY_PYDATA_DIR, 10); + auto read_result1 = mio::osecir::read_input_data_county( + model1, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), 1.0, TEST_GERMANY_PYDATA_DIR, + 10); - auto read_result2 = - mio::osecir::read_input_data(model2, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), - 1.0, TEST_GERMANY_PYDATA_DIR, 10); + auto read_result2 = mio::osecir::read_input_data(model2, {2020, 12, 01}, {1002}, + std::vector(size_t(num_age_groups), 1.0), 1.0, + TEST_GERMANY_PYDATA_DIR, 10); - auto read_result_district = - mio::osecir::read_input_data(model3, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), - 1.0, mio::path_join(TEST_DATA_DIR, "District", "pydata"), 10); + auto read_result_district = mio::osecir::read_input_data( + model3, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), 1.0, + mio::path_join(TEST_DATA_DIR, "District", "pydata"), 10); EXPECT_THAT(read_result1, IsSuccess()); EXPECT_THAT(read_result2, IsSuccess()); @@ -1321,20 +1321,21 @@ TEST_F(ModelTestOdeSecir, export_time_series_init) // Test exporting time series std::vector> models{model}; - EXPECT_THAT(mio::osecir::export_input_data_county_timeseries( + EXPECT_THAT(mio::osecir::export_input_data_county_timeseries( models, tmp_results_dir, {1002}, {2020, 12, 01}, std::vector(size_t(num_age_groups), 1.0), 1.0, 2, mio::path_join(pydata_dir_Germany, "county_divi_ma7.json"), mio::path_join(pydata_dir_Germany, "cases_all_county_age_ma7.json"), mio::path_join(pydata_dir_Germany, "county_current_population.json")), IsSuccess()); - auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); + auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); EXPECT_THAT(data_extrapolated, IsSuccess()); // Values were generated by the tested function export_input_data_county_timeseries; // can only check stability of the implementation, not correctness auto expected_results = - mio::read_result(mio::path_join(TEST_DATA_DIR, "export_time_series_initialization_osecir.h5")).value(); + mio::read_result(mio::path_join(TEST_DATA_DIR, "export_time_series_initialization_osecir.h5")) + .value(); EXPECT_THAT(print_wrap(data_extrapolated.value()[0].get_groups().matrix()), MatrixNear(print_wrap(expected_results[0].get_groups().matrix()), 1e-5, 1e-5)); @@ -1351,14 +1352,14 @@ TEST_F(ModelTestOdeSecir, export_time_series_init_old_date) // Test exporting time series std::vector> models{model}; - EXPECT_THAT(mio::osecir::export_input_data_county_timeseries( + EXPECT_THAT(mio::osecir::export_input_data_county_timeseries( models, tmp_results_dir, {1002}, {1000, 12, 01}, std::vector(size_t(num_age_groups), 1.0), 1.0, 0, mio::path_join(pydata_dir_Germany, "county_divi_ma7.json"), mio::path_join(pydata_dir_Germany, "cases_all_county_age_ma7.json"), mio::path_join(pydata_dir_Germany, "county_current_population.json")), IsSuccess()); - auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); + auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); EXPECT_THAT(data_extrapolated, IsSuccess()); auto results_extrapolated = data_extrapolated.value()[0].get_groups().get_value(0); @@ -1366,7 +1367,7 @@ TEST_F(ModelTestOdeSecir, export_time_series_init_old_date) // read population data std::string path = mio::path_join(pydata_dir_Germany, "county_current_population.json"); const std::vector region{1002}; - auto population_data = mio::read_population_data(path, region).value(); + auto population_data = mio::read_population_data(path, region).value(); // So, the expected values are the population data in the susceptible compartments and zeros in the other compartments. for (size_t i = 0; i < num_age_groups; i++) { @@ -1385,9 +1386,9 @@ TEST_F(ModelTestOdeSecir, model_initialization) auto model_vector = std::vector>{model}; const auto pydata_dir_Germany = mio::path_join(TEST_DATA_DIR, "Germany", "pydata"); - EXPECT_THAT(mio::osecir::read_input_data_county(model_vector, {2020, 12, 01}, {1002}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - pydata_dir_Germany, 2, false), + EXPECT_THAT(mio::osecir::read_input_data_county(model_vector, {2020, 12, 01}, {1002}, + std::vector(size_t(num_age_groups), 1.0), 1.0, + pydata_dir_Germany, 2, false), IsSuccess()); // Values from data/export_time_series_init_osecir.h5, for reading in comparison @@ -1412,16 +1413,16 @@ TEST_F(ModelTestOdeSecir, model_initialization_old_date) auto model_vector = std::vector>{model}; const auto pydata_dir_Germany = mio::path_join(TEST_DATA_DIR, "Germany", "pydata"); - EXPECT_THAT(mio::osecir::read_input_data_county(model_vector, {1000, 12, 01}, {1002}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - pydata_dir_Germany, 0, false), + EXPECT_THAT(mio::osecir::read_input_data_county(model_vector, {1000, 12, 01}, {1002}, + std::vector(size_t(num_age_groups), 1.0), 1.0, + pydata_dir_Germany, 0, false), IsSuccess()); // if we enter an old date, the model only should be initialized with the population data. // read population data std::string path = mio::path_join(pydata_dir_Germany, "county_current_population.json"); const std::vector region{1002}; - auto population_data = mio::read_population_data(path, region).value(); + auto population_data = mio::read_population_data(path, region).value(); // So, the expected values are the population data in the susceptible compartments and zeros in the other compartments. auto expected_values = @@ -1449,7 +1450,8 @@ TEST(TestOdeSecir, set_divi_data_invalid_dates) auto model_vector = std::vector>{model}; // Test with date before DIVI dataset was available. - EXPECT_THAT(mio::osecir::details::set_divi_data(model_vector, "", {1001}, {2019, 12, 01}, 1.0), IsSuccess()); + EXPECT_THAT(mio::osecir::details::set_divi_data(model_vector, "", {1001}, {2019, 12, 01}, 1.0), + IsSuccess()); // Assure that populations is the same as before. EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), MatrixNear(print_wrap(model.populations.array().cast()), 1e-10, 1e-10)); @@ -1465,7 +1467,7 @@ TEST(TestOdeSecir, set_divi_data_empty_data) std::vector num_icu(regions.size(), 0.0); mio::Date date(2020, 4, 1); - auto result = mio::compute_divi_data(empty_data, regions, date, num_icu); + auto result = mio::compute_divi_data(empty_data, regions, date, num_icu); // Expect failure due to empty DIVI data EXPECT_FALSE(result); @@ -1504,9 +1506,9 @@ TEST_F(ModelTestOdeSecir, set_confirmed_cases_data_with_ICU) // calculate ICU values using set_confirmed_cases_data auto model_vector = std::vector>{model}; auto scaling_factor_inf = std::vector(size_t(model.parameters.get_num_groups()), 1.0); - EXPECT_THAT( - mio::osecir::details::set_confirmed_cases_data(model_vector, case_data, {1002}, mid_day, scaling_factor_inf), - IsSuccess()); + EXPECT_THAT(mio::osecir::details::set_confirmed_cases_data(model_vector, case_data, {1002}, mid_day, + scaling_factor_inf), + IsSuccess()); // Since, TimeInfectedCritical is 1, the number of ICU cases is the difference of confirmed cases between two days, which is 1. // We only have an entry for age group 2. All other age groups should be zero. @@ -1526,7 +1528,7 @@ TEST(TestOdeSecir, read_population_data_failure) invalid_data.push_back(invalid_entry); // Test that read_population_data returns failure with correct message - auto result = mio::read_population_data(invalid_data, {1001}); + auto result = mio::read_population_data(invalid_data, {1001}); EXPECT_FALSE(result); EXPECT_EQ(result.error().code(), mio::StatusCode::InvalidFileFormat); EXPECT_EQ(result.error().message(), "File with county population expected."); @@ -1557,10 +1559,12 @@ TEST(TestOdeSecirIO, read_input_data_county_aggregates_one_group) std::vector scale1{1.0}; // Initialize both models - ASSERT_THAT(mio::osecir::read_input_data_county(models6, date, counties, scale6, 1.0, pydata_dir_Germany), - IsSuccess()); - ASSERT_THAT(mio::osecir::read_input_data_county(models1, date, counties, scale1, 1.0, pydata_dir_Germany), - IsSuccess()); + ASSERT_THAT( + mio::osecir::read_input_data_county(models6, date, counties, scale6, 1.0, pydata_dir_Germany), + IsSuccess()); + ASSERT_THAT( + mio::osecir::read_input_data_county(models1, date, counties, scale1, 1.0, pydata_dir_Germany), + IsSuccess()); // Aggreagate the results from the model with 6 age groups and compare with the model with 1 age group const auto& m6 = models6[0]; @@ -1593,8 +1597,8 @@ TEST(TestOdeSecirIO, set_population_data_single_age_group) std::vector regions = {1002}; // Set population data for both models - EXPECT_THAT(mio::osecir::details::set_population_data(models6, population_data6, regions), IsSuccess()); - EXPECT_THAT(mio::osecir::details::set_population_data(models1, population_data1, regions), IsSuccess()); + EXPECT_THAT(mio::osecir::details::set_population_data(models6, population_data6, regions), IsSuccess()); + EXPECT_THAT(mio::osecir::details::set_population_data(models1, population_data1, regions), IsSuccess()); // Sum all compartments across age groups in 6-group model and compare 1-group model const double tol = 1e-10; @@ -1647,11 +1651,11 @@ TEST(TestOdeSecirIO, set_confirmed_cases_data_single_age_group) std::vector scaling_factors = {1.0}; // Set confirmed cases data for both models - EXPECT_THAT(mio::osecir::details::set_confirmed_cases_data(models6, case_data, regions, mio::Date(2020, 12, 1), - scaling_factors), + EXPECT_THAT(mio::osecir::details::set_confirmed_cases_data(models6, case_data, regions, + mio::Date(2020, 12, 1), scaling_factors), IsSuccess()); - EXPECT_THAT(mio::osecir::details::set_confirmed_cases_data(models1, case_data, regions, mio::Date(2020, 12, 1), - scaling_factors), + EXPECT_THAT(mio::osecir::details::set_confirmed_cases_data(models1, case_data, regions, + mio::Date(2020, 12, 1), scaling_factors), IsSuccess()); // Sum all compartments across age groups in 6-group model should be equal to 1-group model @@ -1691,10 +1695,10 @@ TEST(TestOdeSecirIO, set_divi_data_single_age_group) double scaling_factor_icu = 1.0; mio::Date date(2020, 12, 1); std::string divi_data_path = mio::path_join(TEST_DATA_DIR, "Germany", "pydata", "county_divi_ma7.json"); - auto result_6_groups = - mio::osecir::details::set_divi_data(models_6_groups, divi_data_path, regions, date, scaling_factor_icu); - auto result_1_group = - mio::osecir::details::set_divi_data(models_1_group, divi_data_path, regions, date, scaling_factor_icu); + auto result_6_groups = mio::osecir::details::set_divi_data(models_6_groups, divi_data_path, regions, + date, scaling_factor_icu); + auto result_1_group = mio::osecir::details::set_divi_data(models_1_group, divi_data_path, regions, date, + scaling_factor_icu); EXPECT_THAT(result_6_groups, IsSuccess()); EXPECT_THAT(result_1_group, IsSuccess()); diff --git a/cpp/tests/test_odesecirts.cpp b/cpp/tests/test_odesecirts.cpp index 4bf0d0e1a8..33cc80a694 100755 --- a/cpp/tests/test_odesecirts.cpp +++ b/cpp/tests/test_odesecirts.cpp @@ -141,7 +141,7 @@ TEST(TestOdeSECIRTS, get_flows) TEST(TestOdeSECIRTS, Simulation) { - auto sim = mio::osecirts::Simulation(osecirts_testing_model(), 0, 1); + auto sim = mio::osecirts::Simulation(osecirts_testing_model(), 0, 1); sim.set_integrator_core(std::make_unique>()); sim.advance(1); @@ -258,7 +258,8 @@ TEST(TestOdeSECIRTS, overflow_vaccinations) } // simulate one step with explicit Euler - auto result = mio::simulate_flows>(t0, tmax, dt, model, std::make_unique>()); + auto result = mio::simulate_flows>( + t0, tmax, dt, model, std::make_unique>()); // get the flow indices for each type of vaccination and also the indices of the susceptible compartments auto flow_indx_partial_vaccination = @@ -731,14 +732,14 @@ TEST(TestOdeSECIRTS, read_confirmed_cases) num_InfectedSevere[0] = std::vector(num_age_groups, 0.0); num_icu[0] = std::vector(num_age_groups, 0.0); - ASSERT_THAT(mio::osecirts::details::read_confirmed_cases_data( + ASSERT_THAT(mio::osecirts::details::read_confirmed_cases_data( path, region, {2020, 12, 01}, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_timm, model, std::vector(size_t(num_age_groups), 1.0), 0), IsSuccess()); // read again with invalid date - ASSERT_THAT(mio::osecirts::details::read_confirmed_cases_data( + ASSERT_THAT(mio::osecirts::details::read_confirmed_cases_data( path, region, {3020, 12, 01}, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_timm, model, std::vector(size_t(num_age_groups), 1.0), 0), @@ -746,7 +747,7 @@ TEST(TestOdeSECIRTS, read_confirmed_cases) // call the compute function with empty case data const std::vector empty_case_data; - ASSERT_THAT(mio::osecirts::details::compute_confirmed_cases_data( + ASSERT_THAT(mio::osecirts::details::compute_confirmed_cases_data( empty_case_data, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_timm, region, {2020, 12, 01}, model, std::vector(size_t(num_age_groups), 1.0), 0), @@ -761,7 +762,8 @@ TEST(TestOdeSECIRTS, set_divi_data_invalid_dates) auto model_vector = std::vector>{model}; // Test with date before DIVI dataset was available. - EXPECT_THAT(mio::osecirts::details::set_divi_data(model_vector, "", {1001}, {2019, 12, 01}, 1.0), IsSuccess()); + EXPECT_THAT(mio::osecirts::details::set_divi_data(model_vector, "", {1001}, {2019, 12, 01}, 1.0), + IsSuccess()); // Assure that populations is the same as before. EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), MatrixNear(print_wrap(model.populations.array().cast()), 1e-10, 1e-10)); @@ -808,8 +810,8 @@ TEST(TestOdeSECIRTS, set_confirmed_cases_data_with_ICU) auto model_vector = std::vector>{model}; auto scaling_factor_inf = std::vector(size_t(model.parameters.get_num_groups()), 1.0); - EXPECT_THAT(mio::osecirts::details::set_confirmed_cases_data(model_vector, case_data, {1002}, mid_day, - scaling_factor_inf, immunity_population), + EXPECT_THAT(mio::osecirts::details::set_confirmed_cases_data(model_vector, case_data, {1002}, mid_day, + scaling_factor_inf, immunity_population), IsSuccess()); // Since, TimeInfectedCritical is 1, the number of ICU cases is the difference of confirmed cases between two days, which is 1. @@ -848,17 +850,17 @@ TEST(TestOdeSECIRTS, read_data) {0.61, 0.61, 0.62, 0.62, 0.58, 0.41}, {0.35, 0.35, 0.305, 0.3, 0.385, 0.58}}; - auto read_result1 = mio::osecirts::read_input_data_county(model1, {2020, 12, 01}, {1002}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - TEST_GERMANY_PYDATA_DIR, 10, immunity_population); + auto read_result1 = mio::osecirts::read_input_data_county( + model1, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), 1.0, TEST_GERMANY_PYDATA_DIR, + 10, immunity_population); - auto read_result2 = - mio::osecirts::read_input_data(model2, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), - 1.0, TEST_GERMANY_PYDATA_DIR, 10, immunity_population); + auto read_result2 = mio::osecirts::read_input_data( + model2, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), 1.0, TEST_GERMANY_PYDATA_DIR, + 10, immunity_population); - auto read_result_district = - mio::osecirts::read_input_data(model3, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), - 1.0, mio::path_join(TEST_DATA_DIR, "District/pydata"), 10, immunity_population); + auto read_result_district = mio::osecirts::read_input_data( + model3, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), 1.0, + mio::path_join(TEST_DATA_DIR, "District/pydata"), 10, immunity_population); ASSERT_THAT(read_result1, IsSuccess()); ASSERT_THAT(read_result2, IsSuccess()); @@ -1010,7 +1012,7 @@ TEST(TestOdeSECIRTS, export_time_series_init) {0.35, 0.35, 0.305, 0.3, 0.385, 0.58}}; // Test exporting time series - ASSERT_THAT(mio::osecirts::export_input_data_county_timeseries( + ASSERT_THAT(mio::osecirts::export_input_data_county_timeseries( std::vector>{model}, tmp_results_dir, {0}, {2020, 12, 01}, std::vector(size_t(num_age_groups), 1.0), 1.0, 2, mio::path_join(TEST_DATA_DIR, "county_divi_ma7.json"), @@ -1019,13 +1021,14 @@ TEST(TestOdeSECIRTS, export_time_series_init) mio::path_join(TEST_DATA_DIR, "vacc_county_ageinf_ma7.json")), IsSuccess()); - auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); + auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); ASSERT_THAT(data_extrapolated, IsSuccess()); // Values were generated by the tested function export_input_data_county_timeseries; // can only check stability of the implementation, not correctness auto expected_results = - mio::read_result(mio::path_join(TEST_DATA_DIR, "export_time_series_initialization_osecirts.h5")).value(); + mio::read_result(mio::path_join(TEST_DATA_DIR, "export_time_series_initialization_osecirts.h5")) + .value(); ASSERT_THAT(print_wrap(data_extrapolated.value()[0].get_groups().matrix()), MatrixNear(print_wrap(expected_results[0].get_groups().matrix()), 1e-5, 1e-5)); @@ -1043,9 +1046,9 @@ TEST(TestOdeSECIRTS, model_initialization) {0.61, 0.61, 0.62, 0.62, 0.58, 0.41}, {0.35, 0.35, 0.305, 0.3, 0.385, 0.58}}; - ASSERT_THAT(mio::osecirts::read_input_data_county(model_vector, {2020, 12, 01}, {0}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - TEST_DATA_DIR, 2, immunity_population, false), + ASSERT_THAT(mio::osecirts::read_input_data_county(model_vector, {2020, 12, 01}, {0}, + std::vector(size_t(num_age_groups), 1.0), 1.0, + TEST_DATA_DIR, 2, immunity_population, false), IsSuccess()); // Values from data/export_time_series_init_osecirts.h5, for reading in comparison diff --git a/cpp/tests/test_odesecirvvs.cpp b/cpp/tests/test_odesecirvvs.cpp index 838b15ed2b..c7fc249d7c 100755 --- a/cpp/tests/test_odesecirvvs.cpp +++ b/cpp/tests/test_odesecirvvs.cpp @@ -563,7 +563,7 @@ TEST(TestOdeSECIRVVS, read_confirmed_cases) model[0].parameters.template get>()[(mio::AgeGroup)group]); } - auto read = mio::osecirvvs::details::read_confirmed_cases_data( + auto read = mio::osecirvvs::details::read_confirmed_cases_data( path, region, {2020, 12, 01}, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec, t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, std::vector(size_t(num_age_groups), 1.0)); @@ -579,7 +579,8 @@ TEST(TestOdeSECIRVVS, set_divi_data_invalid_dates) auto model_vector = std::vector>{model}; // Test with date before DIVI dataset was available. - EXPECT_THAT(mio::osecirvvs::details::set_divi_data(model_vector, "", {1001}, {2019, 12, 01}, 1.0), IsSuccess()); + EXPECT_THAT(mio::osecirvvs::details::set_divi_data(model_vector, "", {1001}, {2019, 12, 01}, 1.0), + IsSuccess()); // Assure that populations is the same as before. EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), MatrixNear(print_wrap(model.populations.array().cast()), 1e-10, 1e-10)); @@ -622,9 +623,9 @@ TEST(TestOdeSECIRVVS, set_confirmed_cases_data_with_ICU) // calculate ICU values using set_confirmed_cases_data auto model_vector = std::vector>{model}; auto scaling_factor_inf = std::vector(size_t(model.parameters.get_num_groups()), 1.0); - EXPECT_THAT( - mio::osecirvvs::details::set_confirmed_cases_data(model_vector, case_data, {1002}, mid_day, scaling_factor_inf), - IsSuccess()); + EXPECT_THAT(mio::osecirvvs::details::set_confirmed_cases_data(model_vector, case_data, {1002}, mid_day, + scaling_factor_inf), + IsSuccess()); // Since, TimeInfectedCritical is 1, the number of ICU cases is the difference of confirmed cases between two days, which is 1. // We only have an entry for age group 2. All other age groups should be zero. @@ -660,15 +661,15 @@ TEST(TestOdeSECIRVVS, read_data) const auto pydata_dir_District = mio::path_join(TEST_DATA_DIR, "District", "pydata"); - auto read_result1 = mio::osecirvvs::read_input_data_county(model1, {2020, 12, 01}, {1002}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - TEST_GERMANY_PYDATA_DIR, 30); + auto read_result1 = mio::osecirvvs::read_input_data_county( + model1, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), 1.0, TEST_GERMANY_PYDATA_DIR, + 30); - auto read_result2 = mio::osecirvvs::read_input_data(model2, {2020, 12, 01}, {1002}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - TEST_GERMANY_PYDATA_DIR, 30); + auto read_result2 = mio::osecirvvs::read_input_data(model2, {2020, 12, 01}, {1002}, + std::vector(size_t(num_age_groups), 1.0), + 1.0, TEST_GERMANY_PYDATA_DIR, 30); - auto read_result_district = mio::osecirvvs::read_input_data( + auto read_result_district = mio::osecirvvs::read_input_data( model3, {2020, 12, 01}, {1002}, std::vector(size_t(num_age_groups), 1.0), 1.0, pydata_dir_District, 30); ASSERT_THAT(read_result1, IsSuccess()); @@ -821,7 +822,7 @@ TEST(TestOdeSECIRVVS, export_time_series_init) auto model = make_model(num_age_groups); // Test exporting time series - ASSERT_THAT(mio::osecirvvs::export_input_data_county_timeseries( + ASSERT_THAT(mio::osecirvvs::export_input_data_county_timeseries( std::vector>{model}, tmp_results_dir, {0}, {2020, 12, 01}, std::vector(size_t(num_age_groups), 1.0), 1.0, 2, mio::path_join(TEST_DATA_DIR, "county_divi_ma7.json"), @@ -830,13 +831,14 @@ TEST(TestOdeSECIRVVS, export_time_series_init) mio::path_join(TEST_DATA_DIR, "vacc_county_ageinf_ma7.json")), IsSuccess()); - auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); + auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); ASSERT_THAT(data_extrapolated, IsSuccess()); // Values were generated by the tested function export_input_data_county_timeseries; // can only check stability of the implementation, not correctness auto expected_results = - mio::read_result(mio::path_join(TEST_DATA_DIR, "export_time_series_initialization_osecirvvs.h5")).value(); + mio::read_result(mio::path_join(TEST_DATA_DIR, "export_time_series_initialization_osecirvvs.h5")) + .value(); ASSERT_THAT(print_wrap(data_extrapolated.value()[0].get_groups().matrix()), MatrixNear(print_wrap(expected_results[0].get_groups().matrix()), 1e-5, 1e-5)); @@ -859,7 +861,7 @@ TEST(TestOdeSECIRVVS, export_time_series_init_old_date) model.populations.array().setConstant(0.0); // Test exporting time series - ASSERT_THAT(mio::osecirvvs::export_input_data_county_timeseries( + ASSERT_THAT(mio::osecirvvs::export_input_data_county_timeseries( std::vector>{model}, tmp_results_dir, {0}, {20, 12, 01}, std::vector(size_t(num_age_groups), 1.0), 1.0, 0, mio::path_join(TEST_DATA_DIR, "county_divi_ma7.json"), @@ -868,7 +870,7 @@ TEST(TestOdeSECIRVVS, export_time_series_init_old_date) mio::path_join(TEST_DATA_DIR, "vacc_county_ageinf_ma7.json")), IsSuccess()); - auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); + auto data_extrapolated = mio::read_result(mio::path_join(tmp_results_dir, "Results_rki.h5")); ASSERT_THAT(data_extrapolated, IsSuccess()); auto results_extrapolated = data_extrapolated.value()[0].get_groups().get_value(0); @@ -876,7 +878,7 @@ TEST(TestOdeSECIRVVS, export_time_series_init_old_date) // read population data std::string path = mio::path_join(TEST_DATA_DIR, "county_current_population.json"); const std::vector region{0}; - auto population_data = mio::read_population_data(path, region).value(); + auto population_data = mio::read_population_data(path, region).value(); // So, the expected values are the population data in the susceptible compartments and zeros in the other compartments. for (auto i = 0; i < num_age_groups; i++) { @@ -897,9 +899,9 @@ TEST(TestOdeSECIRVVS, model_initialization) // Vector assignment necessary as read_input_data_county changes model auto model_vector = std::vector>{model}; - ASSERT_THAT(mio::osecirvvs::read_input_data_county(model_vector, {2020, 12, 01}, {0}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - TEST_GERMANY_PYDATA_DIR, 30, false), + ASSERT_THAT(mio::osecirvvs::read_input_data_county(model_vector, {2020, 12, 01}, {0}, + std::vector(size_t(num_age_groups), 1.0), + 1.0, TEST_GERMANY_PYDATA_DIR, 30, false), IsSuccess()); // Values from data/export_time_series_init_osecirvvs.h5, for reading in comparison @@ -938,16 +940,16 @@ TEST(TestOdeSECIRVVS, model_initialization_old_date) auto model_vector = std::vector>{model}; - ASSERT_THAT(mio::osecirvvs::read_input_data(model_vector, {100, 12, 01}, {0}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - TEST_GERMANY_PYDATA_DIR, 0, false), + ASSERT_THAT(mio::osecirvvs::read_input_data(model_vector, {100, 12, 01}, {0}, + std::vector(size_t(num_age_groups), 1.0), 1.0, + TEST_GERMANY_PYDATA_DIR, 0, false), IsSuccess()); // if we enter an old date, the model only should be initialized with the population data. // read population data std::string path = mio::path_join(TEST_DATA_DIR, "county_current_population.json"); const std::vector region{0}; - auto population_data = mio::read_population_data(path, region).value(); + auto population_data = mio::read_population_data(path, region).value(); // So, the expected values are the population data in the susceptible compartments and zeros in the other compartments. for (auto i = 0; i < num_age_groups; i++) { @@ -975,16 +977,16 @@ TEST(TestOdeSECIRVVS, model_initialization_old_date_county) auto model_vector = std::vector>{model}; - ASSERT_THAT(mio::osecirvvs::read_input_data_county(model_vector, {100, 12, 01}, {0}, - std::vector(size_t(num_age_groups), 1.0), 1.0, - TEST_GERMANY_PYDATA_DIR, 0, false), + ASSERT_THAT(mio::osecirvvs::read_input_data_county(model_vector, {100, 12, 01}, {0}, + std::vector(size_t(num_age_groups), 1.0), + 1.0, TEST_GERMANY_PYDATA_DIR, 0, false), IsSuccess()); // if we enter an old date, the model only should be initialized with the population data. // read population data std::string path = mio::path_join(TEST_DATA_DIR, "county_current_population.json"); const std::vector region{0}; - auto population_data = mio::read_population_data(path, region).value(); + auto population_data = mio::read_population_data(path, region).value(); // So, the expected values are the population data in the susceptible compartments and zeros in the other compartments. for (auto i = 0; i < num_age_groups; i++) { @@ -1017,10 +1019,10 @@ TEST(TestOdeSECIRVVS, set_population_data_overflow_vacc) std::string path_pop_data = mio::path_join(TEST_DATA_DIR, "county_current_population.json"); const std::vector region{0}; - auto population_data = mio::read_population_data(path_pop_data, region).value(); + auto population_data = mio::read_population_data(path_pop_data, region).value(); // we choose the date so that no case data is available - ASSERT_THAT(mio::osecirvvs::details::set_population_data( + ASSERT_THAT(mio::osecirvvs::details::set_population_data( model_vector, path_pop_data, mio::path_join(TEST_DATA_DIR, "cases_all_county_age_ma7.json"), {0}, {1000, 12, 01}), IsSuccess()); @@ -1057,10 +1059,10 @@ TEST(TestOdeSECIRVVS, set_population_data_no_data_avail) std::string path_pop_data = mio::path_join(TEST_DATA_DIR, "county_current_population.json"); const std::vector region{0}; - auto population_data = mio::read_population_data(path_pop_data, region).value(); + auto population_data = mio::read_population_data(path_pop_data, region).value(); // we choose the date so that no case data is available - ASSERT_THAT(mio::osecirvvs::details::set_population_data( + ASSERT_THAT(mio::osecirvvs::details::set_population_data( model_vector, path_pop_data, mio::path_join(TEST_DATA_DIR, "cases_all_county_age_ma7.json"), {200}, {1000, 12, 01}), IsSuccess()); @@ -1079,7 +1081,7 @@ TEST(TestOdeSECIRVVS, run_simulation) // Load result of a previous run; only tests stability, not correctness. auto expected_result = - mio::read_result(mio::path_join(TEST_DATA_DIR, "results_osecirvvs.h5")).value()[0].get_groups(); + mio::read_result(mio::path_join(TEST_DATA_DIR, "results_osecirvvs.h5")).value()[0].get_groups(); ASSERT_THAT(print_wrap(result.matrix()), MatrixNear(print_wrap(expected_result.matrix()), 1e-5, 1e-5)); } diff --git a/cpp/tests/test_save_parameters.cpp b/cpp/tests/test_save_parameters.cpp index e08324a9ed..86b3c9f7ef 100644 --- a/cpp/tests/test_save_parameters.cpp +++ b/cpp/tests/test_save_parameters.cpp @@ -256,7 +256,7 @@ TEST(TestSaveParameters, read_graph_without_edges) std::vector ids = {0, 1}; auto graph_no_edges = mio::create_graph_without_edges, mio::MobilityParameters>(models, ids); - auto write_status = mio::write_graph(graph_no_edges, tmp_results_dir); + auto write_status = mio::write_graph(graph_no_edges, tmp_results_dir); ASSERT_THAT(print_wrap(write_status), IsSuccess()); auto read_graph = @@ -433,7 +433,7 @@ TEST(TestSaveParameters, json_graphs_write_read_compare) TempFileRegister file_register; auto graph_dir = file_register.get_unique_path("graph_parameters-%%%%-%%%%"); - auto write_status = mio::write_graph(graph, graph_dir); + auto write_status = mio::write_graph(graph, graph_dir); ASSERT_THAT(print_wrap(write_status), IsSuccess()); auto read_result = mio::read_graph>(graph_dir); @@ -642,7 +642,8 @@ TEST(TestSaveParameters, ReadPopulationDataRKIAges) model[0].parameters.get>()[group] = 0.12 * ((size_t)group + 1); } - auto read_result = mio::osecir::read_input_data_germany(model, date, scaling_factor_inf, scaling_factor_icu, path); + auto read_result = + mio::osecir::read_input_data_germany(model, date, scaling_factor_inf, scaling_factor_icu, path); ASSERT_THAT(print_wrap(read_result), IsSuccess()); std::vector sus = {3444279.56, 7666866.07, 18801541.51, 29518513.27, 16321649.37, 6072737.25}; @@ -752,11 +753,11 @@ TEST(TestSaveParameters, ReadPopulationDataCountyAllAges) model3[0].parameters.get>()[group] = 0.12 * ((size_t)group + 1); } - auto read_result1 = mio::osecir::read_input_data_county(model1, date, county, scaling_factor_inf, - scaling_factor_icu, TEST_GERMANY_PYDATA_DIR); - auto read_result2 = mio::osecir::read_input_data(model2, date, county, scaling_factor_inf, scaling_factor_icu, - TEST_GERMANY_PYDATA_DIR); - auto read_result_district = mio::osecir::read_input_data( + auto read_result1 = mio::osecir::read_input_data_county(model1, date, county, scaling_factor_inf, + scaling_factor_icu, TEST_GERMANY_PYDATA_DIR); + auto read_result2 = mio::osecir::read_input_data(model2, date, county, scaling_factor_inf, + scaling_factor_icu, TEST_GERMANY_PYDATA_DIR); + auto read_result_district = mio::osecir::read_input_data( model3, date, county, scaling_factor_inf, scaling_factor_icu, mio::path_join(TEST_DATA_DIR, "District/pydata")); ASSERT_THAT(print_wrap(read_result1), IsSuccess()); ASSERT_THAT(print_wrap(read_result2), IsSuccess()); @@ -839,14 +840,14 @@ TEST(TestSaveParameters, ExtrapolateRKI) TempFileRegister file_register; auto results_dir = file_register.get_unique_path("ExtrapolateRKI-%%%%-%%%%"); boost::filesystem::create_directory(results_dir); - auto extrapolate_result = mio::osecir::export_input_data_county_timeseries( + auto extrapolate_result = mio::osecir::export_input_data_county_timeseries( model, results_dir, county, date, scaling_factor_inf, scaling_factor_icu, 1, mio::path_join(TEST_DATA_DIR, "county_divi_ma7.json"), mio::path_join(TEST_DATA_DIR, "cases_all_county_age_ma7.json"), mio::path_join(TEST_DATA_DIR, "county_current_population.json")); ASSERT_THAT(print_wrap(extrapolate_result), IsSuccess()); - auto read_result = mio::read_result(mio::path_join(results_dir, "Results_rki.h5")); + auto read_result = mio::read_result(mio::path_join(results_dir, "Results_rki.h5")); ASSERT_THAT(print_wrap(read_result), IsSuccess()); auto& file_results = read_result.value(); auto results = file_results[0].get_groups(); diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index a15c62d419..804a735f1a 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -81,11 +81,12 @@ TEST(TestSaveResult, save_result) std::vector ids = {1, 2}; TempFileRegister file_register; - auto results_file_path = file_register.get_unique_path("test_result-%%%%-%%%%.h5"); - auto save_result_status = mio::save_result(results_from_sim, ids, (int)(size_t)nb_groups, results_file_path); + auto results_file_path = file_register.get_unique_path("test_result-%%%%-%%%%.h5"); + auto save_result_status = + mio::save_result(results_from_sim, ids, (int)(size_t)nb_groups, results_file_path); ASSERT_TRUE(save_result_status); - auto results_from_file = mio::read_result(results_file_path); + auto results_from_file = mio::read_result(results_file_path); ASSERT_TRUE(results_from_file); auto result_from_file = results_from_file.value()[0]; @@ -196,13 +197,13 @@ TEST(TestSaveResult, save_result_with_params) return node.property.get_simulation().get_model(); }); - save_result_status = save_result_with_params(ensemble_results.back(), ensemble_params.back(), {0, 1}, - tmp_results_dir, run_idx); + save_result_status = save_result_with_params(ensemble_results.back(), ensemble_params.back(), + {0, 1}, tmp_results_dir, run_idx); return 0; //function needs to return something }); ASSERT_TRUE(save_result_status); - auto results_from_file = mio::read_result(tmp_results_dir + "/run0/Result.h5"); + auto results_from_file = mio::read_result(tmp_results_dir + "/run0/Result.h5"); ASSERT_TRUE(results_from_file); auto result_from_file = results_from_file.value()[0]; EXPECT_EQ(ensemble_results.back().back().get_num_elements(), result_from_file.get_groups().get_num_elements()); @@ -340,19 +341,19 @@ TEST(TestSaveResult, save_percentiles_and_sums) return 0; //function needs to return something }); - auto save_results_status = save_results(ensemble_results, ensemble_params, {0, 1}, tmp_results_dir); + auto save_results_status = save_results(ensemble_results, ensemble_params, {0, 1}, tmp_results_dir); ASSERT_TRUE(save_results_status); // test percentiles - auto results_from_file_p05 = mio::read_result(tmp_results_dir + "/p05/Results.h5"); + auto results_from_file_p05 = mio::read_result(tmp_results_dir + "/p05/Results.h5"); ASSERT_TRUE(results_from_file_p05); - auto results_from_file_p25 = mio::read_result(tmp_results_dir + "/p25/Results.h5"); + auto results_from_file_p25 = mio::read_result(tmp_results_dir + "/p25/Results.h5"); ASSERT_TRUE(results_from_file_p25); - auto results_from_file_p50 = mio::read_result(tmp_results_dir + "/p50/Results.h5"); + auto results_from_file_p50 = mio::read_result(tmp_results_dir + "/p50/Results.h5"); ASSERT_TRUE(results_from_file_p50); - auto results_from_file_p75 = mio::read_result(tmp_results_dir + "/p75/Results.h5"); + auto results_from_file_p75 = mio::read_result(tmp_results_dir + "/p75/Results.h5"); ASSERT_TRUE(results_from_file_p75); - auto results_from_file_p95 = mio::read_result(tmp_results_dir + "/p95/Results.h5"); + auto results_from_file_p95 = mio::read_result(tmp_results_dir + "/p95/Results.h5"); ASSERT_TRUE(results_from_file_p95); auto result_from_file = results_from_file_p25.value()[0]; @@ -361,35 +362,35 @@ TEST(TestSaveResult, save_percentiles_and_sums) result_from_file.get_groups().get_num_time_points()); // results_run - auto results_run0 = mio::read_result(tmp_results_dir + "/results_run0.h5"); + auto results_run0 = mio::read_result(tmp_results_dir + "/results_run0.h5"); ASSERT_TRUE(results_run0); - auto results_run0_sum = mio::read_result(tmp_results_dir + "/results_run0_sum.h5"); + auto results_run0_sum = mio::read_result(tmp_results_dir + "/results_run0_sum.h5"); ASSERT_TRUE(results_run0_sum); - auto results_run1 = mio::read_result(tmp_results_dir + "/results_run1.h5"); + auto results_run1 = mio::read_result(tmp_results_dir + "/results_run1.h5"); ASSERT_TRUE(results_run1); - auto results_run1_sum = mio::read_result(tmp_results_dir + "/results_run1_sum.h5"); + auto results_run1_sum = mio::read_result(tmp_results_dir + "/results_run1_sum.h5"); ASSERT_TRUE(results_run1_sum); - auto results_run2 = mio::read_result(tmp_results_dir + "/results_run2.h5"); + auto results_run2 = mio::read_result(tmp_results_dir + "/results_run2.h5"); ASSERT_TRUE(results_run2); - auto results_run2_sum = mio::read_result(tmp_results_dir + "/results_run2_sum.h5"); + auto results_run2_sum = mio::read_result(tmp_results_dir + "/results_run2_sum.h5"); ASSERT_TRUE(results_run2_sum); // test save edges (percentiles and results from single runs) std::vector> pairs_edges = {{0, 1}}; - auto save_edges_status = save_edges(ensemble_edges, pairs_edges, tmp_results_dir, true, true); + auto save_edges_status = save_edges(ensemble_edges, pairs_edges, tmp_results_dir, true, true); ASSERT_TRUE(save_edges_status); // percentiles - auto results_edges_from_file_p05 = mio::read_result(tmp_results_dir + "/p05/Edges.h5"); + auto results_edges_from_file_p05 = mio::read_result(tmp_results_dir + "/p05/Edges.h5"); ASSERT_TRUE(results_edges_from_file_p05); - auto results_edges_from_file_p25 = mio::read_result(tmp_results_dir + "/p25/Edges.h5"); + auto results_edges_from_file_p25 = mio::read_result(tmp_results_dir + "/p25/Edges.h5"); ASSERT_TRUE(results_edges_from_file_p25); - auto results_edges_from_file_p50 = mio::read_result(tmp_results_dir + "/p50/Edges.h5"); + auto results_edges_from_file_p50 = mio::read_result(tmp_results_dir + "/p50/Edges.h5"); ASSERT_TRUE(results_edges_from_file_p50); - auto results_edges_from_file_p75 = mio::read_result(tmp_results_dir + "/p75/Edges.h5"); + auto results_edges_from_file_p75 = mio::read_result(tmp_results_dir + "/p75/Edges.h5"); ASSERT_TRUE(results_edges_from_file_p75); - auto results_edges_from_file_p95 = mio::read_result(tmp_results_dir + "/p95/Edges.h5"); + auto results_edges_from_file_p95 = mio::read_result(tmp_results_dir + "/p95/Edges.h5"); ASSERT_TRUE(results_edges_from_file_p95); auto result_edges_from_file = results_edges_from_file_p25.value()[0]; @@ -398,11 +399,11 @@ TEST(TestSaveResult, save_percentiles_and_sums) result_edges_from_file.get_groups().get_num_time_points()); // single runs - auto results_edges_run0 = mio::read_result(tmp_results_dir + "/Edges_run0.h5"); + auto results_edges_run0 = mio::read_result(tmp_results_dir + "/Edges_run0.h5"); ASSERT_TRUE(results_edges_run0); - auto results_edges_run1 = mio::read_result(tmp_results_dir + "/Edges_run1.h5"); + auto results_edges_run1 = mio::read_result(tmp_results_dir + "/Edges_run1.h5"); ASSERT_TRUE(results_edges_run1); - auto results_edges_run2 = mio::read_result(tmp_results_dir + "/Edges_run2.h5"); + auto results_edges_run2 = mio::read_result(tmp_results_dir + "/Edges_run2.h5"); ASSERT_TRUE(results_edges_run2); } @@ -428,11 +429,11 @@ TEST(TestSaveResult, save_edges) // save the results to a file TempFileRegister file_register; auto results_file_path = file_register.get_unique_path("test_result-%%%%-%%%%.h5"); - auto save_edges_status = mio::save_edges(results_edges, pairs_edges, results_file_path); + auto save_edges_status = mio::save_edges(results_edges, pairs_edges, results_file_path); ASSERT_TRUE(save_edges_status); // read the results back in and check that they are correct. - auto results_from_file = mio::read_result(results_file_path); + auto results_from_file = mio::read_result(results_file_path); ASSERT_TRUE(results_from_file); // group 0 @@ -495,5 +496,6 @@ TEST(TestSaveEdges, save_edges_empty_ts) auto results_file_path = file_register.get_unique_path("TestEdges-%%%%-%%%%.h5"); // Call the save_edges function and check if it returns a failure - ASSERT_THAT(save_edges(results, pairs_edges, results_file_path), IsFailure(mio::StatusCode::InvalidValue)); + ASSERT_THAT(save_edges(results, pairs_edges, results_file_path), + IsFailure(mio::StatusCode::InvalidValue)); } diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/io/result_io.h b/pycode/memilio-simulation/memilio/simulation/bindings/io/result_io.h index b34a06c255..73e79a1967 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/io/result_io.h +++ b/pycode/memilio-simulation/memilio/simulation/bindings/io/result_io.h @@ -39,8 +39,8 @@ void bind_save_results(pybind11::module_& m) const std::vector>& ensemble_params, const std::vector& county_ids, const std::string& result_dir, bool save_single_runs, bool save_percentiles) { boost::filesystem::path dir(result_dir); - auto ioresult = mio::save_results(ensemble_results, ensemble_params, county_ids, dir, save_single_runs, - save_percentiles); + auto ioresult = mio::save_results(ensemble_results, ensemble_params, county_ids, dir, + save_single_runs, save_percentiles); return NULL; }); } diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp index ed3ad0bf97..0d1d5afcae 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp @@ -259,14 +259,16 @@ PYBIND11_MODULE(_simulation_osecir, m) mio::Graph, mio::MobilityParameters>& params_graph, const std::vector& scaling_factor_inf, double scaling_factor_icu, double tnt_capacity_factor, int num_days = 0, bool export_time_series = false) { - auto result = mio::set_nodes< - double, // FP - mio::osecir::TestAndTraceCapacity, mio::osecir::ContactPatterns, - mio::osecir::Model, mio::MobilityParameters, mio::osecir::Parameters, - decltype(mio::osecir::read_input_data_county>), decltype(mio::get_node_ids)>( - params, start_date, end_date, data_dir, population_data_path, is_node_for_county, params_graph, - mio::osecir::read_input_data_county>, mio::get_node_ids, scaling_factor_inf, - scaling_factor_icu, tnt_capacity_factor, num_days, export_time_series); + auto result = + mio::set_nodes, mio::osecir::ContactPatterns, + mio::osecir::Model, mio::MobilityParameters, + mio::osecir::Parameters, + decltype(mio::osecir::read_input_data_county>), + decltype(mio::get_node_ids)>( + params, start_date, end_date, data_dir, population_data_path, is_node_for_county, params_graph, + mio::osecir::read_input_data_county>, mio::get_node_ids, + scaling_factor_inf, scaling_factor_icu, tnt_capacity_factor, num_days, export_time_series); return pymio::check_and_throw(result); }, py::return_value_policy::move); @@ -307,7 +309,7 @@ PYBIND11_MODULE(_simulation_osecir, m) [](std::vector>& model, mio::Date date, const std::vector& county, const std::vector& scaling_factor_inf, double scaling_factor_icu, const std::string& dir, int num_days = 0, bool export_time_series = false) { - auto result = mio::osecir::read_input_data_county>( + auto result = mio::osecir::read_input_data_county>( model, date, county, scaling_factor_inf, scaling_factor_icu, dir, num_days, export_time_series); return pymio::check_and_throw(result); }, diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp index e14e8b34f4..2b671d349d 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp @@ -312,10 +312,10 @@ PYBIND11_MODULE(_simulation_osecirvvs, m) mio::osecirvvs::TestAndTraceCapacity, mio::osecirvvs::ContactPatterns, mio::osecirvvs::Model, mio::MobilityParameters, mio::osecirvvs::Parameters, - decltype(mio::osecirvvs::read_input_data_county>), + decltype(mio::osecirvvs::read_input_data_county>), decltype(mio::get_node_ids)>( params, start_date, end_date, data_dir, population_data_path, is_node_for_county, params_graph, - mio::osecirvvs::read_input_data_county>, mio::get_node_ids, + mio::osecirvvs::read_input_data_county>, mio::get_node_ids, scaling_factor_inf, scaling_factor_icu, tnt_capacity_factor, num_days, export_time_series); return pymio::check_and_throw(result); }, @@ -360,7 +360,7 @@ PYBIND11_MODULE(_simulation_osecirvvs, m) [](std::vector>& model, mio::Date date, const std::vector& county, const std::vector& scaling_factor_inf, double scaling_factor_icu, const std::string& dir, int num_days = 0, bool export_time_series = false) { - auto result = mio::osecirvvs::read_input_data_county>( + auto result = mio::osecirvvs::read_input_data_county>( model, date, county, scaling_factor_inf, scaling_factor_icu, dir, num_days, export_time_series); return pymio::check_and_throw(result); }, diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/simulation.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/simulation.cpp index a98012dce8..3a9e24f4c4 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/simulation.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/simulation.cpp @@ -59,7 +59,7 @@ PYBIND11_MODULE(_simulation, m) pymio::bind_CustomIndexArray, mio::AgeGroup>(m, "AgeGroupArray"); pymio::bind_class>(m, "AgeGroup") .def(py::init()); - + pymio::bind_CustomIndexArray(m, "AgeGroupSimulationDayArray"); pymio::bind_class>(m, "SimulationDay") @@ -70,7 +70,8 @@ PYBIND11_MODULE(_simulation, m) auto damping_class = pymio::bind_class, pymio::EnablePickling::Required>(m, "Damping"); pymio::bind_damping_members(damping_class); - auto dampings_class = pymio::bind_class, pymio::EnablePickling::Required>(m, "Dampings"); + auto dampings_class = + pymio::bind_class, pymio::EnablePickling::Required>(m, "Dampings"); pymio::bind_dampings_members(dampings_class); pymio::bind_time_series(m, "TimeSeries"); @@ -99,12 +100,13 @@ PYBIND11_MODULE(_simulation, m) pymio::bind_class, pymio::EnablePickling::Required>(m, "MobilityDampings"); pymio::bind_dampings_members(mobility_dampings_class); - auto mobility_coeffs_class = - pymio::bind_class, pymio::EnablePickling::Required>(m, "MobilityCoefficients"); + auto mobility_coeffs_class = pymio::bind_class, pymio::EnablePickling::Required>( + m, "MobilityCoefficients"); pymio::bind_damping_expression_members(mobility_coeffs_class); - auto mobility_coeff_group_class = pymio::bind_class, pymio::EnablePickling::Required>( - m, "MobilityCoefficientGroup"); + auto mobility_coeff_group_class = + pymio::bind_class, pymio::EnablePickling::Required>( + m, "MobilityCoefficientGroup"); pymio::bind_damping_expression_group_members(mobility_coeff_group_class); pymio::bind_dynamicNPI_members(m, "DynamicNPIs"); @@ -127,12 +129,12 @@ PYBIND11_MODULE(_simulation, m) return std::vector>(h.begin(), h.end()); }, py::arg("state_id"), py::arg("start_date") = mio::Date(std::numeric_limits::min(), 1, 1), - py::arg("end_date") = mio::Date(std::numeric_limits::max(), 1, 1)); + py::arg("end_date") = mio::Date(std::numeric_limits::max(), 1, 1)); m.def( "read_mobility_plain", [](const std::string& filename) { - auto result = mio::read_mobility_plain(filename); + auto result = mio::read_mobility_plain(filename); return pymio::check_and_throw(result); }, py::return_value_policy::move); From ccddcd1963398a672256e6fd0855af0728e53c29 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Tue, 16 Sep 2025 08:57:02 +0000 Subject: [PATCH 2/7] Safer FP to int conversion --- cpp/models/ode_secirts/parameters_io.h | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/cpp/models/ode_secirts/parameters_io.h b/cpp/models/ode_secirts/parameters_io.h index 5af585c396..c026f56414 100644 --- a/cpp/models/ode_secirts/parameters_io.h +++ b/cpp/models/ode_secirts/parameters_io.h @@ -743,12 +743,13 @@ IOResult set_vaccination_data(std::vector>& model, const std::ve auto num_groups = model[0].parameters.get_num_groups(); - auto days_until_effective_n = - (int)(floor(model[0].parameters.template get>()[AgeGroup(0)])); - auto days_until_effective_pi = - (int)(floor(model[0].parameters.template get>()[AgeGroup(0)])); - auto days_until_effective_ii = - (int)(floor(model[0].parameters.template get>()[AgeGroup(0)])); + int days_until_effective_n = static_cast(floor( + static_cast(model[0].parameters.template get>()[AgeGroup(0)]))); + int days_until_effective_pi = static_cast(floor( + static_cast(model[0].parameters.template get>()[AgeGroup(0)]))); + int days_until_effective_ii = static_cast( + floor(static_cast(model[0].parameters.template get>()[AgeGroup(0)]))); + // iterate over regions (e.g., counties) for (size_t i = 0; i < model.size(); ++i) { // iterate over age groups in region From 5641fbd14d222c274384ad23d744170bb7432b6a Mon Sep 17 00:00:00 2001 From: julianlitz Date: Tue, 16 Sep 2025 09:49:24 +0000 Subject: [PATCH 3/7] std::round and static_cast --- cpp/models/ode_secir/parameters_io.h | 22 ++++++++++-------- cpp/models/ode_secirts/parameters_io.h | 32 +++++++++++++------------- 2 files changed, 28 insertions(+), 26 deletions(-) diff --git a/cpp/models/ode_secir/parameters_io.h b/cpp/models/ode_secir/parameters_io.h index fc71756de9..af5b2c41b5 100644 --- a/cpp/models/ode_secir/parameters_io.h +++ b/cpp/models/ode_secir/parameters_io.h @@ -219,6 +219,8 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect const std::vector& region, Date date, const std::vector& scaling_factor_inf) { + using std::round; + const size_t num_age_groups = ConfirmedCasesDataEntry::age_group_names.size(); // allow single scalar scaling that is broadcast to all age groups assert(scaling_factor_inf.size() == 1 || scaling_factor_inf.size() == num_age_groups); @@ -251,16 +253,16 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect // reuse group 0 parameters for all RKI age groups const size_t group = (model_groups == num_age_groups) ? ag : 0; - t_Exposed[node].push_back( - static_cast(std::round(model[node].parameters.template get>()[(AgeGroup)group]))); - t_InfectedNoSymptoms[node].push_back(static_cast( - std::round(model[node].parameters.template get>()[(AgeGroup)group]))); - t_InfectedSymptoms[node].push_back(static_cast( - std::round(model[node].parameters.template get>()[(AgeGroup)group]))); - t_InfectedSevere[node].push_back(static_cast( - std::round(model[node].parameters.template get>()[(AgeGroup)group]))); - t_InfectedCritical[node].push_back(static_cast( - std::round(model[node].parameters.template get>()[(AgeGroup)group]))); + t_Exposed[node].push_back(static_cast( + round(static_cast(model[node].parameters.template get>()[(AgeGroup)group])))); + t_InfectedNoSymptoms[node].push_back(static_cast(round( + static_cast(model[node].parameters.template get>()[(AgeGroup)group])))); + t_InfectedSymptoms[node].push_back(static_cast(round( + static_cast(model[node].parameters.template get>()[(AgeGroup)group])))); + t_InfectedSevere[node].push_back(static_cast(round( + static_cast(model[node].parameters.template get>()[(AgeGroup)group])))); + t_InfectedCritical[node].push_back(static_cast(round( + static_cast(model[node].parameters.template get>()[(AgeGroup)group])))); mu_C_R[node].push_back( model[node].parameters.template get>()[(AgeGroup)group]); diff --git a/cpp/models/ode_secirts/parameters_io.h b/cpp/models/ode_secirts/parameters_io.h index c026f56414..d9295f36a7 100644 --- a/cpp/models/ode_secirts/parameters_io.h +++ b/cpp/models/ode_secirts/parameters_io.h @@ -120,18 +120,18 @@ IOResult compute_confirmed_cases_data( auto age = (size_t)entry.age_group; // (rounded) transition times - const int t_exposed = - static_cast(round(params_region.template get>()[entry.age_group])); - int t_InfectedNoSymptoms = - static_cast(round(params_region.template get>()[entry.age_group])); - int t_InfectedSymptoms = - static_cast(round(params_region.template get>()[entry.age_group])); - const int t_InfectedSevere = - static_cast(round(params_region.template get>()[entry.age_group])); - const int t_InfectedCritical = - static_cast(round(params_region.template get>()[entry.age_group])); - const int t_imm_interval_i = - static_cast(round(params_region.template get>()[entry.age_group])); + const int t_exposed = static_cast( + round(static_cast(params_region.template get>()[entry.age_group]))); + int t_InfectedNoSymptoms = static_cast( + round(static_cast(params_region.template get>()[entry.age_group]))); + int t_InfectedSymptoms = static_cast( + round(static_cast(params_region.template get>()[entry.age_group]))); + const int t_InfectedSevere = static_cast( + round(static_cast(params_region.template get>()[entry.age_group]))); + const int t_InfectedCritical = static_cast( + round(static_cast(params_region.template get>()[entry.age_group]))); + const int t_imm_interval_i = static_cast( + round(static_cast(params_region.template get>()[entry.age_group]))); // transition probabilities FP recoveredPerInfectedNoSymptoms = @@ -141,10 +141,10 @@ IOResult compute_confirmed_cases_data( // if we select a layer with better immunity (layer > 0), we need to adjust the times and transition rates if (layer > 0) { - t_InfectedNoSymptoms = static_cast(round( - t_InfectedNoSymptoms * params_region.template get>()[entry.age_group])); - t_InfectedSymptoms = static_cast(round( - t_InfectedSymptoms * params_region.template get>()[entry.age_group])); + t_InfectedNoSymptoms = static_cast(round(static_cast( + t_InfectedNoSymptoms * params_region.template get>()[entry.age_group]))); + t_InfectedSymptoms = static_cast(round(static_cast( + t_InfectedSymptoms * params_region.template get>()[entry.age_group]))); const FP red_fact_exp = layer == 1 ? params_region.template get>()[entry.age_group] From 052e31f242296554652f625bc5a0878f7d89330f Mon Sep 17 00:00:00 2001 From: julianlitz Date: Fri, 26 Sep 2025 14:17:57 +0000 Subject: [PATCH 4/7] Fix FP metaprogramming madness. --- .../mobility/metapopulation_mobility_instant.h | 14 +++++++------- cpp/models/ode_secir/model.h | 4 ++-- cpp/models/ode_secirts/model.h | 4 ++-- cpp/models/ode_secirvvs/model.h | 4 ++-- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 4a01e1392e..5be92d3f50 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -478,7 +478,7 @@ void calculate_mobility_returns(Eigen::Ref::Vector> mobi * detect a get_infections_relative function for the Model type. */ template -using get_infections_relative_expr_t = decltype(get_infections_relative( +using get_infections_relative_expr_t = decltype(get_infections_relative( std::declval(), std::declval(), std::declval>&>())); /** @@ -509,7 +509,7 @@ FP get_infections_relative(const SimulationNode& node, FP t, const Eige * detect a get_mobility_factors function for the Model type. */ template -using get_mobility_factors_expr_t = decltype(get_mobility_factors( +using get_mobility_factors_expr_t = decltype(get_mobility_factors( std::declval(), std::declval(), std::declval>&>())); /** @@ -541,8 +541,8 @@ auto get_mobility_factors(const SimulationNode& node, FP t, const Eigen * detect a get_mobility_factors function for the Model type. */ template -using test_commuters_expr_t = decltype(test_commuters( - std::declval(), std::declval>&>(), std::declval())); +using test_commuters_expr_t = decltype(test_commuters( + std::declval(), std::declval>>(), std::declval())); /** * Test persons when moving from their source node. @@ -581,7 +581,7 @@ void mio::MobilityEdge::apply_mobility(FP t, FP dt, SimulationNode& if (dyn_npis.get_thresholds().size() > 0 && floating_point_greater_equal(t, m_t_last_dynamic_npi_check + dyn_npis.get_interval().get())) { auto inf_rel = - get_infections_relative(node_from, t, node_from.get_last_state()) * dyn_npis.get_base_value(); + get_infections_relative(node_from, t, node_from.get_last_state()) * dyn_npis.get_base_value(); auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); if (exceeded_threshold != dyn_npis.get_thresholds().end() && (exceeded_threshold->first > m_dynamic_npi.first || @@ -602,8 +602,8 @@ void mio::MobilityEdge::apply_mobility(FP t, FP dt, SimulationNode& if (m_return_times.get_time(i) <= t) { auto v0 = find_value_reverse(node_to.get_result(), m_mobile_population.get_time(i), 1e-10, 1e-10); assert(v0 != node_to.get_result().rend() && "unexpected error."); - calculate_mobility_returns(m_mobile_population[i], node_to.get_simulation(), *v0, - m_mobile_population.get_time(i), dt); + calculate_mobility_returns(m_mobile_population[i], node_to.get_simulation(), *v0, + m_mobile_population.get_time(i), dt); //the lower-order return calculation may in rare cases produce negative compartments, //especially at the beginning of the simulation. diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h index 132fbd545e..f1541acc93 100644 --- a/cpp/models/ode_secir/model.h +++ b/cpp/models/ode_secir/model.h @@ -659,7 +659,7 @@ IOResult get_reproduction_number(FP t_value, const Simulation& sim * @tparam FP floating point type, e.g., double. * @tparam Base simulation type that uses a secir compartment model; see Simulation. */ -template , FP>> +template >> auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) { auto& params = sim.get_model().parameters; @@ -689,7 +689,7 @@ auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen return factors; } -template , FP>> +template >> auto test_commuters(Simulation& sim, Eigen::Ref> mobile_population, FP time) { auto& model = sim.get_model(); diff --git a/cpp/models/ode_secirts/model.h b/cpp/models/ode_secirts/model.h index 2d66b9c165..5ff833bc1f 100644 --- a/cpp/models/ode_secirts/model.h +++ b/cpp/models/ode_secirts/model.h @@ -896,7 +896,7 @@ FP get_infections_relative(const Simulation& sim, FP /*t*/, const Eige * @return vector expression, same size as y, with migration factors per compartment. * @tparam Base simulation type that uses a SECIRS-type compartment model. see Simulation. */ -template , FP>> +template >> auto get_migration_factors(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) { auto& params = sim.get_model().parameters; @@ -943,7 +943,7 @@ auto get_migration_factors(const Simulation& sim, FP /*t*/, const Eigen::R * @param[in,out] migrated Vector representing the number of commuters in each state. * @param[in] time Current simulation time, used to determine the commuter detection period. */ -template , FP>> +template >> auto test_commuters(Simulation& sim, Eigen::Ref> migrated, FP time) { auto& model = sim.get_model(); diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h index 5a30fe330f..780b850c88 100644 --- a/cpp/models/ode_secirvvs/model.h +++ b/cpp/models/ode_secirvvs/model.h @@ -825,7 +825,7 @@ FP get_infections_relative(const Simulation& sim, FP /*t*/, const Eige * @return vector expression, same size as y, with mobility factors per compartment. * @tparam Base simulation type that uses a secir compartment model. see Simulation. */ -template , FP>> +template >> auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) { @@ -866,7 +866,7 @@ auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen::Re return factors; } -template , FP>> +template >> auto test_commuters(Simulation& sim, Eigen::Ref> mobile_population, FP time) { auto& model = sim.get_model(); From 31e299a08cea08bf0304f054103285ea1e4485da Mon Sep 17 00:00:00 2001 From: julianlitz Date: Fri, 26 Sep 2025 19:40:25 +0000 Subject: [PATCH 5/7] FP fixes in Graph_Simulation.h --- cpp/memilio/mobility/graph_simulation.h | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/cpp/memilio/mobility/graph_simulation.h b/cpp/memilio/mobility/graph_simulation.h index 33770f8a47..a88efcfd7f 100644 --- a/cpp/memilio/mobility/graph_simulation.h +++ b/cpp/memilio/mobility/graph_simulation.h @@ -149,7 +149,7 @@ class GraphSimulationStochastic { } - void advance(FP t_max) + void advance(ScalarType t_max) { //draw normalized waiting time ScalarType normalized_waiting_time = ExponentialDistribution::get_instance()(m_rng, 1.0); @@ -296,6 +296,8 @@ class FeedbackGraphSimulation void advance(FP t_max) { + using std::min; + // Initialize global and regional ICU occupancy if not done yet if (!m_initialized) { calculate_global_icu_occupancy(); @@ -304,7 +306,7 @@ class FeedbackGraphSimulation } while (m_t < t_max) { - FP dt_eff = std::min(m_dt, t_max - m_t); + FP dt_eff = min(m_dt, t_max - m_t); // assign the regional and global ICU occupancy to each node distribute_icu_data(); @@ -344,11 +346,11 @@ class FeedbackGraphSimulation auto& first_node_sim = m_graph.nodes()[0].property.get_simulation(); auto& first_node_icu = first_node_sim.get_parameters().template get>(); m_global_icu_occupancy = mio::TimeSeries(first_node_icu.get_num_elements()); - m_global_icu_occupancy.add_time_point(0.0, Eigen::VectorXd::Zero(first_node_icu.get_num_elements())); + m_global_icu_occupancy.add_time_point(0.0, Eigen::VectorX::Zero(first_node_icu.get_num_elements())); FP total_population = 0; for (Eigen::Index i = 0; i < first_node_icu.get_num_time_points(); ++i) { - Eigen::VectorXd sum = Eigen::VectorXd::Zero(first_node_icu.get_num_elements()); + Eigen::VectorX sum = Eigen::VectorX::Zero(first_node_icu.get_num_elements()); for (auto& node : m_graph.nodes()) { auto& sim = node.property.get_simulation(); auto& icu_history = sim.get_parameters().template get>(); @@ -378,10 +380,10 @@ class FeedbackGraphSimulation Eigen::Index num_elements = first_node_icu.get_num_elements(); // For each region: vector of summed ICU values for all time points - std::unordered_map> regional_sums; + std::unordered_map>> regional_sums; for (auto const& [region_id, _] : regional_population) { regional_sums[region_id] = - std::vector(num_time_points, Eigen::VectorXd::Zero(num_elements)); + std::vector>(num_time_points, Eigen::VectorX::Zero(num_elements)); } // Sum up @@ -412,8 +414,8 @@ class FeedbackGraphSimulation */ void update_global_icu_occupancy() { - FP total_population = 0; - Eigen::VectorXd sum = Eigen::VectorXd::Zero(m_global_icu_occupancy.get_num_elements()); + FP total_population = 0; + Eigen::VectorX sum = Eigen::VectorX::Zero(m_global_icu_occupancy.get_num_elements()); for (auto& node : m_graph.nodes()) { auto& sim = node.property.get_simulation(); auto& icu_history = sim.get_parameters().template get>(); @@ -432,7 +434,7 @@ class FeedbackGraphSimulation std::unordered_map regional_population; for (auto& [region_id, regional_data] : m_regional_icu_occupancy) { - Eigen::VectorXd sum = Eigen::VectorXd::Zero(regional_data.get_num_elements()); + Eigen::VectorX sum = Eigen::VectorX::Zero(regional_data.get_num_elements()); for (auto& node : m_graph.nodes()) { if (mio::regions::get_state_id(node.id).get() == region_id) { auto& sim = node.property.get_simulation(); From c20d181bb644a9ddca61e7db72160c04a06fd2d4 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 29 Sep 2025 07:46:31 +0000 Subject: [PATCH 6/7] Resolve Issue 1382, 1383, 1384 --- .../mobility/metapopulation_mobility_instant.h | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 5be92d3f50..58fbb73a29 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -491,7 +491,7 @@ using get_infections_relative_expr_t = decltype(get_infections_relative( * @param t the current simulation time */ template ::value, void*> = nullptr> + std::enable_if_t::value, void*> = nullptr> FP get_infections_relative(const SimulationNode& /*node*/, FP /*t*/, const Eigen::Ref>& /*y*/) { @@ -499,10 +499,10 @@ FP get_infections_relative(const SimulationNode& /*node*/, FP /*t*/, return 0; } template ::value, void*> = nullptr> + std::enable_if_t::value, void*> = nullptr> FP get_infections_relative(const SimulationNode& node, FP t, const Eigen::Ref>& y) { - return get_infections_relative(node.get_simulation(), t, y); + return get_infections_relative(node.get_simulation(), t, y); } /** @@ -524,17 +524,17 @@ using get_mobility_factors_expr_t = decltype(get_mobility_factors( * @return a vector expression, same size as y, with the factor for each compartment. */ template ::value, void*> = nullptr> + std::enable_if_t::value, void*> = nullptr> auto get_mobility_factors(const SimulationNode& /*node*/, FP /*t*/, const Eigen::Ref>& y) { return Eigen::VectorX::Ones(y.rows()); } template ::value, void*> = nullptr> + std::enable_if_t::value, void*> = nullptr> auto get_mobility_factors(const SimulationNode& node, FP t, const Eigen::Ref>& y) { - return get_mobility_factors(node.get_simulation(), t, y); + return get_mobility_factors(node.get_simulation(), t, y); } /** @@ -555,16 +555,16 @@ using test_commuters_expr_t = decltype(test_commuters( * @param t the current simulation time. */ template ::value, void*> = nullptr> + std::enable_if_t::value, void*> = nullptr> void test_commuters(SimulationNode& /*node*/, Eigen::Ref> /*mobile_population*/, FP /*time*/) { } template ::value, void*> = nullptr> + std::enable_if_t::value, void*> = nullptr> void test_commuters(SimulationNode& node, Eigen::Ref> mobile_population, FP time) { - return test_commuters(node.get_simulation(), mobile_population, time); + return test_commuters(node.get_simulation(), mobile_population, time); } template From 11bb56c7d08d0d9b45b413afdb6d836457f55e0a Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:15:10 +0200 Subject: [PATCH 7/7] Implement suggestion to merge round and trunc namespace region --- cpp/memilio/ad/ad.h | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/cpp/memilio/ad/ad.h b/cpp/memilio/ad/ad.h index 92cb90e98c..8ea41afc6f 100644 --- a/cpp/memilio/ad/ad.h +++ b/cpp/memilio/ad/ad.h @@ -27,7 +27,7 @@ #include #include -// Extend automatic differentiation (AD) library to support std::round. +// Extend automatic differentiation (AD) library to support std::round and std::trunc. namespace ad { namespace internal @@ -58,14 +58,7 @@ static inline double round(const ad::internal::unary_intermediate static inline double trunc(const ad::internal::active_type& x)