Skip to content

1348 Fix initialization recovered ode secir #1353

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Aug 20, 2025
Merged

1348 Fix initialization recovered ode secir #1353

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
99b0ccf
fix recovered ?
HenrZu Aug 15, 2025
5196125
typo ....
HenrZu Aug 15, 2025
cd39863
update tests
HenrZu Aug 15, 2025
e733f84
unscale to init rec
HenrZu Aug 19, 2025
37949ed
[ci skip] comment on scaling/ not scaling
HenrZu Aug 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 14 additions & 0 deletions cpp/models/ode_secir/parameters_io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,9 @@ IOResult<void> read_confirmed_cases_data(

if (date_df == offset_date_by_days(date, 0)) {
num_InfectedSymptoms[age] += scaling_factor_inf[age] * region_entry.num_confirmed;
// We intentionally do NOT multiply recovered with the scaling_factor_inf here.
// If we apply the scaling factor to recovered as well, we would implicitly
// assume that this factor holds for the entire historical period up to t0, which is not valid.
num_rec[age] += region_entry.num_confirmed;
}
if (date_df == offset_date_by_days(date, days_surplus)) {
Expand Down Expand Up @@ -159,6 +162,17 @@ IOResult<void> read_confirmed_cases_data(
auto& num_icu = vnum_icu[region_idx];

for (size_t i = 0; i < ConfirmedCasesDataEntry::age_group_names.size(); i++) {
// R(t0) = ΣC(t0) − I(t0) − H(t0) − U(t0) − D(t0)
// subtract currently infectious/hospitalized/ICU/dead
// Note: We divide I/H/U/D by scaling_factor_inf to "unscale" these contributions back to the
// reported level before subtracting from recovered. If we also applied the scaling factor to
// recovered, we would implicitly assume that the same underreporting applies to the entire
// history up to t0, which would be wrong. The scaling factor should reflect underreporting
// around t0 only.
num_rec[i] -=
(num_InfectedSymptoms[i] / scaling_factor_inf[i] + num_InfectedSevere[i] / scaling_factor_inf[i] +
num_icu[i] / scaling_factor_inf[i] + num_death[i] / scaling_factor_inf[i]);

auto try_fix_constraints = [region, i](double& value, double error, auto str) {
if (value < error) {
//this should probably return a failure
Expand Down
Binary file modified cpp/tests/data/export_time_series_initialization_osecir.h5
View file
Open in desktop
Binary file not shown.
20 changes: 10 additions & 10 deletions cpp/tests/test_odesecir.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1282,11 +1282,11 @@ TEST_F(ModelTestOdeSecir, read_input_data)

// 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)
(Eigen::ArrayXd(num_age_groups * Eigen::Index(mio::osecir::InfectionState::Count)) << 10287.4, 2.82575, 1.25589,
0, 5.85714, 0, 1.42857, 1.33333, 20.9357, 0, 19105.3, 6.59341, 4.23862, 0, 11, 0, 2.54286, 1.33333, 73.9571, 0,
73901.1, 41.9152, 21.6641, 0, 54.5714, 0, 22, 1.33333, 444.393, 0, 82602.5, 39.5604, 21.0361, 0, 51.4286, 0,
15.1714, 1.33333, 399.686, 0, 43757.4, 8.32025, 5.18053, 0, 11.4286, 0, 4.22857, 1.33333, 134.264, 8, 15654.2,
10.5181, 3.2967, 0, 3.28571, 0, 0.657143, 1.33333, 38.3, 7.57143)
.finished();

EXPECT_THAT(print_wrap(model1[0].populations.array().cast<double>()),
Expand Down Expand Up @@ -1379,11 +1379,11 @@ TEST_F(ModelTestOdeSecir, model_initialization)
// 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)
(Eigen::ArrayXd(num_age_groups * Eigen::Index(mio::osecir::InfectionState::Count)) << 10287.4, 2.82575, 1.25589,
0, 5.85714, 0, 1.42857, 1.33333, 20.9357, 0, 19105.3, 6.59341, 4.23862, 0, 11, 0, 2.54286, 1.33333, 73.9571, 0,
73901.1, 41.9152, 21.6641, 0, 54.5714, 0, 22, 1.33333, 444.393, 0, 82602.5, 39.5604, 21.0361, 0, 51.4286, 0,
15.1714, 1.33333, 399.686, 0, 43757.4, 8.32025, 5.18053, 0, 11.4286, 0, 4.22857, 1.33333, 134.264, 8, 15654.2,
10.5181, 3.2967, 0, 3.28571, 0, 0.657143, 1.33333, 38.3, 7.57143)
.finished();

EXPECT_THAT(print_wrap(model_vector[0].populations.array().cast<double>()),
Expand Down
16 changes: 8 additions & 8 deletions cpp/tests/test_save_parameters.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -644,13 +644,13 @@ TEST(TestSaveParameters, ReadPopulationDataRKIAges)
auto read_result = mio::osecir::read_input_data_germany(model, date, scaling_factor_inf, scaling_factor_icu, path);
ASSERT_THAT(print_wrap(read_result), IsSuccess());

std::vector<double> sus = {3443857.42, 7665093.95, 18792870.93, 29503629.76, 16307262.45, 6049150.54};
std::vector<double> sus = {3444279.56, 7666866.07, 18801541.51, 29518513.27, 16321649.37, 6072737.25};
std::vector<double> exp = {433.015, 1771.61, 8856.33, 14757.62, 7222.86, 6626.07};
std::vector<double> car = {434.444, 1772.14, 8724.49, 14386.90, 6995.14, 6307.14};
std::vector<double> inf = {375.429, 1393.43, 6007.14, 8438.71, 3377.57, 2421.57};
std::vector<double> hosp = {39.9614, 303.191, 1934.84, 3621.2, 1793.39, 1557.03};
std::vector<double> icu = {47.6813, 190.725, 429.132, 762.901, 1192.03, 1716.53};
std::vector<double> rec = {23557.7, 78946.3, 398585.142, 487273.71, 178660.14, 96021.9};
std::vector<double> rec = {23135.58, 77174.17, 389914.56, 472390.20, 164273.23, 72435.15};
std::vector<double> death = {2, 4, 48, 1137.86, 8174.14, 18528.9};

for (size_t i = 0; i < 6; i++) {
Expand Down Expand Up @@ -694,13 +694,13 @@ TEST(TestSaveParameters, ReadPopulationDataStateAllAges)
mio::osecir::read_input_data_state(model, date, state, scaling_factor_inf, scaling_factor_icu, path);
ASSERT_THAT(print_wrap(read_result), IsSuccess());

std::vector<double> sus = {116692.2, 283912.8, 622795.86, 1042178.3, 606450.7, 212836.9};
std::vector<double> sus = {116699.98, 283936.40, 622935.07, 1042390.64, 606622.90, 213104.13};
std::vector<double> exp = {8.57143, 30.5357, 149.388, 228.809, 87.1429, 99.2857};
std::vector<double> car = {7.77778, 26.0714, 143.061, 217.143, 84.8571, 92.1429};
std::vector<double> inf = {7.00000, 18.7143, 97.7143, 122.000, 40.8571, 36.1429};
std::vector<double> hosp = {0.707143, 3.92857, 30.6429, 50.5371, 20.35, 19.9886};
std::vector<double> icu = {0.274725, 1.0989, 2.47253, 4.3956, 6.86813, 9.89011};
std::vector<double> rec = {393.143, 1216.14, 5467.86, 6543.57, 2281.29, 1045.71};
std::vector<double> rec = {385.36, 1192.59, 5328.66, 6331.18, 2109.05, 778.53};
std::vector<double> death = {0, 0, 0, 16.2857, 99.5714, 198.286};

for (size_t i = 0; i < 6; i++) {
Expand Down Expand Up @@ -761,13 +761,13 @@ TEST(TestSaveParameters, ReadPopulationDataCountyAllAges)
ASSERT_THAT(print_wrap(read_result2), IsSuccess());
ASSERT_THAT(print_wrap(read_result_district), IsSuccess());

std::vector<double> sus = {10284.13, 19082.86, 73783.12, 82494.81, 43725.08, 15612.70};
std::vector<double> sus = {10284.66, 19086.96, 73805.00, 82515.81, 43739.68, 15633.34};
std::vector<double> exp = {0.571429, 4.82143, 20.8163, 22.1429, 4.57143, 4.64286};
std::vector<double> car = {0.557143, 4.46429, 22.0408, 20.7143, 4.28571, 4.64286};
std::vector<double> inf = {0.42857, 3.285714, 15.2857, 13.0000, 2.42857, 2.00000};
std::vector<double> hosp = {0.0942857, 0.691429, 4.90286, 5.34286, 1.41429, 2.45143};
std::vector<double> icu = {0.0769231, 0.307692, 0.692308, 1.23077, 1.92308, 2.76923};
std::vector<double> rec = {35, 108.571, 640.143, 573.429, 180.429, 75.5714};
std::vector<double> rec = {34.47, 104.47, 618.26, 552.43, 165.83, 54.93};
std::vector<double> death = {0, 0, 0, 0, 10, 14.4286};

for (size_t i = 0; i < 6; i++) {
Expand Down Expand Up @@ -850,13 +850,13 @@ TEST(TestSaveParameters, ExtrapolateRKI)
auto& file_results = read_result.value();
auto results = file_results[0].get_groups();

std::vector<double> sus = {10284.1, 19082.9, 73783.1, 82494.8, 43725.1, 15612.7};
std::vector<double> sus = {10284.66, 19086.96, 73805.00, 82515.81, 43739.68, 15633.34};
std::vector<double> exp = {0.571429, 4.82143, 20.8163, 22.1429, 4.57143, 4.64286};
std::vector<double> car = {0.557143, 4.46429, 22.0408, 20.7143, 4.28571, 4.64286};
std::vector<double> inf = {0.428571, 3.28571, 15.2857, 13.0000, 2.42857, 2.00000};
std::vector<double> hosp = {0.0942857, 0.691429, 4.90286, 5.34286, 1.41429, 2.45143};
std::vector<double> icu = {0.0769231, 0.307692, 0.692308, 1.23077, 1.92308, 2.76923};
std::vector<double> rec = {35, 108.571, 640.143, 573.429, 180.429, 75.5714};
std::vector<double> rec = {34.47, 104.47, 618.26, 552.43, 165.83, 54.93};
std::vector<double> death = {0, 0, 0, 0, 10, 14.4286};

for (size_t i = 0; i < 6; i++) {
Expand Down