diff --git a/cpp/models/ode_secir/parameters_io.cpp b/cpp/models/ode_secir/parameters_io.cpp index 7fe6fa6fa6..68923adb90 100644 --- a/cpp/models/ode_secir/parameters_io.cpp +++ b/cpp/models/ode_secir/parameters_io.cpp @@ -58,7 +58,7 @@ int get_region_id(int id) } IOResult read_confirmed_cases_data( - std::string const& path, std::vector& rki_data, std::vector const& vregion, Date date, + 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, @@ -74,12 +74,12 @@ IOResult read_confirmed_cases_data( }); if (max_date_entry == rki_data.end()) { log_error("RKI data file is empty."); - return failure(StatusCode::InvalidFileFormat, path + ", file is empty."); + return failure(StatusCode::InvalidFileFormat, "RKI file is empty."); } auto max_date = max_date_entry->date; if (max_date < date) { log_error("Specified date does not exist in RKI data"); - return failure(StatusCode::OutOfRange, path + ", 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); @@ -277,4 +277,4 @@ read_population_data(const std::string& path, const std::vector& vregion, b #endif // MEMILIO_HAS_JSONCPP -GCC_CLANG_DIAGNOSTIC(pop) \ No newline at end of file +GCC_CLANG_DIAGNOSTIC(pop) diff --git a/cpp/models/ode_secir/parameters_io.h b/cpp/models/ode_secir/parameters_io.h index 3998837548..7713232ce4 100644 --- a/cpp/models/ode_secir/parameters_io.h +++ b/cpp/models/ode_secir/parameters_io.h @@ -21,6 +21,7 @@ #define ODESECIR_PARAMETERS_IO_H #include "memilio/config.h" +#include #ifdef MEMILIO_HAS_JSONCPP @@ -37,28 +38,17 @@ namespace osecir namespace details { /** - * @brief interpolates age_ranges to param_ranges and saves ratios in interpolation - * @param age_ranges original age ranges of the data - * @param interpolation vector of ratios that are aplied to the data of age_ranges - * @param carry_over boolean vector which indicates whether there is an overflow from one age group to the next while interpolating data - */ -void interpolate_ages(const std::vector& age_ranges, std::vector>& interpolation, - std::vector& carry_over); - -/** - * @brief reads populations data from RKI - * @param[in] path Path to RKI file - * @param[in] rki_data Vector of ConfirmedCasesDataEntry%s - * @param[in] region vector of keys of the region of interest - * @param[in] year Specifies year at which the data is read - * @param[in] month Specifies month at which the data is read - * @param[in] day Specifies day at which the data is read - * @param[in, out] num_* output vector for number of people in the corresponding compartement - * @param[in] t_* vector average time it takes to get from one compartement to another for each age group - * @param[in] mu_* vector probabilities to get from one compartement to another for each age group + * @brief reads populations data from RKI. + * @param[in] rki_data Vector of ConfirmedCasesDataEntry%s. + * @param[in] vregion Vector of keys of the region of interest. + * @param[in] date Date at which the data is read. + * @param[in, out] vnum_* Output vector for number of people in the corresponding compartement. + * @param[in] vt_* vector Average time it takes to get from one compartement to another for each age group. + * @param[in] vmu_* vector Probabilities to get from one compartement to another for each age group. + * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. */ IOResult read_confirmed_cases_data( - std::string const& path, std::vector& rki_data, std::vector const& vregion, Date date, + 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, @@ -70,19 +60,17 @@ IOResult read_confirmed_cases_data( const std::vector& scaling_factor_inf); /** - * @brief sets populations data from json file with multiple age groups into a Model with one age group - * @tparam FP floating point data type, e.g., double - * @param[in, out] model vector of objects in which the data is set - * @param[in] path Path to confirmed cases file - * @param[in] region vector of keys of the region of interest - * @param[in] year Specifies year at which the data is read - * @param[in] month Specifies month at which the data is read - * @param[in] day Specifies day at which the data is read - * @param[in] scaling_factor_inf factors by which to scale the confirmed cases of rki data - */ + * @brief Sets populations data from already read case data with multiple age groups into a Model with one age group. + * @tparam FP Floating point data type, e.g., double. + * @param[in, out] model Vector of models in which the data is set. + * @param[in] case_data List of confirmed cases data entries. + * @param[in] region Vector of keys of the region of interest. + * @param[in] date Date at which the data is read. + * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. + */ template -IOResult set_confirmed_cases_data(std::vector>& model, const std::string& path, - std::vector const& region, Date date, +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 size_t num_age_groups = ConfirmedCasesDataEntry::age_group_names.size(); @@ -99,8 +87,6 @@ IOResult set_confirmed_cases_data(std::vector>& model, const std std::vector> mu_H_U{model.size()}; std::vector> mu_U_D{model.size()}; - BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path)); - for (size_t node = 0; node < model.size(); ++node) { for (size_t group = 0; group < num_age_groups; group++) { @@ -131,7 +117,7 @@ IOResult set_confirmed_cases_data(std::vector>& model, const std 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(path, case_data, region, date, num_Exposed, num_InfectedNoSymptoms, + 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)); @@ -162,27 +148,43 @@ IOResult set_confirmed_cases_data(std::vector>& model, const std } /** - * @brief reads number of ICU patients from DIVI register into Parameters - * @param path Path to DIVI file - * @param vregion Keys of the region of interest - * @param year Specifies year at which the data is read - * @param month Specifies month at which the data is read - * @param day Specifies day at which the data is read - * @param vnum_icu number of ICU patients + * @brief Sets the infected population for a given model based on confirmed cases data. Here, we + * read the case data from a file. + * @tparam FP Floating point data type, e.g., double. + * @param[in, out] Vector of models for which the confirmed cases data will be set. + * @param[in] path Path to the confirmed cases data file. + * @param[in] region Vector of keys of the region of interest. + * @param[in] date Date at which the data is read. + * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. + */ +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) +{ + 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)); + return success(); +} + +/** + * @brief Reads number of ICU patients from DIVI register into Parameters. + * @param[in] path Path to DIVI file. + * @param[in] vregion Keys of the region of interest. + * @param date Date for which we initialize. + * @param vnum_icu Number of ICU patients. */ IOResult read_divi_data(const std::string& path, const std::vector& vregion, Date date, std::vector& vnum_icu); /** - * @brief sets populations data from DIVI register into Model - * @tparam FP floating point data type, e.g., double - * @param model vector of objects in which the data is set - * @param path Path to DIVI file - * @param vregion vector of keys of the regions of interest - * @param year Specifies year at which the data is read - * @param month Specifies month at which the data is read - * @param day Specifies day at which the data is read - * @param 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, @@ -214,43 +216,54 @@ IOResult set_divi_data(std::vector>& model, const std::string& p } /** - * @brief Reads population data from census data - * @tparam FP floating point data type, e.g., double - * @param[in] path Path to RKI file - * @param[in] vregion Vector of keys of the regions of interest - * @param[in] accumulate_age_groups Specifies whether population data sould be accumulated to one age group + * @brief Reads population data from census data. + * @tparam FP floating point data type, e.g., double. + * @param[in] path Path to RKI file. + * @param[in] vregion Vector of keys of the regions of interest. + * @param[in] accumulate_age_groups Specifies whether population data should be accumulated to one age group. */ IOResult>> read_population_data(const std::string& path, const std::vector& vregion, bool accumulate_age_groups = false); /** - * @brief sets population data from census data - * @tparam FP floating point data type, e.g., double - * @param[in, out] model vector of objects in which the data is set - * @param[in] path Path to RKI file - * @param[in] vregion vector of keys of the regions of interest - * @param[in] accumulate_age_groups specifies whether population data sould be accumulated to one age group - */ +* @brief Sets population data from census data which has been read into num_population. +* @tparam FP floating point data type, e.g., double. +* @param[in, out] model Vector of models in which the data is set. There should be one model per region. +* @param[in] num_population Vector of population data. The size should be the same as vregion and model. +* @param[in] vregion Vector of keys of the regions of interest. +*/ template -IOResult set_population_data(std::vector>& model, const std::string& path, - const std::vector& vregion, bool accumulate_age_groups = false) +IOResult set_population_data(std::vector>& model, + const std::vector>& num_population, + const std::vector& vregion) { - BOOST_OUTCOME_TRY(auto&& num_population, read_population_data(path, vregion, accumulate_age_groups)); - + assert(num_population.size() == vregion.size()); + assert(model.size() == vregion.size()); for (size_t region = 0; region < vregion.size(); region++) { - if (std::accumulate(num_population[region].begin(), num_population[region].end(), 0.0) > 0) { - auto num_groups = model[region].parameters.get_num_groups(); - for (auto i = AgeGroup(0); i < num_groups; i++) { - model[region].populations.template set_difference_from_group_total( - {i, InfectionState::Susceptible}, num_population[region][size_t(i)]); - } - } - else { - log_warning("No population data available for region " + std::to_string(region) + - ". Population data has not been set."); + auto num_groups = model[region].parameters.get_num_groups(); + for (auto i = AgeGroup(0); i < num_groups; i++) { + model[region].populations.template set_difference_from_group_total( + {i, InfectionState::Susceptible}, num_population[region][size_t(i)]); } } + return success(); +} +/** +* @brief Sets population data from census data into a 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 RKI file containing population data. +* @param[in] vregion Vector of keys of the regions of interest. +*/ +template +IOResult set_population_data(std::vector>& model, const std::string& path, + const std::vector& vregion) +{ + // Specifies whether population data should be accumulated to one age group. + const bool is_single_age_group = static_cast(model[0].parameters.get_num_groups()) == 1; + BOOST_OUTCOME_TRY(const auto&& num_population, read_population_data(path, vregion, is_single_age_group)); + BOOST_OUTCOME_TRY(set_population_data(model, num_population, vregion)); return success(); } @@ -259,145 +272,70 @@ IOResult set_population_data(std::vector>& model, const std::str #ifdef MEMILIO_HAS_HDF5 /** -* @brief sets populations data from RKI into a Model with one age group -* @param model vector of objects in which the data is set -* @param data_dir Path to RKI files -* @param results_dir Path to result files -* @param region vector of keys of the region of interest -* @param year Specifies year at which the data is read -* @param month Specifies month at which the data is read -* @param day Specifies day at which the data is read -* @param scaling_factor_inf factors by which to scale the confirmed cases of rki data +* @brief Uses the initialisation method, which uses the reported data to set the initial conditions for the model for a given day. +* The initialisation is applied for a predefined number of days and finally saved in a timeseries for each region. In the end, +* we save the files "Results_rki.h5" and "Results_rki_sum.h5" in the results_dir. +* Results_rki.h5 contains a time series for each region and Results_rki_sum.h5 contains the sum of all regions. +* @param[in] models Vector of models in which the data is set. Copy is made to avoid changing the original model. +* @param[in] results_dir Path to result files. +* @param[in] region Vector of keys of the region of interest. +* @param[in] date Date for which the data should be read. +* @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. +* @param[in] scaling_factor_icu Factor by which to scale the icu cases of divi data. +* @param[in] num_days Number of days to be simulated/initialized. +* @param[in] divi_data_path Path to DIVI file. +* @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& model, const std::string& 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) +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) { - std::vector> t_InfectedNoSymptoms{model.size()}; - std::vector> t_Exposed{model.size()}; - std::vector> t_InfectedSymptoms{model.size()}; - 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 sum_mu_I_U(region.size(), 0); - std::vector> mu_I_U{model.size()}; - - const size_t num_age_groups = ConfirmedCasesDataEntry::age_group_names.size(); + const auto num_age_groups = (size_t)models[0].parameters.get_num_groups(); + assert(scaling_factor_inf.size() == num_age_groups); + assert(num_age_groups == ConfirmedCasesDataEntry::age_group_names.size()); + assert(models.size() == region.size()); + std::vector> extrapolated_data( + region.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_age_groups)); + BOOST_OUTCOME_TRY(auto&& num_population, details::read_population_data(population_data_path, region)); BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(confirmed_cases_path)); - for (size_t node = 0; node < model.size(); node++) { - for (size_t group = 0; group < ConfirmedCasesDataEntry::age_group_names.size(); group++) { - 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)]))); + for (int t = 0; t <= num_days; ++t) { + auto offset_day = offset_date_by_days(date, t); - mu_C_R[node].push_back( - model[node].parameters.template get>()[AgeGroup(group)]); - mu_I_H[node].push_back( - model[node].parameters.template get>()[AgeGroup(group)]); - mu_H_U[node].push_back(model[node].parameters.template get>()[AgeGroup(group)]); - mu_U_D[node].push_back(model[node].parameters.template get>()[AgeGroup(group)]); - - sum_mu_I_U[node] += - model[node].parameters.template get>()[AgeGroup(group)] * - model[node].parameters.template get>()[AgeGroup(group)]; - mu_I_U[node].push_back( - model[node].parameters.template get>()[AgeGroup(group)] * - model[node].parameters.template get>()[AgeGroup(group)]); + if (offset_day > Date(2020, 4, 23)) { + BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, region, offset_day, scaling_factor_icu)); + } + else { + log_warning("No DIVI data available for date: {}-{}-{}", offset_day.day, offset_day.month, offset_day.year); } - } - std::vector> rki_data( - region.size(), TimeSeries::zero(num_days, (size_t)InfectionState::Count * num_age_groups)); - - for (size_t t = 0; t < static_cast(num_days); ++t) { - std::vector> num_InfectedSymptoms( - model.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - std::vector> num_InfectedSymptomsConfirmed( - model.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - std::vector> num_death( - model.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - std::vector> num_rec( - model.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - std::vector> num_Exposed( - model.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - std::vector> num_InfectedNoSymptoms( - model.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - std::vector> num_InfectedNoSymptomsConfirmed( - model.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - std::vector> num_InfectedSevere( - model.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - std::vector> dummy_icu( - model.size(), std::vector(ConfirmedCasesDataEntry::age_group_names.size(), 0.0)); - std::vector num_icu(model.size(), 0.0); - - BOOST_OUTCOME_TRY(details::read_confirmed_cases_data( - confirmed_cases_path, case_data, region, date, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, - num_InfectedSevere, dummy_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(details::read_divi_data(divi_data_path, region, date, num_icu)); - bool interpolate_rki_age_groups = num_age_groups == 1 ? true : false; - BOOST_OUTCOME_TRY(auto&& num_population, - details::read_population_data(population_data_path, region, interpolate_rki_age_groups)); - - for (size_t i = 0; i < region.size(); i++) { - for (size_t age = 0; age < num_age_groups; age++) { - rki_data[i][t]((size_t)InfectionState::Exposed + (size_t)InfectionState::Count * age) = - num_Exposed[i][age]; - rki_data[i][t]((size_t)InfectionState::InfectedNoSymptoms + (size_t)InfectionState::Count * age) = - num_InfectedNoSymptoms[i][age]; - rki_data[i][t]((size_t)InfectionState::InfectedNoSymptomsConfirmed + - (size_t)InfectionState::Count * age) = num_InfectedNoSymptomsConfirmed[i][age]; - rki_data[i][t]((size_t)InfectionState::InfectedSymptoms + (size_t)InfectionState::Count * age) = - num_InfectedSymptoms[i][age]; - rki_data[i][t]((size_t)InfectionState::InfectedSymptomsConfirmed + - (size_t)InfectionState::Count * age) = num_InfectedSymptomsConfirmed[i][age]; - rki_data[i][t]((size_t)InfectionState::InfectedSevere + (size_t)InfectionState::Count * age) = - num_InfectedSevere[i][age]; - rki_data[i][t]((size_t)InfectionState::InfectedCritical + (size_t)InfectionState::Count * age) = - scaling_factor_icu * num_icu[i] * mu_I_U[i][age] / sum_mu_I_U[i]; - rki_data[i][t]((size_t)InfectionState::Recovered + (size_t)InfectionState::Count * age) = - num_rec[i][age]; - rki_data[i][t]((size_t)InfectionState::Dead + (size_t)InfectionState::Count * age) = num_death[i][age]; - rki_data[i][t]((size_t)InfectionState::Susceptible + (size_t)InfectionState::Count * age) = - num_population[i][age] - num_Exposed[i][age] - num_InfectedNoSymptoms[i][age] - - num_InfectedNoSymptomsConfirmed[i][age] - num_InfectedSymptoms[i][age] - - num_InfectedSymptomsConfirmed[i][age] - num_InfectedSevere[i][age] - num_rec[i][age] - - num_death[i][age] - - rki_data[i][t]((size_t)InfectionState::InfectedCritical + (size_t)InfectionState::Count * age); - } + 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(); } - date = offset_date_by_days(date, 1); } - auto num_groups = (int)(size_t)model[0].parameters.get_num_groups(); - BOOST_OUTCOME_TRY(save_result(rki_data, region, num_groups, path_join(dir, "Results_rki.h5"))); - auto rki_data_sum = sum_nodes(std::vector>>{rki_data}); - BOOST_OUTCOME_TRY(save_result({rki_data_sum[0][0]}, {0}, num_groups, path_join(dir, "Results_rki_sum.h5"))); + BOOST_OUTCOME_TRY( + 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}); + BOOST_OUTCOME_TRY( + 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&, const std::string&, std::vector const&, - Date, const std::vector&, double, int, const std::string&, - const std::string&, const std::string&) +IOResult +export_input_data_county_timeseries(std::vector models, const std::string& 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) { mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data."); return success(); @@ -405,12 +343,12 @@ IOResult export_input_data_county_timeseries(std::vector&, const st #endif // MEMILIO_HAS_HDF5 /** - * @brief reads population data from population files for the whole country - * @param model vector of model in which the data is set - * @param date Date for which the data should be read - * @param scaling_factor_inf factors by which to scale the confirmed cases of rki data - * @param scaling_factor_icu factor by which to scale the icu cases of divi data - * @param dir directory of files + * @brief Reads population data from population files for the whole country. + * @param[in, out] model Vector of model in which the data is set. + * @param[in] date Date for which the data should be read. + * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. + * @param[in] scaling_factor_icu Factor by which to scale the icu cases of divi data. + * @param[in] dir Directory of files. */ template IOResult read_input_data_germany(std::vector& model, Date date, @@ -426,19 +364,18 @@ IOResult read_input_data_germany(std::vector& model, Date date, } BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(dir, "cases_all_age_ma7.json"), {0}, date, scaling_factor_inf)); - BOOST_OUTCOME_TRY( - details::set_population_data(model, path_join(dir, "county_current_population.json"), {0}, false)); + BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(dir, "county_current_population.json"), {0})); return success(); } /** - * @brief reads population data from population files for the specefied state - * @param model vector of model in which the data is set - * @param date Date for which the data should be read - * @param state vector of region keys of states of interest - * @param scaling_factor_inf factors by which to scale the confirmed cases of rki data - * @param scaling_factor_icu factor by which to scale the icu cases of divi data - * @param dir directory of files + * @brief Reads population data from population files for the specefied state. + * @param[in, out] model Vector of model in which the data is set. + * @param[in] date Date for which the data should be read. + * @param[in] state Vector of region keys of states of interest. + * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. + * @param[in] scaling_factor_icu Factor by which to scale the icu cases of divi data. + * @param[in] dir Directory of files. */ template IOResult read_input_data_state(std::vector& model, Date date, std::vector& state, @@ -455,21 +392,20 @@ IOResult read_input_data_state(std::vector& model, Date date, std:: BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(dir, "cases_all_state_age_ma7.json"), state, date, scaling_factor_inf)); - BOOST_OUTCOME_TRY( - details::set_population_data(model, path_join(dir, "county_current_population.json"), state, false)); + BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(dir, "county_current_population.json"), state)); return success(); } /** - * @brief reads population data from population files for the specefied county - * @param model vector of model in which the data is set - * @param date Date for which the data should be read - * @param county vector of region keys of counties of interest - * @param scaling_factor_inf factors by which to scale the confirmed cases of rki data - * @param scaling_factor_icu factor by which to scale the icu cases of divi data - * @param dir directory of files - * @param num_days [Default: 0] Number of days to be simulated; required to extrapolate real data - * @param 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. + * @brief Reads population data from population files for the specefied county. + * @param[in, out] model Vector of model in which the data is set. + * @param[in] date Date for which the data should be read. + * @param[in] county Vector of region keys of counties of interest. + * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. + * @param[in] scaling_factor_icu Factor by which to scale the icu cases of divi data. + * @param[in] dir Directory of files. + * @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 IOResult read_input_data_county(std::vector& model, Date date, const std::vector& county, @@ -486,7 +422,7 @@ IOResult read_input_data_county(std::vector& model, Date date, cons BOOST_OUTCOME_TRY(details::set_confirmed_cases_data( model, path_join(dir, "pydata/Germany", "cases_all_county_age_ma7.json"), county, date, scaling_factor_inf)); BOOST_OUTCOME_TRY(details::set_population_data( - model, path_join(dir, "pydata/Germany", "county_current_population.json"), county, false)); + model, path_join(dir, "pydata/Germany", "county_current_population.json"), county)); if (export_time_series) { // Use only if extrapolated real data is needed for comparison. EXPENSIVE ! @@ -505,13 +441,13 @@ IOResult read_input_data_county(std::vector& model, Date date, cons /** * @brief reads population data from population files for the specified nodes - * @param model vector of model in which the data is set - * @param date Date for which the data should be read - * @param county vector of region keys of interest - * @param scaling_factor_inf factors by which to scale the confirmed cases of rki data - * @param scaling_factor_icu factor by which to scale the icu cases of divi data - * @param dir directory of files - * @param age_group_names strings specifying age group names + * @param[in, out] model vector of model in which the data is set + * @param[in] date Date for which the data should be read + * @param[in] county vector of region keys of interest + * @param[in] scaling_factor_inf factors by which to scale the confirmed cases of rki data + * @param[in] scaling_factor_icu factor by which to scale the icu cases of divi data + * @param[in] dir directory of files + * @param[in] age_group_names strings specifying age group names */ template IOResult read_input_data(std::vector& model, Date date, const std::vector& node_ids, @@ -527,9 +463,7 @@ IOResult read_input_data(std::vector& model, Date date, const std:: } BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(data_dir, "confirmed_cases.json"), node_ids, date, scaling_factor_inf)); - bool single_age_group = scaling_factor_inf.size() == 1 ? true : false; - BOOST_OUTCOME_TRY( - details::set_population_data(model, path_join(data_dir, "population_data.json"), node_ids, single_age_group)); + BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(data_dir, "population_data.json"), node_ids)); if (export_time_series) { // Use only if extrapolated real data is needed for comparison. EXPENSIVE ! diff --git a/cpp/models/ode_secirvvs/parameters_io.h b/cpp/models/ode_secirvvs/parameters_io.h index c6233ae368..89cfd4d3d7 100644 --- a/cpp/models/ode_secirvvs/parameters_io.h +++ b/cpp/models/ode_secirvvs/parameters_io.h @@ -40,13 +40,13 @@ namespace details { /** * @brief Reads subpopulations of infection states from transformed RKI cases file. - * @param path Path to transformed RKI cases file. - * @param vregion vector of keys of the region of interest - * @param date Date for which the arrays are initialized - * @param num_* output vector for number of people in the corresponding compartement - * @param t_* average time it takes to get from one compartement to another (vector with one element per age group) - * @param mu_* probabilities to get from one compartement to another (vector with one element per age group) - * @param scaling_factor_inf Factor for scaling the confirmed cases to account for estimated undetected cases. + * @param[in] path Path to transformed RKI cases file. + * @param[in] vregion Vector of keys of the region of interest. + * @param[in] date Date for which the arrays are initialized. + * @param[in, out] num_* Output vector for number of people in the corresponding compartement. + * @param[in] vt_* Average time it takes to get from one compartement to another (vector with one element per age group). + * @param[in] vmu_* Probabilities to get from one compartement to another (vector with one element per age group). + * @param[in] scaling_factor_inf Factor for scaling the confirmed cases to account for estimated undetected cases. * @see mio::read_confirmed_cases_data * @{ */ @@ -74,12 +74,12 @@ IOResult read_confirmed_cases_data( /**@}*/ /** - * @brief Reads confirmed cases data and translates data of day t0-delay to recovered compartment, - * @param path Path to RKI confirmed cases file. - * @param vregion vector of keys of the region of interest - * @param date Date for which the arrays are initialized - * @param num_rec output vector for number of people in the compartement recovered - * @param delay number of days in the past the are used to set recovered compartment. + * @brief Reads confirmed cases data and translates data of day t0-delay to recovered compartment. + * @param[in] path Path to RKI confirmed cases file. + * @param[in] vregion Vector of keys of the region of interest. + * @param[in] date Date for which the arrays are initialized. + * @param[in, out] num_rec Output vector for number of people in the compartement recovered. + * @param[in] delay Number of days in the past the are used to set recovered compartment. * @see mio::read_confirmed_cases_data * @{ */ @@ -92,25 +92,24 @@ IOResult read_confirmed_cases_data_fix_recovered(std::string const& path, /**@}*/ /** - * @brief sets populations data from a transformed RKI cases file into a Model. - * @param model vector of objects in which the data is set - * @param path Path to transformed RKI cases file - * @param region vector of keys of the region of interest - * @param date Date for which the arrays are initialized - * @param scaling_factor_inf factors by which to scale the confirmed cases of - * rki data + * @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::string& path, +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) + const std::vector& scaling_factor_inf, bool set_death = false) { 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); - BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path)); - std::vector> t_Exposed{model.size()}; std::vector> t_InfectedNoSymptoms{model.size()}; std::vector> t_InfectedSymptoms{model.size()}; @@ -160,7 +159,7 @@ IOResult set_confirmed_cases_data(std::vector& model, const std::st } } - BOOST_OUTCOME_TRY(read_confirmed_cases_data(rki_data, region, date, num_Exposed, num_InfectedNoSymptoms, + 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)); @@ -179,6 +178,16 @@ IOResult set_confirmed_cases_data(std::vector& model, const std::st model[county].populations[{AgeGroup(i), InfectionState::InfectedSevereNaive}] = num_InfectedSevere[county][i]; model[county].populations[{AgeGroup(i), InfectionState::SusceptibleImprovedImmunity}] = num_rec[county][i]; + if (set_death) { + // in set_confirmed_cases_data initilization, deaths are now set to 0. In order to visualize + // the extrapolated real number of deaths, they have to be set here. In the comparison of data + // it has to be paid attention to the fact, the the simulation starts with deaths=0 + // while this method starts with deaths=number of reported deaths so far... + // Additionally, we set the number of reported deaths to DeadNaive since no information on that is + // available here. + // Do only add deaths after substraction. + model[county].populations[{AgeGroup(i), InfectionState::DeadNaive}] = num_death[county][i]; + } } // } @@ -245,7 +254,7 @@ IOResult set_confirmed_cases_data(std::vector& model, const std::st } } - BOOST_OUTCOME_TRY(read_confirmed_cases_data(rki_data, region, date, num_Exposed, num_InfectedNoSymptoms, + 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)); @@ -328,7 +337,7 @@ IOResult set_confirmed_cases_data(std::vector& model, const std::st } } - BOOST_OUTCOME_TRY(read_confirmed_cases_data(rki_data, region, date, num_Exposed, num_InfectedNoSymptoms, + 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)); @@ -354,15 +363,36 @@ IOResult set_confirmed_cases_data(std::vector& model, const std::st std::to_string(region[county]) + ". Population data has not been set."); } } + + return success(); +} + +/** + * @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) +{ + 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)); return success(); } /** * @brief reads number of ICU patients from DIVI register into Parameters - * @param path Path to transformed DIVI file - * @param vregion Keys of the region of interest - * @param date Date for which the arrays are initialized - * @param vnum_icu number of ICU patients + * @param[in] path Path to transformed DIVI file + * @param[in] vregion Keys of the region of interest + * @param[in] date Date for which the arrays are initialized + * @param[in, out] vnum_icu number of ICU patients * @see mio::read_divi_data * @{ */ @@ -374,11 +404,11 @@ IOResult read_divi_data(const std::vector& divi_data, const std /** * @brief sets populations data from DIVI register into Model - * @param model vector of objects in which the data is set - * @param path Path to transformed DIVI file - * @param vregion vector of keys of the regions of interest - * @param date Date for which the arrays are initialized - * @param scaling_factor_icu factor by which to scale the icu cases of divi data + * @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, @@ -411,8 +441,8 @@ IOResult set_divi_data(std::vector& model, const std::string& path, /** * @brief reads population data from census data. - * @param path Path to population data file. - * @param vregion vector of keys of the regions of interest + * @param[in] path Path to population data file. + * @param[in] vregion vector of keys of the regions of interest * @see mio::read_population_data * @{ */ @@ -422,16 +452,23 @@ IOResult>> read_population_data(const std::vecto const std::vector& vregion); /**@}*/ +/** +* @brief sets population data from census data which has been read into num_population +* @param[in, out] model vector of objects in which the data is set +* @param[in] num_population vector of population data +* @param[in] case_data vector of case data. Each inner vector represents a different region +* @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::string& path, const std::string& path_rki, +IOResult set_population_data(std::vector& model, const std::vector>& num_population, + const std::vector& case_data, const std::vector& vregion, Date date) { - BOOST_OUTCOME_TRY(auto&& num_population, read_population_data(path, vregion)); - auto num_age_groups = ConfirmedCasesDataEntry::age_group_names.size(); - std::vector> num_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(path_rki, vregion, date, num_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(), 0.0) > 0) { @@ -440,7 +477,7 @@ IOResult set_population_data(std::vector& model, const std::string& double S_v = std::min( model[region].parameters.template get>()[{i, SimulationDay(0)}] + - num_rec[region][size_t(i)], + 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)}] - @@ -606,19 +643,48 @@ IOResult set_population_data(std::vector& model, const std::string& return success(); } -template -IOResult set_vaccination_data(std::vector>& model, const std::string& path, Date date, - const std::vector& vregion, int num_days) +/** +* @brief Sets population data from census data which has been read into num_population. +* @param[in, out] model Vector of objects in which the data is set. +* @param[in] path Path to population data file. +* @param[in] path_rki Path to RKI cases data file. +* @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::string& path, const std::string& path_rki, + const std::vector& vregion, Date date) { - BOOST_OUTCOME_TRY(auto&& vacc_data, read_vaccination_data(path)); + BOOST_OUTCOME_TRY(auto&& num_population, details::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)); + return success(); +} +/** + * @brief Sets vaccination data for models stored in a vector. + * + * @tparam FP Floating point type used in the Model objects. + * @param[in, out] model Vector of Model objects in which the vaccination data is set. + * @param[in] vacc_data Vector of VaccinationDataEntry objects containing the vaccination data. + * @param[in] date Start date for the simulation. + * @param[in] vregion Vector of region identifiers. + * @param[in] num_days Number of days for which the simulation is run. + */ +template +IOResult set_vaccination_data(std::vector>& model, const std::vector& vacc_data, + Date date, const std::vector& vregion, int num_days) +{ auto num_groups = model[0].parameters.get_num_groups(); - auto days_until_effective1 = - (int)(double)model[0].parameters.template get>()[AgeGroup(0)]; - auto days_until_effective2 = - (int)(double)model[0].parameters.template get>()[AgeGroup(0)]; - auto vaccination_distance = (int)(double)model[0].parameters.template get>()[AgeGroup(0)]; + // 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)])); + auto vaccination_distance = + static_cast(static_cast(model[0].parameters.template get>()[AgeGroup(0)])); // iterate over regions (e.g., counties) for (size_t i = 0; i < model.size(); ++i) { @@ -733,656 +799,116 @@ IOResult set_vaccination_data(std::vector>& model, const std::st return success(); } +/** + * @brief Reads vaccination data from a file and sets it for each model. + * + * @tparam FP Floating point type used in the Model objects. + * @param[in, out] model Vector of Model objects in which the vaccination data is set. + * @param[in] path Path to vaccination data file. + * @param[in] date Start date for the simulation. + * @param[in] vregion Vector of region identifiers. + * @param[in] num_days Number of days for which the simulation is run. + */ +template +IOResult set_vaccination_data(std::vector>& model, const std::string& path, Date date, + const std::vector& vregion, int num_days) +{ + BOOST_OUTCOME_TRY(auto&& vacc_data, read_vaccination_data(path)); + BOOST_OUTCOME_TRY(set_vaccination_data(model, vacc_data, date, vregion, num_days)); + return success(); +} + } // namespace details #ifdef MEMILIO_HAS_HDF5 /** - * @brief Exports the time series of extrapolated real data according to - * the extrapolation / approximation method used to initialize the model from - * real world data. - (This is the vector-valued functionality of set_confirmed_cases_data()) - * @param model vector of objects in which the data is set - * @param data_dir Path to transformed RKI cases files - * @param results_dir Path to result files - * @param start_date Start date of the time series to be exported. - * @param region vector of keys of the region of interest - * @param scaling_factor_inf Factor for scaling the confirmed cases to account for an estimated number of undetected cases. - * @param scaling_factor_icu Factor for scaling the reported ICU cases to account for possibly unreported ICU cases. - * @param num_days Number of days for which the time series is exported. - * @param divi_data_path path to divi data file - * @param confirmed_cases_path path to confirmed cases file - * @param population_data_path path to population data file - */ +* @brief Uses the initialisation method, which uses the reported data to set the initial conditions for the model for a given day. +* The initialisation is applied for a predefined number of days and finally saved in a timeseries for each region. In the end, +* we save the files "Results_rki.h5" and "Results_rki_sum.h5" in the results_dir. +* Results_rki.h5 contains a time series for each region and Results_rki_sum.h5 contains the sum of all regions. +* @param[in] models Vector of models in which the data is set. Copy is made to avoid changing the original model. +* @param[in] results_dir Path to result files. +* @param[in] counties Vector of keys of the counties of interest. +* @param[in] date Date for which the data should be read. +* @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. +* @param[in] scaling_factor_icu Factor by which to scale the icu cases of divi data. +* @param[in] num_days Number of days to be simulated/initialized. +* @param[in] divi_data_path Path to DIVI file. +* @param[in] confirmed_cases_path Path to confirmed cases file. +* @param[in] population_data_path Path to population data file. +* @param[in] vaccination_data_path Path to vaccination data file. +*/ template IOResult export_input_data_county_timeseries( - const std::vector& model, const std::string& dir, std::vector const& region, Date start_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) + 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::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path, + const std::string& vaccination_data_path = "") { - auto num_age_groups = (size_t)model[0].parameters.get_num_groups(); - assert(scaling_factor_inf.size() == num_age_groups); - assert(num_age_groups == ConfirmedCasesDataEntry::age_group_names.size()); - assert(model.size() == region.size()); - - BOOST_OUTCOME_TRY(auto&& rki_data, read_confirmed_cases_data(confirmed_cases_path)); - BOOST_OUTCOME_TRY(auto&& population_data, read_population_data(population_data_path)); - BOOST_OUTCOME_TRY(auto&& divi_data, read_divi_data(divi_data_path)); - - /* functionality copy from set_confirmed_cases_data() here splitted in params */ - /* which do not need to be reset for each day and compartments sizes that are */ - /* set later for each day */ - /*----------- UNVACCINATED -----------*/ - // data needs to be int, because access to data-specific confirmed cases - // is done with these parameters. TODO: Rounding instead - // of casting to int should be introduced in the future. - std::vector> t_InfectedNoSymptoms_uv{model.size()}; - std::vector> t_Exposed_uv{model.size()}; - std::vector> t_InfectedSymptoms_uv{model.size()}; - std::vector> t_InfectedSevere_uv{model.size()}; - std::vector> t_InfectedCritical_uv{model.size()}; - - std::vector> mu_C_R_uv{model.size()}; - std::vector> mu_I_H_uv{model.size()}; - std::vector> mu_H_U_uv{model.size()}; - // ICU data is not age-resolved. Use a partition of unity defined by - // the age-dependent probability I->H->U divided by the sum over all - // age groups of all of these probabilities. - std::vector sum_mu_I_U_uv(model.size(), 0); - std::vector> mu_I_U_uv{model.size()}; - - for (size_t county = 0; county < model.size(); county++) { - for (size_t group = 0; group < num_age_groups; group++) { - t_Exposed_uv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedNoSymptoms_uv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedSymptoms_uv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedSevere_uv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedCritical_uv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - - mu_C_R_uv[county].push_back( - model[county].parameters.template get>()[(AgeGroup)group]); - mu_I_H_uv[county].push_back( - model[county].parameters.template get>()[(AgeGroup)group]); - mu_H_U_uv[county].push_back( - model[county].parameters.template get>()[(AgeGroup)group]); - - /* begin: NOT in set_confirmed_cases_data() */ - sum_mu_I_U_uv[county] += - model[county].parameters.template get>()[AgeGroup(group)] * - model[county].parameters.template get>()[AgeGroup(group)]; - mu_I_U_uv[county].push_back( - model[county].parameters.template get>()[AgeGroup(group)] * - model[county].parameters.template get>()[AgeGroup(group)]); - /* end: NOT in set_confirmed_cases_data() */ - } + 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)); + + BOOST_OUTCOME_TRY(auto&& case_data, read_confirmed_cases_data(confirmed_cases_path)); + BOOST_OUTCOME_TRY(auto&& population_data, details::read_population_data(population_data_path, counties)); + + // empty vector if set_vaccination_data is not set + std::vector vacc_data; + if (!vaccination_data_path.empty()) { + BOOST_OUTCOME_TRY(vacc_data, read_vaccination_data(vaccination_data_path)); } - /*----------- PARTIALLY VACCINATED -----------*/ - // data needs to be int, because access to data-specific confirmed cases - // is done with these parameters. TODO: Rounding instead - // of casting to int should be introduced in the future. - std::vector> t_InfectedNoSymptoms_pv{model.size()}; - std::vector> t_Exposed_pv{model.size()}; - std::vector> t_InfectedSymptoms_pv{model.size()}; - std::vector> t_InfectedSevere_pv{model.size()}; - std::vector> t_InfectedCritical_pv{model.size()}; - - std::vector> mu_C_R_pv{model.size()}; - std::vector> mu_I_H_pv{model.size()}; - std::vector> mu_H_U_pv{model.size()}; - - // ICU data is not age-resolved. Use a partition of unity defined by - // the age-dependent probability I->H->U divided by the sum over all - // age groups of all of these probabilities. - std::vector sum_mu_I_U_pv(model.size(), 0); - std::vector> mu_I_U_pv{model.size()}; - for (size_t county = 0; county < model.size(); county++) { - for (size_t group = 0; group < num_age_groups; group++) { - double reduc_t = model[0].parameters.template get>()[(AgeGroup)group]; - t_Exposed_pv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedNoSymptoms_pv[county].push_back(static_cast(std::round( - model[county].parameters.template get>()[(AgeGroup)group] * reduc_t))); - t_InfectedSymptoms_pv[county].push_back(static_cast(std::round( - model[county].parameters.template get>()[(AgeGroup)group] * reduc_t))); - t_InfectedSevere_pv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedCritical_pv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - - double exp_fact_part_immune = - model[county].parameters.template get>()[(AgeGroup)group]; - double inf_fact_part_immune = - model[county].parameters.template get>()[(AgeGroup)group]; - double hosp_fact_part_immune = - model[county] - .parameters.template get>()[(AgeGroup)group]; - double icu_fact_part_immune = - model[county] - .parameters.template get>()[(AgeGroup)group]; - mu_C_R_pv[county].push_back(( - 1 - inf_fact_part_immune / exp_fact_part_immune * - (1 - model[county] - .parameters.template get>()[(AgeGroup)group]))); - mu_I_H_pv[county].push_back( - hosp_fact_part_immune / inf_fact_part_immune * - model[county].parameters.template get>()[(AgeGroup)group]); - // transfer from H to U, D unchanged. - mu_H_U_pv[county].push_back( - icu_fact_part_immune / hosp_fact_part_immune * - model[county].parameters.template get>()[(AgeGroup)group]); + for (int t = 0; t <= num_days; ++t) { + auto offset_day = offset_date_by_days(date, t); - sum_mu_I_U_pv[county] += - icu_fact_part_immune / hosp_fact_part_immune * - model[county].parameters.template get>()[AgeGroup(group)] * - hosp_fact_part_immune / inf_fact_part_immune * - model[county].parameters.template get>()[AgeGroup(group)]; - mu_I_U_pv[county].push_back( - icu_fact_part_immune / hosp_fact_part_immune * - model[county].parameters.template get>()[AgeGroup(group)] * - hosp_fact_part_immune / inf_fact_part_immune * - model[county].parameters.template get>()[AgeGroup(group)]); + if (!vaccination_data_path.empty()) { + BOOST_OUTCOME_TRY(details::set_vaccination_data(models, vacc_data, offset_day, counties, num_days)); } - } - /*----------- FULLY VACCINATED -----------*/ - // data needs to be int, because access to data-specific confirmed cases - // is done with these parameters. TODO: Rounding instead - // of casting to int should be introduced in the future. - std::vector> t_InfectedNoSymptoms_fv{model.size()}; - std::vector> t_Exposed_fv{model.size()}; - std::vector> t_InfectedSymptoms_fv{model.size()}; - std::vector> t_InfectedSevere_fv{model.size()}; - std::vector> t_InfectedCritical_fv{model.size()}; - - std::vector> mu_C_R_fv{model.size()}; - std::vector> mu_I_H_fv{model.size()}; - std::vector> mu_H_U_fv{model.size()}; - // ICU data is not age-resolved. Use a partition of unity defined by - // the age-dependent probability I->H->U divided by the sum over all - // age groups of all of these probabilities. - std::vector sum_mu_I_U_fv(model.size(), 0); - std::vector> mu_I_U_fv{model.size()}; - for (size_t county = 0; county < model.size(); county++) { - for (size_t group = 0; group < num_age_groups; group++) { - double reduc_t = model[0].parameters.template get>()[(AgeGroup)group]; - t_Exposed_fv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedNoSymptoms_fv[county].push_back(static_cast(std::round( - model[county].parameters.template get>()[(AgeGroup)group] * reduc_t))); - t_InfectedSymptoms_fv[county].push_back(static_cast(std::round( - model[county].parameters.template get>()[(AgeGroup)group] * reduc_t))); - t_InfectedSevere_fv[county].push_back(static_cast( - std::round(model[county].parameters.template get>()[(AgeGroup)group]))); - t_InfectedCritical_fv[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_fv[county].push_back(( - 1 - reduc_immune_inf / reduc_immune_exp * - (1 - model[county] - .parameters.template get>()[(AgeGroup)group]))); - mu_I_H_fv[county].push_back( - reduc_immune_hosp / reduc_immune_inf * - model[county].parameters.template get>()[(AgeGroup)group]); - // transfer from H to U, D unchanged. - mu_H_U_fv[county].push_back( - reduc_immune_icu / reduc_immune_hosp * - model[county].parameters.template get>()[(AgeGroup)group]); - - sum_mu_I_U_fv[county] += - reduc_immune_icu / reduc_immune_hosp * - model[county].parameters.template get>()[AgeGroup(group)] * - reduc_immune_hosp / reduc_immune_inf * - model[county].parameters.template get>()[AgeGroup(group)]; - mu_I_U_fv[county].push_back( - reduc_immune_icu / reduc_immune_hosp * - model[county].parameters.template get>()[AgeGroup(group)] * - reduc_immune_hosp / reduc_immune_inf * - model[county].parameters.template get>()[AgeGroup(group)]); - } - } - std::vector> extrapolated_rki( - model.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_age_groups)); - - for (size_t day = 0; day <= static_cast(num_days); day++) { - auto date = offset_date_by_days(start_date, int(day)); - - // unvaccinated - std::vector> num_Exposed_uv(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedNoSymptoms_uv(model.size(), - std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedSymptoms_uv(model.size(), - std::vector(num_age_groups, 0.0)); - // potential TODO: these confirmed are only confirmed by commuting, set to zero here. Adapt if generalized! - std::vector> num_InfectedNoSymptomsConfirmed_uv(model.size(), - std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedSymptomsConfirmed_uv(model.size(), - std::vector(num_age_groups, 0.0)); - // end TODO - std::vector> num_rec_uv(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedSevere_uv(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_death_uv(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> dummy_icu(model.size(), std::vector(num_age_groups, 0.0)); - BOOST_OUTCOME_TRY(details::read_confirmed_cases_data( - rki_data, region, date, num_Exposed_uv, num_InfectedNoSymptoms_uv, num_InfectedSymptoms_uv, - num_InfectedSevere_uv, dummy_icu, num_death_uv, num_rec_uv, t_Exposed_uv, t_InfectedNoSymptoms_uv, - t_InfectedSymptoms_uv, t_InfectedSevere_uv, t_InfectedCritical_uv, mu_C_R_uv, mu_I_H_uv, mu_H_U_uv, - scaling_factor_inf)); - - // partially vaccinated - std::vector> num_Exposed_pv(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedNoSymptoms_pv(model.size(), - std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedSymptoms_pv(model.size(), - std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedSevere_pv(model.size(), std::vector(num_age_groups, 0.0)); - // potential TODO: these confirmed are only confirmed by commuting, set to zero here. Adapt if generalized! - std::vector> num_InfectedNoSymptomsConfirmed_pv(model.size(), - std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedSymptomsConfirmed_pv(model.size(), - std::vector(num_age_groups, 0.0)); - // end TODO - std::vector> dummy_death(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> dummy_rec(model.size(), std::vector(num_age_groups, 0.0)); - for (size_t county = 0; county < model.size(); county++) { - dummy_death[county] = std::vector(num_age_groups, 0.0); - dummy_icu[county] = std::vector(num_age_groups, 0.0); + // 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. + if (offset_day > Date(2020, 4, 23)) { + BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu)); } - BOOST_OUTCOME_TRY(details::read_confirmed_cases_data( - rki_data, region, date, num_Exposed_pv, num_InfectedNoSymptoms_pv, num_InfectedSymptoms_pv, - num_InfectedSevere_pv, dummy_icu, dummy_death, dummy_rec, t_Exposed_pv, t_InfectedNoSymptoms_pv, - t_InfectedSymptoms_pv, t_InfectedSevere_pv, t_InfectedCritical_pv, mu_C_R_pv, mu_I_H_pv, mu_H_U_pv, - scaling_factor_inf)); - - // fully vaccinated - std::vector> num_Exposed_fv(model.size(), std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedNoSymptoms_fv(model.size(), - std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedSymptoms_fv(model.size(), - std::vector(num_age_groups, 0.0)); - // potential TODO: these confirmed are only confirmed by commuting, set to zero here. Adapt if generalized! - std::vector> num_InfectedNoSymptomsConfirmed_fv(model.size(), - std::vector(num_age_groups, 0.0)); - std::vector> num_InfectedSymptomsConfirmed_fv(model.size(), - std::vector(num_age_groups, 0.0)); - // end TODO - std::vector> num_InfectedSevere_fv(model.size(), std::vector(num_age_groups, 0.0)); - for (size_t county = 0; county < model.size(); county++) { - dummy_rec[county] = std::vector(num_age_groups, 0.0); - dummy_death[county] = std::vector(num_age_groups, 0.0); - dummy_icu[county] = std::vector(num_age_groups, 0.0); - } - BOOST_OUTCOME_TRY(details::read_confirmed_cases_data( - rki_data, region, date, num_Exposed_fv, num_InfectedNoSymptoms_fv, num_InfectedSymptoms_fv, - num_InfectedSevere_fv, dummy_icu, dummy_death, dummy_rec, t_Exposed_fv, t_InfectedNoSymptoms_fv, - t_InfectedSymptoms_fv, t_InfectedSevere_fv, t_InfectedCritical_fv, mu_C_R_fv, mu_I_H_fv, mu_H_U_fv, - scaling_factor_inf)); - - // ICU only read for compartment InfectionState::InfectedCritical and then distributed later - std::vector dummy_icu2(model.size(), 0.0); - BOOST_OUTCOME_TRY(details::read_divi_data(divi_data, region, date, dummy_icu2)); - - std::vector> num_icu(model.size(), std::vector(num_age_groups, 0.0)); - for (size_t county = 0; county < region.size(); county++) { - for (size_t age = 0; age < num_age_groups; age++) { - num_icu[county][age] = - scaling_factor_icu * dummy_icu2[county] * mu_I_U_uv[county][age] / sum_mu_I_U_uv[county]; - } + else { + log_warning("No DIVI data available for for date: {}-{}-{}", offset_day.day, offset_day.month, + offset_day.year); } - // read population basics - BOOST_OUTCOME_TRY(auto&& num_population, details::read_population_data(population_data, region)); - - std::vector> num_rec(model.size(), std::vector(num_age_groups, 0.0)); - BOOST_OUTCOME_TRY(details::read_confirmed_cases_data_fix_recovered(rki_data, region, date, num_rec, 14.)); - - for (size_t county = 0; county < region.size(); county++) { - if (std::accumulate(num_population[county].begin(), num_population[county].end(), 0.0) > 0) { - for (size_t age = 0; age < num_age_groups; age++) { - - auto age_group_offset = age * (size_t)InfectionState::Count; - double S_v = std::min(model[county].parameters.template get>()[{ - AgeGroup(age), SimulationDay(day)}] + - num_rec[county][age], - num_population[county][age]); - double S_pv = std::max(model[county].parameters.template get>()[{ - AgeGroup(age), SimulationDay(day)}] - - model[county].parameters.template get>()[{ - AgeGroup(age), SimulationDay(day)}], - 0.0); // use std::max with 0 - double S; - if (num_population[county][age] - S_pv - S_v < 0.0) { - log_warning("Number of vaccinated greater than population at county {}, age group {}: {} + " - "{} > {}.", - county, age, S_pv, S_v, num_population[county][age]); - S = 0.0; - S_v = num_population[county][age] - S_pv; - } - else { - S = num_population[county][age] - S_pv - S_v; - } - - double denom_E = - 1 / (S + - S_pv * model[county] - .parameters.template get>()[AgeGroup(age)] + - S_v * model[county] - .parameters.template get>()[AgeGroup(age)]); - double denom_C = - 1 / (S + - S_pv * model[county] - .parameters.template get>()[AgeGroup(age)] + - S_v * model[county] - .parameters.template get>()[AgeGroup(age)]); - double denom_I = - 1 / (S + - S_pv * model[county] - .parameters - .template get>()[AgeGroup(age)] + - S_v * model[county] - .parameters - .template get>()[AgeGroup(age)]); - double denom_HU = - 1 / (S + - S_pv * model[county] - .parameters.template get< - ReducInfectedSevereCriticalDeadPartialImmunity>()[AgeGroup(age)] + - S_v * model[county] - .parameters.template get< - ReducInfectedSevereCriticalDeadImprovedImmunity>()[AgeGroup(age)]); - - extrapolated_rki[county][day]((size_t)InfectionState::ExposedNaive + age_group_offset) = - S * denom_E * num_Exposed_uv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::ExposedPartialImmunity + age_group_offset) = - S_pv * - model[county].parameters.template get>()[AgeGroup(age)] * - denom_E * num_Exposed_pv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::ExposedImprovedImmunity + age_group_offset) = - S_v * - model[county].parameters.template get>()[AgeGroup(age)] * - denom_E * num_Exposed_fv[county][age]; - - extrapolated_rki[county][day]((size_t)InfectionState::InfectedNoSymptomsNaive + age_group_offset) = - S * denom_C * num_InfectedNoSymptoms_uv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedNoSymptomsPartialImmunity + - age_group_offset) = - S_pv * denom_C * num_InfectedNoSymptoms_pv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedNoSymptomsImprovedImmunity + - age_group_offset) = - S_v * denom_C * num_InfectedNoSymptoms_fv[county][age]; - - extrapolated_rki[county][day]((size_t)InfectionState::InfectedNoSymptomsNaiveConfirmed + - age_group_offset) = - S * denom_C * num_InfectedNoSymptomsConfirmed_uv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedNoSymptomsPartialImmunityConfirmed + - age_group_offset) = - S_pv * denom_C * num_InfectedNoSymptomsConfirmed_pv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed + - age_group_offset) = - S_v * denom_C * num_InfectedNoSymptomsConfirmed_fv[county][age]; - - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsNaive + age_group_offset) = - S * denom_I * num_InfectedSymptoms_uv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsPartialImmunity + - age_group_offset) = - S_pv * - model[county] - .parameters.template get>()[AgeGroup(age)] * - denom_I * num_InfectedSymptoms_pv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsImprovedImmunity + - age_group_offset) = - S_v * - model[county] - .parameters.template get>()[AgeGroup(age)] * - denom_I * num_InfectedSymptoms_fv[county][age]; - - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsNaiveConfirmed + - age_group_offset) = - S * denom_I * num_InfectedSymptomsConfirmed_uv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsPartialImmunityConfirmed + - age_group_offset) = - S_pv * - model[county] - .parameters.template get>()[AgeGroup(age)] * - denom_I * num_InfectedSymptomsConfirmed_pv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsImprovedImmunityConfirmed + - age_group_offset) = - S_v * - model[county] - .parameters.template get>()[AgeGroup(age)] * - denom_I * num_InfectedSymptomsConfirmed_fv[county][age]; - - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSevereNaive + age_group_offset) = - S * denom_HU * num_InfectedSevere_uv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSeverePartialImmunity + - age_group_offset) = - S_pv * - model[county] - .parameters - .template get>()[AgeGroup(age)] * - denom_HU * num_InfectedSevere_pv[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSevereImprovedImmunity + - age_group_offset) = - S_v * - model[county] - .parameters - .template get>()[AgeGroup(age)] * - denom_HU * num_InfectedSevere_fv[county][age]; - - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalNaive + age_group_offset) = - S * denom_HU * num_icu[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalPartialImmunity + - age_group_offset) = - S_pv * - model[county] - .parameters - .template get>()[AgeGroup(age)] * - denom_HU * num_icu[county][age]; - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalImprovedImmunity + - age_group_offset) = - S_v * - model[county] - .parameters - .template get>()[AgeGroup(age)] * - denom_HU * num_icu[county][age]; - - extrapolated_rki[county][day]((size_t)InfectionState::SusceptibleImprovedImmunity + - age_group_offset) = - model[county].parameters.template get>()[{AgeGroup(age), - SimulationDay(day)}] + - num_rec_uv[county][age] - - (extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsNaive + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsPartialImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsImprovedImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsNaiveConfirmed + - age_group_offset) + - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedSymptomsPartialImmunityConfirmed + age_group_offset) + - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedSymptomsImprovedImmunityConfirmed + age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSevereNaive + age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSeverePartialImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSevereImprovedImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalNaive + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalPartialImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalImprovedImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::DeadNaive + age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::DeadPartialImmunity + age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::DeadImprovedImmunity + - age_group_offset)); - - extrapolated_rki[county][day]((size_t)InfectionState::SusceptibleImprovedImmunity + - age_group_offset) = - std::min( - S_v - - extrapolated_rki[county][day]((size_t)InfectionState::ExposedImprovedImmunity + - age_group_offset) - - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedNoSymptomsImprovedImmunity + age_group_offset) - - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed + - age_group_offset) - - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsImprovedImmunity + - age_group_offset) - - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedSymptomsImprovedImmunityConfirmed + - age_group_offset) - - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSevereImprovedImmunity + - age_group_offset) - - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalImprovedImmunity + - age_group_offset) - - extrapolated_rki[county][day]((size_t)InfectionState::DeadImprovedImmunity + - age_group_offset), - std::max(0.0, - double(extrapolated_rki[county][day]( - (size_t)InfectionState::SusceptibleImprovedImmunity + age_group_offset)))); - - extrapolated_rki[county][day]((size_t)InfectionState::SusceptiblePartialImmunity + - age_group_offset) = - std::max(0.0, - S_pv - - extrapolated_rki[county][day]((size_t)InfectionState::ExposedPartialImmunity + - age_group_offset) - - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedNoSymptomsPartialImmunity + age_group_offset) - - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedNoSymptomsPartialImmunityConfirmed + - age_group_offset) - - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedSymptomsPartialImmunity + age_group_offset) - - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedSymptomsPartialImmunityConfirmed + - age_group_offset) - - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedSeverePartialImmunity + age_group_offset) - - extrapolated_rki[county][day]( - (size_t)InfectionState::InfectedCriticalPartialImmunity + age_group_offset) - - extrapolated_rki[county][day]((size_t)InfectionState::DeadPartialImmunity + - age_group_offset)); - - extrapolated_rki[county][day]((size_t)InfectionState::SusceptibleNaive + age_group_offset) = - num_population[county][age] - - (extrapolated_rki[county][day]((size_t)InfectionState::SusceptiblePartialImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::SusceptibleImprovedImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::ExposedNaive + age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::ExposedPartialImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::ExposedImprovedImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedNoSymptomsNaive + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedNoSymptomsPartialImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedNoSymptomsImprovedImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsNaive + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsPartialImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSymptomsImprovedImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSevereNaive + age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSeverePartialImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedSevereImprovedImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalNaive + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalPartialImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::InfectedCriticalImprovedImmunity + - age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::DeadNaive + age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::DeadPartialImmunity + age_group_offset) + - extrapolated_rki[county][day]((size_t)InfectionState::DeadImprovedImmunity + - age_group_offset)); - - // in set_confirmed_cases_data initilization, deaths are now set to 0. In order to visualize - // the extrapolated real number of deaths, they have to be set here. In the comparison of data - // it has to be paid attention to the fact, the the simulation starts with deaths=0 - // while this method starts with deaths=number of reported deaths so far... - // Additionally, we set the number of reported deaths to DeadNaive since no information on that is - // available here. - // Do only add deaths after substraction. - extrapolated_rki[county][day]((size_t)InfectionState::DeadNaive + age_group_offset) = - num_death_uv[county][age]; - } - } - else { - log_warning("No population data available for region " + std::to_string(county) + - ". Population data has not been set."); + 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)); + + for (size_t r = 0; r < counties.size(); r++) { + extrapolated_data[r][t] = models[r].get_initial_values(); + // in set_population_data the number of death individuals is subtracted from the SusceptibleImprovedImmunity compartment. + // Since we should be independent whether we consider them or not, we add them back here before we save the data. + for (size_t age = 0; age < num_groups; age++) { + extrapolated_data[r][t][(size_t)InfectionState::SusceptibleImprovedImmunity + + age * (size_t)InfectionState::Count] += + extrapolated_data[r][t][(size_t)InfectionState::DeadNaive + age * (size_t)InfectionState::Count]; } } - log_info("extrapolated real data for date: {}-{}-{}", date.day, date.month, date.year); - } - /* end: similar functionality in set_confirmed_cases_data(), here only for vector of TimeSeries */ - auto num_groups = (int)(size_t)model[0].parameters.get_num_groups(); - BOOST_OUTCOME_TRY(save_result(extrapolated_rki, region, num_groups, path_join(dir, "Results_rki.h5"))); - - auto extrapolated_rki_data_sum = sum_nodes(std::vector>>{extrapolated_rki}); - BOOST_OUTCOME_TRY( - save_result({extrapolated_rki_data_sum[0][0]}, {0}, num_groups, path_join(dir, "Results_rki_sum.h5"))); - - return success(); -} - -template -IOResult -export_input_data_county_timeseries(std::vector&& model, const std::string& 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, - bool set_vaccination_data, const std::string& vaccination_data_path) -{ - if (set_vaccination_data) { - BOOST_OUTCOME_TRY(details::set_vaccination_data(model, vaccination_data_path, date, region, num_days)); } + BOOST_OUTCOME_TRY(save_result(extrapolated_data, counties, static_cast(num_groups), + path_join(results_dir, "Results_rki.h5"))); - BOOST_OUTCOME_TRY(export_input_data_county_timeseries(model, dir, region, date, scaling_factor_inf, - scaling_factor_icu, num_days, divi_data_path, - confirmed_cases_path, population_data_path)); + 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 -IOResult export_input_data_county_timeseries(std::vector&, const std::string&, std::vector const&, - Date, const std::vector&, double, int, const std::string&, - const std::string&, const std::string&) -{ - mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data."); - return success(); -} - -template -IOResult export_input_data_county_timeseries(std::vector&&, const std::string&, std::vector const&, - Date, const std::vector&, double, int, const std::string&, - const std::string&, const std::string&, bool, const std::string&) +IOResult export_input_data_county_timeseries(std::vector, const std::string&, const std::vector&, + Date, const std::vector&, const double, const int, + const std::string&, const std::string&, const std::string&, + const std::string&) { mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data."); return success(); @@ -1395,14 +921,14 @@ IOResult export_input_data_county_timeseries(std::vector&&, const s * 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 model Vector of SECIRVVS models, one per county. - * @param date Date for which the data should be read. - * @param county Ids of the counties. - * @param scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county. - * @param scaling_factor_icu Factor of ICU cases to account for underreporting. - * @param dir Directory that contains the data files. - * @param num_days Number of days to be simulated; required to load data for vaccinations during the simulation. - * @param export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files. + * @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] 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, @@ -1438,7 +964,8 @@ IOResult read_input_data_county(std::vector& model, Date date, cons export_input_data_county_timeseries(model, dir, county, date, scaling_factor_inf, scaling_factor_icu, num_days, path_join(dir, "pydata/Germany", "county_divi_ma7.json"), path_join(dir, "pydata/Germany", "cases_all_county_age_ma7.json"), - path_join(dir, "pydata/Germany", "county_current_population.json"))); + path_join(dir, "pydata/Germany", "county_current_population.json"), + path_join(dir, "pydata/Germany", "all_county_ageinf_vacc_ma7.json"))); } return success(); @@ -1449,14 +976,14 @@ IOResult read_input_data_county(std::vector& model, Date date, cons * 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 model Vector of SECIRVVS models, one per county. - * @param date Date for which the data should be read. - * @param county Ids of the counties. - * @param scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county. - * @param scaling_factor_icu Factor of ICU cases to account for underreporting. - * @param dir Directory that contains the data files. - * @param num_days Number of days to be simulated; required to load data for vaccinations during the simulation. - * @param export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files. + * @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] 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, @@ -1490,8 +1017,8 @@ IOResult read_input_data(std::vector& model, Date date, const std:: "For simulation runs over the same time period, deactivate it."); BOOST_OUTCOME_TRY(export_input_data_county_timeseries( model, data_dir, node_ids, date, scaling_factor_inf, scaling_factor_icu, num_days, - path_join(data_dir, "critical_cases.json"), path_join(data_dir, "confirmed_cases.json"), - path_join(data_dir, "population_data.json"))); + path_join(data_dir, "divi_data.json"), path_join(data_dir, "confirmed_cases.json"), + path_join(data_dir, "population_data.json"), path_join(data_dir, "vaccination_data.json"))); } return success(); diff --git a/cpp/tests/data/export_time_series_initialization_osecir.h5 b/cpp/tests/data/export_time_series_initialization_osecir.h5 new file mode 100644 index 0000000000..34fe8bb25a Binary files /dev/null and b/cpp/tests/data/export_time_series_initialization_osecir.h5 differ diff --git a/cpp/tests/test_odesecir.cpp b/cpp/tests/test_odesecir.cpp index 90c8334a45..cb133d06bb 100644 --- a/cpp/tests/test_odesecir.cpp +++ b/cpp/tests/test_odesecir.cpp @@ -20,6 +20,7 @@ #include "distributions_helpers.h" #include "load_test_data.h" #include "matchers.h" +#include "temp_file_register.h" #include "ode_secir/model.h" #include "ode_secir/parameter_space.h" @@ -730,7 +731,7 @@ TEST(TestOdeSecir, testModelConstraints) } } -TEST(Secir, testAndTraceCapacity) +TEST(TestOdeSecir, testAndTraceCapacity) { double cont_freq = 10; @@ -780,7 +781,7 @@ TEST(Secir, testAndTraceCapacity) dydt_default[(size_t)mio::osecir::InfectionState::Exposed]); } -TEST(Secir, getInfectionsRelative) +TEST(TestOdeSecir, getInfectionsRelative) { size_t num_groups = 3; mio::osecir::Model model((int)num_groups); @@ -799,7 +800,7 @@ TEST(Secir, getInfectionsRelative) (100. + 50. + 25.) / (10'000 + 20'000 + 40'000)); } -TEST(Secir, get_reproduction_number) +TEST(TestOdeSecir, get_reproduction_number) { const size_t num_groups = 3; mio::osecir::Model model((int)num_groups); @@ -1000,7 +1001,7 @@ TEST(Secir, get_mobility_factors) } } -TEST(Secir, test_commuters) +TEST(TestOdeSecir, test_commuters) { auto model = mio::osecir::Model(2); auto mobility_factor = 0.1; @@ -1034,7 +1035,7 @@ TEST(Secir, test_commuters) 1e-5); } -TEST(Secir, check_constraints_parameters) +TEST(TestOdeSecir, check_constraints_parameters) { auto model = mio::osecir::Model(1); EXPECT_EQ(model.parameters.check_constraints(), 0); @@ -1108,7 +1109,7 @@ TEST(Secir, check_constraints_parameters) mio::set_log_level(mio::LogLevel::warn); } -TEST(Secir, apply_constraints_parameters) +TEST(TestOdeSecir, apply_constraints_parameters) { auto model = mio::osecir::Model(1); auto indx_agegroup = mio::AgeGroup(0); @@ -1188,7 +1189,7 @@ TEST(Secir, apply_constraints_parameters) #if defined(MEMILIO_HAS_JSONCPP) -TEST(Secir, read_population_data_one_age_group) +TEST(TestOdeSecir, read_population_data_one_age_group) { std::string path = mio::path_join(TEST_DATA_DIR, "county_current_population.json"); const std::vector region{1001}; @@ -1202,5 +1203,208 @@ TEST(Secir, read_population_data_one_age_group) EXPECT_EQ(result_multiple_age_groups[0].size(), 6); EXPECT_EQ(result_multiple_age_groups[0][0], 3433.0); } +#if defined(MEMILIO_HAS_HDF5) +class ModelTestOdeSecir : public testing::Test +{ + +public: + ModelTestOdeSecir() + : model(6) + , num_age_groups(6) + { + } + mio::osecir::Model model; + size_t num_age_groups; + +protected: + void SetUp() override + { + model.parameters.set(60); + model.parameters.set>(0.2); + + for (auto i = mio::AgeGroup(0); i < (mio::AgeGroup)num_age_groups; i++) { + model.parameters.get>()[i] = 3.2; + model.parameters.get>()[i] = 2.0; + model.parameters.get>()[i] = 5.8; + model.parameters.get>()[i] = 9.5; + model.parameters.get>()[i] = 7.1; + + model.parameters.get>()[i] = 0.05; + model.parameters.get>()[i] = 0.7; + model.parameters.get>()[i] = 0.09; + model.parameters.get>()[i] = 0.25; + model.parameters.get>()[i] = 0.45; + model.parameters.get>()[i] = 0.2; + model.parameters.get>()[i] = 0.25; + model.parameters.get>()[i] = 0.3; + } + model.apply_constraints(); + } +}; + +// Initialize model with data frrom external reported data. Check if the model is initialized correctly. +TEST_F(ModelTestOdeSecir, read_input_data) +{ + auto model1 = std::vector>{model}; + 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_DATA_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_DATA_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, "pydata/District"), 10); + + EXPECT_THAT(read_result1, IsSuccess()); + EXPECT_THAT(read_result2, IsSuccess()); + EXPECT_THAT(read_result_district, IsSuccess()); + + // values were generated by the tested function; can only check stability of the implementation, not correctness + auto expected_values = + (Eigen::ArrayXd(num_age_groups * Eigen::Index(mio::osecir::InfectionState::Count)) << 10280, 2.82575, 1.25589, + 0, 5.85714, 0, 1.42857, 1.33333, 28.2857, 0, 19091.3, 6.59341, 4.23862, 0, 11, 0, 2.54286, 1.33333, 88, 0, + 73821.7, 41.9152, 21.6641, 0, 54.5714, 0, 22, 1.33333, 523.857, 0, 82533.6, 39.5604, 21.0361, 0, 51.4286, 0, + 15.1714, 1.33333, 468.571, 0, 43733.1, 8.32025, 5.18053, 0, 11.4286, 0, 4.22857, 1.33333, 158.571, 8, 15642.7, + 10.5181, 3.2967, 0, 3.28571, 0, 0.657143, 1.33333, 49.8571, 7.57143) + .finished(); + + EXPECT_THAT(print_wrap(model1[0].populations.array().cast()), + MatrixNear(print_wrap(expected_values), 1e-5, 1e-5)); + EXPECT_THAT(print_wrap(model2[0].populations.array().cast()), + MatrixNear(print_wrap(expected_values), 1e-5, 1e-5)); + EXPECT_THAT(print_wrap(model3[0].populations.array().cast()), + MatrixNear(print_wrap(expected_values), 1e-5, 1e-5)); +} + +// Test if the export of the time series is called not failing and results are the same as the initialization. +TEST_F(ModelTestOdeSecir, export_time_series_init) +{ + TempFileRegister temp_file_register; + auto tmp_results_dir = temp_file_register.get_unique_path(); + EXPECT_THAT(mio::create_directory(tmp_results_dir), IsSuccess()); + + // Test exporting time series + std::vector> models{model}; + 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(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")), + IsSuccess()); + + 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(); + + EXPECT_THAT(print_wrap(data_extrapolated.value()[0].get_groups().matrix()), + MatrixNear(print_wrap(expected_results[0].get_groups().matrix()), 1e-5, 1e-5)); +} + +// Test the output of the function for a day way in the past. The model should be initialized with the population data since no Case data is available there. +TEST_F(ModelTestOdeSecir, export_time_series_init_old_date) +{ + mio::set_log_level(mio::LogLevel::off); + TempFileRegister temp_file_register; + auto tmp_results_dir = temp_file_register.get_unique_path(); + EXPECT_THAT(mio::create_directory(tmp_results_dir), IsSuccess()); + + // Test exporting time series + std::vector> models{model}; + 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(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")), + IsSuccess()); + + 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); + + // 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{1002}; + auto population_data = mio::osecir::details::read_population_data(path, region, false).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++) { + EXPECT_EQ(results_extrapolated(i * Eigen::Index(mio::osecir::InfectionState::Count)), population_data[0][i]); + } + // sum of all compartments should be equal to the population + EXPECT_EQ(results_extrapolated.sum(), std::accumulate(population_data[0].begin(), population_data[0].end(), 0.0)); + mio::set_log_level(mio::LogLevel::warn); +} + +// // Model initialization should return same start values as export time series on that day +TEST_F(ModelTestOdeSecir, model_initialization) +{ + // Vector assignment necessary as read_input_data_county changes model + auto model_vector = std::vector>{model}; + + 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, + TEST_DATA_DIR, 2, false), + IsSuccess()); + + // Values from data/export_time_series_init_osecir.h5, for reading in comparison + // operator for return of mio::read_result and model population needed. + auto expected_values = + (Eigen::ArrayXd(num_age_groups * Eigen::Index(mio::osecir::InfectionState::Count)) << 10280, 2.82575, 1.25589, + 0, 5.85714, 0, 1.42857, 1.33333, 28.2857, 0, 19091.3, 6.59341, 4.23862, 0, 11, 0, 2.54286, 1.33333, 88, 0, + 73821.7, 41.9152, 21.6641, 0, 54.5714, 0, 22, 1.33333, 523.857, 0, 82533.6, 39.5604, 21.0361, 0, 51.4286, 0, + 15.1714, 1.33333, 468.571, 0, 43733.1, 8.32025, 5.18053, 0, 11.4286, 0, 4.22857, 1.33333, 158.571, 8, 15642.7, + 10.5181, 3.2967, 0, 3.28571, 0, 0.657143, 1.33333, 49.8571, 7.57143) + .finished(); + + EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), + MatrixNear(print_wrap(expected_values), 1e-5, 1e-5)); +} + +// Calling the model initialization with a date way in the past should only initialize the model with the population data. +TEST_F(ModelTestOdeSecir, model_initialization_old_date) +{ + mio::set_log_level(mio::LogLevel::off); + // Vector assignment necessary as read_input_data_county changes model + auto model_vector = std::vector>{model}; + + 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, + TEST_DATA_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{1002}; + auto population_data = mio::osecir::details::read_population_data(path, region, false).value(); + + // So, the expected values are the population data in the susceptible compartments and zeros in the other compartments. + auto expected_values = + (Eigen::ArrayXd(num_age_groups * Eigen::Index(mio::osecir::InfectionState::Count)) << population_data[0][0], 0, + 0, 0, 0, 0, 0, 0, 0, 0, population_data[0][1], 0, 0, 0, 0, 0, 0, 0, 0, 0, population_data[0][2], 0, 0, 0, 0, 0, + 0, 0, 0, 0, population_data[0][3], 0, 0, 0, 0, 0, 0, 0, 0, 0, population_data[0][4], 0, 0, 0, 0, 0, 0, 0, 0, 0, + population_data[0][5], 0, 0, 0, 0, 0, 0, 0, 0, 0) + .finished(); + + auto results_extrapolated = model_vector[0].populations.array().cast(); + + for (size_t i = 0; i < num_age_groups; i++) { + EXPECT_EQ(results_extrapolated(i * Eigen::Index(mio::osecir::InfectionState::Count)), population_data[0][i]); + } + // sum of all compartments should be equal to the population + EXPECT_EQ(results_extrapolated.sum(), std::accumulate(population_data[0].begin(), population_data[0].end(), 0.0)); + mio::set_log_level(mio::LogLevel::warn); +} + +#endif #endif diff --git a/cpp/tests/test_odesecirvvs.cpp b/cpp/tests/test_odesecirvvs.cpp index 65f7a10f90..b5148a4ac3 100644 --- a/cpp/tests/test_odesecirvvs.cpp +++ b/cpp/tests/test_odesecirvvs.cpp @@ -45,7 +45,7 @@ #include "ode_secirvvs/analyze_result.h" #include "gtest/gtest.h" -#include "gmock/gmock-matchers.h" +#include "gmock/gmock.h" #include #include #include @@ -743,7 +743,7 @@ TEST(TestOdeSECIRVVS, export_time_series_init) std::vector(size_t(num_age_groups), 1.0), 1.0, 2, 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"), true, + mio::path_join(TEST_DATA_DIR, "county_current_population.json"), mio::path_join(TEST_DATA_DIR, "vacc_county_ageinf_ma7.json")), IsSuccess()); @@ -759,6 +759,53 @@ TEST(TestOdeSECIRVVS, export_time_series_init) MatrixNear(print_wrap(expected_results[0].get_groups().matrix()), 1e-5, 1e-5)); } +TEST(TestOdeSECIRVVS, export_time_series_init_old_date) +{ + mio::set_log_level(mio::LogLevel::off); + TempFileRegister temp_file_register; + auto tmp_results_dir = temp_file_register.get_unique_path(); + ASSERT_THAT(mio::create_directory(tmp_results_dir), IsSuccess()); + + auto num_age_groups = 6; // Data to be read requires RKI confirmed cases data age groups + auto model = make_model(num_age_groups); + + // set vaccinations to zero + model.parameters.get>().array().setConstant(0); + model.parameters.get>().array().setConstant(0); + // set all compartments to zero + model.populations.array().setConstant(0.0); + + // Test exporting time series + 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"), + mio::path_join(TEST_DATA_DIR, "cases_all_county_age_ma7.json"), + mio::path_join(TEST_DATA_DIR, "county_current_population.json"), + 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")); + ASSERT_THAT(data_extrapolated, IsSuccess()); + auto results_extrapolated = data_extrapolated.value()[0].get_groups().get_value(0); + + // 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::osecirvvs::details::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++) { + EXPECT_NEAR(results_extrapolated(i * Eigen::Index(mio::osecirvvs::InfectionState::Count)), + population_data[0][i], 1e-10); + } + // sum of all compartments should be equal to the population + EXPECT_NEAR(results_extrapolated.sum(), std::accumulate(population_data[0].begin(), population_data[0].end(), 0.0), + 1e-5); + mio::set_log_level(mio::LogLevel::warn); +} + // Model initialization should return same start values as export time series on that day TEST(TestOdeSECIRVVS, model_initialization) { @@ -795,6 +842,148 @@ TEST(TestOdeSECIRVVS, model_initialization) MatrixNear(print_wrap(expected_values), 1e-5, 1e-5)); } +TEST(TestOdeSECIRVVS, model_initialization_old_date) +{ + mio::set_log_level(mio::LogLevel::off); + constexpr auto num_age_groups = 6; // Data to be read requires RKI confirmed cases data age groups + auto model = make_model(num_age_groups); + // set vaccinations to zero + model.parameters.get>().array().setConstant(0); + model.parameters.get>().array().setConstant(0); + // set all compartments to zero + model.populations.array().setConstant(0.0); + + 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_DATA_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::osecirvvs::details::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++) { + EXPECT_NEAR( + model_vector[0].populations.array().cast()(i * Eigen::Index(mio::osecirvvs::InfectionState::Count)), + population_data[0][i], 1e-5); + } + + // sum of all compartments should be equal to the population + EXPECT_NEAR(model_vector[0].populations.array().cast().sum(), + std::accumulate(population_data[0].begin(), population_data[0].end(), 0.0), 1e-5); + mio::set_log_level(mio::LogLevel::warn); +} + +TEST(TestOdeSECIRVVS, model_initialization_old_date_county) +{ + mio::set_log_level(mio::LogLevel::off); + constexpr auto num_age_groups = 6; // Data to be read requires RKI confirmed cases data age groups + auto model = make_model(num_age_groups); + // set vaccinations to zero + model.parameters.get>().array().setConstant(0); + model.parameters.get>().array().setConstant(0); + // set all compartments to zero + model.populations.array().setConstant(0.0); + + 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_DATA_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::osecirvvs::details::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++) { + EXPECT_NEAR( + model_vector[0].populations.array().cast()(i * Eigen::Index(mio::osecirvvs::InfectionState::Count)), + population_data[0][i], 1e-5); + } + + // sum of all compartments should be equal to the population + EXPECT_NEAR(model_vector[0].populations.array().cast().sum(), + std::accumulate(population_data[0].begin(), population_data[0].end(), 0.0), 1e-5); + mio::set_log_level(mio::LogLevel::warn); +} + +TEST(TestOdeSECIRVVS, set_population_data_overflow_vacc) +{ + auto num_age_groups = 6; // Data to be read requires RKI confirmed cases data age groups + auto model = make_model(num_age_groups); + // set all compartments to zero + model.populations.array().setConstant(0.0); + + model.parameters + .template get>()[{mio::AgeGroup(0), mio::SimulationDay(0)}] = + 1e7 + 1; + + model.parameters + .template get>()[{mio::AgeGroup(0), mio::SimulationDay(0)}] = 1e7; + + auto model_vector = std::vector>{model}; + + std::string path_pop_data = mio::path_join(TEST_DATA_DIR, "county_current_population.json"); + const std::vector region{0}; + auto population_data = mio::osecirvvs::details::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( + model_vector, path_pop_data, mio::path_join(TEST_DATA_DIR, "cases_all_county_age_ma7.json"), {0}, + {1000, 12, 01}), + IsSuccess()); + + EXPECT_NEAR( + double(model_vector[0].populations[{mio::AgeGroup(0), mio::osecirvvs::InfectionState::SusceptibleNaive}]), 0.0, + 1e-10); + EXPECT_NEAR(double(model_vector[0].populations[{mio::AgeGroup(0), + mio::osecirvvs::InfectionState::SusceptiblePartialImmunity}]), + 1, 1e-10); + EXPECT_NEAR(double(model_vector[0].populations[{mio::AgeGroup(0), + mio::osecirvvs::InfectionState::SusceptibleImprovedImmunity}]), + population_data[0][0] - 1, 1e-10); + + EXPECT_NEAR(model_vector[0].populations.get_group_total(mio::AgeGroup(0)), population_data[0][0], 1e-9); +} + +TEST(TestOdeSECIRVVS, set_population_data_no_data_avail) +{ + auto num_age_groups = 6; // Data to be read requires RKI confirmed cases data age groups + auto model = make_model(num_age_groups); + // set all compartments to zero + model.populations.array().setConstant(0.0); + + // if the number of vaccinated individuals is greater than the population, we must limit the number of vaccinated. + model.parameters + .template get>()[{mio::AgeGroup(0), mio::SimulationDay(0)}] = 0; + + model.parameters + .template get>()[{mio::AgeGroup(0), mio::SimulationDay(0)}] = 0; + + auto model_vector = std::vector>{model}; + + std::string path_pop_data = mio::path_join(TEST_DATA_DIR, "county_current_population.json"); + const std::vector region{0}; + auto population_data = mio::osecirvvs::details::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( + model_vector, path_pop_data, mio::path_join(TEST_DATA_DIR, "cases_all_county_age_ma7.json"), {200}, + {1000, 12, 01}), + IsSuccess()); + + EXPECT_NEAR(model.populations.get_total(), 0.0, 1e-10); +} + TEST(TestOdeSECIRVVS, run_simulation) { auto num_age_groups = 3;