diff --git a/cpp/models/ode_secir/parameters_io.h b/cpp/models/ode_secir/parameters_io.h index ab334d922d..8471e25a00 100644 --- a/cpp/models/ode_secir/parameters_io.h +++ b/cpp/models/ode_secir/parameters_io.h @@ -75,7 +75,17 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect const std::vector& scaling_factor_inf) { const size_t num_age_groups = ConfirmedCasesDataEntry::age_group_names.size(); - assert(scaling_factor_inf.size() == num_age_groups); + // allow single scalar scaling that is broadcast to all age groups + assert(scaling_factor_inf.size() == 1 || scaling_factor_inf.size() == num_age_groups); + + // Set scaling factors to match num age groups + std::vector scaling_factor_inf_full; + if (scaling_factor_inf.size() == 1) { + scaling_factor_inf_full.assign(num_age_groups, scaling_factor_inf[0]); + } + else { + scaling_factor_inf_full = scaling_factor_inf; + } std::vector> t_InfectedNoSymptoms{model.size()}; std::vector> t_Exposed{model.size()}; @@ -89,7 +99,12 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect std::vector> mu_U_D{model.size()}; for (size_t node = 0; node < model.size(); ++node) { - for (size_t group = 0; group < num_age_groups; group++) { + const size_t model_groups = (size_t)model[node].parameters.get_num_groups(); + assert(model_groups == 1 || model_groups == num_age_groups); + for (size_t ag = 0; ag < num_age_groups; ag++) { + // If the model has fewer groups than casedata entries available, + // reuse group 0 parameters for all RKI age groups + const size_t group = (model_groups == num_age_groups) ? ag : 0; t_Exposed[node].push_back( static_cast(std::round(model[node].parameters.template get>()[(AgeGroup)group]))); @@ -121,26 +136,49 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect 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)); + t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf_full)); for (size_t node = 0; node < model.size(); node++) { if (std::accumulate(num_InfectedSymptoms[node].begin(), num_InfectedSymptoms[node].end(), 0.0) > 0) { size_t num_groups = (size_t)model[node].parameters.get_num_groups(); - for (size_t i = 0; i < num_groups; i++) { - model[node].populations[{AgeGroup(i), InfectionState::Exposed}] = num_Exposed[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptoms}] = - num_InfectedNoSymptoms[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsConfirmed}] = 0; - model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptoms}] = - num_InfectedSymptoms[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptomsConfirmed}] = 0; - model[node].populations[{AgeGroup(i), InfectionState::InfectedSevere}] = num_InfectedSevere[node][i]; - // Only set the number of ICU patients here, if the date is not available in the data. + if (num_groups == num_age_groups) { + for (size_t i = 0; i < num_groups; i++) { + model[node].populations[{AgeGroup(i), InfectionState::Exposed}] = num_Exposed[node][i]; + model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptoms}] = + num_InfectedNoSymptoms[node][i]; + model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptoms}] = + num_InfectedSymptoms[node][i]; + model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptomsConfirmed}] = 0; + model[node].populations[{AgeGroup(i), InfectionState::InfectedSevere}] = + num_InfectedSevere[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 { + const auto sum_vec = [](const std::vector& v) { + return std::accumulate(v.begin(), v.end(), 0.0); + }; + const size_t i0 = 0; + model[node].populations[{AgeGroup(i0), InfectionState::Exposed}] = sum_vec(num_Exposed[node]); + model[node].populations[{AgeGroup(i0), InfectionState::InfectedNoSymptoms}] = + sum_vec(num_InfectedNoSymptoms[node]); + model[node].populations[{AgeGroup(i0), InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model[node].populations[{AgeGroup(i0), InfectionState::InfectedSymptoms}] = + sum_vec(num_InfectedSymptoms[node]); + model[node].populations[{AgeGroup(i0), InfectionState::InfectedSymptomsConfirmed}] = 0; + model[node].populations[{AgeGroup(i0), InfectionState::InfectedSevere}] = + sum_vec(num_InfectedSevere[node]); if (!is_divi_data_available(date)) { - model[node].populations[{AgeGroup(i), InfectionState::InfectedCritical}] = num_icu[node][i]; + model[node].populations[{AgeGroup(i0), InfectionState::InfectedCritical}] = sum_vec(num_icu[node]); } - model[node].populations[{AgeGroup(i), InfectionState::Dead}] = num_death[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::Recovered}] = num_rec[node][i]; + model[node].populations[{AgeGroup(i0), InfectionState::Dead}] = sum_vec(num_death[node]); + model[node].populations[{AgeGroup(i0), InfectionState::Recovered}] = sum_vec(num_rec[node]); } } else { @@ -231,10 +269,20 @@ IOResult set_population_data(std::vector>& model, assert(num_population.size() == vregion.size()); assert(model.size() == vregion.size()); for (size_t region = 0; region < vregion.size(); region++) { - auto num_groups = model[region].parameters.get_num_groups(); - for (auto i = AgeGroup(0); i < num_groups; i++) { + const auto model_groups = (size_t)model[region].parameters.get_num_groups(); + const auto data_groups = num_population[region].size(); + assert(data_groups == model_groups || (model_groups == 1 && data_groups >= 1)); + + if (data_groups == model_groups) { + for (auto i = AgeGroup(0); i < model[region].parameters.get_num_groups(); i++) { + model[region].populations.template set_difference_from_group_total( + {i, InfectionState::Susceptible}, num_population[region][(size_t)i]); + } + } + else if (model_groups == 1 && data_groups >= 1) { + const double total = std::accumulate(num_population[region].begin(), num_population[region].end(), 0.0); model[region].populations.template set_difference_from_group_total( - {i, InfectionState::Susceptible}, num_population[region][size_t(i)]); + {AgeGroup(0), InfectionState::Susceptible}, total); } } return success(); @@ -283,8 +331,8 @@ IOResult export_input_data_county_timeseries( const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path) { const auto num_age_groups = (size_t)models[0].parameters.get_num_groups(); - assert(scaling_factor_inf.size() == num_age_groups); - assert(num_age_groups == ConfirmedCasesDataEntry::age_group_names.size()); + // allow scalar scaling factor as convenience for 1-group models + assert(scaling_factor_inf.size() == 1 || scaling_factor_inf.size() == num_age_groups); assert(models.size() == region.size()); std::vector> extrapolated_data( region.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_age_groups)); diff --git a/cpp/tests/test_odesecir.cpp b/cpp/tests/test_odesecir.cpp index 023c56fc6b..010fc84202 100644 --- a/cpp/tests/test_odesecir.cpp +++ b/cpp/tests/test_odesecir.cpp @@ -29,6 +29,7 @@ #include "memilio/io/parameters_io.h" #include "memilio/data/analyze_result.h" #include "memilio/math/adapt_rk.h" +#include "memilio/geography/regions.h" #include @@ -1518,5 +1519,184 @@ TEST(TestOdeSecir, read_population_data_failure) EXPECT_EQ(result.error().message(), "File with county population expected."); } +TEST(TestOdeSecirIO, read_input_data_county_aggregates_one_group) +{ + // Set up two models with different number of age groups. + const size_t num_age_groups = 6; + std::vector> models6{mio::osecir::Model((int)num_age_groups)}; + std::vector> models1{mio::osecir::Model(1)}; + + // Relevant parameters for model with 6 age groups + for (auto i = mio::AgeGroup(0); i < (mio::AgeGroup)num_age_groups; ++i) { + models6[0].parameters.get>()[i] = 0.2; + models6[0].parameters.get>()[i] = 0.25; + } + + // Relevant parameters for model with 1 age group + models1[0].parameters.get>()[mio::AgeGroup(0)] = 0.2; + models1[0].parameters.get>()[mio::AgeGroup(0)] = 0.25; + + const auto pydata_dir_Germany = mio::path_join(TEST_DATA_DIR, "Germany", "pydata"); + const std::vector counties{1002}; + const auto date = mio::Date(2020, 12, 1); + + std::vector scale6(num_age_groups, 1.0); + std::vector scale1{1.0}; + + // Initialize both models + ASSERT_THAT(mio::osecir::read_input_data_county(models6, date, counties, scale6, 1.0, pydata_dir_Germany), + IsSuccess()); + ASSERT_THAT(mio::osecir::read_input_data_county(models1, date, counties, scale1, 1.0, pydata_dir_Germany), + IsSuccess()); + + // Aggreagate the results from the model with 6 age groups and compare with the model with 1 age group + const auto& m6 = models6[0]; + const auto& m1 = models1[0]; + const double tol = 1e-10; + for (int s = 0; s < (int)mio::osecir::InfectionState::Count; ++s) { + double sum6 = 0.0; + for (size_t ag = 0; ag < num_age_groups; ++ag) { + sum6 += m6.populations[{mio::AgeGroup(ag), (mio::osecir::InfectionState)s}].value(); + } + const double v1 = m1.populations[{mio::AgeGroup(0), (mio::osecir::InfectionState)s}].value(); + EXPECT_NEAR(sum6, v1, tol); + } + + // Total population + EXPECT_NEAR(m6.populations.get_total(), m1.populations.get_total(), tol); +} + +TEST(TestOdeSecirIO, set_population_data_single_age_group) +{ + const size_t num_age_groups = 6; + + // Create two models: one with 6 age groups, one with 1 age group + std::vector> models6{mio::osecir::Model((int)num_age_groups)}; + std::vector> models1{mio::osecir::Model(1)}; + + // Test population data with 6 different values for age groups + std::vector> population_data6 = {{10000.0, 20000.0, 30000.0, 25000.0, 15000.0, 8000.0}}; + std::vector> population_data1 = {{108000.0}}; // sum of all age groups + std::vector regions = {1002}; + + // Set population data for both models + EXPECT_THAT(mio::osecir::details::set_population_data(models6, population_data6, regions), IsSuccess()); + EXPECT_THAT(mio::osecir::details::set_population_data(models1, population_data1, regions), IsSuccess()); + + // Sum all compartments across age groups in 6-group model and compare 1-group model + const double tol = 1e-10; + for (int s = 0; s < (int)mio::osecir::InfectionState::Count; ++s) { + double sum6 = 0.0; + for (size_t ag = 0; ag < num_age_groups; ++ag) { + sum6 += models6[0].populations[{mio::AgeGroup(ag), (mio::osecir::InfectionState)s}].value(); + } + double val1 = models1[0].populations[{mio::AgeGroup(0), (mio::osecir::InfectionState)s}].value(); + + EXPECT_NEAR(sum6, val1, tol); + } + + // Total population should also match + EXPECT_NEAR(models6[0].populations.get_total(), models1[0].populations.get_total(), tol); +} + +TEST(TestOdeSecirIO, set_confirmed_cases_data_single_age_group) +{ + const size_t num_age_groups = 6; + + // Create two models: one with 6 age groups, one with 1 age group + std::vector> models6{mio::osecir::Model((int)num_age_groups)}; + std::vector> models1{mio::osecir::Model(1)}; + + // Create case data for all 6 age groups over multiple days (current day + 6 days back) + std::vector case_data; + + for (int day_offset = -6; day_offset <= 0; ++day_offset) { + mio::Date current_date = mio::offset_date_by_days(mio::Date(2020, 12, 1), day_offset); + + for (int age_group = 0; age_group < 6; ++age_group) { + double base_confirmed = 80.0 + age_group * 8.0 + (day_offset + 6) * 5.0; + double base_recovered = 40.0 + age_group * 4.0 + (day_offset + 6) * 3.0; + double base_deaths = 3.0 + age_group * 0.5 + (day_offset + 6) * 0.5; + + mio::ConfirmedCasesDataEntry entry{base_confirmed, + base_recovered, + base_deaths, + current_date, + mio::AgeGroup(age_group), + {}, + mio::regions::CountyId(1002), + {}}; + case_data.push_back(entry); + } + } + + std::vector regions = {1002}; + std::vector scaling_factors = {1.0}; + + // Set confirmed cases data for both models + EXPECT_THAT(mio::osecir::details::set_confirmed_cases_data(models6, case_data, regions, mio::Date(2020, 12, 1), + scaling_factors), + IsSuccess()); + EXPECT_THAT(mio::osecir::details::set_confirmed_cases_data(models1, case_data, regions, mio::Date(2020, 12, 1), + scaling_factors), + IsSuccess()); + + // Sum all compartments across age groups in 6-group model should be equal to 1-group model + for (int s = 0; s < (int)mio::osecir::InfectionState::Count; ++s) { + double sum6 = 0.0; + for (size_t ag = 0; ag < num_age_groups; ++ag) { + sum6 += models6[0].populations[{mio::AgeGroup(ag), (mio::osecir::InfectionState)s}].value(); + } + + double val1 = models1[0].populations[{mio::AgeGroup(0), (mio::osecir::InfectionState)s}].value(); + + EXPECT_NEAR(sum6, val1, 1e-10); + } + + // Total population + EXPECT_NEAR(models6[0].populations.get_total(), models1[0].populations.get_total(), 1e-10); +} + +TEST(TestOdeSecirIO, set_divi_data_single_age_group) +{ + // Create models with 6 age groups and 1 age group + std::vector> models_6_groups{mio::osecir::Model(6)}; + std::vector> models_1_group{mio::osecir::Model(1)}; + + // Set relevant parameters for all age groups + for (int i = 0; i < 6; i++) { + models_6_groups[0].parameters.get>()[mio::AgeGroup(i)] = 0.2; + models_6_groups[0].parameters.get>()[mio::AgeGroup(i)] = 0.25; + } + + // Set relevant parameters for 1 age group model + models_1_group[0].parameters.get>()[mio::AgeGroup(0)] = 0.2; + models_1_group[0].parameters.get>()[mio::AgeGroup(0)] = 0.25; + + // Apply DIVI data to both models + std::vector regions = {1002}; + double scaling_factor_icu = 1.0; + mio::Date date(2020, 12, 1); + std::string divi_data_path = mio::path_join(TEST_DATA_DIR, "Germany", "pydata", "county_divi_ma7.json"); + auto result_6_groups = + mio::osecir::details::set_divi_data(models_6_groups, divi_data_path, regions, date, scaling_factor_icu); + auto result_1_group = + mio::osecir::details::set_divi_data(models_1_group, divi_data_path, regions, date, scaling_factor_icu); + + EXPECT_THAT(result_6_groups, IsSuccess()); + EXPECT_THAT(result_1_group, IsSuccess()); + + // Calculate totals after applying DIVI data + double total_icu_6_groups_after = 0.0; + for (int i = 0; i < 6; i++) { + total_icu_6_groups_after += + models_6_groups[0].populations[{mio::AgeGroup(i), mio::osecir::InfectionState::InfectedCritical}].value(); + } + double icu_1_group_after = + models_1_group[0].populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}].value(); + + EXPECT_NEAR(total_icu_6_groups_after, icu_1_group_after, 1e-10); +} + #endif #endif diff --git a/docs/requirements.txt b/docs/requirements.txt index 4de20881c1..724674aef9 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,3 +1,6 @@ +cmake>=3.26 +ninja +scikit-build>=0.18 sphinx==7.1.2 sphinx-rtd-theme==1.3.0rc1 sphinx-copybutton