documentation

cpp/simulations/IDE_paper/ide_convergence_rate.cpp

@@ -505,7 +505,7 @@ mio::IOResult<void> simulate_ode_and_ide(ScalarType t0, ScalarType tmax, ScalarT
                                  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";
+                std::cout << "Successfully saved the IDE simulation results. \n\n";
             }
             else {
                 std::cout << "Error occured while saving the IDE simulation results. \n";
@@ -521,21 +521,23 @@ 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/";
-    // 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;
-    // Decide if we want to run the IDE simulation or just the ODE simulation, e.g., to create a ground truth.
-    bool ide_simulation = true;
 
     // 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 = 1;
-    // The IDE model will be simulated using a fixed step size dt=10^{-ide_exponent}.
-    ScalarType ide_exponent = 1;
+    ScalarType ode_exponent = 6;
+    // The IDE model will be simulated using a fixed step size dt=10^{-ide_exponent} for ide_exponent in ide_exponents.
+    std::vector<ScalarType> ide_exponents = {1, 2, 3, 4};
+    // 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;
+    // Decide if we want to run the IDE simulation or just the ODE simulation, e.g., to create a ground truth.
+    bool ide_simulation = true;
 
-    auto result = simulate_ode_and_ide(t0, tmax, ode_exponent, ide_exponent, save_exponent, result_dir, ide_simulation);
+    mio::IOResult<void> result = mio::success();
+    for (ScalarType ide_exponent : ide_exponents)
+        result = simulate_ode_and_ide(t0, tmax, ode_exponent, ide_exponent, save_exponent, result_dir, ide_simulation);
 
     if (!result) {
         printf("%s\n", result.error().formatted_message().c_str());

cpp/simulations/IDE_paper/ide_covid_scenario.cpp

@@ -45,7 +45,7 @@ using Vector = Eigen::Matrix<ScalarType, Eigen::Dynamic, 1>;
 
 // Used parameters.
 std::map<std::string, ScalarType> simulation_parameter = {{"t0", 0.},
-                                                          {"dt", 0.1},
+                                                          {"dt", 0.01},
                                                           {"total_population", 83155031.},
                                                           {"total_confirmed_cases", 0.}, // set by RKI data
                                                           {"deaths", 0.}, // set by RKI data

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

cpp/simulations/IDE_paper/plot_changepoint.py

@@ -26,16 +26,24 @@
 from memilio.epidata import getDataIntoPandasDataFrame as gd
 
 
-def plot_changepoint(files, legendplot, flows=True, fileending="", save=True, save_dir='plots/'):
+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"])
 
-    # helmholtzdarkblue, helmholtzclaim
+    # Define colors, we use helmholtzdarkblue, helmholtzclaim.
     colors = [(0, 40/255, 100/255), (20/255, 200/255, 255/255)]
     linestyles = ['-', '--']
-    # add results to plot
+    # Add results to plot-
     for file in range(len(files)):
-        # load data
+        # Load data.
         h5file = h5py.File(str(files[file]) + '.h5', 'r')
 
         if (len(list(h5file.keys())) > 1):
@@ -45,53 +53,32 @@ def plot_changepoint(files, legendplot, flows=True, fileending="", save=True, sa
 
         data = h5file[list(h5file.keys())[0]]
 
-        if flows:
-            # As there should be only one Group, total is the simulation result
-            total = data['Total'][:, :]
-        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]]
+
+        # 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
+        # Get indices where dates are >=0.
         indices = np.where(dates >= 0)
-        # plot data
-        if flows:
-            # ODE
-            if file == 0:
-                # transform cumulative flows to flows 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])
-
-                date_idx = -int(tmax/timestep)-2
-                print(
-                    f"New infections at {dates[date_idx]}: {total[date_idx,0]}")
-                date_idx = -int(tmax/timestep)-1
-                print(
-                    f"New infections at {dates[date_idx]}: {total[date_idx,0]}")
-                date_idx = -int(tmax/timestep)
-                print(
-                    f"New infections at {dates[date_idx]}: {total[date_idx,0]}")
-        else:
-            incidence = (total[:-1, 0]-total[1:, 0])/(dates[1:]-dates[:-1])
-            ax.plot(dates[indices[0][1:]], incidence, label=legendplot[file],
+
+        # 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)
@@ -105,24 +92,27 @@ def plot_changepoint(files, legendplot, flows=True, fileending="", save=True, sa
 
     plt.tight_layout()
 
-    # save result
-    if save:
+    # 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__':
-    legendplot = list(["ODE", "IDE"])
-    # Path to simulation results
-    data_dir = os.path.join(os.path.dirname(
-        __file__), "..", "results/fictional/covasim/")
-
-    plot_changepoint([os.path.join(data_dir, f"fictional_ode_covasim_0.5_12_0.0100_flows"),
-                     os.path.join(data_dir, f"fictional_ide_covasim_0.5_12_0.0100_flows")],
-                     legendplot, flows=True, fileending="0.5_12_0.0100", save=True, save_dir='plots/covasim/changepoints/')
-
-    plot_changepoint([os.path.join(data_dir, f"fictional_ode_covasim_2.0_12_0.0100_flows"),
-                     os.path.join(data_dir, f"fictional_ide_covasim_2.0_12_0.0100_flows")],
-                     legendplot, flows=True, fileending="2.0_12_0.0100", save=True, save_dir='plots/covasim/changepoints/')
+    # 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)