1315 support IO functionality with single age group for ODE SECIR

cpp/models/ode_secir/parameters_io.h

@@ -75,7 +75,17 @@ IOResult<void> set_confirmed_cases_data(std::vector<Model<FP>>& model, std::vect
                                         const std::vector<double>& 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<double> 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<std::vector<int>> t_InfectedNoSymptoms{model.size()};
     std::vector<std::vector<int>> t_Exposed{model.size()};
@@ -89,7 +99,12 @@ IOResult<void> set_confirmed_cases_data(std::vector<Model<FP>>& model, std::vect
     std::vector<std::vector<double>> 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<int>(std::round(model[node].parameters.template get<TimeExposed<FP>>()[(AgeGroup)group])));
@@ -121,26 +136,49 @@ IOResult<void> set_confirmed_cases_data(std::vector<Model<FP>>& 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<double>& 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<void> set_population_data(std::vector<Model<FP>>& 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<AgeGroup>(
+                    {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<AgeGroup>(
-                {i, InfectionState::Susceptible}, num_population[region][size_t(i)]);
+                {AgeGroup(0), InfectionState::Susceptible}, total);
         }
     }
     return success();
@@ -283,8 +331,8 @@ IOResult<void> 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<TimeSeries<double>> extrapolated_data(
         region.size(), TimeSeries<double>::zero(num_days + 1, (size_t)InfectionState::Count * num_age_groups));

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 <gtest/gtest.h>
 
@@ -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<mio::osecir::Model<double>> models6{mio::osecir::Model<double>((int)num_age_groups)};
+    std::vector<mio::osecir::Model<double>> models1{mio::osecir::Model<double>(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<mio::osecir::SeverePerInfectedSymptoms<double>>()[i] = 0.2;
+        models6[0].parameters.get<mio::osecir::CriticalPerSevere<double>>()[i]         = 0.25;
+    }
+
+    // Relevant parameters for model with 1 age group
+    models1[0].parameters.get<mio::osecir::SeverePerInfectedSymptoms<double>>()[mio::AgeGroup(0)] = 0.2;
+    models1[0].parameters.get<mio::osecir::CriticalPerSevere<double>>()[mio::AgeGroup(0)]         = 0.25;
+
+    const auto pydata_dir_Germany = mio::path_join(TEST_DATA_DIR, "Germany", "pydata");
+    const std::vector<int> counties{1002};
+    const auto date = mio::Date(2020, 12, 1);
+
+    std::vector<double> scale6(num_age_groups, 1.0);
+    std::vector<double> 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<mio::osecir::Model<double>> models6{mio::osecir::Model<double>((int)num_age_groups)};
+    std::vector<mio::osecir::Model<double>> models1{mio::osecir::Model<double>(1)};
+
+    // Test population data with 6 different values for age groups
+    std::vector<std::vector<double>> population_data6 = {{10000.0, 20000.0, 30000.0, 25000.0, 15000.0, 8000.0}};
+    std::vector<std::vector<double>> population_data1 = {{108000.0}}; // sum of all age groups
+    std::vector<int> 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<mio::osecir::Model<double>> models6{mio::osecir::Model<double>((int)num_age_groups)};
+    std::vector<mio::osecir::Model<double>> models1{mio::osecir::Model<double>(1)};
+
+    // Create case data for all 6 age groups over multiple days (current day + 6 days back)
+    std::vector<mio::ConfirmedCasesDataEntry> 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<int> regions            = {1002};
+    std::vector<double> 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<mio::osecir::Model<double>> models_6_groups{mio::osecir::Model<double>(6)};
+    std::vector<mio::osecir::Model<double>> models_1_group{mio::osecir::Model<double>(1)};
+
+    // Set relevant parameters for all age groups
+    for (int i = 0; i < 6; i++) {
+        models_6_groups[0].parameters.get<mio::osecir::SeverePerInfectedSymptoms<double>>()[mio::AgeGroup(i)] = 0.2;
+        models_6_groups[0].parameters.get<mio::osecir::CriticalPerSevere<double>>()[mio::AgeGroup(i)]         = 0.25;
+    }
+
+    // Set relevant parameters for 1 age group model
+    models_1_group[0].parameters.get<mio::osecir::SeverePerInfectedSymptoms<double>>()[mio::AgeGroup(0)] = 0.2;
+    models_1_group[0].parameters.get<mio::osecir::CriticalPerSevere<double>>()[mio::AgeGroup(0)]         = 0.25;
+
+    // Apply DIVI data to both models
+    std::vector<int> 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

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