1087 adapt parameters_io of ide model (#1095)

- Shifted the deaths data to get data according to the day of death instead of the day of the positive test.
- Corrected bug to allow correct initialization of flow InfectedNoSymptomsToRecovered.

Co-authored-by: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com>

Author: annawendler

cpp/models/ide_secir/parameters_io.cpp

@@ -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);
@@ -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.
@@ -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;
@@ -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.));
@@ -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++) {
@@ -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)]) *
@@ -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

cpp/models/ide_secir/parameters_io.h

@@ -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.
 * 

cpp/tests/data/ide-parameters-io-compare.csv

@@ -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

cpp/tests/test_ide_parameters_io.cpp

@@ -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;
@@ -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);
 
@@ -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);
@@ -125,4 +144,4 @@ TEST(TestIDEParametersIo, ParametersIoRKIFailure)
 
     // Reactive log output.
     mio::set_log_level(mio::LogLevel::warn);
-}
+}