1036 Simulations from paper to investigate behavior of IDE model

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)

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()

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 <typename DataEntry>
     friend IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const& path, Date date,
                                             ScalarType scale_confirmed_cases);
 };

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

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()

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 <[email protected]>
+*
+* 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 <iostream>
+
+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;
+}

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 <[email protected]>
+#
+# 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}")