diff --git a/cpp/memilio/io/epi_data.h b/cpp/memilio/io/epi_data.h index bb71bcc799..dd2fe5d8d6 100644 --- a/cpp/memilio/io/epi_data.h +++ b/cpp/memilio/io/epi_data.h @@ -230,6 +230,18 @@ class DiviEntry } }; +/** + * @brief Checks if DIVI data is available for a given date. + * @param date The date to check. + * @return True if date is within 2020-04-23 and 2024-07-21, false otherwise. + */ +inline bool is_divi_data_available(const Date& date) +{ + static const Date divi_data_start(2020, 4, 23); + static const Date divi_data_end(2024, 7, 21); + return date >= divi_data_start && date <= divi_data_end; +} + /** * Deserialize a list of DiviEntry from json. * @param jsvalue Json value that contains DIVI data. diff --git a/cpp/models/ode_secir/parameters_io.cpp b/cpp/models/ode_secir/parameters_io.cpp index 1e4c4b65ff..0cebd99a04 100644 --- a/cpp/models/ode_secir/parameters_io.cpp +++ b/cpp/models/ode_secir/parameters_io.cpp @@ -123,8 +123,6 @@ IOResult read_confirmed_cases_data( auto date_df = region_entry.date; auto age = size_t(region_entry.age_group); - bool read_icu = false; //params.populations.get({age, SecirCompartments::U}) == 0; - if (date_df == offset_date_by_days(date, 0)) { num_InfectedSymptoms[age] += scaling_factor_inf[age] * region_entry.num_confirmed; num_rec[age] += region_entry.num_confirmed; @@ -147,16 +145,12 @@ IOResult read_confirmed_cases_data( } if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age])) { num_InfectedSevere[age] -= mu_I_H[age] * scaling_factor_inf[age] * region_entry.num_confirmed; - if (read_icu) { - num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed; - } + num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed; } if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age] - t_InfectedCritical[age])) { num_death[age] += region_entry.num_deaths; - if (read_icu) { - num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed; - } + num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed; } } } diff --git a/cpp/models/ode_secir/parameters_io.h b/cpp/models/ode_secir/parameters_io.h index d11f608be5..5e6649866b 100644 --- a/cpp/models/ode_secir/parameters_io.h +++ b/cpp/models/ode_secir/parameters_io.h @@ -134,8 +134,12 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect num_InfectedSymptoms[node][i]; model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptomsConfirmed}] = 0; model[node].populations[{AgeGroup(i), InfectionState::InfectedSevere}] = num_InfectedSevere[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::Dead}] = num_death[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::Recovered}] = num_rec[node][i]; + // Only set the number of ICU patients here, if the date is not available in the data. + if (!is_divi_data_available(date)) { + model[node].populations[{AgeGroup(i), InfectionState::InfectedCritical}] = num_icu[node][i]; + } + model[node].populations[{AgeGroup(i), InfectionState::Dead}] = num_death[node][i]; + model[node].populations[{AgeGroup(i), InfectionState::Recovered}] = num_rec[node][i]; } } else { @@ -190,6 +194,12 @@ template IOResult set_divi_data(std::vector>& model, const std::string& path, const std::vector& vregion, Date date, double scaling_factor_icu) { + // DIVI dataset will no longer be updated from CW29 2024 on. + if (!is_divi_data_available(date)) { + log_warning("No DIVI data available for date: {}-{}-{}", date.year, date.month, date.day, + ". ICU compartment will be set based on Case data."); + return success(); + } std::vector sum_mu_I_U(vregion.size(), 0); std::vector> mu_I_U{model.size()}; for (size_t region = 0; region < vregion.size(); region++) { @@ -306,13 +316,7 @@ IOResult export_input_data_county_timeseries( for (int t = 0; t <= num_days; ++t) { auto offset_day = offset_date_by_days(date, t); - 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); - } - + BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, region, offset_day, scaling_factor_icu)); BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(models, case_data, region, offset_day, scaling_factor_inf)); BOOST_OUTCOME_TRY(details::set_population_data(models, num_population, region)); for (size_t r = 0; r < region.size(); r++) { @@ -355,13 +359,8 @@ IOResult read_input_data_germany(std::vector& model, Date date, const std::vector& scaling_factor_inf, double scaling_factor_icu, const std::string& dir) { - if (date > Date(2020, 4, 23)) { - BOOST_OUTCOME_TRY( - details::set_divi_data(model, path_join(dir, "germany_divi.json"), {0}, date, scaling_factor_icu)); - } - else { - log_warning("No DIVI data available for this date"); - } + BOOST_OUTCOME_TRY( + details::set_divi_data(model, path_join(dir, "germany_divi.json"), {0}, date, scaling_factor_icu)); 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})); @@ -382,14 +381,9 @@ IOResult read_input_data_state(std::vector& model, Date date, std:: const std::vector& scaling_factor_inf, double scaling_factor_icu, const std::string& dir) { - if (date > Date(2020, 4, 23)) { - BOOST_OUTCOME_TRY( - details::set_divi_data(model, path_join(dir, "state_divi.json"), state, date, scaling_factor_icu)); - } - else { - log_warning("No DIVI data available for this date"); - } + BOOST_OUTCOME_TRY( + details::set_divi_data(model, path_join(dir, "state_divi.json"), state, date, scaling_factor_icu)); 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)); @@ -412,13 +406,8 @@ IOResult read_input_data_county(std::vector& model, Date date, cons const std::vector& scaling_factor_inf, double scaling_factor_icu, const std::string& dir, int num_days = 0, bool export_time_series = false) { - if (date > Date(2020, 4, 23)) { - BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(dir, "pydata/Germany", "county_divi_ma7.json"), - county, date, scaling_factor_icu)); - } - else { - log_warning("No DIVI data available for this date"); - } + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(dir, "pydata/Germany", "county_divi_ma7.json"), county, + date, scaling_factor_icu)); 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( @@ -454,13 +443,8 @@ IOResult read_input_data(std::vector& model, Date date, const std:: const std::vector& scaling_factor_inf, double scaling_factor_icu, const std::string& data_dir, int num_days = 0, bool export_time_series = false) { - if (date > Date(2020, 4, 23)) { - BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(data_dir, "critical_cases.json"), node_ids, date, - scaling_factor_icu)); - } - else { - log_warning("No DIVI data available for this date"); - } + BOOST_OUTCOME_TRY( + details::set_divi_data(model, path_join(data_dir, "critical_cases.json"), node_ids, date, scaling_factor_icu)); BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(data_dir, "confirmed_cases.json"), node_ids, date, scaling_factor_inf)); BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(data_dir, "population_data.json"), node_ids)); diff --git a/cpp/models/ode_secirts/parameters_io.h b/cpp/models/ode_secirts/parameters_io.h index 971ab51b5b..5431fb5576 100644 --- a/cpp/models/ode_secirts/parameters_io.h +++ b/cpp/models/ode_secirts/parameters_io.h @@ -100,7 +100,7 @@ IOResult compute_confirmed_cases_data( std::vector>& vnum_InfectedSevere, std::vector>& vnum_icu, std::vector>& vnum_death, std::vector>& vnum_timm_i, std::vector const& vregion, Date date, const std::vector& model, - const std::vector& scaling_factor_inf, const size_t layer, bool read_icu = false) + const std::vector& scaling_factor_inf, const size_t layer) { auto max_date_entry = std::max_element(case_data.begin(), case_data.end(), [](auto&& a, auto&& b) { return a.date < b.date; @@ -211,17 +211,13 @@ IOResult compute_confirmed_cases_data( } if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms - t_InfectedSevere)) { num_InfectedSevere[age] -= severePerInfectedSymptoms * scaling_factor_inf[age] * entry.num_confirmed; - if (read_icu) { - num_icu[age] += - severePerInfectedSymptoms * criticalPerSevere * scaling_factor_inf[age] * entry.num_confirmed; - } + num_icu[age] += + severePerInfectedSymptoms * criticalPerSevere * scaling_factor_inf[age] * entry.num_confirmed; } if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms - t_InfectedSevere - t_InfectedCritical)) { num_death[age] += entry.num_deaths; - if (read_icu) { - num_icu[age] -= - severePerInfectedSymptoms * criticalPerSevere * scaling_factor_inf[age] * entry.num_confirmed; - } + num_icu[age] -= + severePerInfectedSymptoms * criticalPerSevere * scaling_factor_inf[age] * entry.num_confirmed; } if (entry.date == offset_date_by_days(date, 0 - t_imm_interval_i)) { num_imm[age] -= entry.num_confirmed; @@ -305,13 +301,12 @@ IOResult read_confirmed_cases_data( std::vector>& vnum_InfectedNoSymptoms, std::vector>& vnum_InfectedSymptoms, std::vector>& vnum_InfectedSevere, std::vector>& vnum_icu, std::vector>& vnum_death, std::vector>& vnum_timm_i, - const std::vector& model, const std::vector& scaling_factor_inf, const size_t layer, - bool read_icu = false) + const std::vector& model, const std::vector& scaling_factor_inf, const size_t layer) { BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path)); return compute_confirmed_cases_data(case_data, vnum_Exposed, vnum_InfectedNoSymptoms, vnum_InfectedSymptoms, vnum_InfectedSevere, vnum_icu, vnum_death, vnum_timm_i, vregion, date, model, - scaling_factor_inf, layer, read_icu); + scaling_factor_inf, layer); } /** @@ -354,7 +349,7 @@ set_confirmed_cases_data(std::vector& model, const std::vector denom_E(num_age_groups, 0.0); std::vector denom_I_NS(num_age_groups, 0.0); std::vector denom_I_Sy(num_age_groups, 0.0); - std::vector denom_I_Sev(num_age_groups, 0.0); + std::vector denom_I_Sev_Cr(num_age_groups, 0.0); /*----------- Naive immunity -----------*/ for (size_t county = 0; county < model.size(); county++) { @@ -391,7 +386,7 @@ set_confirmed_cases_data(std::vector& model, const std::vector>()[(AgeGroup)group]); - denom_I_Sev[group] = + denom_I_Sev_Cr[group] = 1 / (immunity_population[0][group] + immunity_population[1][group] * model[county].parameters.template get>()[( @@ -418,7 +413,12 @@ set_confirmed_cases_data(std::vector& model, const std::vector& model, const std::vector>()[(AgeGroup)i] * - denom_I_Sev[i] * num_InfectedSevere[county][i]; + denom_I_Sev_Cr[i] * num_InfectedSevere[county][i]; + // Only set the number of ICU patients here, if the date is not available in the data. + if (!is_divi_data_available(date)) { + model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalPartialImmunity}] = + immunity_population[1][i] * + model[county] + .parameters.template get>()[(AgeGroup)i] * + denom_I_Sev_Cr[i] * num_icu[county][i]; + } // the += is necessary because we already set the previous vaccinated individuals model[county].populations[{AgeGroup(i), InfectionState::TemporaryImmunePartialImmunity}] += immunity_population[1][i] * @@ -511,7 +519,16 @@ set_confirmed_cases_data(std::vector& model, const std::vector>()[(AgeGroup)i] * - denom_I_Sev[i] * num_InfectedSevere[county][i]; + denom_I_Sev_Cr[i] * num_InfectedSevere[county][i]; + // Only set the number of ICU patients here, if the date is not available in the data. + if (!is_divi_data_available(date)) { + model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalImprovedImmunity}] = + immunity_population[2][i] * + model[county] + .parameters.template get>()[(AgeGroup)i] * + denom_I_Sev_Cr[i] * num_icu[county][i]; + } + // the += is necessary because we already set the previous vaccinated individuals model[county].populations[{AgeGroup(i), InfectionState::TemporaryImmuneImprovedImmunity}] += immunity_population[2][i] * @@ -644,6 +661,12 @@ template IOResult set_divi_data(std::vector& model, const std::string& path, const std::vector& vregion, Date date, FP scaling_factor_icu) { + // DIVI dataset will no longer be updated from CW29 2024 on. + if (!is_divi_data_available(date)) { + log_warning("No DIVI data available for date: {}-{}-{}", date.year, date.month, date.day, + ". ICU compartment will be set based on Case data."); + return success(); + } std::vector sum_mu_I_U(vregion.size(), 0); std::vector> mu_I_U{model.size()}; for (size_t region = 0; region < vregion.size(); region++) { @@ -993,13 +1016,7 @@ IOResult export_input_data_county_timeseries( // 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)); - } - else { - log_warning("No DIVI data available for for date: {}-{}-{}", offset_day.day, offset_day.month, - offset_day.year); - } + BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu)); BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(models, case_data, counties, offset_day, scaling_factor_inf, immunity_population)); @@ -1073,13 +1090,8 @@ IOResult read_input_data_county(std::vector& model, Date date, cons // 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 (date > Date(2020, 4, 23)) { - BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(dir, "pydata/Germany", "county_divi_ma7.json"), - county, date, scaling_factor_icu)); - } - else { - log_warning("No DIVI data available for this date"); - } + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(dir, "pydata/Germany", "county_divi_ma7.json"), county, + date, scaling_factor_icu)); BOOST_OUTCOME_TRY( details::set_confirmed_cases_data(model, path_join(dir, "pydata/Germany", "cases_all_county_age_ma7.json"), @@ -1136,13 +1148,8 @@ IOResult read_input_data(std::vector& model, Date date, const std:: // 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 (date > Date(2020, 4, 23)) { - BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(data_dir, "critical_cases.json"), node_ids, date, - scaling_factor_icu)); - } - else { - log_warning("No DIVI data available for this date"); - } + BOOST_OUTCOME_TRY( + details::set_divi_data(model, path_join(data_dir, "critical_cases.json"), node_ids, date, scaling_factor_icu)); BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(data_dir, "confirmed_cases.json"), node_ids, date, scaling_factor_inf, immunity_population)); diff --git a/cpp/models/ode_secirvvs/parameters_io.cpp b/cpp/models/ode_secirvvs/parameters_io.cpp index 86d28ff545..32132cbf8e 100644 --- a/cpp/models/ode_secirvvs/parameters_io.cpp +++ b/cpp/models/ode_secirvvs/parameters_io.cpp @@ -115,8 +115,6 @@ IOResult read_confirmed_cases_data( auto& mu_I_H = vmu_I_H[region_idx]; auto& mu_H_U = vmu_H_U[region_idx]; - bool read_icu = false; // params.populations.get({age, SecirCompartments::U}) == 0; - auto age = (size_t)entry.age_group; if (entry.date == offset_date_by_days(date, 0)) { num_InfectedSymptoms[age] += scaling_factor_inf[age] * entry.num_confirmed; @@ -138,16 +136,12 @@ IOResult read_confirmed_cases_data( } if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age])) { num_InfectedSevere[age] -= mu_I_H[age] * scaling_factor_inf[age] * entry.num_confirmed; - if (read_icu) { - num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed; - } + num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed; } if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age] - t_InfectedCritical[age])) { num_death[age] += entry.num_deaths; - if (read_icu) { - num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed; - } + num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed; } } } diff --git a/cpp/models/ode_secirvvs/parameters_io.h b/cpp/models/ode_secirvvs/parameters_io.h index d5f27c2f71..850ffff7de 100644 --- a/cpp/models/ode_secirvvs/parameters_io.h +++ b/cpp/models/ode_secirvvs/parameters_io.h @@ -177,6 +177,10 @@ IOResult set_confirmed_cases_data(std::vector& model, model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsNaiveConfirmed}] = 0; model[county].populations[{AgeGroup(i), InfectionState::InfectedSevereNaive}] = num_InfectedSevere[county][i]; + // Only set the number of ICU patients here, if the date is not available in the data. + if (!is_divi_data_available(date)) { + model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalNaive}] = num_icu[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 @@ -272,6 +276,11 @@ IOResult set_confirmed_cases_data(std::vector& model, model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsPartialImmunityConfirmed}] = 0; model[county].populations[{AgeGroup(i), InfectionState::InfectedSeverePartialImmunity}] = num_InfectedSevere[county][i]; + // Only set the number of ICU patients here, if the date is not available in the data. + if (!is_divi_data_available(date)) { + model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalPartialImmunity}] = + num_icu[county][i]; + } } // } if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) == 0) { @@ -343,7 +352,6 @@ IOResult set_confirmed_cases_data(std::vector& model, t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf)); for (size_t county = 0; county < model.size(); county++) { - // if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) > 0) { size_t num_groups = (size_t)model[county].parameters.get_num_groups(); for (size_t i = 0; i < num_groups; i++) { model[county].populations[{AgeGroup(i), InfectionState::ExposedImprovedImmunity}] = num_Exposed[county][i]; @@ -355,8 +363,12 @@ IOResult set_confirmed_cases_data(std::vector& model, model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] = 0; model[county].populations[{AgeGroup(i), InfectionState::InfectedSevereImprovedImmunity}] = num_InfectedSevere[county][i]; + // Only set the number of ICU patients here, if the date is not available in the data. + if (!is_divi_data_available(date)) { + model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalImprovedImmunity}] = + num_icu[county][i]; + } } - // } if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) == 0) { log_warning("No infections for vaccinated reported on date " + std::to_string(date.year) + "-" + std::to_string(date.month) + "-" + std::to_string(date.day) + " for region " + @@ -414,6 +426,12 @@ template IOResult set_divi_data(std::vector& model, const std::string& path, const std::vector& vregion, Date date, double scaling_factor_icu) { + // DIVI dataset will no longer be updated from CW29 2024 on. + if (!is_divi_data_available(date)) { + log_warning("No DIVI data available for date: {}-{}-{}", date.year, date.month, date.day, + ". ICU compartment will be set based on Case data."); + return success(); + } std::vector sum_mu_I_U(vregion.size(), 0); std::vector> mu_I_U{model.size()}; for (size_t region = 0; region < vregion.size(); region++) { @@ -872,13 +890,7 @@ IOResult export_input_data_county_timeseries( // 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)); - } - else { - log_warning("No DIVI data available for for date: {}-{}-{}", offset_day.day, offset_day.month, - offset_day.year); - } + BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu)); BOOST_OUTCOME_TRY( details::set_confirmed_cases_data(models, case_data, counties, offset_day, scaling_factor_inf, true)); @@ -941,13 +953,8 @@ IOResult read_input_data_county(std::vector& model, Date date, cons // 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 (date > Date(2020, 4, 23)) { - BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(dir, "pydata/Germany", "county_divi_ma7.json"), - county, date, scaling_factor_icu)); - } - else { - log_warning("No DIVI data available for this date"); - } + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(dir, "pydata/Germany", "county_divi_ma7.json"), county, + date, scaling_factor_icu)); 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)); @@ -997,13 +1004,8 @@ IOResult read_input_data(std::vector& model, Date date, const std:: // 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 (date > Date(2020, 4, 23)) { - BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(data_dir, "critical_cases.json"), node_ids, date, - scaling_factor_icu)); - } - else { - log_warning("No DIVI data available for this date"); - } + BOOST_OUTCOME_TRY( + details::set_divi_data(model, path_join(data_dir, "critical_cases.json"), node_ids, date, scaling_factor_icu)); BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(data_dir, "confirmed_cases.json"), node_ids, date, scaling_factor_inf)); diff --git a/cpp/tests/data/export_time_series_initialization_osecirvvs.h5 b/cpp/tests/data/export_time_series_initialization_osecirvvs.h5 index 3eaea56319..9d4ac4d8ed 100644 Binary files a/cpp/tests/data/export_time_series_initialization_osecirvvs.h5 and b/cpp/tests/data/export_time_series_initialization_osecirvvs.h5 differ diff --git a/cpp/tests/test_epi_data_io.cpp b/cpp/tests/test_epi_data_io.cpp index 18934a8603..0797365323 100644 --- a/cpp/tests/test_epi_data_io.cpp +++ b/cpp/tests/test_epi_data_io.cpp @@ -297,6 +297,15 @@ TEST(TestEpiDataIo, read_divi_data) ASSERT_EQ(divi_data[3].num_icu, 0.35437); } +TEST(TestEpiDataIo, is_divi_data_available) +{ + EXPECT_FALSE(mio::is_divi_data_available(mio::Date(2020, 4, 22))); // Before start + EXPECT_FALSE(mio::is_divi_data_available(mio::Date(2024, 7, 22))); // After end + EXPECT_TRUE(mio::is_divi_data_available(mio::Date(2020, 4, 23))); // Start date + EXPECT_TRUE(mio::is_divi_data_available(mio::Date(2022, 1, 1))); // Inside range + EXPECT_TRUE(mio::is_divi_data_available(mio::Date(2024, 7, 21))); // End date +} + TEST(TestEpiDataIo, read_confirmed_cases_data) { auto case_data = mio::read_confirmed_cases_data(mio::path_join(TEST_DATA_DIR, "test_cases_all_age.json")).value(); diff --git a/cpp/tests/test_odesecir.cpp b/cpp/tests/test_odesecir.cpp index 88c6599f9f..bc46600c8f 100644 --- a/cpp/tests/test_odesecir.cpp +++ b/cpp/tests/test_odesecir.cpp @@ -1434,5 +1434,71 @@ TEST_F(ModelTestOdeSecir, model_initialization_old_date) mio::set_log_level(mio::LogLevel::warn); } +TEST(TestOdeSecir, set_divi_data_invalid_dates) +{ + mio::set_log_level(mio::LogLevel::off); + auto model = mio::osecir::Model(1); + model.populations.array().setConstant(1); + auto model_vector = std::vector>{model}; + + // Test with date before DIVI dataset was available. + EXPECT_THAT(mio::osecir::details::set_divi_data(model_vector, "", {1001}, {2019, 12, 01}, 1.0), IsSuccess()); + // Assure that populations is the same as before. + EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), + MatrixNear(print_wrap(model.populations.array().cast()), 1e-10, 1e-10)); + + // Test with data after DIVI dataset was no longer updated. + EXPECT_THAT(mio::osecir::details::set_divi_data(model_vector, "", {1001}, {2025, 12, 01}, 1.0), IsSuccess()); + EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), + MatrixNear(print_wrap(model.populations.array().cast()), 1e-10, 1e-10)); + + mio::set_log_level(mio::LogLevel::warn); +} + +TEST_F(ModelTestOdeSecir, set_confirmed_cases_data_with_ICU) +{ + // set params + for (auto age_group = mio::AgeGroup(0); age_group < (mio::AgeGroup)num_age_groups; age_group++) { + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + } + + // read case data + auto case_data = + mio::read_confirmed_cases_data(mio::path_join(TEST_DATA_DIR, "cases_all_county_age_ma7.json")).value(); + + // Change dates of the case data so that no ICU data is available at that time. + // Also, increase the number of confirmed cases by 1 each day. + const auto t0 = mio::Date(2025, 1, 1); + auto day_add = 0; + for (auto& entry : case_data) { + entry.date = offset_date_by_days(t0, day_add); + entry.num_confirmed = day_add; + day_add++; + } + + // get day in mid of the data + auto mid_day = case_data[(size_t)case_data.size() / 2].date; + + // calculate ICU values using set_confirmed_cases_data + auto model_vector = std::vector>{model}; + auto scaling_factor_inf = std::vector(size_t(model.parameters.get_num_groups()), 1.0); + EXPECT_THAT( + mio::osecir::details::set_confirmed_cases_data(model_vector, case_data, {1002}, mid_day, scaling_factor_inf), + IsSuccess()); + + // Since, TimeInfectedCritical is 1, the number of ICU cases is the difference of confirmed cases between two days, which is 1. + // We only have an entry for age group 2. All other age groups should be zero. + for (size_t i = 0; i < num_age_groups; ++i) { + const auto expected_value = (i == 2) ? 1.0 : 0.0; + const auto actual_value = + model_vector[0].populations[{mio::AgeGroup(i), mio::osecir::InfectionState::InfectedCritical}].value(); + EXPECT_NEAR(actual_value, expected_value, 1e-10); + } +} + #endif #endif diff --git a/cpp/tests/test_odesecirts.cpp b/cpp/tests/test_odesecirts.cpp index a95268cccd..bbdf7976dc 100644 --- a/cpp/tests/test_odesecirts.cpp +++ b/cpp/tests/test_odesecirts.cpp @@ -737,14 +737,14 @@ TEST(TestOdeSECIRTS, read_confirmed_cases) ASSERT_THAT(mio::osecirts::details::read_confirmed_cases_data( path, region, {2020, 12, 01}, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_timm, model, - std::vector(size_t(num_age_groups), 1.0), 0, true), + std::vector(size_t(num_age_groups), 1.0), 0), IsSuccess()); // read again with invalid date ASSERT_THAT(mio::osecirts::details::read_confirmed_cases_data( path, region, {3020, 12, 01}, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_timm, model, - std::vector(size_t(num_age_groups), 1.0), 0, true), + std::vector(size_t(num_age_groups), 1.0), 0), IsFailure(mio::StatusCode::OutOfRange)); // call the compute function with empty case data @@ -752,10 +752,100 @@ TEST(TestOdeSECIRTS, read_confirmed_cases) ASSERT_THAT(mio::osecirts::details::compute_confirmed_cases_data( empty_case_data, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_timm, region, {2020, 12, 01}, model, - std::vector(size_t(num_age_groups), 1.0), 0, true), + std::vector(size_t(num_age_groups), 1.0), 0), IsFailure(mio::StatusCode::InvalidValue)); } +TEST(TestOdeSECIRTS, set_divi_data_invalid_dates) +{ + mio::set_log_level(mio::LogLevel::off); + auto model = mio::osecirts::Model(1); + model.populations.array().setConstant(1); + auto model_vector = std::vector>{model}; + + // Test with date before DIVI dataset was available. + EXPECT_THAT(mio::osecirts::details::set_divi_data(model_vector, "", {1001}, {2019, 12, 01}, 1.0), IsSuccess()); + // Assure that populations is the same as before. + EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), + MatrixNear(print_wrap(model.populations.array().cast()), 1e-10, 1e-10)); + + // Test with data after DIVI dataset was no longer updated. + EXPECT_THAT(mio::osecirts::details::set_divi_data(model_vector, "", {1001}, {2025, 12, 01}, 1.0), IsSuccess()); + EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), + MatrixNear(print_wrap(model.populations.array().cast()), 1e-10, 1e-10)); + + mio::set_log_level(mio::LogLevel::warn); +} + +TEST(TestOdeSECIRTS, set_confirmed_cases_data_with_ICU) +{ + const auto num_age_groups = 6; + auto model = mio::osecirts::Model(num_age_groups); + model.populations.array().setConstant(1); + + const std::vector> immunity_population = {{0.04, 0.04, 0.075, 0.08, 0.035, 0.01}, + {0.61, 0.61, 0.62, 0.62, 0.58, 0.41}, + {0.35, 0.35, 0.305, 0.3, 0.385, 0.58}}; + + // set params + for (auto age_group = mio::AgeGroup(0); age_group < (mio::AgeGroup)num_age_groups; age_group++) { + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + } + + // read case data + auto case_data = + mio::read_confirmed_cases_data(mio::path_join(TEST_DATA_DIR, "cases_all_county_age_ma7.json")).value(); + + // Change dates of the case data so that no ICU data is available at that time. + // Also, increase the number of confirmed cases by 1 each day. + const auto t0 = mio::Date(2025, 1, 1); + auto day_add = 0; + for (auto& entry : case_data) { + entry.date = offset_date_by_days(t0, day_add); + entry.num_confirmed = day_add; + day_add++; + } + + // get day in mid of the data + auto mid_day = case_data[(size_t)case_data.size() / 2].date; + + // calculate ICU values using set_confirmed_cases_data + auto model_vector = std::vector>{model}; + auto scaling_factor_inf = std::vector(size_t(model.parameters.get_num_groups()), 1.0); + + EXPECT_THAT(mio::osecirts::details::set_confirmed_cases_data(model_vector, case_data, {1002}, mid_day, + scaling_factor_inf, immunity_population), + IsSuccess()); + + // Since, TimeInfectedCritical is 1, the number of ICU cases is the difference of confirmed cases between two days, which is 1. + // We only have an entry for age group 2. All other age groups should be zero. + for (int i = 0; i < num_age_groups; ++i) { + const auto expected_value = (i == 2) ? 1.0 : 0.0; + + auto actual_value_naive = + model_vector[0] + .populations[{mio::AgeGroup(i), mio::osecirts::InfectionState::InfectedCriticalNaive}] + .value(); + EXPECT_NEAR(actual_value_naive, expected_value * immunity_population[0][i], 1e-10); + + auto actual_value_pi = + model_vector[0] + .populations[{mio::AgeGroup(i), mio::osecirts::InfectionState::InfectedCriticalPartialImmunity}] + .value(); + EXPECT_NEAR(actual_value_pi, expected_value * immunity_population[1][i], 1e-10); + + auto actual_value_ii = + model_vector[0] + .populations[{mio::AgeGroup(i), mio::osecirts::InfectionState::InfectedCriticalImprovedImmunity}] + .value(); + EXPECT_NEAR(actual_value_ii, expected_value * immunity_population[2][i], 1e-10); + } +} + TEST(TestOdeSECIRTS, read_data) { auto num_age_groups = 6; //reading data requires RKI data age groups diff --git a/cpp/tests/test_odesecirvvs.cpp b/cpp/tests/test_odesecirvvs.cpp index dfdeb95727..2a26ae8aee 100644 --- a/cpp/tests/test_odesecirvvs.cpp +++ b/cpp/tests/test_odesecirvvs.cpp @@ -570,6 +570,91 @@ TEST(TestOdeSECIRVVS, read_confirmed_cases) ASSERT_THAT(read, IsSuccess()); } +TEST(TestOdeSECIRVVS, set_divi_data_invalid_dates) +{ + mio::set_log_level(mio::LogLevel::off); + auto model = mio::osecirvvs::Model(1); + model.populations.array().setConstant(1); + auto model_vector = std::vector>{model}; + + // Test with date before DIVI dataset was available. + EXPECT_THAT(mio::osecirvvs::details::set_divi_data(model_vector, "", {1001}, {2019, 12, 01}, 1.0), IsSuccess()); + // Assure that populations is the same as before. + EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), + MatrixNear(print_wrap(model.populations.array().cast()), 1e-10, 1e-10)); + + // Test with data after DIVI dataset was no longer updated. + EXPECT_THAT(mio::osecirvvs::details::set_divi_data(model_vector, "", {1001}, {2025, 12, 01}, 1.0), IsSuccess()); + EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast()), + MatrixNear(print_wrap(model.populations.array().cast()), 1e-10, 1e-10)); + + mio::set_log_level(mio::LogLevel::warn); +} + +TEST(TestOdeSECIRVVS, set_confirmed_cases_data_with_ICU) +{ + const auto num_age_groups = 6; + auto model = mio::osecirvvs::Model(num_age_groups); + model.populations.array().setConstant(1); + + // set params + for (auto age_group = mio::AgeGroup(0); age_group < (mio::AgeGroup)num_age_groups; age_group++) { + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + model.parameters.get>()[age_group] = 1.0; + } + + // read case data + auto case_data = + mio::read_confirmed_cases_data(mio::path_join(TEST_DATA_DIR, "cases_all_county_age_ma7.json")).value(); + + // Change dates of the case data so that no ICU data is available at that time. + // Also, increase the number of confirmed cases by 1 each day. + const auto t0 = mio::Date(2025, 1, 1); + auto day_add = 0; + for (auto& entry : case_data) { + entry.date = offset_date_by_days(t0, day_add); + entry.num_confirmed = day_add; + day_add++; + } + + // get day in mid of the data + auto mid_day = case_data[(size_t)case_data.size() / 2].date; + + // calculate ICU values using set_confirmed_cases_data + auto model_vector = std::vector>{model}; + auto scaling_factor_inf = std::vector(size_t(model.parameters.get_num_groups()), 1.0); + EXPECT_THAT( + mio::osecirvvs::details::set_confirmed_cases_data(model_vector, case_data, {1002}, mid_day, scaling_factor_inf), + IsSuccess()); + + // Since, TimeInfectedCritical is 1, the number of ICU cases is the difference of confirmed cases between two days, which is 1. + // We only have an entry for age group 2. All other age groups should be zero. + for (int i = 0; i < num_age_groups; ++i) { + const auto expected_value = (i == 2) ? 1.0 : 0.0; + + auto actual_value_naive = + model_vector[0] + .populations[{mio::AgeGroup(i), mio::osecirvvs::InfectionState::InfectedCriticalNaive}] + .value(); + EXPECT_NEAR(actual_value_naive, expected_value, 1e-10); + + auto actual_value_pi = + model_vector[0] + .populations[{mio::AgeGroup(i), mio::osecirvvs::InfectionState::InfectedCriticalPartialImmunity}] + .value(); + EXPECT_NEAR(actual_value_pi, expected_value, 1e-10); + + auto actual_value_ii = + model_vector[0] + .populations[{mio::AgeGroup(i), mio::osecirvvs::InfectionState::InfectedCriticalImprovedImmunity}] + .value(); + EXPECT_NEAR(actual_value_ii, expected_value, 1e-10); + } +} + TEST(TestOdeSECIRVVS, read_data) { auto num_age_groups = 6; //reading data requires RKI data age groups @@ -595,17 +680,17 @@ TEST(TestOdeSECIRVVS, read_data) auto expected_values = (Eigen::ArrayXd(num_age_groups * Eigen::Index(mio::osecirvvs::InfectionState::Count)) << 8792.15, 175.889, 3.21484, 0.0633116, 0.221057, 1.42882, 0.0351731, 0.29682, 0, 0, 0, 6.93838, 0.0725173, 0.206715, 0, 0, 0, - 0.0337498, 1.23324e-05, 0.000208293, 0.0292822, 5.8568e-05, 0.000406386, 1340.42, 0, 0, 0, 17067.6, 220.137, + 0.0337498, 1.23324e-05, 0.000208293, 0.0292822, 5.8568e-05, 0.000406386, 1340.42, 0, 0, 0, 17067.7, 220.137, 7.64078, 0.0970237, 0.381933, 4.91193, 0.0779655, 0.741778, 0, 0, 0, 11.7286, 0.0890643, 0.286235, 0, 0, 0, - 0.0434344, 8.40756e-06, 0.000160098, 0.0294125, 3.7932e-05, 0.000296738, 1891.19, 0, 0, 0, 72501, 176.267, + 0.0434344, 8.40756e-06, 0.000160098, 0.0294125, 3.7932e-05, 0.000296738, 1891.18, 0, 0, 0, 72501, 176.267, 47.227, 0.113013, 0.490073, 24.4094, 0.0730141, 0.765246, 0, 0, 0, 64.6789, 0.0855947, 0.303032, 0, 0, 0, - 1.23754, 4.5968e-05, 0.000964262, 0.0751837, 1.82724e-05, 0.000157466, 1670.26, 0, 0, 0, 80790.1, 184.645, + 1.23754, 4.5968e-05, 0.000964262, 0.0751837, 1.82724e-05, 0.000157466, 1670.26, 0, 0, 0, 80791.1, 184.645, 44.5477, 0.100229, 0.50512, 23.6881, 0.0666206, 0.811467, 0, 0, 0, 58.9805, 0.0758111, 0.31192, 0, 0, 0, - 3.75961, 0.000136175, 0.00331973, 0.486628, 0.000111199, 0.00111367, 2022.58, 0, 0, 0, 41581, 177.478, 9.27393, - 0.0389771, 0.216151, 5.77433, 0.030336, 0.4066, 0, 0, 0, 13.3664, 0.0312302, 0.141394, 0, 0, 0, 3.119, - 0.000209444, 0.00561852, 2.60439, 0.00111169, 0.0122515, 2136.6, 0, 0, 0, 13223.8, 216.037, 11.1838, 0.179986, + 3.75961, 0.000136175, 0.00331973, 0.486628, 0.000111199, 0.00111367, 2021.62, 0, 0, 0, 41582.5, 177.478, + 9.27393, 0.0389771, 0.216151, 5.77433, 0.030336, 0.4066, 0, 0, 0, 13.3664, 0.0312302, 0.141394, 0, 0, 0, 3.119, + 0.000209444, 0.00561852, 2.60439, 0.00111169, 0.0122515, 2135.09, 0, 0, 0, 13224.1, 216.037, 11.1838, 0.179986, 0.863926, 3.50537, 0.0705169, 0.818075, 0, 0, 0, 3.52982, 0.0331744, 0.130002, 0, 0, 0, 0.695168, 0.000190699, - 0.00442784, 4.67895, 0.00764769, 0.0729502, 2253.61, 0, 0, 0) + 0.00442784, 4.67895, 0.00764769, 0.0729502, 2253.25, 0, 0, 0) .finished(); ASSERT_THAT(print_wrap(model1[0].populations.array().cast()), @@ -830,11 +915,11 @@ TEST(TestOdeSECIRVVS, model_initialization) 0.000325876, 0.0011537, 0, 0, 0, 1.24042, 1.74173e-07, 3.65358e-06, 0.0753588, 6.9234e-08, 5.96637e-07, 1671.81, 0, 0, 0, 3.00317e+07, 184.888, 44.9988, 0.000272769, 0.00137466, 23.9279, 0.000181305, 0.00220837, 0, 0, 0, 59.4274, 0.000205796, 0.000846734, 0, 0, 0, 3.76905, 3.67799e-07, 8.9664e-06, 0.48785, 3.00341e-07, - 3.00797e-06, 2022.51, 0, 0, 0, 1.65123e+07, 177.579, 9.4638, 0.000100211, 0.00055573, 5.89255, 7.79946e-05, + 3.00797e-06, 2021.56, 0, 0, 0, 1.65123e+07, 177.579, 9.4638, 0.000100211, 0.00055573, 5.89255, 7.79946e-05, 0.00104538, 0, 0, 0, 13.5709, 7.98864e-05, 0.000361685, 0, 0, 0, 3.13496, 5.30384e-07, 1.4228e-05, 2.61772, - 2.81518e-06, 3.1025e-05, 2136.56, 0, 0, 0, 6.17983e+06, 216.328, 11.9625, 0.000412312, 0.00197908, 3.74944, + 2.81518e-06, 3.1025e-05, 2135.05, 0, 0, 0, 6.17984e+06, 216.328, 11.9625, 0.000412312, 0.00197908, 3.74944, 0.00016154, 0.00187405, 0, 0, 0, 3.71387, 7.47535e-05, 0.000292941, 0, 0, 0, 0.707117, 4.15435e-07, - 9.64602e-06, 4.75937, 1.66604e-05, 0.000158922, 2253.59, 0, 0, 0) + 9.64602e-06, 4.75937, 1.66604e-05, 0.000158922, 2253.23, 0, 0, 0) .finished(); ASSERT_THAT(print_wrap(model_vector[0].populations.array().cast()),