Skip to content

1087 adapt parameters_io of ide model #1095

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 4 commits into from
Aug 9, 2024
Merged

1087 adapt parameters_io of ide model #1095

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
84104be
adapt parameters_io
annawendler Aug 7, 2024
08fceb5
ScalarType to int
annawendler Aug 7, 2024
72ef6c6
add test and refactor code
lenaploetzke Aug 8, 2024
444e0c6
modify test
annawendler Aug 9, 2024
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
78 changes: 66 additions & 12 deletions cpp/models/ide_secir/parameters_io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,10 @@ IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const&
log_error("Specified date does not exist in RKI data.");
return failure(StatusCode::OutOfRange, path + ", specified date does not exist in RKI data.");
}
auto min_date_entry = std::min_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
return a.date < b.date;
});
auto min_date = min_date_entry->date;

// Get (global) support_max to determine how many flows in the past we have to compute.
ScalarType global_support_max = model.get_global_support_max(dt);
Expand All @@ -71,17 +75,33 @@ IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const&
0, TimeSeries<ScalarType>::Vector::Constant((int)InfectionState::Count, 0));
}

// The first time we need is -4 * global_support_max.
Eigen::Index start_shift = 4 * global_support_max_index;
// The last time needed is dependent on the mean stay time in the Exposed compartment and
// the mean stay time of asymptomatic individuals in InfectedNoSymptoms.
// Set the Dead compartment to 0 so that RKI data can be added correctly.
model.m_populations[0][Eigen::Index(InfectionState::Dead)] = 0;

// Define variables for the mean of transitions used.
ScalarType mean_ExposedToInfectedNoSymptoms =
model.parameters.get<TransitionDistributions>()[Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)]
.get_mean(dt);
ScalarType mean_InfectedNoSymptomsToInfectedSymptoms =
model.parameters
.get<TransitionDistributions>()[Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]
.get<TransitionDistributions>()[(Eigen::Index)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]
.get_mean(dt);
ScalarType mean_InfectedSymptomsToInfectedSevere =
model.parameters
.get<TransitionDistributions>()[(Eigen::Index)InfectionTransition::InfectedSymptomsToInfectedSevere]
.get_mean();
ScalarType mean_InfectedSevereToInfectedCritical =
model.parameters
.get<TransitionDistributions>()[(Eigen::Index)InfectionTransition::InfectedSevereToInfectedCritical]
.get_mean();
ScalarType mean_InfectedCriticalToDead =
model.parameters.get<TransitionDistributions>()[(Eigen::Index)InfectionTransition::InfectedCriticalToDead]
.get_mean();

// The first time we need is -4 * global_support_max.
Eigen::Index start_shift = 4 * global_support_max_index;
// The last time needed is dependent on the mean stay time in the Exposed compartment and
// the mean stay time of asymptomatic individuals in InfectedNoSymptoms.
Eigen::Index last_time_index_needed =
Eigen::Index(std::ceil((mean_ExposedToInfectedNoSymptoms + mean_InfectedNoSymptomsToInfectedSymptoms) / dt));
// Create TimeSeries with zeros. The index of time zero is start_shift.
Expand All @@ -99,6 +119,7 @@ IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const&

bool min_offset_needed_avail = false;
bool max_offset_needed_avail = false;

// Go through the entries of rki_data and check if date is needed for calculation. Confirmed cases are scaled.
// Define dummy variables to store the first and the last index of the TimeSeries where the considered entry of rki_data is potentially needed.
Eigen::Index idx_needed_first = 0;
Expand All @@ -110,6 +131,9 @@ IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const&
if (offset == min_offset_needed) {
min_offset_needed_avail = true;
}
if (offset == max_offset_needed) {
max_offset_needed_avail = true;
}
// Smallest index for which the entry is needed.
idx_needed_first =
Eigen::Index(std::max(std::floor((offset - model.m_transitions.get_time(0) - 1) / dt), 0.));
Expand Down Expand Up @@ -137,21 +161,51 @@ IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const&
}
}

if (offset == max_offset_needed) {
max_offset_needed_avail = true;
// Compute Dead by shifting RKI data according to mean stay times.
// This is done because the RKI reports death with the date of positive test instead of the date of deaths.
if (offset == std::floor(-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
mean_InfectedCriticalToDead)) {
model.m_populations[0][Eigen::Index(InfectionState::Dead)] +=
(1 - (-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
mean_InfectedCriticalToDead -
std::floor(-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
mean_InfectedCriticalToDead))) *
entry.num_deaths;
}
if (offset == std::ceil(-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
mean_InfectedCriticalToDead)) {
model.m_populations[0][Eigen::Index(InfectionState::Dead)] +=
(-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
mean_InfectedCriticalToDead -
std::floor(-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
mean_InfectedCriticalToDead)) *
entry.num_deaths;
}

if (offset == 0) {
model.m_populations[0][Eigen::Index(InfectionState::Dead)] = entry.num_deaths;
model.m_total_confirmed_cases = scale_confirmed_cases * entry.num_confirmed;
}
}
}

if (!max_offset_needed_avail || !min_offset_needed_avail) {
if (!max_offset_needed_avail) {
log_error("Necessary range of dates needed to compute initial values does not exist in RKI data.");
return failure(StatusCode::OutOfRange, path + ", necessary range of dates does not exist in RKI data.");
}

if (!min_offset_needed_avail) {
std::string min_date_string =
std::to_string(min_date.day) + "." + std::to_string(min_date.month) + "." + std::to_string(min_date.year);
// Get first date that is needed.
mio::Date min_offset_date = offset_date_by_days(date, int(min_offset_needed));
std::string min_offset_date_string = std::to_string(min_offset_date.day) + "." +
std::to_string(min_offset_date.month) + "." +
std::to_string(min_offset_date.year);
log_warning("RKI data is needed from " + min_offset_date_string +
" to compute initial values. RKI data is only available from " + min_date_string +
". Missing dates were set to 0.");
}

//--- Calculate the flows "after" InfectedNoSymptomsToInfectedSymptoms. ---
// Compute flow InfectedSymptomsToInfectedSevere for -3 * global_support_max, ..., 0.
for (Eigen::Index i = -3 * global_support_max_index; i <= 0; i++) {
Expand Down Expand Up @@ -183,10 +237,10 @@ IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const&
}

//--- Calculate the remaining flows. ---
// Compute flow ExposedToInfectedNoSymptoms for -global_support_max, ..., 0.
// Compute flow ExposedToInfectedNoSymptoms for -2 * global_support_max, ..., 0.
// Use mean value of the TransitionDistribution InfectedNoSymptomsToInfectedSymptoms for the calculation.
Eigen::Index index_shift_mean = Eigen::Index(std::round(mean_InfectedNoSymptomsToInfectedSymptoms / dt));
for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) {
model.m_transitions[i + start_shift][Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)] =
(1 / model.parameters.get<TransitionProbabilities>()[Eigen::Index(
InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) *
Expand Down Expand Up @@ -226,4 +280,4 @@ IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const&
} // namespace isecir
} // namespace mio

#endif // MEMILIO_HAS_JSONCPP
#endif // MEMILIO_HAS_JSONCPP
3 changes: 2 additions & 1 deletion cpp/models/ide_secir/parameters_io.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,8 @@ namespace isecir
* The flow InfectedNoSymptomsToInfectedSymptoms is calculated with the standard formula from the IDE model
* using the results for ExposedToInfectedNoSymptoms.
*
* The number of deaths used in the model is set to the number given in the RKI data.
* The number of deaths used in the model is computed by shifting the reported RKI data according to the mean values of the transitions
* InfectedSymptomsToInfectedSevere, InfectedSevereToInfectedCritical and InfectedCriticalToDead.
* We also set the number of total confirmed cases in the model.
* Therefore the initialization method using the total confirmed cases is used in the model. See also the documentation of the model class.
*
Expand Down
6 changes: 3 additions & 3 deletions cpp/tests/data/ide-parameters-io-compare.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Time S->E E->C C->I C->R I->H I->R H->U H->R U->D U->R
-2.00000000 25.50000000 0.50000000 4.50000000 0.03661165 1.83235622 1.67500000 0.50414276 0.50414276 0.13557638 0.13557638
-1.50000000 30.00000000 25.50000000 0.25000000 1.95558262 1.78897083 1.77038735 0.73448590 0.73448590 0.20627078 0.20627078
-1.00000000 30.00000000 25.50000000 0.25000000 6.46338835 0.89673307 1.18750000 0.80038452 0.80038452 0.29817594 0.29817594
-2.00000000 25.50000000 0.50000000 4.50000000 2.37500000 1.83235622 1.67500000 0.50414276 0.50414276 0.13557638 0.13557638
-1.50000000 30.00000000 25.50000000 0.25000000 2.70298071 1.78897083 1.77038735 0.73448590 0.73448590 0.20627078 0.20627078
-1.00000000 30.00000000 25.50000000 0.25000000 6.50000000 0.89673307 1.18750000 0.80038452 0.80038452 0.29817594 0.29817594
-0.50000000 0.30000000 30.00000000 12.75000000 11.24892225 1.43851502 1.35149035 0.71427386 0.71427386 0.36054581 0.36054581
0.00000000 0.30000000 30.00000000 12.75000000 13.87500000 4.10519684 3.25000000 0.84440788 0.84440788 0.38336812 0.38336812
27 changes: 23 additions & 4 deletions cpp/tests/test_ide_parameters_io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,9 @@ TEST(TestIDEParametersIo, RKIcompareWithPreviousRun)
// Compare with previous run.
auto compare = load_test_data_csv<ScalarType>("ide-parameters-io-compare.csv");

EXPECT_NEAR(model.m_populations.get_value(0)[(Eigen::Index)mio::isecir::InfectionState::Dead], 6.3, 1e-4);
EXPECT_NEAR(model.m_total_confirmed_cases, 70, 1e-4);

ASSERT_EQ(compare.size(), static_cast<size_t>(model.m_transitions.get_num_time_points()));
for (size_t i = 0; i < compare.size(); i++) {
ASSERT_EQ(compare[i].size(), static_cast<size_t>(model.m_transitions.get_num_elements()) + 1) << "at row " << i;
Expand All @@ -91,12 +94,13 @@ TEST(TestIDEParametersIo, RKIcompareWithPreviousRun)
// Check some cases where computation of initial values for an IDE model based on RKI data should fail.
TEST(TestIDEParametersIo, ParametersIoRKIFailure)
{
// Define start date and the total population used for the initialization.
// Set up.
ScalarType total_population = 80 * 1e6;
ScalarType dt = 0.5;

// Initialize model.
// The number of deaths will be overwritten if real data is used for initialization. Therefore, an arbitrary number is used for the number of deaths.
// The number of deaths will be overwritten if real data is used for initialization.
// Therefore, an arbitrary number is used for the number of deaths.
mio::isecir::Model model(mio::TimeSeries<ScalarType>((int)mio::isecir::InfectionTransition::Count),
total_population, -1);

Expand All @@ -110,13 +114,28 @@ TEST(TestIDEParametersIo, ParametersIoRKIFailure)

ASSERT_THAT(print_wrap(status), IsFailure(mio::StatusCode::OutOfRange));

// --- Case where not all needed dates are provided.
// --- Case where not all needed dates from the future are provided.
start_date = mio::Date(2020, 6, 7);
status =
mio::isecir::set_initial_flows(model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date);

ASSERT_THAT(print_wrap(status), IsFailure(mio::StatusCode::OutOfRange));

// --- Case where not all needed dates from the past are provided.
start_date = mio::Date(2020, 5, 24);
status =
mio::isecir::set_initial_flows(model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date);
// Check that status is Success as just a warning is logged.
ASSERT_THAT(print_wrap(status), IsSuccess());
// Check that the flow InfectedNoSymptomsToInfectedSymptoms has actually been set to 0.
// As start_date is the first date in the data file, there is some infection data for all times >-1
// (as we interpolate linearly in between -1 and 0). Therefore, we only expect the flow InfectedNoSymptomsToInfectedSymptoms
// to be zero for times <=-1.
for (Eigen::Index i = 0; i < model.m_transitions.get_num_time_points() - 2; i++) {
EXPECT_EQ(0., model.m_transitions.get_value(
i)[(Eigen::Index)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]);
}

// --- Case with empty RKI data file.
status =
mio::isecir::set_initial_flows(model, dt, mio::path_join(TEST_DATA_DIR, "test_empty_file.json"), start_date);
Expand All @@ -125,4 +144,4 @@ TEST(TestIDEParametersIo, ParametersIoRKIFailure)

// Reactive log output.
mio::set_log_level(mio::LogLevel::warn);
}
}