diff --git a/cpp/simulations/IDE_paper/ide_convergence_rate.cpp b/cpp/simulations/IDE_paper/ide_convergence_rate.cpp index 3cb093eb4c..286ee7057e 100644 --- a/cpp/simulations/IDE_paper/ide_convergence_rate.cpp +++ b/cpp/simulations/IDE_paper/ide_convergence_rate.cpp @@ -505,7 +505,7 @@ mio::IOResult 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 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 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()); diff --git a/cpp/simulations/IDE_paper/ide_covid_scenario.cpp b/cpp/simulations/IDE_paper/ide_covid_scenario.cpp index 2c320cfc81..c8bbdc8403 100644 --- a/cpp/simulations/IDE_paper/ide_covid_scenario.cpp +++ b/cpp/simulations/IDE_paper/ide_covid_scenario.cpp @@ -45,7 +45,7 @@ using Vector = Eigen::Matrix; // Used parameters. std::map 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 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 index 238da87642..c26e3f026e 100644 --- a/cpp/simulations/IDE_paper/plot_changepoint.py +++ b/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,8 +92,8 @@ 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", @@ -114,15 +101,18 @@ def plot_changepoint(files, legendplot, flows=True, fileending="", save=True, sa 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) diff --git a/tools/get_lognormal_parameters.py b/tools/get_lognormal_parameters.py deleted file mode 100644 index f922184934..0000000000 --- a/tools/get_lognormal_parameters.py +++ /dev/null @@ -1,59 +0,0 @@ -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) - - # print("Test mean: ", mean_test) - - # print(mean_lognorm, variance_lognorm) - - 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/tools/investigate_age_distribution.py b/tools/investigate_age_distribution.py deleted file mode 100644 index 01a96358f5..0000000000 --- a/tools/investigate_age_distribution.py +++ /dev/null @@ -1,273 +0,0 @@ -import os -import pandas as pd -import matplotlib.pyplot as plt -import numpy as np - - -def get_df_daily(): - # Read file. - datafile = os.path.join(os.path.dirname( - __file__), "..", "data", "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_population_per_agegroup(): - # 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(start_date, T_IH, T_HU, T_U): - # Get dataframe with daily confirmed cases. - df = get_df_daily() - - # 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(start_dates, T_IH, T_HU, T_U): - - 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( - start_dates[0], T_IH, T_HU, T_U) - proportions_october = get_relevant_confirmed_cases( - start_dates[1], T_IH, T_HU, T_U) - population_per_agegroup = get_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() - - plt.savefig(f"plots/proportions_per_agegroup.png", dpi=500) - - -def covasim_to_rki_agegroups(mu_covasim): - mu_rki = np.zeros(6) - # Covasim age groups are split every 10 years. - # A00-A04 - mu_rki[0] = (0.5*mu_covasim[0])/0.5 - # A05-A14 - mu_rki[1] = (0.5*mu_covasim[0] + 0.5*mu_covasim[1])/1. - # A15-A34 - mu_rki[2] = (0.5*mu_covasim[1] + mu_covasim[2] + 0.5 * mu_covasim[3])/2. - # A35-A59 - mu_rki[3] = (0.5 * mu_covasim[3] + mu_covasim[4] + mu_covasim[5])/2.5 - # A60-A79 - mu_rki[4] = (mu_covasim[6] + mu_covasim[7])/2. - # A80+ - mu_rki[5] = (mu_covasim[8] + mu_covasim[9])/2. - - return mu_rki - - -def compute_covasim_probs_per_rki_agegroup(): - # Vectors with probabilities are directly from Covasim Paper. These are adjusted accordingly to our probabilities by division. - # In Covasim, there are probabilities given from C->I, C->H, C->U and C->D. - # We calculate our needed transition probabilities by dividing accordingly. - - mu_CI_covasim = np.array( - [0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.90]) - - mu_CH_covasim = np.array([0.00050, 0.00165, 0.00720, 0.02080, 0.03430, - 0.07650, 0.13280, 0.20655, 0.24570, 0.24570]) - mu_IH_covasim = mu_CH_covasim/mu_CI_covasim - - mu_CU_covasim = np.array([0.00003, 0.00008, 0.00036, 0.00104, 0.00216, - 0.00933, 0.03639, 0.08923, 0.17420, 0.17420]) - mu_HU_covasim = mu_CU_covasim/(mu_CH_covasim) - - mu_CD_covasim = np.array([0.00002, 0.00002, 0.00010, 0.00032, 0.00098, - 0.00265, 0.00766, 0.02439, 0.08292, 0.16190]) - mu_UD_covasim = mu_CD_covasim/mu_CU_covasim - - # # Compute average by population just to test. - # agegroups_covasim = np.array([7752706.0, 7581868, 9483430, 10871964, 10070748, - # 13304542, 10717241, 7436098, 5092743, 843691]) - # total_pop = agegroups_covasim.sum() - - # mu_CI_average = 0 - # mu_IH_average = 0 - # mu_HU_average = 0 - # mu_UD_average = 0 - # for i in range(len(agegroups_covasim)): - # mu_CI_average += mu_CI_covasim[i]*agegroups_covasim[i]/total_pop - # mu_IH_average += mu_IH_covasim[i]*agegroups_covasim[i]/total_pop - # mu_HU_average += mu_HU_covasim[i]*agegroups_covasim[i]/total_pop - # mu_UD_average += mu_UD_covasim[i]*agegroups_covasim[i]/total_pop - - # print("Covasim probs by pop: ", mu_CI_average, - # mu_IH_average, mu_HU_average, mu_UD_average) - - # Convert from 10 agegroups from Covasim Paper to 6 age groups according to RKI data. - mu_CI_rki = covasim_to_rki_agegroups(mu_CI_covasim) - mu_IH_rki = covasim_to_rki_agegroups(mu_IH_covasim) - mu_HU_rki = covasim_to_rki_agegroups(mu_HU_covasim) - mu_UD_rki = covasim_to_rki_agegroups(mu_UD_covasim) - - return mu_CI_rki, mu_IH_rki, mu_HU_rki, mu_UD_rki - - -def compute_adapted_mu(start_date, T_IH, T_HU, T_U): - mu_CI_age, mu_IH_age, mu_HU_age, mu_UD_age = compute_covasim_probs_per_rki_agegroup() - - population_share = get_relevant_confirmed_cases( - start_date, T_IH, T_HU, T_U) - - mu_CI = 0 - mu_IH = 0 - mu_HU = 0 - mu_UD = 0 - for i in range(len(population_share)): - mu_CI += mu_CI_age[i] * population_share[i] - mu_IH += mu_IH_age[i] * population_share[i] - mu_HU += mu_HU_age[i] * population_share[i] - mu_UD += mu_UD_age[i] * population_share[i] - - return mu_CI, mu_IH, mu_HU, mu_UD - - -def compute_mu_by_population(): - mu_CI_age, mu_IH_age, mu_HU_age, mu_UD_age = compute_covasim_probs_per_rki_agegroup() - population_per_agegroup = get_population_per_agegroup() - - mu_CI = 0 - mu_IH = 0 - mu_HU = 0 - mu_UD = 0 - for i in range(len(population_per_agegroup)): - mu_CI += mu_CI_age[i] * population_per_agegroup[i] - mu_IH += mu_IH_age[i] * population_per_agegroup[i] - mu_HU += mu_HU_age[i] * population_per_agegroup[i] - mu_UD += mu_UD_age[i] * population_per_agegroup[i] - - return mu_CI, mu_IH, mu_HU, mu_UD - - -def mu_assessment_by_cases(start_date, T_IH, T_HU): - population_share = get_relevant_confirmed_cases( - start_date, T_IH, T_HU) - - mu_CI_assessment = np.array([0.75, 0.75, 0.8, 0.8, 0.8, 0.8]) - mu_IH_assessment = np.array([0.0075, 0.0075, 0.019, 0.0615, 0.165, 0.225]) - mu_HU_assessment = np.array([0.075, 0.075, 0.15, 0.15, 0.3, 0.4]) - mu_UD_assessment = np.array([0.05, 0.05, 0.14, 0.14, 0.4, 0.6]) - - mu_CI = 0 - mu_IH = 0 - mu_HU = 0 - mu_UD = 0 - for i in range(len(population_share)): - mu_CI += mu_CI_assessment[i] * population_share[i] - mu_IH += mu_IH_assessment[i] * population_share[i] - mu_HU += mu_HU_assessment[i] * population_share[i] - mu_UD += mu_UD_assessment[i] * population_share[i] - - return mu_CI, mu_IH, mu_HU, mu_UD - - -def main(): - - df_daily = get_df_daily() - - T_IH = 6.6 - T_HU = 1.5 - T_U = 15.230258 - - # # from Assessment Paper - # mu_IH = 0.0786429 - # mu_HU = 0.173176 - - # # from Covasim - # mu_IH = 0.104907 - # mu_HU = 0.369201 - - start_dates = ["2020-06-01", "2020-10-01"] - - plot(start_dates, T_IH, T_HU, T_U) - - mu_CI, mu_IH, mu_HU, mu_UD = compute_adapted_mu( - start_dates[0], T_IH, T_HU, T_U) - print(f"mu {start_dates[0]}: {mu_CI}, {mu_IH}, {mu_HU}, {mu_UD}") - - mu_CI, mu_IH, mu_HU, mu_UD = compute_adapted_mu( - start_dates[1], T_IH, T_HU, T_U) - print(f"mu {start_dates[1]}: {mu_CI}, {mu_IH}, {mu_HU}, {mu_UD}") - - mu_CI, mu_IH, mu_HU, mu_UD = compute_mu_by_population() - print(f"mu by population: {mu_CI}, {mu_IH}, {mu_HU}, {mu_UD}") - - # mu_CI, mu_IH, mu_HU, mu_UD = mu_assessment_by_cases(start_dates[0], T_IH, T_HU) - # print(f"Assessment mu by cases June: {mu_CI}, {mu_IH}, {mu_HU}, {mu_UD}") - - # mu_CI, mu_IH, mu_HU, mu_UD = mu_assessment_by_cases(start_dates[1], T_IH, T_HU) - # print(f"Assessment mu by cases October: {mu_CI}, {mu_IH}, {mu_HU}, {mu_UD}") - - -if __name__ == "__main__": - main() diff --git a/tools/plot_changepoint.py b/tools/plot_changepoint.py deleted file mode 100644 index 2946c6c9c4..0000000000 --- a/tools/plot_changepoint.py +++ /dev/null @@ -1,111 +0,0 @@ -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, legendplot, flows=True, fileending="", save=True, save_dir='plots/'): - - fig, ax = plt.subplots() - - # 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 - 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]] - - dates = data['Time'][:] - - timestep = np.diff(dates)[0] - tmax = dates[-1] - - # 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], - color=colors[file], linestyle=linestyles[file]) - - h5file.close() - - # ax.set_title(secir_dict[i], fontsize=8) - # axs[int(i/2), i % 2].set_ylim(bottom=0) - 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: - 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/') diff --git a/tools/plot_conv_paper.py b/tools/plot_conv_paper.py deleted file mode 100644 index 6cb2dcc08e..0000000000 --- a/tools/plot_conv_paper.py +++ /dev/null @@ -1,323 +0,0 @@ -############################################################################# -# 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, exponent_ode, flows=False): - """ Read groundtruth from data. We define the groundtruth as the results obtained by the ODE model with timestep dt=1e-6.""" - - model = 'ode' - results = {model: []} - savefreq = 4 - if flows: - h5file = h5py.File(os.path.join( - data_dir, f'result_{model}_flows_dt=1e-{exponent_ode:.0f}_savefrequency{savefreq:.0f}.h5'), 'r') - else: - h5file = h5py.File(os.path.join( - data_dir, f'result_{model}_dt=1e-{exponent_ode:.0f}_savefrequency{savefreq:.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, exponent_ode, 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.""" - - 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-{exponent_ode:.0f}.h5'), 'r') - else: - h5file = h5py.File(os.path.join( - data_dir, f'result_{model}_dt=1e-{exponent:.0f}_init_dt_ode=1e-{exponent_ode:.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): - """ Compute norm of a time series.""" - - norm = np.sqrt(timestep * np.sum(timeseries**2)) - return norm - - -def compute_relerror_norm_l2(groundtruth, results, timestep_ode, timesteps_ide, flows=False): - """ Compute norm of the difference between time series from ODE and time series - from IDE for all compartments/flows.""" - - 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/timestep_ode - 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, compartment=None, save=False): - """ Plot errors against timesteps.""" - - # 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: - if flows: - if not os.path.isdir('plots/flows'): - os.makedirs('plots/flows') - plt.savefig(f'plots/flows/convergence_all_flows.png', format='png', bbox_extra_artists=(lgd,), bbox_inches='tight', - dpi=500) - else: - if not os.path.isdir('plots/compartments'): - os.makedirs('plots/compartments') - plt.savefig(f'plots/compartments/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=False): - """ Plot errors against timesteps. This function creates one plot to depict all compartments/flows, respectively.""" - - 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: - if flows: - if not os.path.isdir('plots/flows'): - os.makedirs('plots/flows') - plt.savefig(f'plots/flows/convergence_flows.png', format='png', - dpi=500) - else: - if not os.path.isdir('plots/compartments'): - os.makedirs('plots/compartments') - plt.savefig(f'plots/compartments/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.""" - 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(): - data_dir = os.path.join(os.path.dirname( - __file__), "..", "results") - - exponents_ide = [1, 2, 3, 4] - timesteps_ide = [] - for x in exponents_ide: - timesteps_ide.append(pow(10, -x)) - - flows = True - - groundtruth = read_groundtruth(data_dir, 6, flows=flows) - timestep_ode = 1e-4 - - results = read_data(data_dir, 6, exponents_ide, flows=flows) - - relerrors_l2 = compute_relerror_norm_l2( - groundtruth, results, timestep_ode, timesteps_ide, flows=flows) - - plot_convergence_oneplot( - relerrors_l2, timesteps_ide, flows=flows, save=True) - plot_convergence(relerrors_l2, timesteps_ide, flows=flows, save=True) - - """order = compute_order_of_convergence( - relerrors_l2, timesteps_ide, flows=flows) - - print('Orders of convergence: ', order)""" - - -if __name__ == '__main__': - main() diff --git a/tools/plot_real_scenario.py b/tools/plot_real_scenario.py deleted file mode 100644 index 0f47933611..0000000000 --- a/tools/plot_real_scenario.py +++ /dev/null @@ -1,514 +0,0 @@ -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, used for plotting real data. -# Probabilities from Assessment paper -parameters = { - 'TimeExposed': 4.5, - 'TimeInfectedNoSymptoms': 2.527617, # Covasim: 3.18163, ## Assessment: 2.52762 - 'TimeInfectedSymptoms': 7.889900, # 7.85313, ' 7.8899 - 'TimeInfectedSevere': 15.225278, # 11.9713, # 15.2253 - 'TimeInfectedNoSymptomsToInfectedSymptoms': 1.1, - 'TimeInfectedSymptomsToInfectedSevere': 6.6, - 'TimeInfectedSymptomsToRecovered': 8.0, - 'TimeInfectedSevereToInfectedCritical': 1.5, - 'TimeInfectedCriticalToDead': 10.7, - 'InfectedSymptomsPerInfectedNoSymptoms': 0.793099, # 0.698315 #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, mu_IH): - parameters["TimeInfectedCriticalToDead"] = T_UD - parameters["SeverePerInfectedSymptoms"] = mu_IH - return parameters - - -def load_data(file, start_date, simulation_time, T_UD, mu_IH): - """ 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, T_UD, mu_IH): - - datafile = os.path.join(os.path.dirname( - __file__), "..", "data", "pydata", "Germany", "cases_all_germany.json") - data_rki = load_data(datafile, start_date, simulation_time, T_UD, mu_IH) - - # datafile_ma7 = os.path.join(os.path.dirname( - # __file__), "..", "data", "pydata", "Germany", "cases_all_germany_ma7.json") - # data_rki = load_data(datafile_ma7, start_date, simulation_time, T_UD) - - # 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, T_UD, mu_IH, 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, T_UD, mu_IH) - - # 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, T_UD) - - 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") - # ax.plot(np.linspace(0, simulation_time, simulation_time + 1), - # data_rki_ma7["NewInfectionsDay"], color='gray', label="Extrapolated RKI data MA7") - - 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, T_UD, mu_IH, 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, T_UD, mu_IH) - - 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, T_UD, mu_IH) - - # print("Infectedsymptoms from RKI (ma7) on 1.10.2020: ", data_rki_ma7[data_rki_ma7["Date"]=="2020-10-01"]["InfectedSymptoms"].values[0]) - - 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/") # different_UD/ - - start_date = '2020-10-1' - simulation_time = 45 - timestep = "0.1000" - probs = "assessment_probs" - # probs = "covasim_probs" - # probs = "different_UD" - - # if probs == "assessment_probs": - # parameters = parameters_assessment_probs - # elif probs == "different_UD": - # parameters = parameters_assessment_probs_differentUD - # else: - # parameters = parameters_covasim_probs - - T_UD = parameters["TimeInfectedCriticalToDead"] - - 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, T_UD, - fileending=f"{start_date}_{simulation_time}_{timestep}", save=True, save_dir=f"plots/real/{start_date}/{simulation_time}/{probs}/") - - 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, T_UD, - fileending=f"{start_date}_{simulation_time}_{timestep}", save=True, save_dir=f"plots/real/{start_date}/{simulation_time}/{probs}/") - - 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/real/{start_date}/{simulation_time}/{probs}/') diff --git a/tools/run_and_plot_scenario.py b/tools/run_and_plot_scenario.py deleted file mode 100644 index 948919b1ae..0000000000 --- a/tools/run_and_plot_scenario.py +++ /dev/null @@ -1,189 +0,0 @@ -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 -from investigate_age_distribution import compute_adapted_mu, mu_assessment_by_cases - - -def run_real_scenario(data_dir, save_dir, start_date, simulation_time, timestep, mu_IH, mu_HU, mu_UD, T_UD, T_UR, std_UD, std_UR, scale_contacts, scale_confirmed_cases): - - year = start_date.split("-")[0] - month = start_date.split("-")[1] - day = start_date.split("-")[2] - - shape_UD, scale_UD = get_lognormal_parameters(T_UD, std_UD) - shape_UR, scale_UR = get_lognormal_parameters(T_UR, std_UR) - - subprocess.call([f"./build/bin/ide_real_scenario", data_dir, save_dir, f"{mu_UD}", f"{T_UD}", f"{T_UR}", f"{shape_UD}", f"{scale_UD}", f"{shape_UR}", f"{scale_UR}", - f"{year}", f"{month}", f"{day}", f"{simulation_time}", f"{timestep}", f"{scale_contacts}", f"{mu_IH}", f"{mu_HU}"]) - - -def contact_scaling(save_dir, start_date, simulation_time, timestep, T_UD, mu_IH): - scale_contacts = get_scale_contacts([os.path.join( - save_dir, f"ide_{start_date}_{simulation_time}_{timestep}_flows")], pd.Timestamp(start_date), simulation_time, T_UD, mu_IH) - - 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, T_UD, mu_IH): - 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, T_UD, mu_IH, - 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, T_UD, mu_IH, - 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, mu_IH, mu_HU, mu_UD, T_UD, T_UR, std_UD, std_UR): - # 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, mu_IH, mu_HU, mu_UD, T_UD, T_UR, std_UD, std_UR, scale_contacts, scale_confirmed_cases) - # # 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, T_UD, mu_IH) - # # 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, mu_IH, mu_HU, mu_UD, T_UD, T_UR, std_UD, std_UR, scale_contacts, scale_confirmed_cases) - plot_real_scenario(save_dir, plot_dir, start_date, - simulation_time, timestep, T_UD, mu_IH) - - -def october_scenario(save_folder, timestep, mu_IH, mu_HU, mu_UD, T_UD, std_UD, T_UR, std_UR): - start_date = '2020-10-1' - simulation_time = 45 - # timestep = "0.1000" - - probs = [mu_UD] - T_UDs = [T_UD] - T_URs = [T_UR] - - # # Try differen parameters regarding U. - # # Assessment - # mu_UD_assessment = 0.217177 - # # Covasim - # mu_UD_covasim = 0.387803 - # probs = [mu_UD_covasim, mu_UD_assessment] # mu_UD_assessment, - # T_UD = 10.7 - # # T_UDs = [10.7] - # T_UDs = [T_UD-i*0.5 for i in range(5)] - # std_UD = 4.8 - # T_UR = 18.1 - # # T_URs = [18.1] - # T_URs = [T_UR-i*0.5 for i in range(5)] - # std_UR = 6.3 - - for T_UD in T_UDs: - for T_UR in T_URs: - for mu_UD in probs: - print(f"Parameters: {mu_IH, mu_HU, mu_UD, T_UD, T_UR}") - data_dir = "./data" - save_dir = f"./results/real/{save_folder}/{start_date}/{mu_IH}_{mu_HU}_{mu_UD}_{T_UD}_{T_UR}/" - plot_dir = f"plots/real/{save_folder}/{start_date}/{mu_IH}_{mu_HU}_{mu_UD}_{T_UD}_{T_UR}/" - - run_scenario(data_dir, save_dir, plot_dir, start_date, simulation_time, - timestep, mu_IH, mu_HU, mu_UD, T_UD, T_UR, std_UD, std_UR) - - -def june_scenario(save_folder, timestep, mu_IH, mu_HU, mu_UD, T_UD, std_UD, T_UR, std_UR): - start_date = '2020-6-1' - simulation_time = 45 - # timestep = "0.1000" - - data_dir = "./data" - save_dir = f"./results/real/{save_folder}/{start_date}/{mu_IH}_{mu_HU}_{mu_UD}_{T_UD}_{T_UR}/" - plot_dir = f"plots/real/{save_folder}/{start_date}/{mu_IH}_{mu_HU}_{mu_UD}_{T_UD}_{T_UR}/" - - run_scenario(data_dir, save_dir, plot_dir, start_date, simulation_time, - timestep, mu_IH, mu_HU, mu_UD, T_UD, T_UR, std_UD, std_UR) - - -def main(): - # Folder in plots/real/save_folder where plots will be stored. - - T_IH = 6.6 - T_HU = 1.5 - T_UD = 10.7 - std_UD = 4.8 - T_UR = 18.1 - std_UR = 6.3 - - # Always use mu_CI=0.793099 from Assessment paper - - # save_folder = "mu_weighed_by_cases" - # mu_CI_june, mu_IH_june, mu_HU_june, mu_UD_june = compute_adapted_mu("2020-06-01", T_IH, T_HU) - # june_scenario(save_folder, mu_IH, mu_HU_june, mu_UD_june, T_UD, std_UD, T_UR, std_UR) - - # mu_CI_october, mu_IH_october, mu_HU_october, mu_UD_october = compute_adapted_mu("2020-10-01", T_IH, T_HU) - # october_scenario(save_folder, mu_IH, mu_HU_october, mu_UD_october, T_UD, std_UD, T_UR, std_UR) - - # save_folder = "new_covasim_probs" - # mu_CI_june, mu_IH_june, mu_HU_june, mu_UD_june = compute_adapted_mu( - # "2020-06-01", T_IH, T_HU) - # june_scenario(save_folder, mu_IH_june, mu_HU_june, - # mu_UD_june, T_UD, std_UD, T_UR, std_UR) - - # mu_CI_october, mu_IH_october, mu_HU_october, mu_UD_october = compute_adapted_mu( - # "2020-10-01", T_IH, T_HU) - # october_scenario(save_folder, mu_IH_october, mu_HU_october, - # mu_UD_october, T_UD, std_UD, T_UR, std_UR) - - # save_folder = "assessment_by_cases" - # mu_CI_june, mu_IH_june, mu_HU_june, mu_UD_june =mu_assessment_by_cases("2020-06-01", T_IH, T_HU) - # june_scenario(save_folder, mu_IH_june, mu_HU_june, - # mu_UD_june, T_UD, std_UD, T_UR, std_UR) - - # mu_CI_october, mu_IH_october, mu_HU_october, mu_UD_october = mu_assessment_by_cases("2020-10-01", T_IH, T_HU) - # october_scenario(save_folder, mu_IH_october, mu_HU_october, - # mu_UD_october, T_UD, std_UD, T_UR, std_UR) - - # save_folder = "average_first_by_cases_only_UD" - # # maximum - # mu_CI_june=0.793099 - # mu_IH_june= 0.0786429 - # mu_HU_june= 0.173176 - # mu_UD_june = 0.5475428440089082 - # # mu_CI_june, mu_IH_june, mu_HU_june, mu_UD_june =mu_assessment_by_cases("2020-06-01", T_IH, T_HU) - # june_scenario(save_folder, mu_IH_june, mu_HU_june, - # mu_UD_june, T_UD, std_UD, T_UR, std_UR) - - # mu_CI_october=0.793099 - # mu_IH_october= 0.0786429 - # mu_HU_october= 0.173176 - # mu_UD_october = 0.4376045329412043 - # # mu_CI_october, mu_IH_october, mu_HU_october, mu_UD_october = mu_assessment_by_cases("2020-10-01", T_IH, T_HU) - # october_scenario(save_folder, mu_IH_october, mu_HU_october, - # mu_UD_october, T_UD, std_UD, T_UR, std_UR) - - save_folder = "paper_dt=0.01" - timestep = "0.0100" - mu_IH_october = 0.078643 - mu_HU_october = 0.173176 - mu_UD_october = 0.36245545116366784 - october_scenario(save_folder, timestep, mu_IH_october, mu_HU_october, - mu_UD_october, T_UD, std_UD, T_UR, std_UR) - - -if __name__ == "__main__": - - main()