diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index bee7f41eaa..cd0a36e697 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -174,6 +174,7 @@ endif() if(MEMILIO_BUILD_SIMULATIONS) add_subdirectory(simulations) + add_subdirectory(simulations/IDE_paper) endif() if(MEMILIO_BUILD_BENCHMARKS) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 1d13454b6a..bc4e5ed826 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -164,7 +164,7 @@ if(MEMILIO_HAS_HDF5) endif() if(MEMILIO_HAS_JSONCPP) - add_executable(ide_initialization_example ide_initialization.cpp) - target_link_libraries(ide_initialization_example PRIVATE memilio ide_secir) - target_compile_options(ide_initialization_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(ide_initialization_example ide_initialization.cpp) + target_link_libraries(ide_initialization_example PRIVATE memilio ide_secir) + target_compile_options(ide_initialization_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) endif() diff --git a/cpp/models/ide_secir/model.h b/cpp/models/ide_secir/model.h index 9e02bdd50c..9a04a1a2f6 100644 --- a/cpp/models/ide_secir/model.h +++ b/cpp/models/ide_secir/model.h @@ -25,6 +25,7 @@ #include "ide_secir/infection_state.h" #include "memilio/config.h" #include "memilio/epidemiology/age_group.h" +#include "memilio/io/epi_data.h" #include "memilio/utils/custom_index_array.h" #include "memilio/utils/time_series.h" @@ -358,6 +359,7 @@ class Model friend class Simulation; // In set_initial_flows(), we compute initial flows based on RKI data using the (private) compute_flow() function // which is why it is defined as a friend function. + template friend IOResult set_initial_flows(Model& model, ScalarType dt, std::string const& path, Date date, ScalarType scale_confirmed_cases); }; diff --git a/cpp/models/ide_secir/parameters_io.cpp b/cpp/models/ide_secir/parameters_io.cpp deleted file mode 100644 index 39fcd07dc3..0000000000 --- a/cpp/models/ide_secir/parameters_io.cpp +++ /dev/null @@ -1,343 +0,0 @@ -/* -* Copyright (C) 2020-2024 MEmilio -* -* Authors: Lena Ploetzke, Anna Wendler -* -* Contact: Martin J. Kuehn -* -* Licensed under the Apache License, Version 2.0 (the "License"); -* you may not use this file except in compliance with the License. -* You may obtain a copy of the License at -* -* http://www.apache.org/licenses/LICENSE-2.0 -* -* Unless required by applicable law or agreed to in writing, software -* distributed under the License is distributed on an "AS IS" BASIS, -* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -* See the License for the specific language governing permissions and -* limitations under the License. -*/ - -#include "ide_secir/parameters_io.h" -#include "memilio/config.h" -#include "memilio/epidemiology/age_group.h" -#include -#include - -#ifdef MEMILIO_HAS_JSONCPP - -#include "ide_secir/model.h" -#include "ide_secir/infection_state.h" -#include "memilio/math/eigen.h" -#include "memilio/io/epi_data.h" -#include "memilio/io/io.h" -#include "memilio/utils/date.h" - -#include -#include - -namespace mio -{ -namespace isecir -{ - -IOResult set_initial_flows(Model& model, ScalarType dt, std::string const& path, Date date, - ScalarType scale_confirmed_cases) -{ - //--- Preparations --- - // Try to get RKI data from path. - BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path)); - auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { - return a.date < b.date; - }); - if (max_date_entry == rki_data.end()) { - log_error("RKI data file is empty."); - return failure(StatusCode::InvalidFileFormat, path + ", file is empty."); - } - auto max_date = max_date_entry->date; - if (max_date < date) { - 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); - Eigen::Index global_support_max_index = Eigen::Index(std::ceil(global_support_max / dt)); - - // Get the number of AgeGroups. - const size_t num_age_groups = ConfirmedCasesDataEntry::age_group_names.size(); - - // m_transitions should be empty at the beginning. - - if (model.m_transitions.get_num_time_points() > 0) { - model.m_transitions = TimeSeries(Eigen::Index(InfectionTransition::Count) * num_age_groups); - } - if (model.m_populations.get_time(0) != 0) { - model.m_populations.remove_last_time_point(); - model.m_populations.add_time_point( - 0, TimeSeries::Vector::Constant((int)InfectionState::Count * num_age_groups, 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. - // The mean stay time in a compartment may be dependent on the AgeGroup. - CustomIndexArray mean_ExposedToInfectedNoSymptoms = - CustomIndexArray(AgeGroup(num_age_groups), 0.); - CustomIndexArray mean_InfectedNoSymptomsToInfectedSymptoms = - CustomIndexArray(AgeGroup(num_age_groups), 0.); - CustomIndexArray mean_InfectedSymptomsToInfectedSevere = - CustomIndexArray(AgeGroup(num_age_groups), 0.); - CustomIndexArray mean_InfectedSevereToInfectedCritical = - CustomIndexArray(AgeGroup(num_age_groups), 0.); - CustomIndexArray mean_InfectedCriticalToDead = - CustomIndexArray(AgeGroup(num_age_groups), 0.); - Eigen::Index last_time_index_needed = 0; - - for (AgeGroup group = AgeGroup(0); group < AgeGroup(num_age_groups); group++) { - // Set the Dead compartment to 0 so that RKI data can be added correctly. - int Di = model.get_state_flat_index(Eigen::Index(InfectionState::Dead), group); - model.m_populations[0][Di] = 0; - - mean_ExposedToInfectedNoSymptoms[group] = - model.parameters - .get()[group][Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)] - .get_mean(dt); - mean_InfectedNoSymptomsToInfectedSymptoms[group] = - model.parameters - .get()[group] - [Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] - .get_mean(dt); - mean_InfectedSymptomsToInfectedSevere[group] = - model.parameters - .get()[group] - [Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere)] - .get_mean(dt); - mean_InfectedSevereToInfectedCritical[group] = - model.parameters - .get()[group] - [Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical)] - .get_mean(dt); - mean_InfectedCriticalToDead[group] = - model.parameters - .get()[group][Eigen::Index(InfectionTransition::InfectedCriticalToDead)] - .get_mean(dt); - if (last_time_index_needed < - Eigen::Index(std::ceil( - (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt))) { - last_time_index_needed = Eigen::Index(std::ceil( - (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt)); - } - } - - // Create TimeSeries with zeros. The index of time zero is start_shift. - for (Eigen::Index i = -start_shift; i <= last_time_index_needed; i++) { - // Add time point. - model.m_transitions.add_time_point( - i * dt, TimeSeries::Vector::Constant((size_t)InfectionTransition::Count * num_age_groups, 0.)); - } - - // Setting the entries in m_total_confirmed_cases to zero before overwriting it with the RKI data. - model.m_total_confirmed_cases = CustomIndexArray(AgeGroup(num_age_groups), 0.); - //--- Calculate the flow InfectedNoSymptomsToInfectedSymptoms using the RKI data and store in the m_transitions object.--- - ScalarType min_offset_needed = std::ceil( - model.m_transitions.get_time(0) - - 1); // Need -1 if first time point is integer and just the floor value if not, therefore use ceil and -1 - ScalarType max_offset_needed = std::ceil(model.m_transitions.get_last_time()); - 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; - Eigen::Index idx_needed_last = 0; - ScalarType time_idx = 0; - for (auto&& entry : rki_data) { - int offset = get_offset_in_days(entry.date, date); - AgeGroup group = entry.age_group; - if ((offset >= min_offset_needed) && (offset <= max_offset_needed)) { - 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.)); - // Biggest index for which the entry is needed. - idx_needed_last = Eigen::Index(std::min(std::ceil((offset - model.m_transitions.get_time(0) + 1) / dt), - double(model.m_transitions.get_num_time_points() - 1))); - - int INStISyi = model.get_transition_flat_index( - Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), group); - - for (Eigen::Index i = idx_needed_first; i <= idx_needed_last; i++) { - - time_idx = model.m_transitions.get_time(i); - if (offset == int(std::floor(time_idx))) { - model.m_transitions[i][INStISyi] += - (1 - (time_idx - std::floor(time_idx))) * scale_confirmed_cases * entry.num_confirmed; - } - if (offset == int(std::ceil(time_idx))) { - model.m_transitions[i][INStISyi] += - (time_idx - std::floor(time_idx)) * scale_confirmed_cases * entry.num_confirmed; - } - if (offset == int(std::floor(time_idx - dt))) { - model.m_transitions[i][INStISyi] -= - (1 - (time_idx - dt - std::floor(time_idx - dt))) * scale_confirmed_cases * entry.num_confirmed; - } - if (offset == int(std::ceil(time_idx - dt))) { - model.m_transitions[i][INStISyi] -= - (time_idx - dt - std::floor(time_idx - dt)) * scale_confirmed_cases * entry.num_confirmed; - } - } - - // 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. - int Di = model.get_state_flat_index(Eigen::Index(InfectionState::Dead), group); - if (offset == - std::floor(-mean_InfectedSymptomsToInfectedSevere[group] - - mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group])) { - model.m_populations[0][Di] += - (1 - - (-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] - - mean_InfectedCriticalToDead[group] - - std::floor(-mean_InfectedSymptomsToInfectedSevere[group] - - mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group]))) * - entry.num_deaths; - } - if (offset == - std::ceil(-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] - - mean_InfectedCriticalToDead[group])) { - model.m_populations[0][Di] += - (-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] - - mean_InfectedCriticalToDead[group] - - std::floor(-mean_InfectedSymptomsToInfectedSevere[group] - - mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group])) * - entry.num_deaths; - } - - if (offset == 0) { - model.m_total_confirmed_cases[group] = scale_confirmed_cases * entry.num_confirmed; - } - } - } - - 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. --- - // Set m_support_max_vector and m_derivative_vector in the model which is needed for the following computations. - model.set_transitiondistributions_support_max(dt); - model.set_transitiondistributions_derivative(dt); - - for (AgeGroup group = AgeGroup(0); group < AgeGroup(num_age_groups); group++) { - //--- 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++) { - model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), - Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt, - i + start_shift, group); - } - // Compute flow InfectedSevereToInfectedCritical for -2 * global_support_max, ..., 0. - for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) { - model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), - Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift, - group); - } - // Compute flows from InfectedSymptoms, InfectedSevere and InfectedCritical to Recovered and - // flow InfectedCriticalToDead for -global_support_max, ..., 0. - for (Eigen::Index i = -global_support_max_index; i <= 0; i++) { - // Compute flow InfectedSymptomsToRecovered. - model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered), - Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt, - i + start_shift, group); - // Compute flow InfectedSevereToRecovered. - model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToRecovered), - Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift, - group); - // Compute flow InfectedCriticalToRecovered. - model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered), - Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift, - group); - // Compute flow InfectedCriticalToDead. - model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToDead), - Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift, - group); - } - - //--- Calculate the remaining flows. --- - // 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[group] / dt)); - int EtINSi = - model.get_transition_flat_index(Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), group); - int INStISyi = model.get_transition_flat_index( - Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), group); - for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) { - model.m_transitions[i + start_shift][EtINSi] = - (1 / model.parameters.get()[group][Eigen::Index( - InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) * - model.m_transitions[i + start_shift + index_shift_mean][INStISyi]; - } - - // Compute flow SusceptibleToExposed for -global_support_max, ..., 0. - // Use mean values of the TransitionDistribution ExposedToInfectedNoSymptoms and of the - // TransitionDistribution InfectedNoSymptomsToInfectedSymptoms for the calculation. - index_shift_mean = Eigen::Index(std::round( - (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt)); - int StEi = model.get_transition_flat_index(Eigen::Index(InfectionTransition::SusceptibleToExposed), group); - - for (Eigen::Index i = -global_support_max_index; i <= 0; i++) { - model.m_transitions[i + start_shift][StEi] = - (1 / model.parameters.get()[group][Eigen::Index( - InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) * - model.m_transitions[i + start_shift + index_shift_mean][INStISyi]; - } - - // InfectedNoSymptomsToRecovered for -global_support_max, ..., 0. - // If we previously calculated the transition ExposedToInfectedNoSymptoms, we can calculate this transition - // using the standard formula. - for (Eigen::Index i = -global_support_max_index; i <= 0; i++) { - model.compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered), - Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt, i + start_shift, - group); - } - } - - // At the end of the calculation, delete all time points that are not required for the simulation. - auto transition_copy(model.m_transitions); - model.m_transitions = TimeSeries(Eigen::Index(InfectionTransition::Count) * num_age_groups); - for (Eigen::Index i = -global_support_max_index; i <= 0; i++) { - model.m_transitions.add_time_point(i * dt, transition_copy.get_value(i + start_shift)); - } - - return mio::success(); -} - -} // namespace isecir -} // namespace mio - -#endif // MEMILIO_HAS_JSONCPP diff --git a/cpp/models/ide_secir/parameters_io.h b/cpp/models/ide_secir/parameters_io.h index c909eba2f1..357cbd7adf 100644 --- a/cpp/models/ide_secir/parameters_io.h +++ b/cpp/models/ide_secir/parameters_io.h @@ -21,6 +21,9 @@ #define IDE_INITIALFLOWS_H #include "memilio/config.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/io/epi_data.h" +#include #ifdef MEMILIO_HAS_JSONCPP @@ -34,6 +37,44 @@ namespace mio { namespace isecir { +template ::value, int>* = nullptr> +IOResult> f(std::string const& path) +{ + BOOST_OUTCOME_TRY(auto&& rki_data_read, mio::read_confirmed_cases_noage(path)); + return rki_data_read; +} + +template ::value, int>* = nullptr> +IOResult> f(std::string const& path) +{ + BOOST_OUTCOME_TRY(auto&& rki_data_read, mio::read_confirmed_cases_data(path)); + return rki_data_read; +} + +template ::value, int>* = nullptr> +size_t g() +{ + return 1; +} + +template ::value, int>* = nullptr> +size_t g() +{ + return T::age_group_names.size(); +} + +template ::value, int>* = nullptr> +AgeGroup h(T entry) +{ + unused(entry); + return AgeGroup(0); +} + +template ::value, int>* = nullptr> +AgeGroup h(T entry) +{ + return entry.age_group; +} /** * @brief Computes a TimeSeries of flows to provide initial data for an IDE-SECIR model with data from RKI. @@ -66,12 +107,312 @@ namespace isecir * @param[in] scale_confirmed_cases Factor by which to scale the confirmed cases of rki data to consider unreported cases. * @returns Any io errors that happen during reading of the files. */ + +template IOResult set_initial_flows(Model& model, ScalarType dt, std::string const& path, Date date, - ScalarType scale_confirmed_cases = 1.); + ScalarType scale_confirmed_cases) +{ + + //--- Preparations --- + // Try to get RKI data from path. + std::vector rki_data = f(path).value(); + + auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { + return a.date < b.date; + }); + if (max_date_entry == rki_data.end()) { + log_error("RKI data file is empty."); + return failure(StatusCode::InvalidFileFormat, path + ", file is empty."); + } + auto max_date = max_date_entry->date; + if (max_date < date) { + 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); + Eigen::Index global_support_max_index = Eigen::Index(std::ceil(global_support_max / dt)); + + // Get the number of AgeGroups. + const size_t num_age_groups = g(); + std::cout << "num age groups: " << num_age_groups << std::endl; + + // m_transitions should be empty at the beginning. + + if (model.m_transitions.get_num_time_points() > 0) { + model.m_transitions = TimeSeries(Eigen::Index(InfectionTransition::Count) * num_age_groups); + } + if (model.m_populations.get_time(0) != 0) { + model.m_populations.remove_last_time_point(); + model.m_populations.add_time_point( + 0, TimeSeries::Vector::Constant((int)InfectionState::Count * num_age_groups, 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. + // The mean stay time in a compartment may be dependent on the AgeGroup. + CustomIndexArray mean_ExposedToInfectedNoSymptoms = + CustomIndexArray(AgeGroup(num_age_groups), 0.); + CustomIndexArray mean_InfectedNoSymptomsToInfectedSymptoms = + CustomIndexArray(AgeGroup(num_age_groups), 0.); + CustomIndexArray mean_InfectedSymptomsToInfectedSevere = + CustomIndexArray(AgeGroup(num_age_groups), 0.); + CustomIndexArray mean_InfectedSevereToInfectedCritical = + CustomIndexArray(AgeGroup(num_age_groups), 0.); + CustomIndexArray mean_InfectedCriticalToDead = + CustomIndexArray(AgeGroup(num_age_groups), 0.); + Eigen::Index last_time_index_needed = 0; + + for (AgeGroup group = AgeGroup(0); group < AgeGroup(num_age_groups); group++) { + // Set the Dead compartment to 0 so that RKI data can be added correctly. + int Di = model.get_state_flat_index(Eigen::Index(InfectionState::Dead), group); + model.m_populations[0][Di] = 0; + + mean_ExposedToInfectedNoSymptoms[group] = + model.parameters + .get()[group][Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)] + .get_mean(dt); + mean_InfectedNoSymptomsToInfectedSymptoms[group] = + model.parameters + .get()[group] + [Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] + .get_mean(dt); + mean_InfectedSymptomsToInfectedSevere[group] = + model.parameters + .get()[group] + [Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere)] + .get_mean(dt); + mean_InfectedSevereToInfectedCritical[group] = + model.parameters + .get()[group] + [Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical)] + .get_mean(dt); + mean_InfectedCriticalToDead[group] = + model.parameters + .get()[group][Eigen::Index(InfectionTransition::InfectedCriticalToDead)] + .get_mean(dt); + if (last_time_index_needed < + Eigen::Index(std::ceil( + (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt))) { + last_time_index_needed = Eigen::Index(std::ceil( + (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt)); + } + } + + // Create TimeSeries with zeros. The index of time zero is start_shift. + for (Eigen::Index i = -start_shift; i <= last_time_index_needed; i++) { + // Add time point. + model.m_transitions.add_time_point( + i * dt, TimeSeries::Vector::Constant((size_t)InfectionTransition::Count * num_age_groups, 0.)); + } + + // Setting the entries in m_total_confirmed_cases to zero before overwriting it with the RKI data. + model.m_total_confirmed_cases = CustomIndexArray(AgeGroup(num_age_groups), 0.); + //--- Calculate the flow InfectedNoSymptomsToInfectedSymptoms using the RKI data and store in the m_transitions object.--- + ScalarType min_offset_needed = std::ceil( + model.m_transitions.get_time(0) - + 1); // Need -1 if first time point is integer and just the floor value if not, therefore use ceil and -1 + ScalarType max_offset_needed = std::ceil(model.m_transitions.get_last_time()); + 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; + Eigen::Index idx_needed_last = 0; + ScalarType time_idx = 0; + + for (auto&& entry : rki_data) { + int offset = get_offset_in_days(entry.date, date); + AgeGroup group = h(entry); + + if ((offset >= min_offset_needed) && (offset <= max_offset_needed)) { + 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.)); + // Biggest index for which the entry is needed. + idx_needed_last = Eigen::Index(std::min(std::ceil((offset - model.m_transitions.get_time(0) + 1) / dt), + double(model.m_transitions.get_num_time_points() - 1))); + + int INStISyi = model.get_transition_flat_index( + Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), group); + + for (Eigen::Index i = idx_needed_first; i <= idx_needed_last; i++) { + + time_idx = model.m_transitions.get_time(i); + if (offset == int(std::floor(time_idx))) { + model.m_transitions[i][INStISyi] += + (1 - (time_idx - std::floor(time_idx))) * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == int(std::ceil(time_idx))) { + model.m_transitions[i][INStISyi] += + (time_idx - std::floor(time_idx)) * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == int(std::floor(time_idx - dt))) { + model.m_transitions[i][INStISyi] -= + (1 - (time_idx - dt - std::floor(time_idx - dt))) * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == int(std::ceil(time_idx - dt))) { + model.m_transitions[i][INStISyi] -= + (time_idx - dt - std::floor(time_idx - dt)) * scale_confirmed_cases * entry.num_confirmed; + } + } + + // 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. + int Di = model.get_state_flat_index(Eigen::Index(InfectionState::Dead), group); + if (offset == + std::floor(-mean_InfectedSymptomsToInfectedSevere[group] - + mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group])) { + model.m_populations[0][Di] += + (1 - + (-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] - + mean_InfectedCriticalToDead[group] - + std::floor(-mean_InfectedSymptomsToInfectedSevere[group] - + mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group]))) * + entry.num_deaths; + } + if (offset == + std::ceil(-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] - + mean_InfectedCriticalToDead[group])) { + model.m_populations[0][Di] += + (-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] - + mean_InfectedCriticalToDead[group] - + std::floor(-mean_InfectedSymptomsToInfectedSevere[group] - + mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group])) * + entry.num_deaths; + } + + if (offset == 0) { + model.m_total_confirmed_cases[group] = scale_confirmed_cases * entry.num_confirmed; + } + } + } + + 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. --- + // Set m_support_max_vector and m_derivative_vector in the model which is needed for the following computations. + model.set_transitiondistributions_support_max(dt); + model.set_transitiondistributions_derivative(dt); + + for (AgeGroup group = AgeGroup(0); group < AgeGroup(num_age_groups); group++) { + //--- 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++) { + model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), + Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt, + i + start_shift, group); + } + // Compute flow InfectedSevereToInfectedCritical for -2 * global_support_max, ..., 0. + for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) { + model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), + Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift, + group); + } + // Compute flows from InfectedSymptoms, InfectedSevere and InfectedCritical to Recovered and + // flow InfectedCriticalToDead for -global_support_max, ..., 0. + for (Eigen::Index i = -global_support_max_index; i <= 0; i++) { + // Compute flow InfectedSymptomsToRecovered. + model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered), + Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt, + i + start_shift, group); + // Compute flow InfectedSevereToRecovered. + model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToRecovered), + Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift, + group); + // Compute flow InfectedCriticalToRecovered. + model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered), + Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift, + group); + // Compute flow InfectedCriticalToDead. + model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToDead), + Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift, + group); + } + + //--- Calculate the remaining flows. --- + // 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[group] / dt)); + int EtINSi = + model.get_transition_flat_index(Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), group); + int INStISyi = model.get_transition_flat_index( + Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), group); + for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) { + model.m_transitions[i + start_shift][EtINSi] = + (1 / model.parameters.get()[group][Eigen::Index( + InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) * + model.m_transitions[i + start_shift + index_shift_mean][INStISyi]; + } + + // Compute flow SusceptibleToExposed for -global_support_max, ..., 0. + // Use mean values of the TransitionDistribution ExposedToInfectedNoSymptoms and of the + // TransitionDistribution InfectedNoSymptomsToInfectedSymptoms for the calculation. + index_shift_mean = Eigen::Index(std::round( + (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt)); + int StEi = model.get_transition_flat_index(Eigen::Index(InfectionTransition::SusceptibleToExposed), group); + + for (Eigen::Index i = -global_support_max_index; i <= 0; i++) { + model.m_transitions[i + start_shift][StEi] = + (1 / model.parameters.get()[group][Eigen::Index( + InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) * + model.m_transitions[i + start_shift + index_shift_mean][INStISyi]; + } + + // InfectedNoSymptomsToRecovered for -global_support_max, ..., 0. + // If we previously calculated the transition ExposedToInfectedNoSymptoms, we can calculate this transition + // using the standard formula. + for (Eigen::Index i = -global_support_max_index; i <= 0; i++) { + model.compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered), + Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt, i + start_shift, + group); + } + } + + // At the end of the calculation, delete all time points that are not required for the simulation. + auto transition_copy(model.m_transitions); + model.m_transitions = TimeSeries(Eigen::Index(InfectionTransition::Count) * num_age_groups); + for (Eigen::Index i = -global_support_max_index; i <= 0; i++) { + model.m_transitions.add_time_point(i * dt, transition_copy.get_value(i + start_shift)); + } + + return mio::success(); +} } // namespace isecir } // namespace mio #endif // MEMILIO_HAS_JSONCPP -#endif // IDE_INITIALFLOWS_H \ No newline at end of file +#endif // IDE_INITIALFLOWS_H diff --git a/cpp/models/ide_secir/simulation.h b/cpp/models/ide_secir/simulation.h index 6c070894ef..d90f6ac71a 100644 --- a/cpp/models/ide_secir/simulation.h +++ b/cpp/models/ide_secir/simulation.h @@ -32,14 +32,14 @@ namespace isecir { /** - * run the simulation in discrete steps and report results. + * Run the simulation in discrete steps and report results. */ class Simulation { public: /** - * @brief setup the Simulation for an IDE model. + * @brief Setup the Simulation for an IDE model. * @param[in] model An instance of the IDE model. * @param[in] dt Step size of numerical solver. */ diff --git a/cpp/simulations/IDE_paper/CMakeLists.txt b/cpp/simulations/IDE_paper/CMakeLists.txt new file mode 100644 index 0000000000..01eadb3826 --- /dev/null +++ b/cpp/simulations/IDE_paper/CMakeLists.txt @@ -0,0 +1,16 @@ +if(MEMILIO_HAS_JSONCPP AND MEMILIO_HAS_HDF5) + add_executable(ide_convergence_rate ide_convergence_rate.cpp) + target_link_libraries(ide_convergence_rate PRIVATE memilio ode_secir ide_secir Boost::filesystem ${HDF5_C_LIBRARIES}) + target_compile_options(ide_convergence_rate PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + + add_executable(ide_changepoints ide_changepoints.cpp) + target_link_libraries(ide_changepoints PRIVATE memilio ode_secir ide_secir Boost::filesystem ${HDF5_C_LIBRARIES}) + target_compile_options(ide_changepoints PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + + add_executable(compute_parameters compute_parameters.cpp) + target_compile_options(compute_parameters PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + + add_executable(ide_covid_scenario ide_covid_scenario.cpp) + target_link_libraries(ide_covid_scenario PRIVATE memilio ode_secir ide_secir Boost::filesystem ${HDF5_C_LIBRARIES}) + target_compile_options(ide_covid_scenario PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +endif() diff --git a/cpp/simulations/IDE_paper/compute_parameters.cpp b/cpp/simulations/IDE_paper/compute_parameters.cpp new file mode 100644 index 0000000000..816bdb4207 --- /dev/null +++ b/cpp/simulations/IDE_paper/compute_parameters.cpp @@ -0,0 +1,97 @@ +/* +* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC) +* +* Authors: Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include + +int main() +{ + /** With this file parameters without age distribution can be calculated to match those specified in the + covasim paper (https://doi.org/10.1371/journal.pcbi.1009149). + First, we calculate a weighted average time across the age groups. + If other probabilites than required are given, we calculate the right probabilities. + */ + + // Age group sizes are calculated using table number 12411-04-02-4-B from www.regionalstatistik.de for the date 31.12.2020. + const double age_group_sizes[] = {7752706.0, 7581868, 9483430, 10871964, 10070748, + 13304542, 10717241, 7436098, 5092743, 843691}; + const int total = 83155031; + const int numagegroups = 10; + + // Calculate value for probability InfectedSymptomsPerInfectedNoSymptoms. + const double InfectedSymptomsPerInfectedNoSymptoms[] = {0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.90}; + double resultInfectedSymptomsPerInfectedNoSymptoms = 0; + for (int i = 0; i < numagegroups; i++) { + resultInfectedSymptomsPerInfectedNoSymptoms += age_group_sizes[i] * InfectedSymptomsPerInfectedNoSymptoms[i]; + } + resultInfectedSymptomsPerInfectedNoSymptoms = resultInfectedSymptomsPerInfectedNoSymptoms / total; + + std::cout << "InfectedSymptomsPerInfectedNoSymptoms: " << resultInfectedSymptomsPerInfectedNoSymptoms << std::endl; + + // Calculate value for probability SeverePerInfectedSymptoms. + const double SeverePerInfectedNoSymptoms[] = {0.00050, 0.00165, 0.00720, 0.02080, 0.03430, + 0.07650, 0.13280, 0.20655, 0.24570, 0.24570}; + double average_SeverePerInfectedNoSymptoms = 0; + for (int i = 0; i < numagegroups; i++) { + average_SeverePerInfectedNoSymptoms += age_group_sizes[i] * SeverePerInfectedNoSymptoms[i]; + } + average_SeverePerInfectedNoSymptoms = average_SeverePerInfectedNoSymptoms / total; + double resultSeverePerInfectedSymptoms = + average_SeverePerInfectedNoSymptoms / resultInfectedSymptomsPerInfectedNoSymptoms; + + std::cout << "SeverePerInfectedSymptoms: " << resultSeverePerInfectedSymptoms << std::endl; + + // Calculate value for probability CriticalPerSevere. + const double CriticalPerInfectedNoSymptoms[] = {0.00003, 0.00008, 0.00036, 0.00104, 0.00216, + 0.00933, 0.03639, 0.08923, 0.17420, 0.17420}; + double average_CriticalPerInfectedNoSymptoms = 0; + for (int i = 0; i < numagegroups; i++) { + average_CriticalPerInfectedNoSymptoms += age_group_sizes[i] * CriticalPerInfectedNoSymptoms[i]; + } + average_CriticalPerInfectedNoSymptoms = average_CriticalPerInfectedNoSymptoms / total; + double resultCriticalPerSevere = average_CriticalPerInfectedNoSymptoms / average_SeverePerInfectedNoSymptoms; + + std::cout << "CriticalPerSevere: " << resultCriticalPerSevere << std::endl; + + // Calculate value for probability DeathsPerCritical. + const double DeathsPerInfectedNoSymptoms[] = {0.00002, 0.00002, 0.00010, 0.00032, 0.00098, + 0.00265, 0.00766, 0.02439, 0.08292, 0.16190}; + double average_DeathsPerInfectedNoSymptoms = 0; + for (int i = 0; i < numagegroups; i++) { + average_DeathsPerInfectedNoSymptoms += age_group_sizes[i] * DeathsPerInfectedNoSymptoms[i]; + } + average_DeathsPerInfectedNoSymptoms = average_DeathsPerInfectedNoSymptoms / total; + double resultDeathsPerCritical = average_DeathsPerInfectedNoSymptoms / average_CriticalPerInfectedNoSymptoms; + + std::cout << "DeathsPerCritical: " << resultDeathsPerCritical << std::endl; + + // ---- Stay times for the ODE model. ---- + std::cout << "\nAveraged stay times for the ODE model: " << std::endl; + std::cout << "TimeExposed: " << 4.5 << std::endl; + std::cout << "TimeInfectedNoSymptoms: " + << 1.1 * resultInfectedSymptomsPerInfectedNoSymptoms + + 8. * (1. - resultInfectedSymptomsPerInfectedNoSymptoms) + << std::endl; + std::cout << "TimeInfectedSymptoms: " + << 6.6 * resultSeverePerInfectedSymptoms + 8. * (1. - resultSeverePerInfectedSymptoms) << std::endl; + std::cout << "TimeInfectedSevere: " << 1.5 * resultCriticalPerSevere + 18.1 * (1. - resultCriticalPerSevere) + << std::endl; + std::cout << "TimeInfectedCritical: " << 10.7 * resultDeathsPerCritical + 18.1 * (1. - resultDeathsPerCritical) + << std::endl; +} diff --git a/cpp/simulations/IDE_paper/get_lognormal_parameters.py b/cpp/simulations/IDE_paper/get_lognormal_parameters.py new file mode 100644 index 0000000000..e2f841254d --- /dev/null +++ b/cpp/simulations/IDE_paper/get_lognormal_parameters.py @@ -0,0 +1,74 @@ +############################################################################# +# Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +# +# Authors: Anna Wendler +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# +import numpy as np +from scipy.stats import lognorm + + +def get_lognormal_parameters(mean, std): + """ + Compute shape and scale parameters to use in lognormal distribution for given mean and standard deviation. + The lognormal distribution we consider in state_age_function.h is based on the implementation in scipy and the parameters + shape and scale are defined accordingly. + """ + variance = std**2 + + mean_tmp = np.log(mean**2/np.sqrt(mean**2+variance)) + variance_tmp = np.log(variance/mean**2 + 1) + + shape = np.sqrt(variance_tmp) + scale = np.exp(mean_tmp) + + # Test if mean and std are as expected for computed shape and scale parameters. + mean_lognorm, variance_lognorm = lognorm.stats( + shape, loc=0, scale=scale, moments='mv') + + mean_test = np.exp(scale**2/2) + + if np.abs(mean_lognorm-mean) > 1e-8: + print('Distribution does not have expected mean value.') + + if np.abs(np.sqrt(variance_lognorm)-std) > 1e-8: + print('Distribution does not have expected standard deviation.') + + return round(shape, 8), round(scale, 8) + + +def get_weighted_mean(prob_1, stay_time_1, stay_time_2): + + weighted_mean = prob_1*stay_time_1 + (1-prob_1)*stay_time_2 + + return weighted_mean + + +if __name__ == '__main__': + shape, scale = get_lognormal_parameters(2.183, 1.052) + print(f"{shape:.12f}", f"{scale:.12f}") + + weighted_mean = get_weighted_mean(0.793099, 1.1, 8.0) + print(f"{weighted_mean:.6f}") + + weighted_mean = get_weighted_mean(0.078643, 6.6, 8.0) + print(f"{weighted_mean:.6f}") + + weighted_mean = get_weighted_mean(0.173176, 1.5, 18.1) + print(f"{weighted_mean:.6f}") + + weighted_mean = get_weighted_mean(0.387803, 10.7, 18.1) + print(f"{weighted_mean:.6f}") diff --git a/cpp/simulations/IDE_paper/ide_changepoints.cpp b/cpp/simulations/IDE_paper/ide_changepoints.cpp new file mode 100644 index 0000000000..d8ad99c86a --- /dev/null +++ b/cpp/simulations/IDE_paper/ide_changepoints.cpp @@ -0,0 +1,432 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Anna Wendler, Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#include "memilio/config.h" +#include "memilio/epidemiology/state_age_function.h" +#include "memilio/io/result_io.h" +#include "memilio/io/io.h" +#include "memilio/utils/time_series.h" + +#include "ide_secir/model.h" +#include "ide_secir/infection_state.h" +#include "ide_secir/parameters.h" +#include "ide_secir/simulation.h" + +#include "ode_secir/model.h" +#include "ode_secir/infection_state.h" +#include "ode_secir/parameters.h" + +#include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp" +#include "ode_secir/infection_state.h" +#include +#include + +using Vector = Eigen::Matrix; + +// Used parameters +std::map simulation_parameter = { + {"t0", 0.}, + {"dt", 0.01}, + {"total_population", 83155031.}, + {"total_confirmed_cases", 341223.}, + {"deaths", 0.}, + {"TimeExposed", 4.5}, + {"TimeInfectedNoSymptoms", 2.527617}, + {"TimeInfectedSymptoms", 7.889900}, + {"TimeInfectedSevere", 15.225278}, + {"TimeInfectedCritical", 15.230258}, + {"TransmissionProbabilityOnContact", 0.0733271}, + {"RelativeTransmissionNoSymptoms", 1}, + {"RiskOfInfectionFromSymptomatic", 0.3}, + {"Seasonality", 0.}, + {"InfectedSymptomsPerInfectedNoSymptoms", 0.793099}, + {"SeverePerInfectedSymptoms", 0.078643}, + {"CriticalPerSevere", 0.173176}, + {"DeathsPerCritical", 0.387803}, + {"cont_freq", 3.114219}}; // Computed so that we obtain constant new infections at beginning of simulation. + +/** +* @brief Function to scale the contact matrix according to factor contact_scaling after two days. +* +* @param[in] contact_scaling Factor that is applied to contact matrix after two days. +* @returns Scaled contact matrix. +*/ +mio::UncertainContactMatrix scale_contact_matrix(ScalarType contact_scaling) +{ + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, 1); + if (contact_scaling <= 1.) { + // Perform simulation with a decrease in contacts. + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, simulation_parameter["cont_freq"])); + contact_matrix[0].add_damping(0., mio::SimulationTime(2.)); + contact_matrix[0].add_damping(contact_scaling, mio::SimulationTime(2.1)); + } + else { + // Perform simulation with an increase in contacts. + contact_matrix[0] = + mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, contact_scaling * simulation_parameter["cont_freq"])); + contact_matrix[0].add_damping(1 - 1. / contact_scaling, mio::SimulationTime(-1.)); + contact_matrix[0].add_damping(1 - 1. / contact_scaling, mio::SimulationTime(2.)); + contact_matrix[0].add_damping(0., mio::SimulationTime(2.1)); + } + + return mio::UncertainContactMatrix(contact_matrix); +} + +/** +* @brief Function to compute the initial flows needed for the IDE model where we assume that we have a constant number +* of new transmissions. +* +* @returns TimeSeries containing intitial flows. +*/ +mio::TimeSeries get_initial_flows() +{ + // The initialization vector for the IDE model is calculated by defining transitions. + // Create TimeSeries with num_transitions elements. + int num_transitions = (int)mio::isecir::InfectionTransition::Count; + mio::TimeSeries init(num_transitions); + + // Add time points for initialization of transitions. + /* For this example, the intention is to create nearly constant values for SusceptiblesToExposed flow + at the beginning of the simulation. Therefore we initalize the flows accordingly constant for + SusceptiblesToExposed and derive matching values for the other flows.*/ + // 7-Tage-Inzidenz at 15.10.2020 was 34.1, see https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Situationsberichte/Okt_2020/2020-10-15-de.pdf?__blob=publicationFile. + ScalarType SusceptibleToExposed_const = (34.1 / 7.) * simulation_parameter["total_population"] / 100000.; + Eigen::VectorXd init_transitions(num_transitions); + init_transitions[(int)mio::isecir::InfectionTransition::SusceptibleToExposed] = SusceptibleToExposed_const; + init_transitions[(int)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = SusceptibleToExposed_const; + init_transitions[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = + SusceptibleToExposed_const * simulation_parameter["InfectedSymptomsPerInfectedNoSymptoms"]; + init_transitions[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered] = + SusceptibleToExposed_const * (1 - simulation_parameter["InfectedSymptomsPerInfectedNoSymptoms"]); + init_transitions[(int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] = + init_transitions[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] * + simulation_parameter["SeverePerInfectedSymptoms"]; + init_transitions[(int)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered] = + init_transitions[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] * + (1 - simulation_parameter["SeverePerInfectedSymptoms"]); + init_transitions[(int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] = + init_transitions[(int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] * + simulation_parameter["CriticalPerSevere"]; + init_transitions[(int)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = + init_transitions[(int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] * + (1 - simulation_parameter["CriticalPerSevere"]); + init_transitions[(int)mio::isecir::InfectionTransition::InfectedCriticalToDead] = + init_transitions[(int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] * + simulation_parameter["DeathsPerCritical"]; + init_transitions[(int)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = + init_transitions[(int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] * + (1 - simulation_parameter["DeathsPerCritical"]); + init_transitions = init_transitions * simulation_parameter["dt"]; + + // Add initial time point to time series. + init.add_time_point(-350, init_transitions); + // Add further time points until time 0 with constant values. + while (init.get_last_time() < simulation_parameter["t0"] - 1e-3) { + init.add_time_point(init.get_last_time() + simulation_parameter["dt"], init_transitions); + } + return init; +} + +/** +* @brief Function that simulates from time 0 until tmax using an IDE model where we apply a contact scaling after +* two days. +* +* @param[in] contact_scaling Factor that is applied to contact matrix after two days. +* @param[in] tmax Time up to which we simulate. +* @param[in] save_dir Directory where simulation results will be stored. +* @returns Any io errors that happen. +*/ +mio::IOResult> simulate_ide_model(ScalarType contact_scaling, ScalarType tmax, + std::string save_dir = "") +{ + // Initialize model. + size_t num_agegroups = 1; + mio::CustomIndexArray total_population = + mio::CustomIndexArray(mio::AgeGroup(num_agegroups), + simulation_parameter["total_population"]); + mio::CustomIndexArray deaths = + mio::CustomIndexArray(mio::AgeGroup(num_agegroups), simulation_parameter["deaths"]); + mio::CustomIndexArray total_confirmed_cases = + mio::CustomIndexArray(mio::AgeGroup(num_agegroups), + simulation_parameter["total_confirmed_cases"]); + + mio::isecir::Model model_ide(get_initial_flows(), total_population, deaths, num_agegroups, total_confirmed_cases); + + // Set working parameters. + // Set TransitionDistributions. + mio::ConstantFunction initialfunc(0); + mio::StateAgeFunctionWrapper delaydistributioninit(initialfunc); + std::vector vec_delaydistrib((int)mio::isecir::InfectionTransition::Count, + delaydistributioninit); + // ExposedToInfectedNoSymptoms + mio::LognormSurvivalFunction survivalExposedToInfectedNoSymptoms(0.32459285, 0, 4.26907484); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms].set_state_age_function( + survivalExposedToInfectedNoSymptoms); + // InfectedNoSymptomsToInfectedSymptoms + mio::LognormSurvivalFunction survivalInfectedNoSymptomsToInfectedSymptoms(0.71587510, 0, 0.85135303); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] + .set_state_age_function(survivalInfectedNoSymptomsToInfectedSymptoms); + // InfectedNoSymptomsToRecovered + mio::LognormSurvivalFunction survivalInfectedNoSymptomsToRecovered(0.24622068, 0, 7.7611400); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered].set_state_age_function( + survivalInfectedNoSymptomsToRecovered); + // InfectedSymptomsToInfectedSevere + mio::LognormSurvivalFunction survivalInfectedSymptomsToInfectedSevere(0.66258947, 0, 5.29920733); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere].set_state_age_function( + survivalInfectedSymptomsToInfectedSevere); + // InfectedSymptomsToRecovered + mio::LognormSurvivalFunction survivalInfectedSymptomsToRecovered(0.24622068, 0, 7.76114000); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered].set_state_age_function( + survivalInfectedSymptomsToRecovered); + // InfectedSevereToInfectedCritical + mio::LognormSurvivalFunction survivalInfectedSevereToInfectedCritical(1.01076765, 0, 0.90000000); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical].set_state_age_function( + survivalInfectedSevereToInfectedCritical); + // InfectedSevereToRecovered + mio::LognormSurvivalFunction survivalInfectedSevereToRecovered(0.33816427, 0, 17.09411753); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSevereToRecovered].set_state_age_function( + survivalInfectedSevereToRecovered); + // InfectedCriticalToDead + mio::LognormSurvivalFunction survivalInfectedCriticalToDead(0.42819924, 0, 9.76267505); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedCriticalToDead].set_state_age_function( + survivalInfectedCriticalToDead); + // InfectedCriticalToRecovered + mio::LognormSurvivalFunction survivalInfectedCriticalToRecovered(0.33816427, 0, 17.09411753); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedCriticalToRecovered].set_state_age_function( + survivalInfectedCriticalToRecovered); + + model_ide.parameters.set(vec_delaydistrib); + + // Set other parameters. + std::vector vec_prob((int)mio::isecir::InfectionTransition::Count, 1.); + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] = + simulation_parameter["InfectedSymptomsPerInfectedNoSymptoms"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered)] = + 1 - simulation_parameter["InfectedSymptomsPerInfectedNoSymptoms"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere)] = + simulation_parameter["SeverePerInfectedSymptoms"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToRecovered)] = + 1 - simulation_parameter["SeverePerInfectedSymptoms"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical)] = + simulation_parameter["CriticalPerSevere"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToRecovered)] = + 1 - simulation_parameter["CriticalPerSevere"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedCriticalToDead)] = + simulation_parameter["DeathsPerCritical"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedCriticalToRecovered)] = + 1 - simulation_parameter["DeathsPerCritical"]; + + model_ide.parameters.set(vec_prob); + + model_ide.parameters.get() = scale_contact_matrix(contact_scaling); + + mio::ConstantFunction constfunc(simulation_parameter["TransmissionProbabilityOnContact"]); + mio::StateAgeFunctionWrapper StateAgeFunctionWrapperide(constfunc); + model_ide.parameters.set(StateAgeFunctionWrapperide); + StateAgeFunctionWrapperide.set_distribution_parameter(simulation_parameter["RelativeTransmissionNoSymptoms"]); + model_ide.parameters.set(StateAgeFunctionWrapperide); + StateAgeFunctionWrapperide.set_distribution_parameter(simulation_parameter["RiskOfInfectionFromSymptomatic"]); + model_ide.parameters.set(StateAgeFunctionWrapperide); + + model_ide.set_tol_for_support_max(1e-6); + model_ide.check_constraints(simulation_parameter["dt"]); + + // Simulate. + mio::isecir::Simulation sim(model_ide, simulation_parameter["dt"]); + sim.advance(tmax); + + if (!save_dir.empty()) { + std::string contact_scaling_string = std::to_string(contact_scaling); + std::string tmax_string = std::to_string(tmax); + std::string dt_string = std::to_string(simulation_parameter["dt"]); + + std::string filename_ide = + save_dir + "changepoint_ide_" + contact_scaling_string.substr(0, contact_scaling_string.find(".") + 2) + + "_" + tmax_string.substr(0, tmax_string.find(".")) + "_" + dt_string.substr(0, dt_string.find(".") + 5); + + std::string filename_ide_flows = filename_ide + "_flows.h5"; + mio::IOResult save_result_status_f = + mio::save_result({sim.get_transitions()}, {0}, 1, filename_ide_flows); + + std::string filename_ide_compartments = filename_ide + "_compartments.h5"; + mio::IOResult save_result_status_c = + mio::save_result({sim.get_result()}, {0}, 1, filename_ide_compartments); + } + + // Return vector with initial compartments. + return mio::success(sim.get_result()); +} + +/** +* @brief Function that simulates from time 0 until tmax using an ODE model where we apply a contact scaling after +* two days. +* +* @param[in] init_compartments Vector containing initial values for the compartments. +* @param[in] contact_scaling Factor that is applied to contact matrix after two days. +* @param[in] tmax Time up to which we simulate. +* @param[in] save_dir Directory where simulation results will be stored. +* @returns Any io errors that happen. +*/ +mio::IOResult simulate_ode_model(Vector init_compartments, ScalarType contact_scaling, ScalarType tmax, + std::string save_dir = "") +{ + // Use ODE FlowModel. + mio::osecir::Model model_ode(1); + + // Set initial values for compartments. + // Use mio::isecir::InfectionState when accessing init_compartments since this is computed using the IDE model. + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = + init_compartments[int(mio::isecir::InfectionState::Exposed)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = + init_compartments[int(mio::isecir::InfectionState::InfectedNoSymptoms)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = + init_compartments[int(mio::isecir::InfectionState::InfectedSymptoms)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSevere}] = + init_compartments[int(mio::isecir::InfectionState::InfectedSevere)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}] = + init_compartments[int(mio::isecir::InfectionState::InfectedCritical)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Recovered}] = + init_compartments[int(mio::isecir::InfectionState::Recovered)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Dead}] = + init_compartments[int(mio::isecir::InfectionState::Dead)]; + model_ode.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, + simulation_parameter["total_population"]); + + // Set working parameters. + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeExposed"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedNoSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedSevere"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedCritical"]; + + // Set probabilities that determine proportion between compartments. + model_ode.parameters.get>()[(mio::AgeGroup)0] = + 1 - simulation_parameter["InfectedSymptomsPerInfectedNoSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["SeverePerInfectedSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["CriticalPerSevere"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["DeathsPerCritical"]; + + // Further model parameters. + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TransmissionProbabilityOnContact"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["RelativeTransmissionNoSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["RiskOfInfectionFromSymptomatic"]; + // Choose TestAndTraceCapacity very large so that riskFromInfectedSymptomatic = RiskOfInfectionFromSymptomatic. + model_ode.parameters.get>() = std::numeric_limits::max(); + // Choose ICUCapacity very large so that CriticalPerSevereAdjusted = CriticalPerSevere and deathsPerSevereAdjusted = 0. + model_ode.parameters.get>() = std::numeric_limits::max(); + + // Set Seasonality=0 so that cont_freq_eff is equal to contact_matrix. + model_ode.parameters.set>(simulation_parameter["Seasonality"]); + + model_ode.parameters.get>() = scale_contact_matrix(contact_scaling); + + model_ode.check_constraints(); + + // Set integrator and fix step size. + auto integrator = + std::make_shared>(); + integrator->set_dt_min(simulation_parameter["dt"]); + integrator->set_dt_max(simulation_parameter["dt"]); + + // Simulate. + std::vector> results_ode = mio::osecir::simulate_flows( + simulation_parameter["t0"], tmax, simulation_parameter["dt"], model_ode, integrator); + + // Save results. + if (!save_dir.empty()) { + std::string contact_scaling_string = std::to_string(contact_scaling); + std::string tmax_string = std::to_string(tmax); + std::string dt_string = std::to_string(simulation_parameter["dt"]); + + std::string filename_ode = + save_dir + "changepoint_ode_" + contact_scaling_string.substr(0, contact_scaling_string.find(".") + 2) + + "_" + tmax_string.substr(0, tmax_string.find(".")) + "_" + dt_string.substr(0, dt_string.find(".") + 5); + + std::string filename_ode_flows = filename_ode + "_flows.h5"; + mio::IOResult save_result_status_f = mio::save_result({results_ode[1]}, {0}, 1, filename_ode_flows); + + std::string filename_ode_compartments = filename_ode + "_compartments.h5"; + mio::IOResult save_result_status_c = + mio::save_result({results_ode[0]}, {0}, 1, filename_ode_compartments); + } + + return mio::success(); +} + +int main() +{ + // Paths are valid if file is executed e.g. in memilio/build/bin. + std::string save_dir = "../../data/simulation_results/changepoints/"; + // Make folder if not existent yet. + boost::filesystem::path dir(save_dir); + boost::filesystem::create_directories(dir); + + // Changepoint scenario with halving of contacts after two days. + ScalarType contact_scaling = 0.5; + ScalarType tmax = 12; + + auto result_ide = simulate_ide_model(contact_scaling, tmax, save_dir); + if (!result_ide) { + printf("%s\n", result_ide.error().formatted_message().c_str()); + return -1; + } + + // Use compartments at time 0 from IDE simulation as initial values for ODE model to make results comparable. + Vector compartments = result_ide.value().get_value(0); + + auto result_ode = simulate_ode_model(compartments, contact_scaling, tmax, save_dir); + if (!result_ode) { + printf("%s\n", result_ode.error().formatted_message().c_str()); + return -1; + } + + // Changepoint scenario with doubling of contacts after two days. + contact_scaling = 2.; + tmax = 12; + + result_ide = simulate_ide_model(contact_scaling, tmax, save_dir); + if (!result_ide) { + printf("%s\n", result_ide.error().formatted_message().c_str()); + return -1; + } + + // Use compartments at time 0 from IDE simulation as initial values for ODE model to make results comparable. + compartments = result_ide.value().get_value(0); + + result_ode = simulate_ode_model(compartments, contact_scaling, tmax, save_dir); + if (!result_ode) { + printf("%s\n", result_ode.error().formatted_message().c_str()); + return -1; + } + + return 0; +} diff --git a/cpp/simulations/IDE_paper/ide_convergence_rate.cpp b/cpp/simulations/IDE_paper/ide_convergence_rate.cpp new file mode 100644 index 0000000000..e97eba718c --- /dev/null +++ b/cpp/simulations/IDE_paper/ide_convergence_rate.cpp @@ -0,0 +1,551 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Anna Wendler +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#include "memilio/io/result_io.h" +#include "memilio/utils/time_series.h" +#include "memilio/config.h" +#include "memilio/epidemiology/state_age_function.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/math/floating_point.h" + +#include "ode_secir/infection_state.h" +#include "ode_secir/model.h" + +#include "ide_secir/infection_state.h" +// #include "ide_secir/model.h" +// #include "ide_secir/simulation.h" + +#include +#include + +/** This file can be used to create simulation results in order to compare ODE and IDE models. +* This means that the parameters of the IDE model are set in such a way that the continuous version of the IDE model +* is reduced to the ODE model used. +* The IDE model is initialized with flows from the ODE model so that both models are comparable. +* If these results are generated for different accuracies, the convergence rate of the solution scheme of the IDE model +* can be determined numerically with the ODE model with a high accuracy as ground truth. +*/ + +// Used parameters. +std::map simulation_parameter = {{"t0", 0.}, + {"total_population", 10000.}, + {"TimeExposed", 1.4}, + {"TimeInfectedNoSymptoms", 1.2}, + {"TimeInfectedSymptoms", 0.3}, + {"TimeInfectedSevere", 0.3}, + {"TimeInfectedCritical", 0.3}, + {"TransmissionProbabilityOnContact", 1.}, + {"RelativeTransmissionNoSymptoms", 1.}, + {"RiskOfInfectionFromSymptomatic", 1.}, + {"Seasonality", 0.}, + {"InfectedSymptomsPerInfectedNoSymptoms", 0.5}, + {"SeverePerInfectedSymptoms", 0.5}, + {"CriticalPerSevere", 0.5}, + {"DeathsPerCritical", 0.5}, + {"cont_freq", 1.}}; + +/** +* @brief Takes the compartments of an ODE simulation and computes the respective flows. +* +* With t_max and t_window, it can be determined for which time window the flows will be computed. +* By default, we compute the flows with the same time step size. +* It is also possible to compute the corresponding flows for a bigger time step which can be given by dt_target. +* +* @param[in] model_ode ODE-SECIR model used. +* @param[in] compartments TimeSeries containing compartments from an ODE simulation. +* The time steps must be equidistant. +* @param[out] flows TimeSeries where the computed flows will be stored. +* @param[in] t_max Maximal time for which the flows are computed. +* @param[in] t_window Time window before t_max for which flows will be computed. +* @param[in] dt_target Time step size used for the resulting TimeSeries flows. +* Default is the time step size of the ODE simulation result compartments. +* dt_target should be a multiple of the step size used for the simulation result in compartments. +*/ +void get_flows_from_ode_compartments(mio::osecir::Model& model_ode, + mio::TimeSeries compartments, mio::TimeSeries& flows, + ScalarType t_max, ScalarType t_window, ScalarType dt_target = 0.) +{ + ScalarType dt_ode = compartments.get_time(1) - compartments.get_time(0); + int num_transitions = (int)mio::isecir::InfectionTransition::Count; + // Check that the TimeSeries flows is empty as expected. + if (flows.get_num_time_points() > 0) { + flows = mio::TimeSeries(num_transitions); + } + // If dt_target is not set, use dt_ode. + if (dt_target < 1e-10) { + dt_target = dt_ode; + } + // scale_timesteps is used to get from index wrt ODE timestep to index wrt dt_target. + // Here we assume that the ODE model is solved on a finer (or equal) scale than dt_target. + ScalarType scale_timesteps = dt_target / dt_ode; + + // Compute index variables with respect to dt_target. + Eigen::Index t_window_index = Eigen::Index(std::ceil(t_window / dt_target)); + Eigen::Index t_max_index = Eigen::Index(std::ceil(t_max / dt_target)); + + Eigen::Index flows_start_index = t_max_index - t_window_index + 1; + + // Add time points to TimeSeries flows and set flow Susceptible to Exposed. + for (Eigen::Index i = flows_start_index; i <= t_max_index; i++) { + flows.add_time_point(i * dt_target, mio::TimeSeries::Vector::Constant(num_transitions, 0)); + flows.get_last_value()[Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)] += + compartments[(Eigen::Index)(scale_timesteps * (i - 1))] + [(Eigen::Index)mio::osecir::InfectionState::Susceptible] - + compartments[(Eigen::Index)(scale_timesteps * i)][(Eigen::Index)mio::osecir::InfectionState::Susceptible]; + } + + // --- Compute flows as combination of change in compartments and previously computed flows. + // Flow from Exposed to InfectedNoSymptoms. + for (Eigen::Index i = flows_start_index; i <= t_max_index; i++) { + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = + compartments[(Eigen::Index)(scale_timesteps * (i - 1))] + [(Eigen::Index)mio::osecir::InfectionState::Exposed] - + compartments[(Eigen::Index)(scale_timesteps * i)][(Eigen::Index)mio::osecir::InfectionState::Exposed] + + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::SusceptibleToExposed]; + } + ScalarType out_flow = 0; + // Flow from InfectedNoSymptoms to InfectedSymptoms and from InfectedNoSymptoms to Recovered. + for (Eigen::Index i = flows_start_index; i <= t_max_index; i++) { + out_flow = + compartments[(Eigen::Index)(scale_timesteps * (i - 1))] + [(Eigen::Index)mio::osecir::InfectionState::InfectedNoSymptoms] - + compartments[(Eigen::Index)(scale_timesteps * i)] + [(Eigen::Index)mio::osecir::InfectionState::InfectedNoSymptoms] + + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms]; + flows[i - + flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = + (1 - + model_ode.parameters.get>()[(mio::AgeGroup)0]) * + out_flow; + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered] = + model_ode.parameters.get>()[(mio::AgeGroup)0] * + out_flow; + } + // Flow from InfectedSymptoms to InfectedSevere and from InfectedSymptoms to Recovered. + for (Eigen::Index i = flows_start_index; i <= t_max_index; i++) { + out_flow = compartments[(Eigen::Index)(scale_timesteps * (i - 1))] + [(Eigen::Index)mio::osecir::InfectionState::InfectedSymptoms] - + compartments[(Eigen::Index)(scale_timesteps * i)] + [(Eigen::Index)mio::osecir::InfectionState::InfectedSymptoms] + + flows[i - flows_start_index] + [(Eigen::Index)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]; + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] = + model_ode.parameters.get>()[(mio::AgeGroup)0] * out_flow; + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered] = + (1 - model_ode.parameters.get>()[(mio::AgeGroup)0]) * + out_flow; + } + // Flow from InfectedSevere to InfectedCritical and from InfectedSevere to Recovered. + for (Eigen::Index i = flows_start_index; i <= t_max_index; i++) { + out_flow = compartments[(Eigen::Index)(scale_timesteps * (i - 1))] + [(Eigen::Index)mio::osecir::InfectionState::InfectedSevere] - + compartments[(Eigen::Index)(scale_timesteps * i)] + [(Eigen::Index)mio::osecir::InfectionState::InfectedSevere] + + flows[i - flows_start_index] + [(Eigen::Index)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere]; + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] = + model_ode.parameters.get>()[(mio::AgeGroup)0] * out_flow; + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = + (1 - model_ode.parameters.get>()[(mio::AgeGroup)0]) * out_flow; + } + // Flow from InfectedCritical to Dead and from InfectedCritical to Recovered. + for (Eigen::Index i = flows_start_index; i <= t_max_index; i++) { + out_flow = compartments[(Eigen::Index)(scale_timesteps * (i - 1))] + [(Eigen::Index)mio::osecir::InfectionState::InfectedCritical] - + compartments[(Eigen::Index)(scale_timesteps * i)] + [(Eigen::Index)mio::osecir::InfectionState::InfectedCritical] + + flows[i - flows_start_index] + [(Eigen::Index)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical]; + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::InfectedCriticalToDead] = + model_ode.parameters.get>()[(mio::AgeGroup)0] * out_flow; + flows[i - flows_start_index][(Eigen::Index)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = + (1 - model_ode.parameters.get>()[(mio::AgeGroup)0]) * out_flow; + } +} + +/** +* @brief Computes the inital flows that are needed for an IDE simulation given we have the compartments from an ODE +* simulation for an adequate time window before t0_ide. +* +* Here we assume, that the ODE and the IDE model are matching, i.e. that the parameters of the IDE model are chosen +* such that the continous version reduces to the ODE model. This is achieved by choosing exponentially distributed +* transitions with the corresponding mean stay times etc. +* The time step size of ODE and IDE simulation can be chosen independently. +* However, we assume that the time step size of the IDE model is a multiple of the one of the ODE model. +* +* @param[in] model_ode ODE model that is used. +* @param[in] model_ide IDE model that is used. +* @param[in] compartments TimeSeries containing compartments from a simulation with model_ode. +* @param[in] t0_ide Start time of IDE simulation that we want to compute initial flows for. +* @param[in] dt_ide Time step size of IDE simulation. +*/ +// void compute_initial_flows_for_ide_from_ode(mio::osecir::Model& model_ode, mio::isecir::Model& model_ide, +// mio::TimeSeries compartments, ScalarType t0_ide, +// ScalarType dt_ide) +// { +// std::cout << "Computing initial flows. \n"; + +// // Use t_window=t0_ide to get flows from t0 onwards. +// get_flows_from_ode_compartments(model_ode, compartments, model_ide.m_transitions, t0_ide, t0_ide, dt_ide); +// ScalarType dt_ode = compartments.get_time(1) - compartments.get_time(0); +// // Remove time series from previous run and set initial values in populations. +// if (model_ide.m_populations.get_num_time_points() > 0) { +// model_ide.m_populations = mio::TimeSeries((int)mio::isecir::InfectionState::Count); +// } +// model_ide.m_populations.add_time_point( +// model_ide.m_transitions.get_last_time(), +// mio::TimeSeries::Vector::Constant((int)mio::isecir::InfectionState::Count, 0)); +// model_ide.m_populations[0][Eigen::Index(mio::isecir::InfectionState::Dead)] = +// compartments[(Eigen::Index)compartments.get_num_time_points() - +// (Eigen::Index)((compartments.get_last_time() - t0_ide) / dt_ode) - 1] +// [(Eigen::Index)mio::osecir::InfectionState::Dead]; +// } + +/** +* @brief Function to remove time points from some simulation results so that not every point has to be saved afterwards. +* +* @param[in] simulation_result TimeSeries containing simulation results. Can contain compartments or flows. +* @param[in] saving_dt Step size in between the time points of the TimeSeries with less time points. +* This should be a multiple of the time step size used in simulation_results. +* @param[in] scale Factor by which the TimeSeries values should be scaled. +* @returns TimeSeries with simulation results where some time points have been removed. +*/ +mio::TimeSeries remove_time_points(const mio::TimeSeries& simulation_result, + ScalarType saving_dt, ScalarType scale = 1.) +{ + mio::TimeSeries removed(simulation_result.get_num_elements()); + ScalarType time = simulation_result.get_time(0); + removed.add_time_point(time, scale * simulation_result[0]); + time += saving_dt; + for (int i = 1; i < simulation_result.get_num_time_points(); i++) { + if (mio::floating_point_greater_equal(simulation_result.get_time(i), time, 1e-8)) { + removed.add_time_point(simulation_result.get_time(i), scale * simulation_result[i]); + time += saving_dt; + } + } + return removed; +} + +/** +* @brief Simulates from t0 until tmax using an ODE model and subsequently simulates from +* t0_ide = (tmax-t0)/2 until tmax using a matching IDE model to determine convergence of the IDE solver. +* +* To make both ODE and IDE model comparable, we set the parameters of the IDE model according to the parameters of the +* ODE model. Furthermore, we determine the initial flows for the IDE model based on the ODE results so that we have +* equivalent conditions for both models at t0_ide. + +* The time step size of ODE and IDE simulation can be chosen independently. A vector containing the desired +* ide_exponents is passed for which an IDE simulation is run. If an empty vector is passed, only an ODE simulation is +* run, e.g., to create a ground truth +* However, we assume that the time step size of the IDE model is a multiple of the one of the ODE model. +* +* @param[in] t0 Start time of the ODE simulation. +* @param[in] tmax Maximal time for which we simulate. +* @param[in] ode_exponent The ODE model is simulated using a fixed step size dt=10^{-ode_exponent}. +* @param[in] ide_exponents The IDE model is simulated using fixed step sizes dt=10^{-ide_exponent} for ide_exponent in +* ide_exponents. +* @param[in] save_exponent The results of the ODE model will be saved using the step size 10^{-save_exponent}, should +* not be larger than the maximum ide_exponent. +* @param[in] result_dir Directory where simulation results will be stored. +* @returns Any io errors that happen. +*/ +mio::IOResult simulate_ode_and_ide(ScalarType t0, ScalarType tmax, ScalarType ode_exponent, + ScalarType ide_exponent, ScalarType save_exponent, std::string result_dir) +{ + mio::unused(ide_exponent); + /********************************** + * ODE simulation * + **********************************/ + + // The ODE model is simulated using a fixed step size dt=10^{-ode_exponent}. + ScalarType dt_ode = pow(10, -ode_exponent); + + mio::osecir::Model model_ode(1); + + // Set initial values for compartments. + ScalarType nb_exp_t0 = 20, nb_car_t0 = 20, nb_inf_t0 = 3, nb_hosp_t0 = 1, nb_icu_t0 = 1, nb_rec_t0 = 10, + nb_dead_t0 = 0; + + model_ode.populations.set_total(simulation_parameter["total_population"]); + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = nb_exp_t0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = nb_car_t0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = nb_inf_t0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSevere}] = nb_hosp_t0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}] = nb_icu_t0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Recovered}] = nb_rec_t0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Dead}] = nb_dead_t0; + model_ode.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, + simulation_parameter["total_population"]); + + // Set parameters. + ScalarType cont_freq = simulation_parameter["cont_freq"]; + // Set Seasonality=0 so that cont_freq_eff is equal to contact_matrix. + model_ode.parameters.set>(simulation_parameter["Seasonality"]); + mio::ContactMatrixGroup& contact_matrix = model_ode.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); + model_ode.parameters.get>() = mio::UncertainContactMatrix(contact_matrix); + + // Parameters needed to determine transition rates. + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeExposed"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedNoSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedSevere"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedCritical"]; + + // Set probabilities that determine proportion between compartments. + model_ode.parameters.get>()[(mio::AgeGroup)0] = + 1 - simulation_parameter["InfectedSymptomsPerInfectedNoSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["SeverePerInfectedSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["CriticalPerSevere"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["DeathsPerCritical"]; + + // Further model parameters. + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TransmissionProbabilityOnContact"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["RelativeTransmissionNoSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["RiskOfInfectionFromSymptomatic"]; + // Choose TestAndTraceCapacity very large so that parameter has no effect. + model_ode.parameters.get>() = std::numeric_limits::max(); + // Choose ICUCapacity very large so that parameter has no effect. + model_ode.parameters.get>() = std::numeric_limits::max(); + + auto integrator = + std::make_shared>(); + // Choose dt_min = dt_max so that we have a fixed time step and can compare to IDE. + integrator->set_dt_min(dt_ode); + integrator->set_dt_max(dt_ode); + // Set tolerance as follows so that every time step is only computed once (found out by trying). + integrator->set_rel_tolerance(1e-1); + integrator->set_abs_tolerance(1e-1); + + std::cout << "Starting simulation with ODE model. \n"; + mio::TimeSeries secihurd_ode = + mio::osecir::simulate(t0, tmax, dt_ode, model_ode, integrator); + + if (!result_dir.empty() && save_exponent > 0) { + // Create result directory if not existent yet. + boost::filesystem::path res_dir(result_dir); + boost::filesystem::create_directory(res_dir); + auto save_result_status_ode = + mio::save_result({remove_time_points(secihurd_ode, pow(10, -save_exponent))}, {0}, 1, + result_dir + "result_ode_dt=1e-" + fmt::format("{:.0f}", ode_exponent) + "_savefrequency" + + fmt::format("{:.0f}", save_exponent) + ".h5"); + + // Compute flows from ODE result to store results. + // Note that we are computing \tilde{\sigma} here. To be able to compare flows between different timesteps (of ODE and IDE) + // we need to divide by dt to get \hat{\sigma}. This is done while saving the results. + mio::TimeSeries secihurd_ode_flows((int)mio::isecir::InfectionTransition::Count); + get_flows_from_ode_compartments(model_ode, secihurd_ode, secihurd_ode_flows, tmax, tmax - t0); + auto save_result_status_ode_flows = + mio::save_result({remove_time_points(secihurd_ode_flows, pow(10, -save_exponent), 1. / dt_ode)}, {0}, 1, + result_dir + "result_ode_flows_dt=1e-" + fmt::format("{:.0f}", ode_exponent) + + "_savefrequency" + fmt::format("{:.0f}", save_exponent) + ".h5"); + + if (save_result_status_ode && save_result_status_ode_flows) { + std::cout << "Successfully saved the ODE simulation results. \n\n"; + } + else { + return mio::failure(mio::StatusCode::InvalidValue, + "Error occured while saving the ODE simulation results."); + } + } + + // /********************************** + // * IDE simulation * + // **********************************/ + + // // Start IDE model simulation at half of tmax. + // ScalarType t0_ide = (tmax - t0) / 2.; + // // Number of deaths will be set according to the ODE model later in the function where also the transitions are calculated. + // ScalarType deaths_init_value = 0.; + + // // Initialize model. + // mio::TimeSeries init_transitions((int)mio::isecir::InfectionTransition::Count); + // size_t num_agegroups = 1; + // mio::CustomIndexArray total_population = + // mio::CustomIndexArray(mio::AgeGroup(num_agegroups), + // simulation_parameter["total_population"]); + // mio::CustomIndexArray deaths = + // mio::CustomIndexArray(mio::AgeGroup(num_agegroups), deaths_init_value); + // mio::isecir::Model model_ide(std::move(init_transitions), total_population, deaths, num_agegroups); + + // // Set parameters. + // // Contact matrix; contact_matrix was already defined for ODE. + // model_ide.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + + // // To compare with the ODE model we use ExponentialSurvivalFunctions functions as TransitionDistributions. + // // We set the parameters so that they correspond to the above ODE model. + // mio::ExponentialSurvivalFunction exponential(10.0); + // mio::StateAgeFunctionWrapper delaydistribution(exponential); + // std::vector vec_delaydistrib((int)mio::isecir::InfectionTransition::Count, + // delaydistribution); + // // ExposedToInfectedNoSymptoms + // vec_delaydistrib[(int)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms].set_distribution_parameter( + // 1. / model_ode.parameters.get>()[(mio::AgeGroup)0]); + // // InfectedNoSymptomsToInfectedSymptoms + // vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] + // .set_distribution_parameter( + // 1. / model_ode.parameters.get>()[(mio::AgeGroup)0]); + // // InfectedNoSymptomsToRecovered + // vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered].set_distribution_parameter( + // 1. / model_ode.parameters.get>()[(mio::AgeGroup)0]); + // // InfectedSymptomsToInfectedSevere + // vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] + // .set_distribution_parameter( + // 1. / model_ode.parameters.get>()[(mio::AgeGroup)0]); + // // InfectedSymptomsToRecovered + // vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered].set_distribution_parameter( + // 1. / model_ode.parameters.get>()[(mio::AgeGroup)0]); + // // InfectedSevereToInfectedCritical + // vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] + // .set_distribution_parameter( + // 1. / model_ode.parameters.get>()[(mio::AgeGroup)0]); + // // InfectedSevereToRecovered + // vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSevereToRecovered].set_distribution_parameter( + // 1. / model_ode.parameters.get>()[(mio::AgeGroup)0]); + // // InfectedCriticalToDead + // vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedCriticalToDead].set_distribution_parameter( + // 1. / model_ode.parameters.get>()[(mio::AgeGroup)0]); + // // InfectedCriticalToRecovered + // vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedCriticalToRecovered].set_distribution_parameter( + // 1. / model_ode.parameters.get>()[(mio::AgeGroup)0]); + + // model_ide.parameters.set(vec_delaydistrib); + + // // Set probabilities. + // std::vector vec_prob((int)mio::isecir::InfectionTransition::Count, 0.); + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)] = 1; + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)] = 1; + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] = + // 1 - model_ode.parameters.get>()[(mio::AgeGroup)0]; + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered)] = + // model_ode.parameters.get>()[(mio::AgeGroup)0]; + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere)] = + // model_ode.parameters.get>()[(mio::AgeGroup)0]; + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToRecovered)] = + // 1 - model_ode.parameters.get>()[(mio::AgeGroup)0]; + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical)] = + // model_ode.parameters.get>()[(mio::AgeGroup)0]; + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToRecovered)] = + // 1 - model_ode.parameters.get>()[(mio::AgeGroup)0]; + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedCriticalToDead)] = + // model_ode.parameters.get>()[(mio::AgeGroup)0]; + // vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedCriticalToRecovered)] = + // 1 - model_ode.parameters.get>()[(mio::AgeGroup)0]; + + // model_ide.parameters.set(vec_prob); + + // // Set further parameters. + // mio::ConstantFunction constfunc_proboncontact( + // model_ode.parameters.get>()[(mio::AgeGroup)0]); + // mio::StateAgeFunctionWrapper proboncontact(constfunc_proboncontact); + // model_ide.parameters.set(proboncontact); + + // mio::ConstantFunction constfunc_reltransnosympt( + // model_ode.parameters.get>()[(mio::AgeGroup)0]); + // mio::StateAgeFunctionWrapper reltransnosympt(constfunc_reltransnosympt); + // model_ide.parameters.set(reltransnosympt); + + // mio::ConstantFunction constfunc_riskofinf( + // model_ode.parameters.get>()[(mio::AgeGroup)0]); + // mio::StateAgeFunctionWrapper riskofinf(constfunc_riskofinf); + // model_ide.parameters.set(riskofinf); + + // // for (ScalarType ide_exponent : ide_exponents) { + + // // The IDE model is simulated using a fixed step size dt=10^{-ide_exponent}. + // ScalarType dt_ide = pow(10, -ide_exponent); + + // // Compute initial flows from results of ODE simulation and set initial values for populations. + // compute_initial_flows_for_ide_from_ode(model_ode, model_ide, secihurd_ode, t0_ide, dt_ide); + + // model_ide.check_constraints(dt_ide); + + // // Carry out simulation. + // std::cout << "Starting simulation with IDE model. \n"; + // mio::isecir::Simulation sim(model_ide, dt_ide); + // sim.advance(tmax); + + // std::cout << "Initialization method of the IDE model: " << sim.get_model().get_initialization_method_compartments() + // << "\n"; + // if (!result_dir.empty()) { + // // Save compartments. + // mio::TimeSeries secihurd_ide = sim.get_result(); + // auto save_result_status_ide = + // mio::save_result({secihurd_ide}, {0}, 1, + // result_dir + "result_ide_dt=1e-" + fmt::format("{:.0f}", ide_exponent) + + // "_init_dt_ode=1e-" + fmt::format("{:.0f}", ode_exponent) + ".h5"); + // // Save flows. + // mio::TimeSeries secihurd_ide_flows = sim.get_transitions(); + // auto save_result_status_ide_flows = + // mio::save_result({remove_time_points(secihurd_ide_flows, dt_ide, 1. / dt_ide)}, {0}, 1, + // result_dir + "result_ide_flows_dt=1e-" + fmt::format("{:.0f}", ide_exponent) + + // "_init_dt_ode=1e-" + fmt::format("{:.0f}", ode_exponent) + ".h5"); + // if (save_result_status_ide && save_result_status_ide_flows) { + // std::cout << "Successfully saved the IDE simulation results. \n\n"; + // } + // else { + // std::cout << "Error occured while saving the IDE simulation results. \n"; + // return mio::failure(mio::StatusCode::InvalidValue, + // "Error occured while saving the IDE simulation results."); + // } + // } + return mio::success(); +} + +int main() +{ + // Directory where results will be stored. If this string is empty, results will not be saved. + std::string result_dir = "../../data/simulation_results/convergence/"; + + // General set up. + ScalarType t0 = 0.; + ScalarType tmax = 70.; + // The ODE model will be simulated using a fixed step size dt=10^{-ode_exponent}. + ScalarType ode_exponent = 6; + // The results of the ODE model will be saved using the step size 10^{-save_exponent} + // as for very small step sizes used for the simulation, the number of time points stored gets very big. + ScalarType save_exponent = 4; + // The IDE model will be simulated using a fixed step size dt=10^{-ide_exponent} for ide_exponent in ide_exponents. + std::vector ide_exponents = {1, 2, 3, 4}; + + for (ScalarType ide_exponent : ide_exponents) { + mio::IOResult result = + simulate_ode_and_ide(t0, tmax, ode_exponent, ide_exponent, save_exponent, result_dir); + + if (!result) { + printf("%s\n", result.error().formatted_message().c_str()); + return -1; + } + } + + return 0; +} diff --git a/cpp/simulations/IDE_paper/ide_covid_scenario.cpp b/cpp/simulations/IDE_paper/ide_covid_scenario.cpp new file mode 100644 index 0000000000..7053e5a605 --- /dev/null +++ b/cpp/simulations/IDE_paper/ide_covid_scenario.cpp @@ -0,0 +1,590 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Anna Wendler, Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#include "memilio/config.h" +#include "memilio/epidemiology/state_age_function.h" +#include "memilio/io/epi_data.h" +#include "memilio/io/result_io.h" +#include "memilio/io/io.h" +#include "memilio/utils/time_series.h" + +#include "ide_secir/model.h" +#include "ide_secir/infection_state.h" +#include "ide_secir/parameters.h" +#include "ide_secir/parameters_io.h" +#include "ide_secir/simulation.h" + +#include "memilio/epidemiology/contact_matrix.h" +#include "ode_secir/model.h" +#include "ode_secir/infection_state.h" +#include "ode_secir/parameters.h" + +#include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp" +#include "boost/filesystem.hpp" +#include +#include +#include +#include + +using Vector = Eigen::Matrix; + +// Used parameters. +std::map simulation_parameter = {{"t0", 0.}, + {"dt", 0.01}, + {"total_population", 83155031.}, + {"total_confirmed_cases", 0.}, // set by RKI data + {"deaths", 0.}, // set by RKI data + {"TimeExposed", 4.5}, + {"TimeInfectedNoSymptoms", 2.527617}, + {"TimeInfectedSymptoms", 7.889900}, + {"TimeInfectedSevere", 15.225278}, + {"TimeInfectedCritical", 15.230258}, + {"TransmissionProbabilityOnContact", 0.0733271}, + {"RelativeTransmissionNoSymptoms", 1}, + {"RiskOfInfectionFromSymptomatic", 0.3}, + {"Seasonality", 0.}, + {"InfectedSymptomsPerInfectedNoSymptoms", 0.793099}, + {"SeverePerInfectedSymptoms", 0.078643}, + {"CriticalPerSevere", 0.173176}, + {"DeathsPerCritical", 0.387803}, + {"scale_confirmed_cases", 1.}, + {"scale_contacts", 1.}, + {"lockdown_hard", 371 * 14 / (45 * 401.)}}; + +/** + * Indices of contact matrix corresponding to locations where contacts occur. + */ +enum class ContactLocation +{ + Home = 0, + School, + Work, + Other, + Count, +}; + +/** + * Different types of NPIs, used as DampingType. + */ +enum class Intervention +{ + Home, + SchoolClosure, + HomeOffice, + GatheringBanFacilitiesClosure, + PhysicalDistanceAndMasks, + Count, +}; + +/** + * Different levels of NPIs, used as DampingLevel. + */ +enum class InterventionLevel +{ + Main, + PhysicalDistanceAndMasks, + Holidays, + Count, +}; + +// Map the ContactLocation%s to file names. +static const std::map contact_locations = {{ContactLocation::Home, "home"}, + {ContactLocation::School, "school_pf_eig"}, + {ContactLocation::Work, "work"}, + {ContactLocation::Other, "other"}}; + +/** + * @brief Add NPIs to a given contact matrix from 01/10/2020 on. + * + * NPIs from the Paper "Assessment of effective mitigation ..." (doi: 10.1016/j.mbs.2021.108648) are used with slight + * modifications for a period of 45 days from 01/10/2020 on. + * + * @param[in] contact_matrices The contact matrices to which the NPIs are applied to. + * @param[in] start_date Start date of the simulation. + * @param[in] lockdown_hard Proportion of counties for which a hard lockdown is implemented. + */ +void set_npi_october(mio::ContactMatrixGroup& contact_matrices, mio::Date start_date, ScalarType lockdown_hard) +{ + // ---------------------01/10/2020-------------------------------- + // NPIs from paper for October. + auto offset_npi = mio::SimulationTime(mio::get_offset_in_days(mio::Date(2020, 10, 1), start_date)); + // For the beginning of the time period, we assume only half of the defined proportion of counties is in a hard lockdown. + lockdown_hard = lockdown_hard / 2; + // Contact reduction at home. + ScalarType v = 0.2 * (1 - lockdown_hard) + lockdown_hard * 0.4; + contact_matrices[size_t(ContactLocation::Home)].add_damping(Eigen::MatrixXd::Constant(1, 1, v), + mio::DampingLevel(int(InterventionLevel::Main)), + mio::DampingType(int(Intervention::Home)), offset_npi); + // Home-Office + people stopped working. + v = (0.2 + 0.025) * (1 - lockdown_hard) + lockdown_hard * (0.3 + 0.025); + contact_matrices[size_t(ContactLocation::Work)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::Main)), + mio::DampingType(int(Intervention::HomeOffice)), offset_npi); + // GatheringBanFacilitiesClosure affects ContactLocation Other. + v = 0. * (1 - lockdown_hard) + lockdown_hard * 0.2; + contact_matrices[size_t(ContactLocation::Other)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::Main)), + mio::DampingType(int(Intervention::GatheringBanFacilitiesClosure)), offset_npi); + // PhysicalDistanceAndMasks in all locations. + v = 0.2 * (1 - lockdown_hard) + lockdown_hard * 0.4; + for (auto&& contact_location : contact_locations) { + contact_matrices[size_t(contact_location.first)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::PhysicalDistanceAndMasks)), + mio::DampingType(int(Intervention::PhysicalDistanceAndMasks)), offset_npi); + } + // Remote schooling. + v = lockdown_hard * 0.25; + contact_matrices[size_t(ContactLocation::School)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::Main)), + mio::DampingType(int(Intervention::SchoolClosure)), offset_npi); + + // ---------------------24/10/2020-------------------------------- + /* We assume that the stricter NPIs of November defined in the paper are beginning about a week earlier, + which can be observed in the RKI data. + Moreover the lockdown value of PhysicalDistanceAndMasks in the location school is assumed to apply for all counties.*/ + offset_npi = mio::SimulationTime(mio::get_offset_in_days(mio::Date(2020, 10, 24), start_date)); + // For the second half of the simulation, the proportion of counties in hard lockdown is increased to compensate for the lower proportion before. + lockdown_hard = lockdown_hard * 5; + // Contact reduction at home. + v = 0.4 * (1 - lockdown_hard) + 0.6 * lockdown_hard; + contact_matrices[size_t(ContactLocation::Home)].add_damping(Eigen::MatrixXd::Constant(1, 1, v), + mio::DampingLevel(int(InterventionLevel::Main)), + mio::DampingType(int(Intervention::Home)), offset_npi); + // Home-Office + people stopped working. + v = (0.2 + 0.05) * (1 - lockdown_hard) + lockdown_hard * (0.3 + 0.05); + contact_matrices[size_t(ContactLocation::Work)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::Main)), + mio::DampingType(int(Intervention::HomeOffice)), offset_npi); + // GatheringBanFacilitiesClosure affects ContactLocation Other. + v = 0.6 * (1 - lockdown_hard) + 0.8 * lockdown_hard; + contact_matrices[size_t(ContactLocation::Other)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::Main)), + mio::DampingType(int(Intervention::GatheringBanFacilitiesClosure)), offset_npi); + // PhysicalDistanceAndMasks in ContactLocation%s Home. + v = 0.2 * (1 - lockdown_hard) + lockdown_hard * 0.4; + contact_matrices[size_t(ContactLocation::Home)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::PhysicalDistanceAndMasks)), + mio::DampingType(int(Intervention::PhysicalDistanceAndMasks)), offset_npi); + // PhysicalDistanceAndMasks in ContactLocation%s School. + v = 0.2 * (1 - lockdown_hard) + lockdown_hard * 0.4; + contact_matrices[size_t(ContactLocation::School)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::PhysicalDistanceAndMasks)), + mio::DampingType(int(Intervention::PhysicalDistanceAndMasks)), offset_npi); + // PhysicalDistanceAndMasks in ContactLocation%s Work and Other. + v = 0.4 * (1 - lockdown_hard) + lockdown_hard * 0.6; + contact_matrices[size_t(ContactLocation::Work)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::PhysicalDistanceAndMasks)), + mio::DampingType(int(Intervention::PhysicalDistanceAndMasks)), offset_npi); + contact_matrices[size_t(ContactLocation::Other)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::PhysicalDistanceAndMasks)), + mio::DampingType(int(Intervention::PhysicalDistanceAndMasks)), offset_npi); + // Remote schooling. + v = lockdown_hard * 0.25; + contact_matrices[size_t(ContactLocation::School)].add_damping( + Eigen::MatrixXd::Constant(1, 1, v), mio::DampingLevel(int(InterventionLevel::Main)), + mio::DampingType(int(Intervention::SchoolClosure)), offset_npi); +} + +/** + * @brief Set the contact pattern of parameters for a model without division in age groups. + * + * The contacts are calculated using contact matrices from files in the data directory for different locations. + * Also set nonpharmaceutical interventions influencing the ContactPatterns used for simulation in the timeframe from start_date to end_date. + * + * @param[in] data_dir Directory to files with minimum and baseline contact matrices. + * @param[in] start_date Start date of the simulation. + * @returns Any io errors that happen during reading of the input files. + */ +mio::IOResult define_contact_matrices(const boost::filesystem::path& data_dir, + mio::Date start_date) +{ + // Files in data_dir are containing contact matrices with 6 agegroups. We use this to compute a contact pattern without division of age groups. + // Age group sizes are calculated using table number 12411-04-02-4-B from www.regionalstatistik.de for the date 31.12.2020. + const ScalarType age_group_sizes[] = {3969138.0, 7508662, 18921292, 28666166, 18153339, 5936434}; + const ScalarType total = simulation_parameter["total_population"]; + const int numagegroups = 6; + + auto contact_matrices = mio::ContactMatrixGroup(contact_locations.size(), 1); + // Load and set minimum and baseline contacts for each contact location. + for (auto&& contact_location : contact_locations) { + BOOST_OUTCOME_TRY(auto&& baseline, + mio::read_mobility_plain( + (data_dir / "contacts" / ("baseline_" + contact_location.second + ".txt")).string())); + BOOST_OUTCOME_TRY(auto&& minimum, + mio::read_mobility_plain( + (data_dir / "contacts" / ("minimum_" + contact_location.second + ".txt")).string())); + ScalarType base = 0; + ScalarType min = 0; + for (int i = 0; i < numagegroups; i++) { + for (int j = 0; j < numagegroups; j++) { + // Calculate a weighted average according to the age group sizes of the total contacts. + base += age_group_sizes[i] / total * baseline(i, j); + min += age_group_sizes[i] / total * minimum(i, j); + } + } + contact_matrices[size_t(contact_location.first)].get_baseline() = + simulation_parameter["scale_contacts"] * Eigen::MatrixXd::Constant(1, 1, base); + contact_matrices[size_t(contact_location.first)].get_minimum() = + simulation_parameter["scale_contacts"] * Eigen::MatrixXd::Constant(1, 1, min); + } + + // ----- Add NPIs to the contact matrices. ----- + // Set of NPIs for October. + if (start_date == mio::Date(2020, 10, 1)) { + set_npi_october(contact_matrices, start_date, simulation_parameter["lockdown_hard"]); + } + + return mio::success(contact_matrices); +} + +/** + * @brief Set the contact pattern of parameters for a model without division in age groups without using the + * age-resolved contact_matrices. + * + * In case of only one age group the contact matrix reduces to a 1x1 matrix. + * Instead of using contact matrices from files in the data directory for different locations as in the function + * define_contact_matrices(), we set the contact frequency to the value that we obtained using the function + * define_contact_matrices() above. Accordingly, we set the damping as in the function above to model the implementation + * of NPIs on Oct 24, 2020. + * + * @param[in] start_date Start date of the simulation. + * @returns Any io errors that happen during reading of the input files. + */ +mio::IOResult define_contact_matrices_simplified(mio::Date start_date) +{ + // Set of NPIs for October. + auto contact_matrices = mio::ContactMatrixGroup(1, 1); + auto start_npi_october = start_date; + if (start_npi_october == mio::Date(2020, 10, 1)) { + + contact_matrices[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 7.69129)); + + auto offset_npi = mio::SimulationTime(mio::get_offset_in_days(mio::Date(2020, 10, 24), start_npi_october)); + contact_matrices[0].add_damping(0., mio::SimulationTime(0.1)); + contact_matrices[0].add_damping(1 - 3.51782 / 7.69129, offset_npi); + } + + return mio::success(contact_matrices); +} + +/** +* @brief Simulates using an IDE model. +* +* @param[in] start_date Start date of the simulation +* @param[in] simulation_time Duration of the simulation. +* @param[in] contact_matrices Contact matrices used. +* @param[in] data_dir Directory to files with minimum and baseline contact matrices and reported data from RKI. +* @param[in] save_dir Directory where simulation results will be stored. +* @returns Any IO errros that happen. +*/ +mio::IOResult>> +simulate_ide_model(mio::Date start_date, ScalarType simulation_time, mio::ContactMatrixGroup contact_matrices, + const boost::filesystem::path& data_dir, std::string save_dir = "") + +{ + // Initialize model. + size_t num_agegroups = 1; + mio::CustomIndexArray total_population = + mio::CustomIndexArray(mio::AgeGroup(num_agegroups), + simulation_parameter["total_population"]); + mio::CustomIndexArray deaths = + mio::CustomIndexArray(mio::AgeGroup(num_agegroups), simulation_parameter["deaths"]); + mio::CustomIndexArray total_confirmed_cases = + mio::CustomIndexArray(mio::AgeGroup(num_agegroups), + simulation_parameter["total_confirmed_cases"]); + + mio::isecir::Model model_ide(mio::TimeSeries((int)mio::isecir::InfectionTransition::Count), + total_population, deaths, num_agegroups, total_confirmed_cases); + + // Set working parameters. + // Set TransitionDistributions. + mio::ConstantFunction initialfunc(0); + mio::StateAgeFunctionWrapper delaydistributioninit(initialfunc); + std::vector vec_delaydistrib((int)mio::isecir::InfectionTransition::Count, + delaydistributioninit); + // ExposedToInfectedNoSymptoms + mio::LognormSurvivalFunction survivalExposedToInfectedNoSymptoms(0.32459285, 0, 4.26907484); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms].set_state_age_function( + survivalExposedToInfectedNoSymptoms); + // InfectedNoSymptomsToInfectedSymptoms + mio::LognormSurvivalFunction survivalInfectedNoSymptomsToInfectedSymptoms(0.71587510, 0, 0.85135303); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] + .set_state_age_function(survivalInfectedNoSymptomsToInfectedSymptoms); + // InfectedNoSymptomsToRecovered + mio::LognormSurvivalFunction survivalInfectedNoSymptomsToRecovered(0.24622068, 0, 7.7611400); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered].set_state_age_function( + survivalInfectedNoSymptomsToRecovered); + // InfectedSymptomsToInfectedSevere + mio::LognormSurvivalFunction survivalInfectedSymptomsToInfectedSevere(0.66258947, 0, 5.29920733); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere].set_state_age_function( + survivalInfectedSymptomsToInfectedSevere); + // InfectedSymptomsToRecovered + mio::LognormSurvivalFunction survivalInfectedSymptomsToRecovered(0.24622068, 0, 7.76114000); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered].set_state_age_function( + survivalInfectedSymptomsToRecovered); + // InfectedSevereToInfectedCritical + mio::LognormSurvivalFunction survivalInfectedSevereToInfectedCritical(1.010767652595, 0, 0.90000000); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical].set_state_age_function( + survivalInfectedSevereToInfectedCritical); + // InfectedSevereToRecovered + mio::LognormSurvivalFunction survivalInfectedSevereToRecovered(0.33816427, 0, 17.09411753); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedSevereToRecovered].set_state_age_function( + survivalInfectedSevereToRecovered); + // InfectedCriticalToDead + mio::LognormSurvivalFunction survivalInfectedCriticalToDead(0.42819924, 0, 9.76267505); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedCriticalToDead].set_state_age_function( + survivalInfectedCriticalToDead); + // InfectedCriticalToRecovered + mio::LognormSurvivalFunction survivalInfectedCriticalToRecovered(0.33816427, 0, 17.09411753); + vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedCriticalToRecovered].set_state_age_function( + survivalInfectedCriticalToRecovered); + + model_ide.parameters.set(vec_delaydistrib); + + // Set other parameters. + std::vector vec_prob((int)mio::isecir::InfectionTransition::Count, 1.); + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] = + simulation_parameter["InfectedSymptomsPerInfectedNoSymptoms"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered)] = + 1 - simulation_parameter["InfectedSymptomsPerInfectedNoSymptoms"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere)] = + simulation_parameter["SeverePerInfectedSymptoms"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToRecovered)] = + 1 - simulation_parameter["SeverePerInfectedSymptoms"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical)] = + simulation_parameter["CriticalPerSevere"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToRecovered)] = + 1 - simulation_parameter["CriticalPerSevere"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedCriticalToDead)] = + simulation_parameter["DeathsPerCritical"]; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::InfectedCriticalToRecovered)] = + 1 - simulation_parameter["DeathsPerCritical"]; + + model_ide.parameters.set(vec_prob); + + model_ide.parameters.get() = + mio::UncertainContactMatrix(contact_matrices); + + mio::ConstantFunction constfunc(simulation_parameter["TransmissionProbabilityOnContact"]); + mio::StateAgeFunctionWrapper StateAgeFunctionWrapperide(constfunc); + model_ide.parameters.set(StateAgeFunctionWrapperide); + StateAgeFunctionWrapperide.set_distribution_parameter(simulation_parameter["RelativeTransmissionNoSymptoms"]); + model_ide.parameters.set(StateAgeFunctionWrapperide); + StateAgeFunctionWrapperide.set_distribution_parameter(simulation_parameter["RiskOfInfectionFromSymptomatic"]); + model_ide.parameters.set(StateAgeFunctionWrapperide); + + model_ide.set_tol_for_support_max(1e-6); + + // Set initial flows according to RKI data. + std::string path_rki = mio::path_join((data_dir / "pydata" / "Germany").string(), "cases_all_germany_ma7.json"); + + mio::IOResult init_flows = mio::isecir::set_initial_flows( + model_ide, simulation_parameter["dt"], path_rki, start_date, simulation_parameter["scale_confirmed_cases"]); + + model_ide.check_constraints(simulation_parameter["dt"]); + + // Simulate. + mio::isecir::Simulation sim(model_ide, simulation_parameter["dt"]); + + sim.advance(simulation_time); + + // Save results. + if (!save_dir.empty()) { + std::string tmax_string = std::to_string(simulation_time); + std::string dt_string = std::to_string(simulation_parameter["dt"]); + std::string filename_ide = save_dir + "ide_" + std::to_string(start_date.year) + "-" + + std::to_string(start_date.month) + "-" + std::to_string(start_date.day) + "_" + + tmax_string.substr(0, tmax_string.find(".")) + "_" + + dt_string.substr(0, dt_string.find(".") + 5); + + std::string filename_ide_flows = filename_ide + "_flows.h5"; + mio::IOResult save_result_status_f = + mio::save_result({sim.get_transitions()}, {0}, 1, filename_ide_flows); + std::string filename_ide_compartments = filename_ide + "_compartments.h5"; + mio::IOResult save_result_status_c = + mio::save_result({sim.get_result()}, {0}, 1, filename_ide_compartments); + } + + // Return vector with simulation results. + std::vector> result = {sim.get_result(), sim.get_transitions()}; + return mio::success(result); +} + +/** +* @brief Simulates using an ODE model. +* +* We need intial values for the compartments as input. With this, we can make the starting conditions equivalent to an +* the previously simulated results using an IDE model. +* +* @param[in] start_date Start date of the simulation +* @param[in] simulation_time Duration of the simulation. +* @param[in] init_compartments Vector containing initial values for compartments. +* @param[in] contact_matrices Contact matrices used. +* @param[in] save_dir Directory where simulation results will be stored. +* @returns Any IO errros that happen. +*/ +mio::IOResult simulate_ode_model(mio::Date start_date, ScalarType simulation_time, Vector init_compartments, + mio::ContactMatrixGroup contact_matrices, std::string save_dir = "") +{ + // Use ODE FlowModel. + mio::osecir::Model model_ode(1); + + // Set working parameters. + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeExposed"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedNoSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedSevere"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TimeInfectedCritical"]; + + // Set probabilities that determine proportion between compartments. + model_ode.parameters.get>()[(mio::AgeGroup)0] = + 1 - simulation_parameter["InfectedSymptomsPerInfectedNoSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["SeverePerInfectedSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["CriticalPerSevere"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["DeathsPerCritical"]; + + // Further model parameters. + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["TransmissionProbabilityOnContact"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["RelativeTransmissionNoSymptoms"]; + model_ode.parameters.get>()[(mio::AgeGroup)0] = + simulation_parameter["RiskOfInfectionFromSymptomatic"]; + // Choose TestAndTraceCapacity very large so that riskFromInfectedSymptomatic = RiskOfInfectionFromSymptomatic. + model_ode.parameters.get>() = std::numeric_limits::max(); + // Choose ICUCapacity very large so that CriticalPerSevereAdjusted = CriticalPerSevere and deathsPerSevereAdjusted = 0. + model_ode.parameters.get>() = std::numeric_limits::max(); + + // Set Seasonality=0 so that cont_freq_eff is equal to contact_matrix. + model_ode.parameters.set>(simulation_parameter["Seasonality"]); + + model_ode.parameters.get>() = + mio::UncertainContactMatrix(contact_matrices); + + // Use mio::isecir::InfectionState when accessing init_compartments since this is computed using the IDE model. + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}] = + init_compartments[int(mio::isecir::InfectionState::Susceptible)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = + init_compartments[int(mio::isecir::InfectionState::Exposed)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = + init_compartments[int(mio::isecir::InfectionState::InfectedNoSymptoms)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = + init_compartments[int(mio::isecir::InfectionState::InfectedSymptoms)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSevere}] = + init_compartments[int(mio::isecir::InfectionState::InfectedSevere)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}] = + init_compartments[int(mio::isecir::InfectionState::InfectedCritical)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Recovered}] = + init_compartments[int(mio::isecir::InfectionState::Recovered)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Dead}] = + init_compartments[int(mio::isecir::InfectionState::Dead)]; + + model_ode.check_constraints(); + + // Set integrator and fix step size. + auto integrator = + std::make_shared>(); + integrator->set_dt_min(simulation_parameter["dt"]); + integrator->set_dt_max(simulation_parameter["dt"]); + + // Simulate. + std::vector> results_ode = mio::osecir::simulate_flows( + simulation_parameter["t0"], simulation_time, simulation_parameter["dt"], model_ode, integrator); + + // Save results. + if (!save_dir.empty()) { + std::string simulation_time_string = std::to_string(simulation_time); + std::string dt_string = std::to_string(simulation_parameter["dt"]); + std::string filename_ode = save_dir + "ode_" + std::to_string(start_date.year) + "-" + + std::to_string(start_date.month) + "-" + std::to_string(start_date.day) + "_" + + simulation_time_string.substr(0, simulation_time_string.find(".")) + "_" + + dt_string.substr(0, dt_string.find(".") + 5); + + std::string filename_ode_flows = filename_ode + "_flows.h5"; + mio::IOResult save_result_status_f = mio::save_result({results_ode[1]}, {0}, 1, filename_ode_flows); + std::string filename_ode_compartments = filename_ode + "_compartments.h5"; + mio::IOResult save_result_status_c = + mio::save_result({results_ode[0]}, {0}, 1, filename_ode_compartments); + } + + return mio::success(); +} + +int main(int argc, char** argv) +{ + // Paths are valid if file is executed e.g. in memilio/build/bin. + std::string data_dir_string = "../../data"; + std::string save_dir = "../../data/simulation_results/covid_scenario/"; + + mio::Date start_date(2020, 10, 01); + + ScalarType simulation_time = 45; + + if (argc == 9) { + + data_dir_string = argv[1]; + save_dir = argv[2]; + start_date = mio::Date(std::stoi(argv[3]), std::stoi(argv[4]), std::stoi(argv[5])); + simulation_time = std::stod(argv[6]); + simulation_parameter["dt"] = std::stod(argv[7]); + simulation_parameter["scale_contacts"] = std::stod(argv[8]); + + std::cout << std::setprecision(10) << "Contact scaling: " << simulation_parameter["scale_contacts"] + << std::endl; + } + + // Make folder if not existent yet. + boost::filesystem::path dir(save_dir); + boost::filesystem::create_directories(dir); + + const boost::filesystem::path data_dir = data_dir_string; + + // Set contact matrices. + mio::ContactMatrixGroup contact_matrices = define_contact_matrices(data_dir, start_date).value(); + // mio::ContactMatrixGroup contact_matrices = define_contact_matrices_simplified(start_date).value(); + + // Run IDE simulation. + auto result_ide = simulate_ide_model(start_date, simulation_time, contact_matrices, data_dir, save_dir); + if (!result_ide) { + printf("%s\n", result_ide.error().formatted_message().c_str()); + return -1; + } + + // Use results from IDE simulation as initialization for ODE model. + Vector init_compartments = result_ide.value()[0].get_value(0); + + // Run ODE simulation. + auto result_ode = simulate_ode_model(start_date, simulation_time, init_compartments, contact_matrices, save_dir); + if (!result_ode) { + printf("%s\n", result_ode.error().formatted_message().c_str()); + return -1; + } + + return 0; +} diff --git a/cpp/simulations/IDE_paper/investigate_age_distribution.py b/cpp/simulations/IDE_paper/investigate_age_distribution.py new file mode 100644 index 0000000000..5e2e364dff --- /dev/null +++ b/cpp/simulations/IDE_paper/investigate_age_distribution.py @@ -0,0 +1,159 @@ +import os +import pandas as pd +import matplotlib.pyplot as plt +import numpy as np + + +def get_df_daily(data_dir): + """ Get a dataframe that contains daily new confirmed cases and daily deaths from RKI data. + + @param[in] data_dir Directory where file with RKI data is stored. + @returns Dataframe with daily data on confirmed cases and deaths. + """ + # Read file. + datafile = os.path.join(data_dir, "pydata", "Germany", "cases_all_age_all_dates.json") + df = pd.read_json(datafile) + + # Create df_daily, where daily confirmed (and daily deaths) will be stored. + df_daily = df.copy() + df_daily = df_daily.drop(columns=["Confirmed", "Deaths", "Recovered"]) + df_daily = df_daily[df_daily.Date != "2020-01-01"] + df_daily.insert(2, "DailyConfirmed", value=0.) + df_daily.insert(3, "DailyDeaths", value=0.) + + # Compute daily new confirmed cases from df (where cumulative confirmed cases are stored) for each agegroup. (Same for Deaths.) + agegroups = df["Age_RKI"].unique() + for agegroup in agegroups: + df_daily.loc[df_daily["Age_RKI"] == agegroup, + "DailyConfirmed"] = df.loc[df["Age_RKI"] == agegroup, "Confirmed"].diff().dropna() + df_daily.loc[df_daily["Age_RKI"] == agegroup, + "DailyDeaths"] = df.loc[df["Age_RKI"] == agegroup, "Deaths"].diff().dropna() + + return df_daily + + +def get_proportional_population_per_agegroup(): + """ + Computed the proportion of each age group compared to the total population. + """ + # Population from Table 12411-04-02-4-B from regionalstatistik.de, data from 31.12.2020. + population_per_agegroup = np.array( + [3969138, 7508662, 18921292, 28666166, 18153339, 5936434]) + total_population = population_per_agegroup.sum() + population_per_agegroup = population_per_agegroup/total_population + + return population_per_agegroup + + +def get_relevant_confirmed_cases(data_dir, start_date, T_IH, T_HU, T_U): + """Gets the confirmed cases per age group for the relevant time frame that is determined by T_IH, T_HU and T_U. + + @param[in] data_dir Directory where file with RKI data is stored. + @param[in] start_date Start date of interest. + @param[in] T_IH Mean stay time in Infected compartment before transitioning to Hospitalized. + @param[in] T_HU Mean stay time in Hospitalized compartment before transitioning to ICU. + @param[in] T_U Mean stay time in ICU. + """ + # Get dataframe with daily confirmed cases. + df = get_df_daily(data_dir) + + # Considered age groups. + agegroups = ['A00-A04', 'A05-A14', 'A15-A34', 'A35-A59', 'A60-A79', 'A80+'] + + # Extract relevant dates to be considered. + df_date = df[(df["Date"] >= pd.Timestamp(start_date)-pd.Timedelta(days=T_IH+T_HU+T_U)) + & (df["Date"] <= pd.Timestamp(start_date)-pd.Timedelta(days=T_IH+T_HU))] + + # Get total confirmed cases in considered time frame. + totaldailyconfirmed = df_date.DailyConfirmed.sum() + + # Check the proportion of unknown age group. + unknown_proportion = df_date.loc[df["Age_RKI"] == + "unknown", "DailyConfirmed"].sum()/totaldailyconfirmed + if unknown_proportion > 0.001: + print( + f"A proportion of {unknown_proportion} of the considered cases has unknown age group.") + + # Get total daily confirmed cases in considered time frame per age group. + dailyconfirmed_per_agegroup = [] + for agegroup in agegroups: + dailyconfirmed_per_agegroup.append( + df_date.loc[df["Age_RKI"] == agegroup, "DailyConfirmed"].sum()) + + # Compute proportion of totaldailyconfirmed. + dailyconfirmed_per_agegroup = dailyconfirmed_per_agegroup/totaldailyconfirmed + + return dailyconfirmed_per_agegroup + + +def plot(data_dir, start_dates, T_IH, T_HU, T_U, save_dir=""): + """ Plots proportions per age groups in confirmed cases that we expect in ICU compartment at the start_dates as + well as the proportion per age group of the total population of Germany. + + @param[in] data_dir Directory where file with RKI data is stored. + @param[in] start_dates List of start dates. + @param[in] T_IH Mean stay time in Infected compartment before transitioning to Hospitalized. + @param[in] T_HU Mean stay time in Hospitalized compartment before transitioning to ICU. + @param[in] T_U Mean stay time in ICU. + @param[in] save_dir Directory where plot will be stored. + """ + + agegroups = ['A00-A04', 'A05-A14', 'A15-A34', 'A35-A59', 'A60-A79', 'A80+'] + + # Get proportions per age group for June and Ocotber scenario and corresponding share of total population. + proportions_june = get_relevant_confirmed_cases(data_dir, + start_dates[0], T_IH, T_HU, T_U) + proportions_october = get_relevant_confirmed_cases(data_dir, + start_dates[1], T_IH, T_HU, T_U) + population_per_agegroup = get_proportional_population_per_agegroup() + + proportions = [proportions_june, + proportions_october, population_per_agegroup] + labels = ["June", "October", "Population"] + colors = [plt.cm.viridis(x) for x in np.linspace(0., 0.9, 3)] + + # Plot. + fig, ax = plt.subplots() + x_axis = np.arange(len(agegroups)) + width = 0.2 + shift = -width + for i in range(len(proportions)): + ax.bar(x_axis+shift, proportions[i], width=width, + label=f"{labels[i]}", color=colors[i]) + shift += width + + fig.supxlabel('Age groups') + fig.supylabel('Proportion') + plt.legend() + plt.xticks(x_axis, agegroups) + + plt.tight_layout() + + # Save result. + if save_dir != "": + if not os.path.isdir(save_dir): + os.makedirs(save_dir) + plt.savefig(save_dir + "proportions_per_agegroup.png", + bbox_inches='tight', dpi=500) + + +def main(): + # Path where file with RKI case data is stored. + data_dir = os.path.join(os.path.dirname( + __file__), "../../..", "data/") + + # Path where plots will be stored. + save_dir = os.path.join(os.path.dirname( + __file__), "../../..", "data/plots/") + + # Mean stay times according to Covasim paper. + T_IH = 6.6 + T_HU = 1.5 + T_U = 15.230258 + + start_dates = ["2020-06-01", "2020-10-01"] + + plot(data_dir, start_dates, T_IH, T_HU, T_U, save_dir) + +if __name__ == "__main__": + main() diff --git a/cpp/simulations/IDE_paper/plot_changepoint.py b/cpp/simulations/IDE_paper/plot_changepoint.py new file mode 100644 index 0000000000..c26e3f026e --- /dev/null +++ b/cpp/simulations/IDE_paper/plot_changepoint.py @@ -0,0 +1,118 @@ +############################################################################# +# Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +# +# Authors: Anna Wendler, Lena Ploetzke +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# +import h5py +import os +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt + +from memilio.epidata import getDataIntoPandasDataFrame as gd + + +def plot_changepoint(files, fileending="", save_dir=""): + """ + Plots the result of the changepoint simulation. + + @param[in] files Expects list of two files with ODE and IDE simulation results for flows, respectively, in this order. + @param[in] fileending Determines file ending of saved plot. + @param[in] save_dir Directory where plot will be stored. + """ + + fig, ax = plt.subplots() + legendplot = list(["ODE", "IDE"]) + + # Define colors, we use helmholtzdarkblue, helmholtzclaim. + colors = [(0, 40/255, 100/255), (20/255, 200/255, 255/255)] + linestyles = ['-', '--'] + # Add results to plot- + for file in range(len(files)): + # Load data. + h5file = h5py.File(str(files[file]) + '.h5', 'r') + + if (len(list(h5file.keys())) > 1): + raise gd.DataError("File should contain one dataset.") + if (len(list(h5file[list(h5file.keys())[0]].keys())) > 3): + raise gd.DataError("Expected only one group.") + + data = h5file[list(h5file.keys())[0]] + + + # As there should be only one Group, total is the simulation result. + total = data['Total'][:, :] + + dates = data['Time'][:] + + timestep = np.diff(dates)[0] + tmax = dates[-1] + + # Get indices where dates are >=0. + indices = np.where(dates >= 0) + + # Plot data. + # ODE + if file == 0: + # Transform cumulative flows to absolute flows, + # then transform from flows over time interval to flows at time points. + ax.plot(dates[indices[0][1:]], np.diff(total[indices[0], 0])/np.diff(dates[indices[0]]), label=legendplot[file], + color=colors[file], linestyle=linestyles[file]) + + # IDE + elif file == 1: + # Transform from flows over time interval to flows at time points. + ax.plot(dates[1:], total[1:, 0]/np.diff(dates), label=legendplot[file], + color=colors[file], linestyle=linestyles[file]) + + h5file.close() + + ax.set_xlim(left=0, right=tmax) + ax.grid(True, linestyle='--', alpha=0.5) + ax.legend(fontsize=12) + + fig.supxlabel('Simulation time [days]') + fig.supylabel('Daily new transmissions') + plt.subplots_adjust(left=None, bottom=None, right=None, + top=None, wspace=None, hspace=0.6) + + plt.tight_layout() + + # Save result. + if save_dir != "": + if not os.path.isdir(save_dir): + os.makedirs(save_dir) + plt.savefig(save_dir + f"changepoint_{fileending}.png", + bbox_inches='tight', dpi=500) + + +if __name__ == '__main__': + # Paths are valid if script is executed e.g. in memilio/cpp/simulations/IDE_paper + # Path where simulation results (generated with ide_changepoints.cpp) are stored. + result_dir = os.path.join(os.path.dirname( + __file__), "../../..", "data/simulation_results/changepoints/") + # Path where plots will be stored. + save_dir = os.path.join(os.path.dirname( + __file__), "../../..", "data/plots/changepoints/") + + plot_changepoint([os.path.join(result_dir, f"changepoint_ode_0.5_12_0.0100_flows"), + os.path.join(result_dir, f"changepoint_ide_0.5_12_0.0100_flows")], + fileending="0.5_12_0.0100", save_dir=save_dir) + + plot_changepoint([os.path.join(result_dir, f"changepoint_ode_2.0_12_0.0100_flows"), + os.path.join(result_dir, f"changepoint_ide_2.0_12_0.0100_flows")], + fileending="2.0_12_0.0100", save_dir=save_dir) diff --git a/cpp/simulations/IDE_paper/plot_conv_paper.py b/cpp/simulations/IDE_paper/plot_conv_paper.py new file mode 100644 index 0000000000..2a82c3f50f --- /dev/null +++ b/cpp/simulations/IDE_paper/plot_conv_paper.py @@ -0,0 +1,376 @@ +############################################################################# +# Copyright (C) 2020-2023 German Aerospace Center (DLR-SC) +# +# Authors: Anna Wendler +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# +import h5py +import os +import numpy as np +import matplotlib.pyplot as plt + +import memilio.epidata.getDataIntoPandasDataFrame as gd + +from matplotlib.markers import MarkerStyle +from matplotlib.transforms import Affine2D + + +def read_groundtruth(data_dir, ode_exponent, save_exponent, flows=False): + """ Read groundtruth from data. We define the groundtruth as the results obtained by the ODE model with timestep dt=1e-6. + + @param[in] data_dir Directory where h5 files are stored. + @param[in] ode_exponent Exponent that determines time step size via dt =10^{-ode_exponent}. + @param[in] save_exponent The results of the ODE model were saved using the step size 10^{-save_exponent}. + @param[in] flows Bool that determines whether we consider flows or compartments. Default is False. + @returns Dict with results of ODE model. + """ + model = 'ode' + results = {model: []} + if flows: + h5file = h5py.File(os.path.join( + data_dir, f'result_{model}_flows_dt=1e-{ode_exponent:.0f}_savefrequency{save_exponent:.0f}.h5'), 'r') + else: + h5file = h5py.File(os.path.join( + data_dir, f'result_{model}_dt=1e-{ode_exponent:.0f}_savefrequency{save_exponent:.0f}.h5'), 'r') + + if (len(list(h5file.keys())) > 1): + raise gd.DataError("File should contain one dataset.") + if (len(list(h5file[list(h5file.keys())[0]].keys())) > 3): + raise gd.DataError("Expected only one group.") + + data = h5file[list(h5file.keys())[0]] + + if flows: + # Flows are already scaled to one day. + results[model].append(data['Total'][:, :]) + else: + if len(data['Total'][0]) == 8: + # As there should be only one Group, total is the simulation result + results[model].append(data['Total'][:, :]) + elif len(data['Total'][0]) == 10: + # in ODE there are two compartments we don't use, throw these out + results[model].append( + data['Total'][:, [0, 1, 2, 4, 6, 7, 8, 9]]) + else: + raise gd.DataError( + 'Expected a different size of vector in time series.') + + dates = data['Time'][:] + + h5file.close() + + return results + + +def read_data(data_dir, ode_exponent, exponents_ide, flows=False): + """ Read data into a dict, where the keys correspond to the respective model. + At the moment we are only storing results of the IDE model here. There + we have an array that contains all results for SECIHURD for all time points + for each time step size that is investigated. + + @param[in] data_dir Directory where h5 files are stored. + @param[in] ode_exponent Exponent that determines time step size of ODE simulation via dt =10^{-ode_exponent}. + @param[in] exponents_ide List of considered exponents that determine time step size of IDE simulation via + dt =10^{-exponent_ide}. + @param[in] flows Bool that determines whether we consider flows or compartments. Default is False. + @returns Dict with results of ODE model per time step size. + """ + models = ['ide'] + results = {models[0]: []} + for model in models: + for exponent in exponents_ide: + if flows: + h5file = h5py.File(os.path.join( + data_dir, f'result_{model}_flows_dt=1e-{exponent:.0f}_init_dt_ode=1e-{ode_exponent:.0f}.h5'), 'r') + else: + h5file = h5py.File(os.path.join( + data_dir, f'result_{model}_dt=1e-{exponent:.0f}_init_dt_ode=1e-{ode_exponent:.0f}.h5'), 'r') + + data = h5file[list(h5file.keys())[0]] + + if flows: + # Flows are already scaled to one day. + results[model].append(data['Total'][:, :]) + else: + if len(data['Total'][0]) == 8: + # As there should be only one Group, total is the simulation result. + results[model].append(data['Total'][:, :]) + elif len(data['Total'][0]) == 10: + # In ODE there are two compartments we don't use, throw these out. + results[model].append( + data['Total'][:, [0, 1, 2, 4, 6, 7, 8, 9]]) + else: + raise gd.DataError( + "Expected a different size of vector in time series.") + + h5file.close() + + return results + + +def compute_l2_norm(timeseries, timestep): + """ Computes L2 norm of a time series. + + @param[in] timeseries Considered timeseries. + @param[in] Time step size. + @returns Norm. + """ + norm = np.sqrt(timestep * np.sum(timeseries**2)) + return norm + + +def compute_relerror_norm_l2(groundtruth, results, save_exponent, timesteps_ide, flows=False): + """ Computes relative L2 norm of the difference between time series from ODE and time series + from IDE for all compartments/flows. + + @param[in] groundtruth Result obtained with ODE model. + @param[in] results Results obtained with IDE model fordifferent time step sizes. + @param[in] save_exponent The results of the ODE model were saved using the step size 10^{-save_exponent}. + @param[in] timesteps_ide List of time steps used in IDE simulations. + @param[in] flows Bool that determines whether we consider flows or compartments. Default is False. + @param[in] Array that contains computed errors. + """ + if flows: + num_errors = 10 + else: + num_errors = 8 + errors = [] + + # Compute error. + for i in range(len(results['ide'])): + errors.append([]) + for compartment in range(num_errors): + timestep = timesteps_ide[i] + scale_timesteps = timestep/pow(10, -save_exponent) + num_timepoints = len(results['ide'][i]) + # Only consider compartments of ODE model for t>=t0_IDE. Use that 2*t0_IDE = t_max. + if flows: + difference = groundtruth['ode'][0][int(scale_timesteps*(num_timepoints/2-1))::int( + scale_timesteps)][:, compartment]-results['ide'][i][int(num_timepoints/2-1):][:, compartment] + norm_groundtruth = compute_l2_norm(groundtruth['ode'][0][int(scale_timesteps*(num_timepoints/2-1))::int( + scale_timesteps)][:, compartment], timestep) + errors[i].append(compute_l2_norm( + difference, timestep)/norm_groundtruth) + else: + difference = groundtruth['ode'][0][int(scale_timesteps*(num_timepoints - + 1))::int(scale_timesteps)][:, compartment]-results['ide'][i][:, compartment] + norm_groundtruth = compute_l2_norm(groundtruth['ode'][0][int(scale_timesteps*(num_timepoints - + 1))::int(scale_timesteps)][:, compartment], timestep) + errors[i].append(compute_l2_norm( + difference, timestep)/norm_groundtruth) + + return np.array(errors) + + +def plot_convergence(errors, timesteps_ide, flows=False, save_dir=""): + """ Plots errors against timesteps with a subplot for each compartment /flow. + + @param[in] errors Array that contains computed errors of IDE model compared to groundtruth. + @param[in] timesteps_ide List of time steps used in IDE simulations. + @param[in] flows Bool that determines whether we consider flows or compartments. Default is False. + @param[in] save_dir Directory where plot will be stored. + """ + # Define subplots and labels. + if flows: + fig, axs = plt.subplots( + 5, 2, sharex=True, sharey=True, figsize=(10, 10)) + num_plots = 10 + secir_dict = {0: r"$\sigma_S^E$", 1: r"$\sigma_E^C$", 2: r"$\sigma_C^I$", 3: r"$\sigma_C^R$", 4: r"$\sigma_I^H$", + 5: r"$\sigma_I^R$", 6: r"$\sigma_H^U$", 7: r"$\sigma_H^R$", 8: r"$\sigma_U^D$", 9: r"$\sigma_U^R$"} + labels = [ + r"$\Vert {\widehat{\sigma}}_{\text{IDE}} - {\widehat{\sigma}}_{\text{ODE}}\Vert_{2,\text{rel}}$", r"$\mathcal{O}(\Delta t)$"] + else: + fig, axs = plt.subplots(4, 2, sharex=True, figsize=(10, 8)) + num_plots = 8 + secir_dict = {0: 'Susceptible', 1: 'Exposed', 2: 'Carrier', 3: 'Infected', 4: 'Hospitalized', + 5: 'ICU', 6: 'Recovered', 7: 'Dead'} + labels = [ + r"$\Vert \widehat{Z}_{\text{IDE}} - \widehat{Z}_{\text{ODE}}\Vert_{2,\text{rel}}$", r"$\mathcal{O}(\Delta t)$"] + + # helmholtzdarkblue, helmholtzclaim + colors = [(0, 40/255, 100/255), (20/255, 200/255, 255/255)] + + for i in range(num_plots): + # Plot results. + axs[int(i/2), i % 2].plot(timesteps_ide, + errors[:, i], '-o', color=colors[1]) + + # Plot comparison line for linear convergence. + comparison = [dt for dt in timesteps_ide] + axs[int(i/2), i % 2].plot(timesteps_ide, comparison, + '--', color='gray', linewidth=1.2) + + # Adapt plots. + axs[int(i/2), i % 2].set_xscale("log", base=10) + axs[int(i/2), i % 2].set_yscale("log", base=10) + + axs[int(i/2), i % 2].set_title(secir_dict[i], fontsize=10) + axs[int(i/2), i % 2].grid(True, linestyle='--', alpha=0.6) + + fig.supxlabel(r'Time step $\Delta t$', fontsize=12) + + # Invert x axis only for one plot so that sharex=True and invert_xaxis work as intended. + axs[0, 0].invert_xaxis() + + lgd = fig.legend(labels, ncol=2, loc='outside lower center', + fontsize=14, bbox_to_anchor=(0.5, -0.06), bbox_transform=fig.transFigure) + plt.tight_layout(pad=0, w_pad=0.5, h_pad=0.1) + if save_dir != "": + if not os.path.isdir(save_dir): + os.makedirs(save_dir) + if flows: + plt.savefig(f'{save_dir}/convergence_all_flows.png', format='png', bbox_extra_artists=(lgd,), bbox_inches='tight', + dpi=500) + else: + plt.savefig(f'{save_dir}/convergence_all_compartments.png', format='png', bbox_extra_artists=(lgd,), bbox_inches='tight', + dpi=500) + + +def plot_convergence_oneplot(errors, timesteps_ide, flows=False, save_dir=""): + """ Plot errors against timesteps. This function creates one plot to depict all compartments/flows, respectively. + + @param[in] errors Array that contains computed errors of IDE model compared to groundtruth. + @param[in] timesteps_ide List of time steps used in IDE simulations. + @param[in] flows Bool that determines whether we consider flows or compartments. Default is False. + @param[in] save_dir Directory where plot will be stored. + """ + plt.figure() + + if flows: + secir_dict = {0: r"$\sigma_S^E$", 1: r"$\sigma_E^C$", 2: r"$\sigma_C^I$", 3: r"$\sigma_C^R$", 4: r"$\sigma_I^H$", + 5: r"$\sigma_I^R$", 6: r"$\sigma_H^U$", 7: r"$\sigma_H^R$", 8: r"$\sigma_U^D$", 9: r"$\sigma_U^R$"} + plt.ylabel( + r"$err_{\text{rel}}$", fontsize=10) + else: + secir_dict = {0: 'Susceptible', 1: 'Exposed', 2: 'Carrier', 3: 'Infected', 4: 'Hospitalized', + 5: 'ICU', 6: 'Recovered', 7: 'Dead'} + plt.ylabel( + r"$err_{\text{rel}}$", fontsize=10) + + if flows: + num_lines = 10 + else: + num_lines = 8 + + colors = [plt.cm.viridis(x) for x in np.linspace(0, 1, num_lines)] + + # Angles to rotate markers of the plot. + angles = [0, 45, 0, 45, 0, 45, 0, 45, 0, 45] + + for i in range(num_lines): + # Plot results. + rotation = Affine2D().rotate_deg(angles[i]) + plt.plot(timesteps_ide, + errors[:, i], '-', marker=MarkerStyle('x', 'full', rotation), markersize=5, color=colors[i], label=secir_dict[i]) + + # Plot comparison line for linear convergence. + comparison = [dt for dt in timesteps_ide] + plt.plot(timesteps_ide, comparison, '--', color='gray', + label=r"$\mathcal{O}(\Delta t)$") + + # Adapt plots. + plt.xscale("log", base=10) + plt.yscale("log", base=10) + plt.gca().invert_xaxis() + + plt.xlabel(r'Time step $\Delta t$', fontsize=10) + + plt.legend(fontsize=10, framealpha=0.5, ncol=2) + plt.grid(True, linestyle='--', alpha=0.6) + plt.tight_layout() + + if save_dir != "": + if not os.path.isdir(save_dir): + os.makedirs(save_dir) + if flows: + plt.savefig(f'{save_dir}/convergence_flows.png', format='png', + dpi=500) + else: + plt.savefig(f'{save_dir}/convergence_compartments.png', format='png', + dpi=500) + + +def compute_order_of_convergence(errors, timesteps_ide, flows=False): + """ Compute order of convergence between two consecutive time step sizes. + + @param[in] errors Array that contains computed errors of IDE model compared to groundtruth. + @param[in] timesteps_ide List of time steps used in IDE simulations. + @param[in] flows Bool that determines whether we consider flows or compartments. Default is False. + """ + if flows: + num_orders = 10 + else: + num_orders = 8 + + order = [] + for compartment in range(num_orders): + order.append([]) + for i in range(len(errors)-1): + order[compartment].append(np.log(errors[i+1][compartment]/errors[i][compartment]) / + np.log(timesteps_ide[i+1]/timesteps_ide[i])) + return np.array(order) + + +def main(): + # Path where simulation results (generated with ide_convergence_rate.cpp) are stored. + result_dir = os.path.join(os.path.dirname( + __file__), "../../..", "data/simulation_results/convergence/") + + # Path where plots will be stored. + save_dir = os.path.join(os.path.dirname( + __file__), "../../..", "data/plots/convergence/") + + # The ODE model was simulated using a fixed step size dt=10^{-ode_exponent}. + ode_exponent = 6 + # The results of the ODE model were saved using the step size 10^{-save_exponent} + # as for very small step sizes used for the simulation, the number of time points stored gets very big. + save_exponent = 4 + # The IDE model was simulated using a fixed step size dt=10^{-ide_exponent} for ide_exponent in ide_exponents. + exponents_ide = [1, 2, 3, 4] + # Calculate time steps resulting from exponents_ide. + timesteps_ide = [] + for exp in exponents_ide: + timesteps_ide.append(pow(10, -exp)) + + # Plot compartments and flows. + flow_bools = [False, True] + + for flow_bool in flow_bools: + # Read groundtruth (from ODE model). + groundtruth = read_groundtruth(result_dir, ode_exponent, save_exponent, flow_bool) + + # Read results from IDE simulations. + results = read_data(result_dir, ode_exponent, exponents_ide, flow_bool) + + # Compute relative L2 error norm of IDE results compared to groundtruth. + relerrors_l2 = compute_relerror_norm_l2( + groundtruth, results, save_exponent, timesteps_ide, flow_bool) + + # Plot convergence of all compartments/flows in one plot, respectively. + plot_convergence_oneplot( + relerrors_l2, timesteps_ide, flow_bool, save_dir) + # Plot convergence of all compartments/flows separately. + plot_convergence(relerrors_l2, timesteps_ide, flow_bool, save_dir) + + # # Determine order of convergence + # order = compute_order_of_convergence( + # relerrors_l2, timesteps_ide, flow_bool) + + # print('Orders of convergence: ', order) + + +if __name__ == '__main__': + main() diff --git a/cpp/simulations/IDE_paper/plot_real_scenario.py b/cpp/simulations/IDE_paper/plot_real_scenario.py new file mode 100644 index 0000000000..85b7f24eb7 --- /dev/null +++ b/cpp/simulations/IDE_paper/plot_real_scenario.py @@ -0,0 +1,508 @@ +############################################################################# +# Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +# +# Authors: Anna Wendler, Lena Ploetzke +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# +import h5py +import os +import math +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt + +from memilio.epidata import getDataIntoPandasDataFrame as gd + +# Define parameters used for simulation, also used for plotting reported data. +parameters = { + 'TimeExposed': 4.5, + 'TimeInfectedNoSymptoms': 2.527617, + 'TimeInfectedSymptoms': 7.889900, + 'TimeInfectedSevere': 15.225278, + 'TimeInfectedNoSymptomsToInfectedSymptoms': 1.1, + 'TimeInfectedSymptomsToInfectedSevere': 6.6, + 'TimeInfectedSymptomsToRecovered': 8.0, + 'TimeInfectedSevereToInfectedCritical': 1.5, + 'TimeInfectedCriticalToDead': 10.7, + 'InfectedSymptomsPerInfectedNoSymptoms': 0.793099, + 'SeverePerInfectedSymptoms': 0.078643, + 'start_date': pd.Timestamp('2020.10.01') - pd.DateOffset(days=20), + 'end_date': pd.Timestamp('2020.10.01') + pd.DateOffset(days=30), + 'scaleConfirmed': 1. +} + + +def set_parameters_U(parameters, T_UD=parameters["TimeInfectedCriticalToDead"], mu_IH=parameters["SeverePerInfectedSymptoms"]): + parameters["TimeInfectedCriticalToDead"] = T_UD + parameters["SeverePerInfectedSymptoms"] = mu_IH + return parameters + + +def load_data(file, start_date, simulation_time): + """ Loads RKI data and computes 'InfectedSymptoms', 'Deaths' and 'NewInfectionsDay' using scales, dates etc from the dictionary parameters. + Method matches the method for computing initial values for the LCT model. See also cpp/models/lct_secir/parameters_io.h. + @param[in] file Path to the RKI data file for whole Germany. Can be downloaded eg via pycode/memilio-epidata/memilio/epidata/getCaseData.py. + """ + # set_parameters_U(parameters, T_UD, mu_IH) + parameters['start_date'] = start_date - pd.DateOffset(days=0) + parameters['end_date'] = start_date + pd.DateOffset(days=simulation_time) + # Read data. + df = pd.read_json(file) + df = df.drop(columns=['Recovered']) + + # Remove unnecessary dates. + df = df[(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=-math.ceil(parameters['TimeInfectedSymptomsToInfectedSevere'] + parameters['TimeInfectedSevereToInfectedCritical'] + parameters['TimeInfectedCriticalToDead']))) + & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=math.ceil(parameters['TimeExposed'] + parameters['TimeInfectedNoSymptomsToInfectedSymptoms'])))] + # Scale confirmed cases because of undetected infections. + df['Confirmed'] = parameters['scaleConfirmed'] * df['Confirmed'] + # df2 stores the result of the computation. + df2 = df.copy() + df2 = df2[(df['Date'] >= parameters['start_date']) + & (df['Date'] <= parameters['end_date'])] + df2 = df2.reset_index() + df2 = df2.drop(columns=['index', 'Confirmed', 'Deaths']) + # Calculate individuals in compartment InfectedSymptoms. + help_I = df['Confirmed'][(df['Date'] >= parameters['start_date']) + & (df['Date'] <= parameters['end_date'])].to_numpy() + + # shift according to T_I^H and T_I^R + help_I = help_I - parameters["SeverePerInfectedSymptoms"]*((1 - math.fmod(parameters['TimeInfectedSymptomsToInfectedSevere'], 1)) * df['Confirmed'][(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=-math.floor(parameters['TimeInfectedSymptomsToInfectedSevere']))) + & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=-math.floor(parameters['TimeInfectedSymptomsToInfectedSevere'])))].to_numpy()) \ + - (1-parameters["SeverePerInfectedSymptoms"])*((1 - math.fmod(parameters['TimeInfectedSymptomsToRecovered'], 1)) * df['Confirmed'][(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=-math.floor(parameters['TimeInfectedSymptomsToRecovered']))) + & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=-math.floor(parameters['TimeInfectedSymptomsToRecovered'])))].to_numpy()) + help_I = help_I - parameters["SeverePerInfectedSymptoms"]*(math.fmod(parameters['TimeInfectedSymptomsToInfectedSevere'], 1) * df['Confirmed'][(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=-math.ceil( + parameters['TimeInfectedSymptomsToInfectedSevere']))) & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=-math.ceil(parameters['TimeInfectedSymptomsToInfectedSevere'])))].to_numpy()) \ + - (1-parameters["SeverePerInfectedSymptoms"])*(math.fmod(parameters['TimeInfectedSymptomsToRecovered'], 1) * df['Confirmed'][(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=-math.ceil( + parameters['TimeInfectedSymptomsToRecovered']))) & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=-math.ceil(parameters['TimeInfectedSymptomsToRecovered'])))].to_numpy()) + df2['InfectedSymptoms'] = help_I + # Calculate number of dead individuals. + help_D = (1 - (1 - math.fmod(parameters['TimeInfectedSymptomsToInfectedSevere'] + parameters['TimeInfectedSevereToInfectedCritical'] + parameters['TimeInfectedCriticalToDead'], 1))) * df['Deaths'][(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=-math.ceil(parameters['TimeInfectedSymptomsToInfectedSevere'] + parameters['TimeInfectedSevereToInfectedCritical'] + parameters['TimeInfectedCriticalToDead']))) + & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=-math.ceil(parameters['TimeInfectedSymptomsToInfectedSevere'] + parameters['TimeInfectedSevereToInfectedCritical'] + parameters['TimeInfectedCriticalToDead'])))].to_numpy() + help_D = help_D + (1 - math.fmod(parameters['TimeInfectedSymptomsToInfectedSevere'] + parameters['TimeInfectedSevereToInfectedCritical'] + parameters['TimeInfectedCriticalToDead'], 1)) * df['Deaths'][(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=-math.floor(parameters['TimeInfectedSymptomsToInfectedSevere'] + parameters['TimeInfectedSevereToInfectedCritical'] + parameters['TimeInfectedCriticalToDead']))) + & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=-math.floor(parameters['TimeInfectedSymptomsToInfectedSevere'] + parameters['TimeInfectedSevereToInfectedCritical'] + parameters['TimeInfectedCriticalToDead'])))].to_numpy() + df2['Deaths'] = help_D + # df2['Deaths'] = df['Deaths'][(df['Date'] >= parameters['start_date']) + # & (df['Date'] <= parameters['end_date'])].to_numpy() + # Calculate new infections per day. + fmod = math.fmod( + parameters['TimeInfectedNoSymptomsToInfectedSymptoms'] + parameters['TimeExposed'], 1) + help_newE = fmod * df['Confirmed'][(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=math.ceil(parameters['TimeInfectedNoSymptomsToInfectedSymptoms'] + parameters['TimeExposed']))) + & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=math.ceil(parameters['TimeInfectedNoSymptomsToInfectedSymptoms'] + parameters['TimeExposed'])))].to_numpy() + help_newE = help_newE + (1 - 2 * fmod) * df['Confirmed'][(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=math.floor(parameters['TimeInfectedNoSymptomsToInfectedSymptoms'] + parameters['TimeExposed']))) + & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=math.floor(parameters['TimeInfectedNoSymptomsToInfectedSymptoms'] + parameters['TimeExposed'])))].to_numpy() + help_newE = help_newE - (1 - fmod) * df['Confirmed'][(df['Date'] >= parameters['start_date'] + pd.DateOffset(days=math.floor(parameters['TimeInfectedNoSymptomsToInfectedSymptoms'] + parameters['TimeExposed'] - 1))) + & (df['Date'] <= parameters['end_date'] + pd.DateOffset(days=math.floor(parameters['TimeInfectedNoSymptomsToInfectedSymptoms'] + parameters['TimeExposed'] - 1)))].to_numpy() + df2['NewInfectionsDay'] = help_newE / \ + parameters['InfectedSymptomsPerInfectedNoSymptoms'] + return df2 + + +def get_scale_contacts(files, start_date, simulation_time): + + datafile = os.path.join(os.path.dirname( + __file__), "..", "data", "pydata", "Germany", "cases_all_germany.json") + data_rki = load_data(datafile, start_date, simulation_time) + + # Load IDE data. + for file in range(len(files)): + # load data + h5file = h5py.File(str(files[file]) + '.h5', 'r') + + if (len(list(h5file.keys())) > 1): + raise gd.DataError("File should contain one dataset.") + if (len(list(h5file[list(h5file.keys())[0]].keys())) > 3): + raise gd.DataError("Expected only one group.") + + data = h5file[list(h5file.keys())[0]] + + # As there should be only one Group, total is the simulation result + total = data['Total'][:, :] + + dates = data['Time'][:] + timestep = np.diff(dates)[0] + + date_idx = int(-simulation_time / timestep) + + print( + f"IDE new infections on {dates[date_idx]}: ", total[date_idx, 0] / timestep) + new_infections_ide = total[date_idx, 0] / timestep + + # New infections from RKI at first timestep. + new_infections_rki = data_rki[data_rki["Date"] == start_date]["NewInfectionsDay"].values[0] + timestep * (data_rki[data_rki["Date"] == start_date + pd.DateOffset( + days=1)]["NewInfectionsDay"].values[0] - data_rki[data_rki["Date"] == start_date]["NewInfectionsDay"].values[0]) + print(f"Expected new infections at {timestep}: {new_infections_rki}") + + scale_contacts = new_infections_rki/new_infections_ide + + return scale_contacts + + +def get_scale_confirmed_cases(files, start_date): + for file in range(len(files)): + # load data + h5file = h5py.File(str(files[file]) + '.h5', 'r') + + if (len(list(h5file.keys())) > 1): + raise gd.DataError("File should contain one dataset.") + if (len(list(h5file[list(h5file.keys())[0]].keys())) > 3): + raise gd.DataError("Expected only one group.") + + data = h5file[list(h5file.keys())[0]] + + if len(data['Total'][0]) == 8: + # As there should be only one Group, total is the simulation result + total = data['Total'][:, :] + elif len(data['Total'][0]) == 10: + # in ODE there are two compartments we don't use, throw these out + total = data['Total'][:, [0, 1, 2, 4, 6, 7, 8, 9]] + + dates = data['Time'][:] + + icu_ide = total[:, 5][0] + print(f"ICU in IDE on {start_date}: ", + total[:, 5][0]) + + datafile_icu = os.path.join(os.path.dirname( + __file__), "..", "data", "pydata", "Germany", "germany_divi.json") + + df = pd.read_json(datafile_icu) + + icu_divi = df[df["Date"] == start_date]["ICU"].values[0] + print(f"ICU from DIVI (ma7) on {start_date}: ", + df[df["Date"] == start_date]["ICU"].values[0]) + + scale_confirmed_cases = icu_divi/icu_ide + + return scale_confirmed_cases + + +def plot_new_infections(files, start_date, simulation_time, fileending="", save=True, save_dir='plots/'): + + datafile = os.path.join(os.path.dirname( + __file__), "..", "data", "pydata", "Germany", "cases_all_germany.json") + data_rki = load_data(datafile, start_date, simulation_time) + + fig, ax = plt.subplots() + + ax.scatter(np.linspace(0, simulation_time, simulation_time + 1), + data_rki["NewInfectionsDay"], marker="x", s=20, color='gray', label="Extrapolated RKI data") + + legendplot = list(["ODE", "IDE"]) + # helmholtzdarkblue, helmholtzclaim + colors = [(0, 40 / 255, 100 / 255), (20 / 255, 200 / 255, 255 / 255)] + linestyles = ['-', '-'] + # add results to plot + for file in range(len(files)): + # load data + h5file = h5py.File(str(files[file]) + '.h5', 'r') + + if (len(list(h5file.keys())) > 1): + raise gd.DataError("File should contain one dataset.") + if (len(list(h5file[list(h5file.keys())[0]].keys())) > 3): + raise gd.DataError("Expected only one group.") + + data = h5file[list(h5file.keys())[0]] + + # As there should be only one Group, total is the simulation result + total = data['Total'][:, :] + + dates = data['Time'][:] + timestep = np.diff(dates)[0] + + # get indices where dates are >=0 + # indices = np.where(dates >= 0) + # plot data + # ODE + if file == 0: + print(f"New infections from RKI on {start_date}: ", + data_rki[data_rki["Date"] == start_date]["NewInfectionsDay"].values[0]) + print(f"Expected new infections at {timestep}: ", + data_rki[data_rki["Date"] == start_date]["NewInfectionsDay"].values[0] + timestep * (data_rki[data_rki["Date"] == start_date + pd.DateOffset(days=1)]["NewInfectionsDay"].values[0] - data_rki[data_rki["Date"] == start_date]["NewInfectionsDay"].values[0])) + + print(f"New infections from RKI on {start_date + pd.DateOffset(days=1)}: ", + data_rki[data_rki["Date"] == start_date + pd.DateOffset(days=1)]["NewInfectionsDay"].values[0]) + # transform cumulative flows to flows absolute flows + # then transform from flows over time interval to flows at time + # points + ax.plot(dates[1:], np.diff(total[:, 0]) / np.diff(dates), label=legendplot[file], + color=colors[file], linestyle=linestyles[file]) + + date_idx = 1 + print(f"ODE new infections at dates {dates[date_idx]}: ", (np.diff( + total[:, 0]) / np.diff(dates))[date_idx - 1]) + date_idx = int(1 / np.diff(dates)[0]) + print(f"ODE new infections on {dates[date_idx]}: ", (np.diff( + total[:, 0]) / np.diff(dates))[date_idx - 1]) + + # IDE + elif file == 1: + # transform from flows over time interval to flows at time points + ax.plot(dates[1:], total[1:, 0] / np.diff(dates), label=legendplot[file], + color=colors[file], linestyle=linestyles[file]) + + date_idx = int(-simulation_time / timestep - 1) + print( + f"IDE new infections on {dates[date_idx]}: ", total[date_idx, 0] / timestep) + date_idx = int(-simulation_time / timestep) + print( + f"IDE new infections on {dates[date_idx]}: ", total[date_idx, 0] / timestep) + date_idx = int(-(simulation_time - 1) / timestep - 1) + print( + f"IDE new infections at {dates[date_idx]}: ", total[date_idx, 0] / timestep) + + h5file.close() + + ax.set_xlim(left=0, right=simulation_time) + ax.set_ylim(bottom=0) + if start_date == pd.Timestamp("2020-6-1"): + ax.set_ylim(bottom=0, top=1000) + ax.grid(True, linestyle='--', alpha=0.5) + ax.legend(fontsize=12) + + # Define x-ticks. + datelist = np.array(pd.date_range(parameters["start_date"].date(), + periods=simulation_time+1, freq='D').strftime('%m-%d').tolist()) + tick_range = (np.arange(int((simulation_time) / 5) + 1) * 5) + plt.xticks(tick_range, datelist[tick_range], + rotation=45, fontsize=12) + plt.xticks(np.arange(simulation_time), minor=True) + + fig.supxlabel('Date') + fig.supylabel(r'Daily new transmissions') + plt.subplots_adjust(left=None, bottom=None, right=None, + top=None, wspace=None, hspace=0.6) + + plt.tight_layout() + + # save result + if save: + if not os.path.isdir(save_dir): + os.makedirs(save_dir) + plt.savefig(save_dir + f"NewInfections_{fileending}.png", + bbox_inches='tight', dpi=500) + + +def plot_infectedsymptoms_deaths( + files, start_date, simulation_time, fileending="", save=True, save_dir='plots/'): + + datafile = os.path.join(os.path.dirname( + __file__), "..", "data", "pydata", "Germany", "cases_all_germany.json") + data_rki = load_data(datafile, start_date, simulation_time) + + datafile_ma7 = os.path.join(os.path.dirname( + __file__), "..", "data", "pydata", "Germany", "cases_all_germany_ma7.json") + data_rki_ma7 = load_data(datafile_ma7, start_date, + simulation_time) + + legendplot = list(["ODE", "IDE"]) + # helmholtzdarkblue, helmholtzclaim + colors = [(0, 40 / 255, 100 / 255), (20 / 255, 200 / 255, 255 / 255)] + linestyles = ['-', '-'] + # add results to plot + + compartments = [["InfectedSymptoms", 3], ["Deaths", 7]] + + for compartment in range(len(compartments)): + + print(f"{compartments[compartment][0]} from RKI on {start_date}: ", + data_rki_ma7[data_rki["Date"] == start_date][compartments[compartment][0]].values[0]) + + fig, ax = plt.subplots() + + ax.scatter(np.linspace(0, simulation_time, simulation_time + 1), + data_rki[compartments[compartment][0]], marker="x", s=20, color='gray', label="Extrapolated RKI data") + + for file in range(len(files)): + # load data + h5file = h5py.File(str(files[file]) + '.h5', 'r') + + if (len(list(h5file.keys())) > 1): + raise gd.DataError("File should contain one dataset.") + if (len(list(h5file[list(h5file.keys())[0]].keys())) > 3): + raise gd.DataError("Expected only one group.") + + data = h5file[list(h5file.keys())[0]] + + if len(data['Total'][0]) == 8: + # As there should be only one Group, total is the simulation + # result + total = data['Total'][:, :] + elif len(data['Total'][0]) == 10: + # in ODE there are two compartments we don't use, throw these + # out + total = data['Total'][:, [0, 1, 2, 4, 6, 7, 8, 9]] + + dates = data['Time'][:] + + ax.plot(dates, total[:, compartments[compartment][1]], label=legendplot[file], + color=colors[file], linestyle=linestyles[file]) + + if file == 0: + print(f"{compartments[compartment][0]} in ODE on {start_date}: ", + total[:, compartments[compartment][1]][0]) + + if file == 1: + print(f"{compartments[compartment][0]} in IDE on {start_date}: ", + total[:, compartments[compartment][1]][0]) + + h5file.close() + + ax.set_xlim(left=0, right=simulation_time) + ax.grid(True, linestyle='--', alpha=0.5) + ax.legend(fontsize=12) + + # Define x-ticks. + datelist = np.array(pd.date_range(parameters["start_date"].date(), + periods=simulation_time+1, freq='D').strftime('%m-%d').tolist()) + tick_range = (np.arange(int((simulation_time) / 5) + 1) * 5) + plt.xticks(tick_range, datelist[tick_range], + rotation=45, fontsize=12) + plt.xticks(np.arange(simulation_time), minor=True) + + fig.supxlabel('Date') + if compartment == 0: + fig.supylabel( + f'Mildly symptomatic individuals') + if compartment == 1: + fig.supylabel( + f'Deaths') + plt.subplots_adjust(left=None, bottom=None, right=None, + top=None, wspace=None, hspace=0.6) + plt.tight_layout() + # save result + if save: + if not os.path.isdir(save_dir): + os.makedirs(save_dir) + plt.savefig(save_dir + f"{compartments[compartment][0]}_{fileending}.png", + bbox_inches='tight', dpi=500) + + +def plot_icu( + files, start_date, simulation_time, fileending="", save=True, save_dir='plots/'): + + # datafile_icu_ma7 = os.path.join(os.path.dirname( + # __file__), "..", "data", "pydata", "Germany", "germany_divi_ma7.json") + + datafile_icu = os.path.join(os.path.dirname( + __file__), "..", "data", "pydata", "Germany", "germany_divi.json") + + df = pd.read_json(datafile_icu) + # data_icu_ma7 = load_data(datafile_icu_ma7, start_date, simulation_time) + + # print("Infectedsymptoms from RKI (ma7) on 1.10.2020: ", data_rki_ma7[data_rki_ma7["Date"]=="2020-10-01"]["InfectedSymptoms"].values[0]) + + # helmholtzdarkblue, helmholtzclaim + legendplot = list(["ODE", "IDE"]) + colors = [(0, 40 / 255, 100 / 255), (20 / 255, 200 / 255, 255 / 255)] + linestyles = ['-', '-'] + # add results to plot + + compartments = [["InfectedCritical", 5]] + + for compartment in range(len(compartments)): + + print(f"{compartments[compartment][0]} from DIVI (ma7) on {start_date}: ", + df[df["Date"] == start_date]["ICU"].values[0]) + + fig, ax = plt.subplots() + + ax.scatter(np.linspace(0, simulation_time, simulation_time + 1), + df[(df["Date"] >= start_date) & (df["Date"] <= start_date + pd.DateOffset(days=simulation_time))]["ICU"], marker="x", s=20, color='gray', label="DIVI data") + + for file in range(len(files)): + # load data + h5file = h5py.File(str(files[file]) + '.h5', 'r') + + if (len(list(h5file.keys())) > 1): + raise gd.DataError("File should contain one dataset.") + if (len(list(h5file[list(h5file.keys())[0]].keys())) > 3): + raise gd.DataError("Expected only one group.") + + data = h5file[list(h5file.keys())[0]] + + if len(data['Total'][0]) == 8: + # As there should be only one Group, total is the simulation result + total = data['Total'][:, :] + elif len(data['Total'][0]) == 10: + # in ODE there are two compartments we don't use, throw these out + total = data['Total'][:, [0, 1, 2, 4, 6, 7, 8, 9]] + + dates = data['Time'][:] + + ax.plot(dates, total[:, compartments[compartment][1]], label=legendplot[file], + color=colors[file], linestyle=linestyles[file]) + + if file == 0: + print(f"{compartments[compartment][0]} in ODE on {start_date}: ", + total[:, compartments[compartment][1]][0]) + + if file == 1: + print(f"{compartments[compartment][0]} in IDE on {start_date}: ", + total[:, compartments[compartment][1]][0]) + + h5file.close() + + ax.set_xlim(left=0, right=simulation_time) + ax.grid(True, linestyle='--', alpha=0.5) + ax.legend(fontsize=12) + + # Define x-ticks. + datelist = np.array(pd.date_range(parameters["start_date"].date(), + periods=simulation_time+1, freq='D').strftime('%m-%d').tolist()) + tick_range = (np.arange(int((simulation_time) / 5) + 1) * 5) + plt.xticks(tick_range, datelist[tick_range], + rotation=45, fontsize=12) + plt.xticks(np.arange(simulation_time), minor=True) + + fig.supxlabel('Date') + fig.supylabel( + f'ICU patients') + plt.subplots_adjust(left=None, bottom=None, right=None, + top=None, wspace=None, hspace=0.6) + + plt.tight_layout() + + # save result + if save: + if not os.path.isdir(save_dir): + os.makedirs(save_dir) + plt.savefig(save_dir + f"{compartments[compartment][0]}_{fileending}.png", + bbox_inches='tight', dpi=500) + + +if __name__ == '__main__': + # Path to simulation results + data_dir = os.path.join(os.path.dirname( + __file__), "..", "results/real/") + + start_date = '2020-10-1' + simulation_time = 45 + timestep = "0.0100" + + plot_new_infections([os.path.join(data_dir, f"ode_{start_date}_{simulation_time}_{timestep}_flows"), + os.path.join(data_dir, f"ide_{start_date}_{simulation_time}_{timestep}_flows")], + pd.Timestamp(start_date), simulation_time, + fileending=f"{start_date}_{simulation_time}_{timestep}", save=True, save_dir=f"../plots/covid_scenario_no_contact_scaling/{start_date}/{simulation_time}/") + + plot_infectedsymptoms_deaths([os.path.join(data_dir, f"ode_{start_date}_{simulation_time}_{timestep}_compartments"), + os.path.join(data_dir, f"ide_{start_date}_{simulation_time}_{timestep}_compartments")], + pd.Timestamp( + start_date), simulation_time, + fileending=f"{start_date}_{simulation_time}_{timestep}", save=True, save_dir=f"../plots/covid_scenario_no_contact_scaling/{start_date}/{simulation_time}/") + + plot_icu([os.path.join(data_dir, f"ode_{start_date}_{simulation_time}_{timestep}_compartments"), + os.path.join(data_dir, f"ide_{start_date}_{simulation_time}_{timestep}_compartments")], + pd.Timestamp(start_date), simulation_time, fileending=f"{start_date}_{simulation_time}_{timestep}", save=True, save_dir=f'plots/covid_scenario_no_contact_scaling/{start_date}/{simulation_time}/') diff --git a/cpp/simulations/IDE_paper/run_and_plot_covid_scenario.py b/cpp/simulations/IDE_paper/run_and_plot_covid_scenario.py new file mode 100644 index 0000000000..4df3c173cb --- /dev/null +++ b/cpp/simulations/IDE_paper/run_and_plot_covid_scenario.py @@ -0,0 +1,108 @@ +############################################################################# +# Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +# +# Authors: Anna Wendler +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# +import subprocess +import os +import pandas as pd + +from get_lognormal_parameters import get_lognormal_parameters +from plot_real_scenario import plot_new_infections, plot_infectedsymptoms_deaths, plot_icu, get_scale_contacts, get_scale_confirmed_cases + + +def run_real_scenario(data_dir, save_dir, start_date, simulation_time, timestep, scale_contacts): + + year = start_date.split("-")[0] + month = start_date.split("-")[1] + day = start_date.split("-")[2] + + subprocess.call([f"./build/bin/ide_covid_scenario", data_dir, save_dir, + f"{year}", f"{month}", f"{day}", f"{simulation_time}", f"{timestep}", f"{scale_contacts}"]) + + +def contact_scaling(save_dir, start_date, simulation_time, timestep): + scale_contacts = get_scale_contacts([os.path.join( + save_dir, f"ide_{start_date}_{simulation_time}_{timestep}_flows")], pd.Timestamp(start_date), simulation_time) + + return scale_contacts + + +def confirmed_cases_scaling(save_dir, start_date, simulation_time, timestep): + scale_confirmed_cases = get_scale_confirmed_cases([os.path.join( + save_dir, f"ide_{start_date}_{simulation_time}_{timestep}_compartments")], pd.Timestamp(start_date)) + + return scale_confirmed_cases + + +def plot_real_scenario(save_dir, plot_dir, start_date, simulation_time, timestep): + plot_new_infections([os.path.join(save_dir, f"ode_{start_date}_{simulation_time}_{timestep}_flows"), + os.path.join(save_dir, f"ide_{start_date}_{simulation_time}_{timestep}_flows")], + pd.Timestamp(start_date), simulation_time, + fileending=f"{start_date}_{simulation_time}_{timestep}", save=True, save_dir=plot_dir) + + plot_infectedsymptoms_deaths([os.path.join(save_dir, f"ode_{start_date}_{simulation_time}_{timestep}_compartments"), + os.path.join(save_dir, f"ide_{start_date}_{simulation_time}_{timestep}_compartments")], + pd.Timestamp( + start_date), simulation_time, + fileending=f"{start_date}_{simulation_time}_{timestep}", save=True, save_dir=plot_dir) + + plot_icu([os.path.join(save_dir, f"ode_{start_date}_{simulation_time}_{timestep}_compartments"), + os.path.join(save_dir, f"ide_{start_date}_{simulation_time}_{timestep}_compartments")], + pd.Timestamp(start_date), simulation_time, fileending=f"{start_date}_{simulation_time}_{timestep}", save=True, save_dir=plot_dir) + + +def run_scenario(data_dir, save_dir, plot_dir, start_date, simulation_time, timestep): + # First run the simulation with a contact scaling of 1. + scale_contacts = 1. + scale_confirmed_cases = 1. + run_real_scenario(data_dir, save_dir, start_date, simulation_time, + timestep, scale_contacts) + # Then determine contact scaling such that IDE results and RKI new infections match at first timestep. + scale_contacts = contact_scaling( + save_dir, start_date, simulation_time, timestep) + # scale_confirmed_cases=confirmed_cases_scaling(save_dir, start_date, simulation_time, timestep) + print(scale_confirmed_cases) + # Run simulation with new contact scaling. + run_real_scenario(data_dir, save_dir, start_date, simulation_time, + timestep, scale_contacts) + plot_real_scenario(save_dir, plot_dir, start_date, + simulation_time, timestep) + + +def october_scenario(timestep): + start_date = '2020-10-1' + simulation_time = 45 + + data_dir = "./data" + save_dir = f"./results/real/" + plot_dir = f"./plots/covid_scenario/{start_date}/" + + run_scenario(data_dir, save_dir, plot_dir, start_date, simulation_time, + timestep) + + +def main(): + # Folder in plots/real/save_folder where plots will be stored. + + timestep = "0.0100" + october_scenario(timestep) + + +if __name__ == "__main__": + + main() diff --git a/cpp/tests/test_ide_parameters_io.cpp b/cpp/tests/test_ide_parameters_io.cpp index dcad2031f9..feb6726dbc 100644 --- a/cpp/tests/test_ide_parameters_io.cpp +++ b/cpp/tests/test_ide_parameters_io.cpp @@ -25,6 +25,7 @@ #include "memilio/config.h" #include "memilio/epidemiology/age_group.h" +#include "memilio/io/epi_data.h" #include "memilio/utils/date.h" #include "memilio/io/io.h" @@ -39,6 +40,7 @@ TEST(TestIDEParametersIo, RKIcompareWithPreviousRun) { // Define start date and the total population used for the initialization. size_t num_agegroups = 6; + mio::CustomIndexArray total_population = mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 15 * 1e6); auto start_date = mio::Date(2020, 11, 1); @@ -86,8 +88,8 @@ TEST(TestIDEParametersIo, RKIcompareWithPreviousRun) } // Calculate initialization. - auto status = - mio::isecir::set_initial_flows(model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_age_ma7.json"), start_date); + auto status = mio::isecir::set_initial_flows( + model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_age_ma7.json"), start_date, 1.); ASSERT_THAT(print_wrap(status), IsSuccess()); @@ -147,22 +149,22 @@ TEST(TestIDEParametersIo, ParametersIoRKIFailure) // --- Case where start_date is later than maximal provided date in file. auto start_date = mio::Date(2021, 01, 05); - auto status = - mio::isecir::set_initial_flows(model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_age_ma7.json"), start_date); + auto status = mio::isecir::set_initial_flows( + model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_age_ma7.json"), start_date, 1.); ASSERT_THAT(print_wrap(status), IsFailure(mio::StatusCode::OutOfRange)); // --- Case where not all needed dates from the future are provided. start_date = mio::Date(2020, 12, 31); - status = - mio::isecir::set_initial_flows(model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_age_ma7.json"), start_date); + status = mio::isecir::set_initial_flows( + model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_age_ma7.json"), start_date, 1.); 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, 10, 1); - status = - mio::isecir::set_initial_flows(model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_age_ma7.json"), start_date); + status = mio::isecir::set_initial_flows( + model, dt, mio::path_join(TEST_DATA_DIR, "cases_all_age_ma7.json"), start_date, 1.); // 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. @@ -178,8 +180,8 @@ TEST(TestIDEParametersIo, ParametersIoRKIFailure) } // --- 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); + status = mio::isecir::set_initial_flows( + model, dt, mio::path_join(TEST_DATA_DIR, "test_empty_file.json"), start_date, 1.); ASSERT_THAT(print_wrap(status), IsFailure(mio::StatusCode::InvalidFileFormat)); diff --git a/tools/compare_ide_ode.py b/tools/compare_ide_ode.py new file mode 100644 index 0000000000..b86b70add9 --- /dev/null +++ b/tools/compare_ide_ode.py @@ -0,0 +1,236 @@ +############################################################################# +# Copyright (C) 2020-2023 German Aerospace Center (DLR-SC) +# +# Authors: Lena Ploetzke +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# +"""@plot_results_lct_secir.py +functions to plot results of a simulation with a LCT SECIR model with subcompartments. +There is also a method to compare different results of different models. +The data to be plotted should be stored in a '../data/simulation_lct' folder as .h5 files. +Data could be generated eg by executing the file ./cpp/examples/lct_secir_compare.cpp. +""" + +import h5py +import os +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt + +from memilio.epidata import getDataIntoPandasDataFrame as gd + +# Define compartments +secir_dict = {0: 'Susceptible', 1: 'Exposed', 2: 'Carrier', 3: 'Infected', 4: 'Hospitalized', + 5: 'ICU', 6: 'Recovered', 7: 'Dead'} + + +def compare_results(files, legendplot, flows=True, fileending="", save=True, save_dir=None): + """ Creates a 4x2 Plot with one subplot per compartment and one line per result one wants to compare. + @param[in] files: paths of the files (without file extension .h5) with the simulation results that should be compared. + Results should contain exactly 8 compartments (so use accumulated numbers for LCT models). Names can be given in form of a list. + One could compare results with eg different parameters or different models. + @param[in] legendplot: list with names for the results that should be used for the legend of the plot. + @param[in] save: if save is True, the plot is saved in a folder named Plots. + """ + if flows: + + secir_dict = {0: r"$\sigma_S^E$", 1: r"$\sigma_E^C$", 2: r"$\sigma_C^I$", 3: r"$\sigma_C^R$", 4: r"$\sigma_I^H$", + 5: r"$\sigma_I^R$", 6: r"$\sigma_H^U$", 7: r"$\sigma_H^R$", 8: r"$\sigma_U^D$", 9: r"$\sigma_U^R$"} + + fig, axs = plt.subplots(5, 2, sharex='all', num='Compare files') + num_plots = 10 + + else: + + # Define compartments + secir_dict = {0: 'Susceptible', 1: 'Exposed', 2: 'Carrier', 3: 'Infected', 4: 'Hospitalized', + 5: 'ICU', 6: 'Recovered', 7: 'Dead'} + + fig, axs = plt.subplots(4, 2, sharex='all', num='Compare files') + num_plots = 8 + + # helmholtzdarkblue, helmholtzclaim + colors = [(0, 40/255, 100/255), (20/255, 200/255, 255/255)] + linestyles = ['-', '--'] + # add results to plot + for file in range(len(files)): + # load data + h5file = h5py.File(str(files[file]) + '.h5', 'r') + + if (len(list(h5file.keys())) > 1): + raise gd.DataError("File should contain one dataset.") + if (len(list(h5file[list(h5file.keys())[0]].keys())) > 3): + raise gd.DataError("Expected only one group.") + + data = h5file[list(h5file.keys())[0]] + + if flows: + # As there should be only one Group, total is the simulation result + if len(data['Total'][0]) == 10: + total = data['Total'][:, :] + elif len(data['Total'][0]) == 15: + # in ODE: two compartments more are used which results in more flows; + # throw out additional flows + total = data['Total'][:, [0, 1, 2, 3, 6, 7, 10, 11, 13, 14]] + else: + if len(data['Total'][0]) == 8: + # As there should be only one Group, total is the simulation result + total = data['Total'][:, :] + elif len(data['Total'][0]) == 10: + # in ODE there are two compartments we don't use, throw these out + total = data['Total'][:, [0, 1, 2, 4, 6, 7, 8, 9]] + + dates = data['Time'][:] + + # plot data + if flows: + # ODE + if file == 0: + for i in range(num_plots): + axs[int(i/2), i % 2].plot(dates[1:], + np.diff(total[:, i])/np.diff(dates), label=legendplot[file], color=colors[file], linestyle=linestyles[file]) + # IDE + elif file == 1: + for i in range(num_plots): + axs[int(i/2), i % 2].plot(dates[1:], + total[1:, i]/np.diff(dates), label=legendplot[file], color=colors[file], linestyle=linestyles[file]) + min_date = dates[1] + else: + if file == 0: + for i in range(num_plots): + axs[int(i/2), i % 2].plot(dates, + total[:, i], label=legendplot[file], color=colors[file], linestyle=linestyles[file]) + elif file == 1: + for i in range(num_plots): + axs[int(i/2), i % 2].plot(dates, + total[:, i], label=legendplot[file], color=colors[file], linestyle=linestyles[file]) + min_date = dates[0] + + h5file.close() + + # define some characteristics of the plot + for i in range(num_plots): + axs[int(i/2), i % 2].set_title(secir_dict[i], fontsize=8) + # axs[int(i/2), i % 2].set_ylim(bottom=10**3, top=5*10**3) + axs[int(i/2), i % 2].set_xlim(left=-20, right=dates[-1]) + axs[int(i/2), i % 2].grid(True, linestyle='--') + # axs[int(i/2), i % 2].legend(fontsize=8) + axs[int(i/2), i % 2].ticklabel_format(axis='y', + style='sci', scilimits=(0, 0)) + + labels = ['ODE', 'IDE'] + fig.legend(labels, bbox_to_anchor=(0.1, -0.73, 0.8, 0.8), + fancybox=False, shadow=False, ncol=1) # bbox_to_anchor=(0.1, -0.73, 0.8, 0.8), + + fig.supxlabel(' Time') + fig.supylabel('Number of persons') + plt.subplots_adjust(left=None, bottom=None, right=None, + top=None, wspace=None, hspace=0.6) + # plt.tight_layout(pad=0, w_pad=0.5, h_pad=0.3) + + # save result + if save: + if not os.path.isdir(save_dir): + os.makedirs(save_dir) + plt.savefig(save_dir + f'ide_ode_{fileending}.png', + bbox_inches='tight', dpi=500) + + plt.close() + + # # plot legend separately + # figsize = (2, 0.5) + # fig_leg = plt.figure(figsize=figsize) + # ax_leg = fig_leg.add_subplot(111) + # # add the legend from the previous axes + # ax_leg.legend(*axs[0, 0].get_legend_handles_labels(), loc='center', ncol=2) + # # hide the axes frame and the x/y labels + # ax_leg.axis('off') + # fig_leg.savefig('plots/legend.png', dpi=500) + + +if __name__ == '__main__': + # Path to simulation results + + dt_ode = '1e-1' + dt_ide = '1e-1' + + setting = 2 + + flows = True + + legendplot = list(["ODE", "IDE"]) + + # Plot comparison of ODE and IDE models + + # -------------- Simulations based on LEOSS data --------------- + + # data_dir = os.path.join(os.path.dirname( + # __file__), "../results/fictional/leoss/") + + # compare_results([os.path.join(data_dir, f"fictional_ode_leoss_2.0_12_0.1000_flows"), + # os.path.join(data_dir, f"fictional_ide_leoss_2.0_12_0.1000_flows")], + # legendplot, flows=True, fileending="2.0_12_0.1000_flows", save=True, save_dir='plots/leoss/flows/') + + # # # # # # # # # Simulations based on Covasim data + # data_dir = os.path.join(os.path.dirname( + # __file__), "../results/fictional/covasim/") + + # compare_results([os.path.join(data_dir, f"fictional_ode_covasim_0.5_12_0.1000_flows"), + # os.path.join(data_dir, f"fictional_ide_covasim_0.5_12_0.1000_flows")], + # legendplot, flows=True, fileending="0.5_12_0.1000_flows", save=True, save_dir='plots/covasim/flows/') + + # compare_results([os.path.join(data_dir, f"fictional_ode_covasim_0.0_12_0.1000_flows"), + # os.path.join(data_dir, f"fictional_ide_covasim_0.0_12_0.1000_flows")], + # legendplot, flows=True, fileending="0.0_12_0.1000_flows", save=True, save_dir='plots/covasim/flows/') + + # compare_results([os.path.join(data_dir, f"fictional_ode_covasim_0.0_12_0.1000_compartments"), + # os.path.join(data_dir, f"fictional_ide_covasim_0.0_12_0.1000_compartments")], + # legendplot, flows=False, fileending="0.0_12_0.1000_compartments", save=True, save_dir='plots/covasim/compartments/') + + # # # # Real scenario + data_dir = os.path.join(os.path.dirname( + __file__), "../results/real/") + + # compare_results([os.path.join(data_dir, f"ode_2020-11-01_30_0.1000_flows"), + # os.path.join(data_dir, f"ide_2020-11-01_30_0.1000_flows")], + # legendplot, flows=True, fileending="2020-11-01_30_0.1000_flows", save=True, save_dir='plots/real/') + + # compare_results([os.path.join(data_dir, f"ode_2020-11-01_30_0.1000_compartments"), + # os.path.join(data_dir, f"ide_2020-11-01_30_0.1000_compartments")], + # legendplot, flows=False, fileending="2020-11-01_30_0.1000_compartments", save=True, save_dir='plots/real/') + + start_date = '2020-10-1' + simulation_time = 45 + timestep = '0.1000' + compare_results([os.path.join(data_dir, f"ode_{start_date}_{simulation_time}_{timestep}_flows"), + os.path.join(data_dir, f"ide_{start_date}_{simulation_time}_{timestep}_flows")], + legendplot, flows=True, fileending=f"{start_date}_{simulation_time}_{timestep}_flows", save=True, save_dir=f"plots/real/{start_date}/{simulation_time}/") + + compare_results([os.path.join(data_dir, f"ode_{start_date}_{simulation_time}_{timestep}_compartments"), + os.path.join(data_dir, f"ide_{start_date}_{simulation_time}_{timestep}_compartments")], + legendplot, flows=False, fileending=f"{start_date}_{simulation_time}_{timestep}_compartments", save=True, save_dir=f'plots/real/{start_date}/{simulation_time}/') + + # # # # # Constant init + # data_dir = os.path.join(os.path.dirname( + # __file__), "../results/fictional/constant_init/") + + # compare_results([os.path.join(data_dir, f"ode_constant_init_30_0.1000_flows"), + # os.path.join(data_dir, f"ide_constant_init_30_0.1000_flows")], + # legendplot, flows=True, fileending="30_0.1000_flows", save=True, save_dir='plots/constant_init/') + + # compare_results([os.path.join(data_dir, f"ode_constant_init_30_0.1000_compartments"), + # os.path.join(data_dir, f"ide_constant_init_30_0.1000_compartments")], + # legendplot, flows=False, fileending="30_0.1000_compartments", save=True, save_dir='plots/constant_init/')