From 840303dde091a6544a6609a1e23a1cc4a95bf297 Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Tue, 27 May 2025 21:02:44 +0200 Subject: [PATCH 01/38] Add scripts for analyzing infection states and ICU data from simulation runs --- .../memilio/plot/plotAbmICUAndDeadComp.py | 268 ++++++++++++++++++ .../memilio/plot/plotAbmInfectionStates.py | 258 +++++++++++++++++ 2 files changed, 526 insertions(+) create mode 100644 pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py create mode 100644 pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py diff --git a/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py b/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py new file mode 100644 index 0000000000..6a1af97ed5 --- /dev/null +++ b/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py @@ -0,0 +1,268 @@ +# Python script to analyze bs runs +# input is a bs run folder with the following structure: +# bs_run_folder has a txt file for each bs run +# each txt file has a line for each time step +# each line has a column for each compartment as well as the timestep +# each column has the number of individuals in that compartment +# the first line of each txt file is the header + +import sys +import argparse +import os +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +import matplotlib.colors as colors +import matplotlib.cm as cmx +import matplotlib.patches as mpatches +import matplotlib.lines as mlines +import h5py +from datetime import datetime +from matplotlib.dates import DateFormatter +from scipy.ndimage import gaussian_filter1d +from scipy.signal import savgol_filter + +fontsize = 20 + + +def plot_dead(path): + # we will have a seperate plot the cumulative infected individuals, cumulative symptomatic individuals and cumulative dead individual + # we need to load the data + f_p50 = h5py.File( + path+"/infection_state_per_age_group/0/p50/Results.h5", 'r') + p50_bs = f_p50['0'] + + # do the same for 25 and 75 percentile + f_p25 = h5py.File( + path+"/infection_state_per_age_group/0/p25/Results.h5", 'r') + p25_bs = f_p25['0'] + + f_p75 = h5py.File( + path+"/infection_state_per_age_group/0/p75/Results.h5", 'r') + p75_bs = f_p75['0'] + + # do the same for 05 and 95 percentile + f_p05 = h5py.File( + path+"/infection_state_per_age_group/0/p05/Results.h5", 'r') + p05_bs = f_p05['0'] + + f_p95 = h5py.File( + path+"/infection_state_per_age_group/0/p95/Results.h5", 'r') + p95_bs = f_p95['0'] + + age_group_access = ['Group1', 'Group2', 'Group3', + 'Group4', 'Group5', 'Group6', 'Total'] + + # we need the real data json file cases_all_county_age + df_abb = pd.read_json( + path+"/../../../pydata/Germany/cases_all_county_age_ma1.json") + + # we just need the columns cases and date + # we need to offset the dates by 19 day + df_abb['Date'] = df_abb['Date'] + pd.DateOffset(days=18) + # we need just the dates bewteen 2021-03-01 and 2021-06-01 + df_abb = df_abb[(df_abb['Date'] >= '2021-03-01') & + (df_abb['Date'] <= '2021-06-01')] + # we just need the cases with id 3101 + df_abb = df_abb[df_abb['ID_County'] == 3101] + # df_abb['Deaths'] = np.round(df_abb[['Deaths']].to_numpy()) + + # we need the amount of dead persons for each age group: These are A00-A04, A05-A14, A15-A34, A35-A59, A60-A79, A80+ + age_groups = ['A00-A04', 'A05-A14', 'A15-A34', 'A35-A59', 'A60-A79', 'A80+'] + age_grous_string = ['Age 0-4', 'Age 5-14', 'Age 15-34', 'Age 35-59', 'Age 60-79', 'Age 80+'] + # we need to sum up the amount of dead persons for each age group + + # we want the deaths for the age groups + df_abb = df_abb[['Date', 'Deaths', 'Age_RKI']] + # we want a plot with 2 rows. Second row has a plot with each age group and the simulated and real dead persons + # First row has the cumulative dead persons + fig = plt.figure('Deaths') + fig.set_figwidth(20) + fig.set_figheight(9) + gs = fig.add_gridspec(2,6) + + # we need the cumulative dead persons + ax = fig.add_subplot(gs[0, :]) + df_total_dead = df_abb.groupby('Date').sum()[0:90] + y_real = df_total_dead['Deaths'].to_numpy() + # we need to substract the first value from the rest + y_real = y_real - y_real[0] + + y_sim = p50_bs['Total'][()][:, 7][::24][0:90] + y_sim = y_sim - y_sim[0] + + y_sim25 = p25_bs['Total'][()][:, 7][::24][0:90] + y_sim25 = y_sim25 - y_sim25[0] + + y_sim75 = p75_bs['Total'][()][:, 7][::24][0:90] + y_sim75 = y_sim75 - y_sim75[0] + + y_sim05 = p05_bs['Total'][()][:, 7][::24][0:90] + y_sim05 = y_sim05 - y_sim05[0] + + y_sim95 = p95_bs['Total'][()][:, 7][::24][0:90] + y_sim95 = y_sim95 - y_sim95[0] + + + + # we calculate the RMSE + rmse_dead = np.sqrt(((y_real- y_sim)**2).mean()) + # we need to plot the cumulative dead persons from the real world and from the simulation + + ax.plot(df_total_dead.index, y_sim, color='tab:blue',label='Simulated deaths') + ax.plot(df_total_dead.index, y_real, 'v',color='tab:red', linewidth=4, label='Extrapolated deaths from reported infection case data') + ax.fill_between(df_total_dead.index, y_sim75, y_sim25, + alpha=0.5, color='tab:blue', label='50% Confidence interval') + ax.fill_between(df_total_dead.index,y_sim95, y_sim05, + alpha=0.25, color='tab:blue', label='90% Confidence interval') + # ax.text(0.25, 0.8, 'RMSE: '+str(float("{:.2f}".format(rmse_dead))), horizontalalignment='center', + # verticalalignment='center', transform=plt.gca().transAxes, color='pink', fontsize=15) + ax.set_label('Number of individuals') + ax.set_title('Cumulative Deaths', fontsize=fontsize) + ax.set_ylabel('Number of individuals', fontsize=fontsize-8) + ax.legend(fontsize=fontsize-8) + + # now for each age group + for i, age_group in zip(range(6), age_group_access): + ax = fig.add_subplot(gs[1, i]) + # we need the amount of dead persons for each age group + df_abb_age_group = df_abb[df_abb['Age_RKI'] == age_groups[i]][0:90] + y_real = np.round(df_abb_age_group['Deaths'].to_numpy()) + # we need to plot the dead persons from the real world and from the simulation + ax.plot(df_abb_age_group['Date'], y_real-y_real[0], color='tab:red') + ax.plot(df_abb_age_group['Date'], p50_bs[age_group_access[i]][()][:, 7][::24][0:90]-p50_bs[age_group_access[i]][()][:, 7][::24][0], color='tab:blue') + ax.fill_between(df_abb_age_group['Date'], p75_bs[age_group_access[i]][()][:, 7][::24][0:90]-p75_bs[age_group_access[i]][()][:, 7][::24][0], p25_bs[age_group_access[i]][()][:, 7][::24][0:90]-p25_bs[age_group_access[i]][()][:, 7][::24][0], + alpha=0.5, color='tab:blue') + ax.set_title('Deaths, '+age_grous_string[i]) + ax.set_ybound(lower=0) + ax.set_xticks(df_abb_age_group['Date'][::50]) + ax.tick_params(axis='both', which='major', labelsize=fontsize-10) + ax.tick_params(axis='both', which='minor', labelsize=fontsize-10) + if i == 0: + ax.set_ylabel('Number of individuals',fontsize=fontsize-8) + ax.set_ybound(upper=1) + + plt.show() + +def plot_icu(path): + + df_abb = pd.read_json(path+"/../../../pydata/Germany/county_divi.json") + + perc_of_critical_in_icu_age = [0.55,0.55,0.55,0.56,0.54,0.46] + perc_of_critical_in_icu=0.55 + + age_group_access = ['Group1', 'Group2', 'Group3', + 'Group4', 'Group5', 'Group6', 'Total'] + + + # we just need the columns ICU_low and ICU_hig + df_abb = df_abb[['ID_County', 'ICU', 'Date']] + + df_abb = df_abb[df_abb['ID_County'] == 3101] + # we need just the dates bewteen 2021-03-01 and 2021-06-01 + df_abb = df_abb[(df_abb['Date'] >= '2021-03-01') & + (df_abb['Date'] <= '2021-06-01')] + + # we plot this against this the Amount of persons in the ICU from our model + f_p50 = h5py.File( + path+"/infection_state_per_age_group/0/p50/Results.h5", 'r') + total_50 = f_p50['0']['Total'][()][::24][0:90] + + total_50_age = f_p50['0'][age_group_access[0]][()] + for i in range(6): + total_50_age += f_p50['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] + total_50_age = total_50_age[::24][0:90] + + + # we plot this against this the Amount of persons in the ICU from our model + f_p75 = h5py.File( + path+"/infection_state_per_age_group/0/p75/Results.h5", 'r') + # total_75 = f_p75['0']['Total'][()][::24][0:90] + total_75_age = f_p75['0'][age_group_access[0]][()] + for i in range(6): + total_75_age += f_p75['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] + total_75_age = total_75_age[::24][0:90] + + # same with 25 percentile + f_p25 = h5py.File( + path+"/infection_state_per_age_group/0/p25/Results.h5", 'r') + # total_25 = f_p25['0']['Total'][()][::24][0:90] + total_25_age = f_p25['0'][age_group_access[0]][()] + for i in range(6): + total_25_age += f_p25['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] + total_25_age = total_25_age[::24][0:90] + + # same with 05 and 95 percentile + f_p05 = h5py.File( + path+"/infection_state_per_age_group/0/p05/Results.h5", 'r') + # total_05 = f_p05['0']['Total'][()][::24][0:90] + total_05_age = f_p05['0'][age_group_access[0]][()] + for i in range(6): + total_05_age += f_p05['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] + total_05_age = total_05_age[::24][0:90] + + f_p95 = h5py.File( + path+"/infection_state_per_age_group/0/p95/Results.h5", 'r') + # total_95 = f_p95['0']['Total'][()][::24][0:90] + total_95_age = f_p95['0'][age_group_access[0]][()] + for i in range(6): + total_95_age += f_p95['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] + total_95_age = total_95_age[::24][0:90] + + + ICU_Simulation_one_percentile = np.floor(total_50[:, 5]*perc_of_critical_in_icu) + ICU_Simulation = np.round(total_50_age[:, 5]) + ICU_Simulation75 = np.round(total_75_age[:, 5]) + ICU_Simulation25 = np.round(total_25_age[:, 5]) + ICU_Simulation05 = np.round(total_05_age[:, 5]) + ICU_Simulation95 = np.round(total_95_age[:, 5]) + ICU_Real = df_abb['ICU'][0:90] + + #smooth the data + # ICU_Real = gaussian_filter1d(ICU_Real, sigma=1, mode='nearest') + # ICU_Simulation = gaussian_filter1d(ICU_Simulation, sigma=1, mode='nearest') + + + + # we calculate the RMSE + rmse_ICU = np.sqrt(((ICU_Real - ICU_Simulation_one_percentile)**2).mean()) + + # plot the ICU beds and the ICU beds taken + fig, ax = plt.subplots(1, 1, constrained_layout=True) + fig.set_figwidth(12) + fig.set_figheight(9) + # we plot the ICU_low and the ICU_high + ax.plot(df_abb['Date'][0:90], ICU_Real,'x', color='tab:red', linewidth=10, label='Data') + ax.plot(df_abb['Date'][0:90], ICU_Simulation, color='tab:blue', label='Simulation') + # ax.plot(df_abb['Date'][0:90], ICU_Simulation_one_percentile, color='tab:green', label='Simulated ICU beds') + ax.fill_between(df_abb['Date'][0:90],ICU_Simulation75, ICU_Simulation25, + alpha=0.5, color='tab:blue', label='50% Confidence interval') + ax.fill_between(df_abb['Date'][0:90],ICU_Simulation05, ICU_Simulation95, + alpha=0.25, color='tab:blue', label='90% Confidence interval') + + + # we also write the rmse + # ax.text(0.25, 0.8, 'RMSE: '+str(float("{:.2f}".format(rmse_ICU))), horizontalalignment='center', + # verticalalignment='center', transform=plt.gca().transAxes, color='pink', fontsize=15) + ax.tick_params(axis='both', which='major', labelsize=fontsize-4) + ax.tick_params(axis='both', which='minor', labelsize=fontsize-4) + ax.set_ylabel('Occupied ICU beds', fontsize=fontsize) + ax.set_title('ICU beds', fontsize=fontsize+4) + ax.legend(fontsize=fontsize-4) + plt.show() + + + + +if __name__ == "__main__": + path = "" + + if (len(sys.argv) > 1): + n_runs = sys.argv[1] + else: + n_runs = len([entry for entry in os.listdir(path) + if os.path.isfile(os.path.join(path, entry))]) + + # plot_icu(path) + # plot_dead(path) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py new file mode 100644 index 0000000000..cf33fb8d3d --- /dev/null +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -0,0 +1,258 @@ +# Python script to analyze bs runs +# input is a bs run folder with the following structure: +# bs_run_folder has a txt file for each bs run +# each txt file has a line for each time step +# each line has a column for each compartment as well as the timestep +# each column has the number of individuals in that compartment +# the first line of each txt file is the header + +import sys +import argparse +import os +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +import matplotlib.colors as colors +import matplotlib.cm as cmx +import matplotlib.patches as mpatches +import matplotlib.lines as mlines +import h5py +from datetime import datetime +from matplotlib.dates import DateFormatter +from scipy.ndimage import gaussian_filter1d +from scipy.signal import savgol_filter + +fontsize = 20 + + +def plot_infections_loc_types_avarage(path): + # 50-percentile + f_p50 = h5py.File( + path+"/infection_per_location_type_per_age_group/0/p50/Results.h5", 'r') + p50_bs = f_p50['0'] + total_50 = p50_bs['Total'][()] + + # 25-percentile + f_p25 = h5py.File( + path+"/infection_per_location_type_per_age_group/0/p05/Results.h5", 'r') + p25_bs = f_p25['0'] + total_25 = p25_bs['Total'][()] + + # 75-percentile + f_p75 = h5py.File( + path + "/infection_per_location_type_per_age_group/0/p95/Results.h5", 'r') + p75_bs = f_p75['0'] + total_75 = p75_bs['Total'][()] + + time = p50_bs['Time'][()] + + plot_infection_per_location_type_mean( + time, total_50, total_25, total_75) + +def plot_infection_per_location_type_mean(x, y50, y25, y75): + + plt.figure('Infection_location_types') + plt.title('At which location type an infection happened, avaraged over all runs') + + color_plot = matplotlib.colormaps.get_cmap('Set1').colors + + states_plot = [0, 1, 2, 3, 4, 10] + legend_plot = ['Home', 'School', 'Work', + 'SocialEvent', 'BasicsShop','Event'] + + for i in states_plot: + # rolling average# + ac_color = color_plot[i%len(color_plot)] + if(i > len(color_plot)): + ac_color = "black" + + # we need to sum up every 24 hours + indexer = pd.api.indexers.FixedForwardWindowIndexer(window_size=24) + np_y50 = pd.DataFrame(y50[:, i]).rolling(window=indexer, min_periods=1).sum().to_numpy() + np_y50=np_y50[0::24].flatten() + # now smoothen this with a gaussian filter + np_y50 = gaussian_filter1d(np_y50, sigma=1, mode='nearest') + + plt.plot(x[0::24], np_y50, color=ac_color) + + plt.legend(legend_plot) + + # currently the x axis has the values of the time steps, we need to convert them to dates and set the x axis to dates + start_date = datetime.strptime('2021-03-01', '%Y-%m-%d') + xx = [start_date + pd.Timedelta(days=int(i)) for i in x] + xx = [xx[i].strftime('%Y-%m-%d') for i in range(len(xx))] + # but just take every 10th date to make it more readable + plt.gca().set_xticks(x[::150]) + plt.gca().set_xticklabels(xx[::150]) + plt.gcf().autofmt_xdate() + + plt.xlabel('Date') + plt.ylabel('Number of individuals') + plt.show() + +def plot_infection_states_results(path): + # 50-percentile + f_p50 = h5py.File( + path+"/infection_state_per_age_group/0/p50/Results.h5", 'r') + p50_bs = f_p50['0'] + total_50 = p50_bs['Total'][()] + # 25-percentile + f_p25 = h5py.File( + path+"/infection_state_per_age_group/0/p05/Results.h5", 'r') + p25_bs = f_p25['0'] + total_25 = p25_bs['Total'][()] + # 75-percentile + f_p75 = h5py.File( + path + "/infection_state_per_age_group/0/p95/Results.h5", 'r') + p75_bs = f_p75['0'] + total_75 = p75_bs['Total'][()] + + time = p50_bs['Time'][()] + + plot_infection_states_individual( + time, p50_bs, p25_bs, p75_bs) + plot_infection_states(time, total_50, total_25, total_75) + +def plot_infection_states(x, y50, y25, y75, y_real=None): + plt.figure('Infection_states') + plt.title('Infection states') + + color_plot = matplotlib.colormaps.get_cmap('Set1').colors + + states_plot = [1, 2, 3, 4, 5, 7] + legend_plot = ['E', 'I_NSymp', 'I_Symp', + 'I_Sev', 'I_Crit', 'Dead', 'Sm. re. pos.'] + + for i in states_plot: + plt.plot(x, y50[:, i], color=color_plot[i]) + + + plt.legend(legend_plot) + for i in states_plot: + plt.fill_between(x, y50[:, i], y25[:, i], + alpha=0.5, color=color_plot[i]) + plt.fill_between(x, y50[:, i], y75[:, i], + alpha=0.5, color=color_plot[i]) + + # currently the x axis has the values of the time steps, we need to convert them to dates and set the x axis to dates + start_date = datetime.strptime('2021-03-01', '%Y-%m-%d') + xx = [start_date + pd.Timedelta(days=int(i)) for i in x] + xx = [xx[i].strftime('%Y-%m-%d') for i in range(len(xx))] + # but just take every 10th date to make it more readable + plt.gca().set_xticks(x[::150]) + plt.gca().set_xticklabels(xx[::150]) + plt.gcf().autofmt_xdate() + + plt.xlabel('Time') + plt.ylabel('Number of individuals') + plt.show() + +def plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs): + + + age_group_access = ['Group1', 'Group2', 'Group3', + 'Group4', 'Group5', 'Group6', 'Total'] + + color_plot = matplotlib.colormaps.get_cmap('Set1').colors + + fig, ax = plt.subplots(6, len(age_group_access), constrained_layout=True) + fig.set_figwidth(20) + fig.set_figheight(9) + for j, count in zip(age_group_access, range(len(age_group_access))): + y50 = p50_bs[j][()] + y25 = p25_bs[j][()] + y75 = p75_bs[j][()] + + + # infeced no symptoms + ax_infected_no_symptoms = ax[0, count] + ax_infected_no_symptoms.set_xlabel('time (days)') + ax_infected_no_symptoms.plot( + x, y50[:, 1], color=color_plot[count], label='y50') + ax_infected_no_symptoms.fill_between( + x, y50[:, 1], y25[:, 1], alpha=0.5, color=color_plot[count]) + ax_infected_no_symptoms.fill_between( + x, y50[:, 1], y75[:, 1], alpha=0.5, color=color_plot[count]) + ax_infected_no_symptoms.tick_params(axis='y') + ax_infected_no_symptoms.title.set_text( + '#Infected_no_symptoms, Age{}'.format(j)) + ax_infected_no_symptoms.legend(['Simulation']) + + # Infected_symptoms + ax_infected_symptoms = ax[1, count] + ax_infected_symptoms.set_xlabel('time (days)') + ax_infected_symptoms.plot( + x, y50[:, 2], color=color_plot[count], label='y50') + ax_infected_symptoms.fill_between( + x, y50[:, 2], y25[:, 2], alpha=0.5, color=color_plot[count]) + ax_infected_symptoms.fill_between( + x, y50[:, 2], y75[:, 2], alpha=0.5, color=color_plot[count]) + ax_infected_symptoms.tick_params(axis='y') + ax_infected_symptoms.title.set_text( + '#Infected_symptoms, Age{}'.format(j)) + ax_infected_symptoms.legend(['Simulation']) + + # Severe + ax_severe = ax[2, count] + ax_severe.set_xlabel('time (days)') + ax_severe.plot(x, y50[:, 4], color=color_plot[count], label='y50') + ax_severe.fill_between( + x, y50[:, 4], y25[:, 4], alpha=0.5, color=color_plot[count]) + ax_severe.fill_between( + x, y50[:, 4], y75[:, 4], alpha=0.5, color=color_plot[count]) + ax_severe.tick_params(axis='y') + ax_severe.title.set_text('#Severe, Age{}'.format(j)) + ax_severe.legend(['Simulation']) + + # Critical + ax_critical = ax[3, count] + ax_critical.set_xlabel('time (days)') + ax_critical.plot(x, y50[:, [5]], color=color_plot[count], label='y50') + ax_critical.fill_between( + x, y50[:, 5], y25[:, 5], alpha=0.5, color=color_plot[count]) + ax_critical.fill_between( + x, y50[:, 5], y75[:, 5], alpha=0.5, color=color_plot[count]) + ax_critical.tick_params(axis='y') + ax_critical.title.set_text('#Critical, Age{}'.format(j)) + ax_critical.legend(['Simulation']) + + # Dead + ax_dead = ax[4, count] + ax_dead.set_xlabel('time (days)') + ax_dead.plot(x, y50[:, [7]], color=color_plot[count], label='y50') + ax_dead.fill_between(x, y50[:, 7], y25[:, 7], + alpha=0.5, color=color_plot[count]) + ax_dead.fill_between(x, y50[:, 7], y75[:, 7], + alpha=0.5, color=color_plot[count]) + ax_dead.tick_params(axis='y') + ax_dead.title.set_text('#Dead, Age{}'.format(j)) + ax_dead.legend(['Simulation']) + + # Recovered + ax_dead = ax[5, count] + ax_dead.set_xlabel('time (days)') + ax_dead.plot(x, y50[:, [6]], color=color_plot[count], label='y50') + ax_dead.fill_between(x, y50[:, 6], y25[:, 6], + alpha=0.5, color=color_plot[count]) + ax_dead.fill_between(x, y50[:, 6], y75[:, 6], + alpha=0.5, color=color_plot[count]) + ax_dead.tick_params(axis='y') + ax_dead.title.set_text('#Recovered, Age{}'.format(j)) + ax_dead.legend(['Simulation']) + + # fig.tight_layout() # otherwise the right y-label is slightly clipped + plt.show() + + +if __name__ == "__main__": + path = "" + + if (len(sys.argv) > 1): + n_runs = sys.argv[1] + else: + n_runs = len([entry for entry in os.listdir(path) + if os.path.isfile(os.path.join(path, entry))]) + + # plot_infection_states_results(path) + # plot_infections_loc_types_avarage(path) \ No newline at end of file From c7ffab3663b49f2d14c8f44773aeba2df7769402 Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Tue, 27 May 2025 22:32:04 +0200 Subject: [PATCH 02/38] Remove unused imports and enable plotting functions in infection states script --- pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py | 2 -- pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py | 6 ++---- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py b/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py index 6a1af97ed5..510a88104e 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py @@ -20,8 +20,6 @@ import h5py from datetime import datetime from matplotlib.dates import DateFormatter -from scipy.ndimage import gaussian_filter1d -from scipy.signal import savgol_filter fontsize = 20 diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index cf33fb8d3d..3f92bf4c0b 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -20,8 +20,6 @@ import h5py from datetime import datetime from matplotlib.dates import DateFormatter -from scipy.ndimage import gaussian_filter1d -from scipy.signal import savgol_filter fontsize = 20 @@ -254,5 +252,5 @@ def plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs): n_runs = len([entry for entry in os.listdir(path) if os.path.isfile(os.path.join(path, entry))]) - # plot_infection_states_results(path) - # plot_infections_loc_types_avarage(path) \ No newline at end of file + plot_infection_states_results(path) + plot_infections_loc_types_avarage(path) \ No newline at end of file From e986d96c25db8cbe0bc8d63f1d9a9070d94f7800 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Tue, 27 May 2025 23:49:21 +0200 Subject: [PATCH 03/38] Refactor infection states plotting script: enhance documentation, improve function structure, and add new plotting capabilities for location types. --- .../memilio/plot/plotAbmInfectionStates.py | 472 ++++++++++-------- 1 file changed, 253 insertions(+), 219 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 3f92bf4c0b..352c5f5998 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -1,10 +1,22 @@ -# Python script to analyze bs runs -# input is a bs run folder with the following structure: -# bs_run_folder has a txt file for each bs run -# each txt file has a line for each time step -# each line has a column for each compartment as well as the timestep -# each column has the number of individuals in that compartment -# the first line of each txt file is the header +############################################################################# +# Copyright (C) 2020-2024 MEmilio +# +# Authors: Daniel Abele, Martin J. Kuehn +# +# 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 sys import argparse @@ -13,244 +25,266 @@ import numpy as np import matplotlib.pyplot as plt import matplotlib -import matplotlib.colors as colors -import matplotlib.cm as cmx -import matplotlib.patches as mpatches -import matplotlib.lines as mlines import h5py from datetime import datetime -from matplotlib.dates import DateFormatter - -fontsize = 20 - - -def plot_infections_loc_types_avarage(path): - # 50-percentile - f_p50 = h5py.File( - path+"/infection_per_location_type_per_age_group/0/p50/Results.h5", 'r') - p50_bs = f_p50['0'] - total_50 = p50_bs['Total'][()] - - # 25-percentile - f_p25 = h5py.File( - path+"/infection_per_location_type_per_age_group/0/p05/Results.h5", 'r') - p25_bs = f_p25['0'] - total_25 = p25_bs['Total'][()] - - # 75-percentile - f_p75 = h5py.File( - path + "/infection_per_location_type_per_age_group/0/p95/Results.h5", 'r') - p75_bs = f_p75['0'] - total_75 = p75_bs['Total'][()] - - time = p50_bs['Time'][()] - - plot_infection_per_location_type_mean( - time, total_50, total_25, total_75) - -def plot_infection_per_location_type_mean(x, y50, y25, y75): +from scipy.ndimage import gaussian_filter1d + + +""" Module for plotting infection states and location types from ABM results. +This module provides functions to load and visualize infection states and +location types from simulation results stored in HDF5 format and are output +by the MEmilio agent-based model (ABM). +The used Loggers are: +struct LogInfectionStatePerAgeGroup : mio::LogAlways { + using Type = std::pair; + /** + * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. + * @param[in] sim The simulation of the abm. + * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. + */ + static Type log(const mio::abm::Simulation& sim) + { + + Eigen::VectorXd sum = Eigen::VectorXd::Zero( + Eigen::Index((size_t)mio::abm::InfectionState::Count * sim.get_world().parameters.get_num_groups())); + const auto curr_time = sim.get_time(); + const auto persons = sim.get_world().get_persons(); + + // PRAGMA_OMP(parallel for) + for (auto i = size_t(0); i < persons.size(); ++i) { + auto& p = persons[i]; + if (p.get_should_be_logged()) { + auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) + + ((uint32_t)p.get_infection_state(curr_time)); + // PRAGMA_OMP(atomic) + sum[index] += 1; + } + } + return std::make_pair(curr_time, sum); + } +}; + +struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { + using Type = std::pair; + /** + * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. + * @param[in] sim The simulation of the abm. + * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. + */ + static Type log(const mio::abm::Simulation& sim) + { + + Eigen::VectorXd sum = Eigen::VectorXd::Zero( + Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_world().parameters.get_num_groups())); + auto curr_time = sim.get_time(); + auto prev_time = sim.get_prev_time(); + const auto persons = sim.get_world().get_persons(); + + // PRAGMA_OMP(parallel for) + for (auto i = size_t(0); i < persons.size(); ++i) { + auto& p = persons[i]; + if (p.get_should_be_logged()) { + // PRAGMA_OMP(atomic) + if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) && + (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { + auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) + + ((uint32_t)p.get_location().get_type()); + sum[index] += 1; + } + } + } + return std::make_pair(curr_time, sum); + } +}; + +The output of several Loggers is stored in HDF5 files, with the memilio funciton mio::save_results in mio/io/result_io.h. +""" + + + + +def load_h5_results(base_path, percentile): + """ Reads HDF5 results for a given group and percentile. + + @param[in] base_path Path to results directory. + @param[in] percentile Subdirectory for percentile (e.g. 'p50'). + @return Dictionary with data arrays. + """ + file_path = os.path.join(base_path, percentile, "Results.h5") + with h5py.File(file_path, 'r') as f: + data = {k: v[()] for k, v in f['0'].items()} + return data + +def plot_infections_loc_types_average( + path_to_loc_types, + start_date='2021-03-01', + colormap='Set1', + smooth_sigma=1, + rolling_window=24, + xtick_step=150): + """ Plots average infections per location type. + + @param[in] base_path Path to results directory. + @param[in] start_date Start date as string. + @param[in] colormap Matplotlib colormap. + @param[in] smooth_sigma Sigma for Gaussian smoothing. + @param[in] rolling_window Window size for rolling sum. + @param[in] xtick_step Step size for x-axis ticks. + """ + # Load data + p50 = load_h5_results(path_to_loc_types, "p50") + p25 = load_h5_results(path_to_loc_types, "p05") + p75 = load_h5_results(path_to_loc_types, "p95") + time = p50['Time'] + total_50 = p50['Total'] plt.figure('Infection_location_types') - plt.title('At which location type an infection happened, avaraged over all runs') - - color_plot = matplotlib.colormaps.get_cmap('Set1').colors - - states_plot = [0, 1, 2, 3, 4, 10] - legend_plot = ['Home', 'School', 'Work', - 'SocialEvent', 'BasicsShop','Event'] - - for i in states_plot: - # rolling average# - ac_color = color_plot[i%len(color_plot)] - if(i > len(color_plot)): - ac_color = "black" - - # we need to sum up every 24 hours - indexer = pd.api.indexers.FixedForwardWindowIndexer(window_size=24) - np_y50 = pd.DataFrame(y50[:, i]).rolling(window=indexer, min_periods=1).sum().to_numpy() - np_y50=np_y50[0::24].flatten() - # now smoothen this with a gaussian filter - np_y50 = gaussian_filter1d(np_y50, sigma=1, mode='nearest') - - plt.plot(x[0::24], np_y50, color=ac_color) + plt.title('At which location type an infection happened, averaged over all runs') + color_plot = matplotlib.colormaps.get_cmap(colormap).colors + # If you define further location types, you need to adjust this list + states_plot = [0, 1, 2, 3, 4, 5 , 6] + legend_plot = ['Home', 'School', 'Work', 'SocialEvent', 'BasicsShop', 'Hospital', 'ICU'] + + for idx, i in enumerate(states_plot): + color = color_plot[i % len(color_plot)] if i < len(color_plot) else "black" + # Sum up every 24 hours, then smooth + indexer = pd.api.indexers.FixedForwardWindowIndexer(window_size=rolling_window) + y = pd.DataFrame(total_50[:, i]).rolling(window=indexer, min_periods=1).sum().to_numpy() + y = y[0::rolling_window].flatten() + y = gaussian_filter1d(y, sigma=smooth_sigma, mode='nearest') + plt.plot(time[0::rolling_window], y, color=color) plt.legend(legend_plot) - - # currently the x axis has the values of the time steps, we need to convert them to dates and set the x axis to dates - start_date = datetime.strptime('2021-03-01', '%Y-%m-%d') - xx = [start_date + pd.Timedelta(days=int(i)) for i in x] - xx = [xx[i].strftime('%Y-%m-%d') for i in range(len(xx))] - # but just take every 10th date to make it more readable - plt.gca().set_xticks(x[::150]) - plt.gca().set_xticklabels(xx[::150]) - plt.gcf().autofmt_xdate() - + _format_x_axis(time, start_date, xtick_step) plt.xlabel('Date') plt.ylabel('Number of individuals') plt.show() -def plot_infection_states_results(path): - # 50-percentile - f_p50 = h5py.File( - path+"/infection_state_per_age_group/0/p50/Results.h5", 'r') - p50_bs = f_p50['0'] - total_50 = p50_bs['Total'][()] - # 25-percentile - f_p25 = h5py.File( - path+"/infection_state_per_age_group/0/p05/Results.h5", 'r') - p25_bs = f_p25['0'] - total_25 = p25_bs['Total'][()] - # 75-percentile - f_p75 = h5py.File( - path + "/infection_state_per_age_group/0/p95/Results.h5", 'r') - p75_bs = f_p75['0'] - total_75 = p75_bs['Total'][()] - - time = p50_bs['Time'][()] - - plot_infection_states_individual( - time, p50_bs, p25_bs, p75_bs) - plot_infection_states(time, total_50, total_25, total_75) - -def plot_infection_states(x, y50, y25, y75, y_real=None): - plt.figure('Infection_states') +def plot_infection_states_results( + path_to_infection_states, + start_date='2021-03-01', + colormap='Set1', + xtick_step=150): + """ Loads and plots infection state results. + + @param[in] base_path Path to results directory. + @param[in] start_date Start date as string. + @param[in] colormap Matplotlib colormap. + @param[in] xtick_step Step size for x-axis ticks. + """ + # Load data + p50 = load_h5_results(path_to_infection_states, "p50") + p25 = load_h5_results(path_to_infection_states, "p25") + p75 = load_h5_results(path_to_infection_states, "p75") + time = p50['Time'] + total_50 = p50['Total'] + total_25 = p25['Total'] + total_75 = p75['Total'] + + plot_infection_states_individual(time, p50, p25, p75, colormap) + plot_infection_states(time, total_50, total_25, total_75, start_date, colormap, xtick_step) + +def plot_infection_states( + x, y50, y25, y75, + start_date='2021-03-01', + colormap='Set1', + xtick_step=150): + """ Plots infection states with percentiles. + + @param[in] x Time array. + @param[in] y50 Median values. + @param[in] y25 25th percentile values. + @param[in] y75 75th percentile values. + @param[in] start_date Start date as string. + @param[in] colormap Matplotlib colormap. + @param[in] xtick_step Step size for x-axis ticks. + """ + plt.figure('Infection_states with 50% percentile') plt.title('Infection states') - - color_plot = matplotlib.colormaps.get_cmap('Set1').colors - + color_plot = matplotlib.colormaps.get_cmap(colormap).colors states_plot = [1, 2, 3, 4, 5, 7] - legend_plot = ['E', 'I_NSymp', 'I_Symp', - 'I_Sev', 'I_Crit', 'Dead', 'Sm. re. pos.'] + legend_plot = ['E', 'I_NSymp', 'I_Symp', 'I_Sev', 'I_Crit', 'Dead'] for i in states_plot: plt.plot(x, y50[:, i], color=color_plot[i]) - + plt.fill_between(x, y50[:, i], y25[:, i], alpha=0.5, color=color_plot[i]) + plt.fill_between(x, y50[:, i], y75[:, i], alpha=0.5, color=color_plot[i]) plt.legend(legend_plot) - for i in states_plot: - plt.fill_between(x, y50[:, i], y25[:, i], - alpha=0.5, color=color_plot[i]) - plt.fill_between(x, y50[:, i], y75[:, i], - alpha=0.5, color=color_plot[i]) - - # currently the x axis has the values of the time steps, we need to convert them to dates and set the x axis to dates - start_date = datetime.strptime('2021-03-01', '%Y-%m-%d') - xx = [start_date + pd.Timedelta(days=int(i)) for i in x] - xx = [xx[i].strftime('%Y-%m-%d') for i in range(len(xx))] - # but just take every 10th date to make it more readable - plt.gca().set_xticks(x[::150]) - plt.gca().set_xticklabels(xx[::150]) - plt.gcf().autofmt_xdate() - + _format_x_axis(x, start_date, xtick_step) plt.xlabel('Time') plt.ylabel('Number of individuals') plt.show() -def plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs): - - - age_group_access = ['Group1', 'Group2', 'Group3', - 'Group4', 'Group5', 'Group6', 'Total'] - - color_plot = matplotlib.colormaps.get_cmap('Set1').colors - - fig, ax = plt.subplots(6, len(age_group_access), constrained_layout=True) - fig.set_figwidth(20) - fig.set_figheight(9) - for j, count in zip(age_group_access, range(len(age_group_access))): - y50 = p50_bs[j][()] - y25 = p25_bs[j][()] - y75 = p75_bs[j][()] - - - # infeced no symptoms - ax_infected_no_symptoms = ax[0, count] - ax_infected_no_symptoms.set_xlabel('time (days)') - ax_infected_no_symptoms.plot( - x, y50[:, 1], color=color_plot[count], label='y50') - ax_infected_no_symptoms.fill_between( - x, y50[:, 1], y25[:, 1], alpha=0.5, color=color_plot[count]) - ax_infected_no_symptoms.fill_between( - x, y50[:, 1], y75[:, 1], alpha=0.5, color=color_plot[count]) - ax_infected_no_symptoms.tick_params(axis='y') - ax_infected_no_symptoms.title.set_text( - '#Infected_no_symptoms, Age{}'.format(j)) - ax_infected_no_symptoms.legend(['Simulation']) - - # Infected_symptoms - ax_infected_symptoms = ax[1, count] - ax_infected_symptoms.set_xlabel('time (days)') - ax_infected_symptoms.plot( - x, y50[:, 2], color=color_plot[count], label='y50') - ax_infected_symptoms.fill_between( - x, y50[:, 2], y25[:, 2], alpha=0.5, color=color_plot[count]) - ax_infected_symptoms.fill_between( - x, y50[:, 2], y75[:, 2], alpha=0.5, color=color_plot[count]) - ax_infected_symptoms.tick_params(axis='y') - ax_infected_symptoms.title.set_text( - '#Infected_symptoms, Age{}'.format(j)) - ax_infected_symptoms.legend(['Simulation']) - +def plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs, colormap='Set1'): + """ Plots infection states for each age group. + + @param[in] x Time array. + @param[in] p50_bs Median values by group. + @param[in] p25_bs 25th percentile values by group. + @param[in] p75_bs 75th percentile values by group. + @param[in] colormap Matplotlib colormap. + """ + age_groups = ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total'] # Adjust as needed + color_plot = matplotlib.colormaps.get_cmap(colormap).colors + fig, ax = plt.subplots(6, len(age_groups), constrained_layout=True, figsize=(20, 9)) + + for col_idx, group in enumerate(age_groups): + y50 = p50_bs[group] + y25 = p25_bs[group] + y75 = p75_bs[group] + # Infected no symptoms + _plot_state(ax[0, col_idx], x, y50[:, 1], y25[:, 1], y75[:, 1], color_plot[col_idx], '#Infected_no_symptoms, Age' + str(group)) + # Infected symptoms + _plot_state(ax[1, col_idx], x, y50[:, 2], y25[:, 2], y75[:, 2], color_plot[col_idx], '#Infected_symptoms, Age' + str(group)) # Severe - ax_severe = ax[2, count] - ax_severe.set_xlabel('time (days)') - ax_severe.plot(x, y50[:, 4], color=color_plot[count], label='y50') - ax_severe.fill_between( - x, y50[:, 4], y25[:, 4], alpha=0.5, color=color_plot[count]) - ax_severe.fill_between( - x, y50[:, 4], y75[:, 4], alpha=0.5, color=color_plot[count]) - ax_severe.tick_params(axis='y') - ax_severe.title.set_text('#Severe, Age{}'.format(j)) - ax_severe.legend(['Simulation']) - + _plot_state(ax[2, col_idx], x, y50[:, 4], y25[:, 4], y75[:, 4], color_plot[col_idx], '#Severe, Age' + str(group)) # Critical - ax_critical = ax[3, count] - ax_critical.set_xlabel('time (days)') - ax_critical.plot(x, y50[:, [5]], color=color_plot[count], label='y50') - ax_critical.fill_between( - x, y50[:, 5], y25[:, 5], alpha=0.5, color=color_plot[count]) - ax_critical.fill_between( - x, y50[:, 5], y75[:, 5], alpha=0.5, color=color_plot[count]) - ax_critical.tick_params(axis='y') - ax_critical.title.set_text('#Critical, Age{}'.format(j)) - ax_critical.legend(['Simulation']) - + _plot_state(ax[3, col_idx], x, y50[:, 5], y25[:, 5], y75[:, 5], color_plot[col_idx], '#Critical, Age' + str(group)) # Dead - ax_dead = ax[4, count] - ax_dead.set_xlabel('time (days)') - ax_dead.plot(x, y50[:, [7]], color=color_plot[count], label='y50') - ax_dead.fill_between(x, y50[:, 7], y25[:, 7], - alpha=0.5, color=color_plot[count]) - ax_dead.fill_between(x, y50[:, 7], y75[:, 7], - alpha=0.5, color=color_plot[count]) - ax_dead.tick_params(axis='y') - ax_dead.title.set_text('#Dead, Age{}'.format(j)) - ax_dead.legend(['Simulation']) - + _plot_state(ax[4, col_idx], x, y50[:, 7], y25[:, 7], y75[:, 7], color_plot[col_idx], '#Dead, Age' + str(group)) # Recovered - ax_dead = ax[5, count] - ax_dead.set_xlabel('time (days)') - ax_dead.plot(x, y50[:, [6]], color=color_plot[count], label='y50') - ax_dead.fill_between(x, y50[:, 6], y25[:, 6], - alpha=0.5, color=color_plot[count]) - ax_dead.fill_between(x, y50[:, 6], y75[:, 6], - alpha=0.5, color=color_plot[count]) - ax_dead.tick_params(axis='y') - ax_dead.title.set_text('#Recovered, Age{}'.format(j)) - ax_dead.legend(['Simulation']) - - # fig.tight_layout() # otherwise the right y-label is slightly clipped + _plot_state(ax[5, col_idx], x, y50[:, 6], y25[:, 6], y75[:, 6], color_plot[col_idx], '#Recovered, Age' + str(group)) + + fig.suptitle('Infection states per age group with 50% percentile', fontsize=16) + + # We hide the Legend for the individual plots as it is too cluttered + for ax_row in ax: + for ax_col in ax_row: + ax_col.legend().set_visible(False) + plt.show() +def _plot_state(ax, x, y50, y25, y75, color, title): + """ Helper to plot a single state with fill_between. """ + ax.set_xlabel('time (days)') + ax.plot(x, y50, color=color, label='Median') + ax.fill_between(x, y50, y25, alpha=0.5, color=color) + ax.fill_between(x, y50, y75, alpha=0.5, color=color) + ax.tick_params(axis='y') + ax.set_title(title) + ax.legend(['Simulation']) + +def _format_x_axis(x, start_date, xtick_step): + """ Helper to format x-axis as dates. """ + start = datetime.strptime(start_date, '%Y-%m-%d') + xx = [start + pd.Timedelta(days=int(i)) for i in x] + xx_str = [dt.strftime('%Y-%m-%d') for dt in xx] + plt.gca().set_xticks(x[::xtick_step]) + plt.gca().set_xticklabels(xx_str[::xtick_step]) + plt.gcf().autofmt_xdate() + +def main(): + """ Main function for CLI usage. """ + parser = argparse.ArgumentParser(description="Plot infection state and location type results.") + # parser.add_argument("path", type=str, help="Path to the results folder") + parser.add_argument("--start-date", type=str, default='2021-03-01', help="Simulation start date (YYYY-MM-DD)") + parser.add_argument("--colormap", type=str, default='Set1', help="Matplotlib colormap") + parser.add_argument("--xtick-step", type=int, default=150, help="Step for x-axis ticks") + args = parser.parse_args() -if __name__ == "__main__": - path = "" + plot_infection_states_results("/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_state_per_age_group/0", args.start_date, args.colormap, args.xtick_step) + # plot_infections_loc_types_average("/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_per_location_type_per_age_group/0", args.start_date, args.colormap, xtick_step=args.xtick_step) - if (len(sys.argv) > 1): - n_runs = sys.argv[1] - else: - n_runs = len([entry for entry in os.listdir(path) - if os.path.isfile(os.path.join(path, entry))]) - - plot_infection_states_results(path) - plot_infections_loc_types_avarage(path) \ No newline at end of file +if __name__ == "__main__": + main() From 6e4c6fe7c02a54b3471ae9c77b1480b84661512f Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Wed, 28 May 2025 00:23:11 +0200 Subject: [PATCH 04/38] Refactor plotAbmInfectionStates.py: update authorship, enhance module documentation, and improve function argument handling for better usability --- .../memilio/plot/plotAbmInfectionStates.py | 175 ++++++++++-------- 1 file changed, 95 insertions(+), 80 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 352c5f5998..91a6926c5e 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -1,7 +1,7 @@ ############################################################################# # Copyright (C) 2020-2024 MEmilio # -# Authors: Daniel Abele, Martin J. Kuehn +# Authors: Sascha Korf # # Contact: Martin J. Kuehn # @@ -30,77 +30,74 @@ from scipy.ndimage import gaussian_filter1d -""" Module for plotting infection states and location types from ABM results. -This module provides functions to load and visualize infection states and -location types from simulation results stored in HDF5 format and are output -by the MEmilio agent-based model (ABM). -The used Loggers are: -struct LogInfectionStatePerAgeGroup : mio::LogAlways { - using Type = std::pair; - /** - * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. - * @param[in] sim The simulation of the abm. - * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. - */ - static Type log(const mio::abm::Simulation& sim) - { - - Eigen::VectorXd sum = Eigen::VectorXd::Zero( - Eigen::Index((size_t)mio::abm::InfectionState::Count * sim.get_world().parameters.get_num_groups())); - const auto curr_time = sim.get_time(); - const auto persons = sim.get_world().get_persons(); - - // PRAGMA_OMP(parallel for) - for (auto i = size_t(0); i < persons.size(); ++i) { - auto& p = persons[i]; - if (p.get_should_be_logged()) { - auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) + - ((uint32_t)p.get_infection_state(curr_time)); - // PRAGMA_OMP(atomic) - sum[index] += 1; - } - } - return std::make_pair(curr_time, sum); - } -}; - -struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { - using Type = std::pair; - /** - * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. - * @param[in] sim The simulation of the abm. - * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. - */ - static Type log(const mio::abm::Simulation& sim) - { - - Eigen::VectorXd sum = Eigen::VectorXd::Zero( - Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_world().parameters.get_num_groups())); - auto curr_time = sim.get_time(); - auto prev_time = sim.get_prev_time(); - const auto persons = sim.get_world().get_persons(); - - // PRAGMA_OMP(parallel for) - for (auto i = size_t(0); i < persons.size(); ++i) { - auto& p = persons[i]; - if (p.get_should_be_logged()) { - // PRAGMA_OMP(atomic) - if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) && - (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { - auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) + - ((uint32_t)p.get_location().get_type()); - sum[index] += 1; - } - } - } - return std::make_pair(curr_time, sum); - } -}; - -The output of several Loggers is stored in HDF5 files, with the memilio funciton mio::save_results in mio/io/result_io.h. -""" - - +# Module for plotting infection states and location types from ABM results. +# This module provides functions to load and visualize infection states and +# location types from simulation results stored in HDF5 format and are output +# by the MEmilio agent-based model (ABM). +# The used Loggers are: +# struct LogInfectionStatePerAgeGroup : mio::LogAlways { +# using Type = std::pair; +# /** +# * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. +# * @param[in] sim The simulation of the abm. +# * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. +# */ +# static Type log(const mio::abm::Simulation& sim) +# { +# +# Eigen::VectorXd sum = Eigen::VectorXd::Zero( +# Eigen::Index((size_t)mio::abm::InfectionState::Count * sim.get_world().parameters.get_num_groups())); +# const auto curr_time = sim.get_time(); +# const auto persons = sim.get_world().get_persons(); +# +# // PRAGMA_OMP(parallel for) +# for (auto i = size_t(0); i < persons.size(); ++i) { +# auto& p = persons[i]; +# if (p.get_should_be_logged()) { +# auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) + +# ((uint32_t)p.get_infection_state(curr_time)); +# // PRAGMA_OMP(atomic) +# sum[index] += 1; +# } +# } +# return std::make_pair(curr_time, sum); +# } +# }; +# +# struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { +# using Type = std::pair; +# /** +# * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. +# * @param[in] sim The simulation of the abm. +# * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. +# */ +# static Type log(const mio::abm::Simulation& sim) +# { +# +# Eigen::VectorXd sum = Eigen::VectorXd::Zero( +# Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_world().parameters.get_num_groups())); +# auto curr_time = sim.get_time(); +# auto prev_time = sim.get_prev_time(); +# const auto persons = sim.get_world().get_persons(); +# +# // PRAGMA_OMP(parallel for) +# for (auto i = size_t(0); i < persons.size(); ++i) { +# auto& p = persons[i]; +# if (p.get_should_be_logged()) { +# // PRAGMA_OMP(atomic) +# if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) && +# (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { +# auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) + +# ((uint32_t)p.get_location().get_type()); +# sum[index] += 1; +# } +# } +# } +# return std::make_pair(curr_time, sum); +# } +# }; +# +# The output of the loggers of several runs is stored in HDF5 files, with the memilio funciton mio::save_results in mio/io/result_io.h. def load_h5_results(base_path, percentile): @@ -122,7 +119,7 @@ def plot_infections_loc_types_average( smooth_sigma=1, rolling_window=24, xtick_step=150): - """ Plots average infections per location type. + """ Plots rolling average infections per location type for the median run. @param[in] base_path Path to results directory. @param[in] start_date Start date as string. @@ -133,13 +130,11 @@ def plot_infections_loc_types_average( """ # Load data p50 = load_h5_results(path_to_loc_types, "p50") - p25 = load_h5_results(path_to_loc_types, "p05") - p75 = load_h5_results(path_to_loc_types, "p95") time = p50['Time'] total_50 = p50['Total'] plt.figure('Infection_location_types') - plt.title('At which location type an infection happened, averaged over all runs') + plt.title('Infection per location type for the median run, rolling sum over 24 hours') color_plot = matplotlib.colormaps.get_cmap(colormap).colors # If you define further location types, you need to adjust this list states_plot = [0, 1, 2, 3, 4, 5 , 6] @@ -200,13 +195,16 @@ def plot_infection_states( @param[in] xtick_step Step size for x-axis ticks. """ plt.figure('Infection_states with 50% percentile') - plt.title('Infection states') + plt.title('Infection states with 50% percentile') color_plot = matplotlib.colormaps.get_cmap(colormap).colors states_plot = [1, 2, 3, 4, 5, 7] legend_plot = ['E', 'I_NSymp', 'I_Symp', 'I_Sev', 'I_Crit', 'Dead'] for i in states_plot: plt.plot(x, y50[:, i], color=color_plot[i]) + plt.legend(legend_plot) # Needs to be done here, otherwise the percentage fill_between will not work correctly + + for i in states_plot: plt.fill_between(x, y50[:, i], y25[:, i], alpha=0.5, color=color_plot[i]) plt.fill_between(x, y50[:, i], y75[:, i], alpha=0.5, color=color_plot[i]) @@ -277,14 +275,31 @@ def _format_x_axis(x, start_date, xtick_step): def main(): """ Main function for CLI usage. """ parser = argparse.ArgumentParser(description="Plot infection state and location type results.") - # parser.add_argument("path", type=str, help="Path to the results folder") + parser.add_argument("--path-to-infection-states", help="Path to infection states results") + parser.add_argument("--path-to-loc-types", help="Path to location types results") parser.add_argument("--start-date", type=str, default='2021-03-01', help="Simulation start date (YYYY-MM-DD)") parser.add_argument("--colormap", type=str, default='Set1', help="Matplotlib colormap") parser.add_argument("--xtick-step", type=int, default=150, help="Step for x-axis ticks") args = parser.parse_args() - plot_infection_states_results("/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_state_per_age_group/0", args.start_date, args.colormap, args.xtick_step) - # plot_infections_loc_types_average("/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_per_location_type_per_age_group/0", args.start_date, args.colormap, xtick_step=args.xtick_step) + if args.path_to_infection_states: + plot_infection_states_results( + args.path_to_infection_states, + start_date=args.start_date, + colormap=args.colormap, + xtick_step=args.xtick_step) + if args.path_to_loc_types: + plot_infections_loc_types_average( + args.path_to_loc_types, + start_date=args.start_date, + colormap=args.colormap, + xtick_step=args.xtick_step) + + if not args.path_to_infection_states and not args.path_to_loc_types: + print("Please provide a path to infection states or location types results.") + sys.exit(1) + plt.show() + if __name__ == "__main__": main() From 663a189fa2ef24c6a784b0db11258a1fe35edbd0 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Wed, 28 May 2025 00:24:18 +0200 Subject: [PATCH 05/38] Add unit tests for plotAbmInfectionStates: implement comprehensive test cases for loading H5 results and plotting functions --- .../test_plot_plotAbmInfectionStates.py | 122 ++++++++++++++++++ 1 file changed, 122 insertions(+) create mode 100644 pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py diff --git a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py new file mode 100644 index 0000000000..41d5a31da5 --- /dev/null +++ b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py @@ -0,0 +1,122 @@ +############################################################################# +# Copyright (C) 2020-2024 MEmilio +# +# Authors: Sascha Korf +# +# 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 unittest +from unittest.mock import patch, MagicMock +import numpy as np +import pandas as pd + +import memilio.plot.plotAbmInfectionStates as abm + +class TestPlotAbmInfectionStates(unittest.TestCase): + + @patch('memilio.plot.plotAbmInfectionStates.h5py.File') + def test_load_h5_results(self, mock_h5file): + mock_group = {'Time': np.arange(10), 'Total': np.ones((10, 8))} + mock_h5file().__enter__().get.return_value = {'0': mock_group} + mock_h5file().__enter__().items.return_value = [('0', mock_group)] + mock_h5file().__enter__().__getitem__.return_value = mock_group + mock_h5file().__enter__().items.return_value = [('Time', np.arange(10)), ('Total', np.ones((10, 8)))] + with patch('memilio.plot.plotAbmInfectionStates.h5py.File', mock_h5file): + result = abm.load_h5_results('dummy_path', 'p50') + assert 'Time' in result + assert 'Total' in result + np.testing.assert_array_equal(result['Time'], np.arange(10)) + np.testing.assert_array_equal(result['Total'], np.ones((10, 8))) + + @patch('memilio.plot.plotAbmInfectionStates.load_h5_results') + @patch('memilio.plot.plotAbmInfectionStates.matplotlib') + @patch('memilio.plot.plotAbmInfectionStates.gaussian_filter1d', side_effect=lambda x, sigma, mode: x) + @patch('memilio.plot.plotAbmInfectionStates.pd.DataFrame') + def test_plot_infections_loc_types_average(self, mock_df, mock_gauss, mock_matplotlib, mock_load): + mock_load.return_value = {'Time': np.arange(48), 'Total': np.ones((48, 7))} + mock_df.return_value.rolling.return_value.sum.return_value.to_numpy.return_value = np.ones((48, 1)) + mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1,0,0)]*7 + + # Patch plt.gca().plot to a MagicMock + with patch.object(abm.plt, 'gca') as mock_gca: + mock_ax = MagicMock() + mock_gca.return_value = mock_ax + abm.plot_infections_loc_types_average('dummy_path') + assert mock_ax.plot.called + assert mock_ax.set_xticks.called + assert mock_ax.set_xticklabels.called + + @patch('memilio.plot.plotAbmInfectionStates.load_h5_results') + @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states') + @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states_individual') + def test_plot_infection_states_results(self, mock_indiv, mock_states, mock_load): + mock_load.side_effect = [ + {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones((10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))}, + {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones((10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))}, + {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones((10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))} + ] + abm.plot_infection_states_results('dummy_path') + assert mock_indiv.called + assert mock_states.called + + @patch('memilio.plot.plotAbmInfectionStates.matplotlib') + def test_plot_infection_states(self, mock_matplotlib): + x = np.arange(10) + y50 = np.ones((10, 8)) + y25 = np.zeros((10, 8)) + y75 = np.ones((10, 8))*2 + mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1,0,0)]*8 + + # Patch plt.gca().plot and fill_between + with patch.object(abm.plt, 'gca') as mock_gca: + mock_ax = MagicMock() + mock_gca.return_value = mock_ax + abm.plot_infection_states(x, y50, y25, y75) + assert mock_ax.plot.called + assert mock_ax.fill_between.called + assert mock_ax.set_xticks.called + assert mock_ax.set_xticklabels.called + + @patch('memilio.plot.plotAbmInfectionStates.matplotlib') + def test_plot_infection_states_individual(self, mock_matplotlib): + x = np.arange(10) + group_data = np.ones((10, 8)) + p50_bs = {g: group_data for g in ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} + p25_bs = {g: group_data for g in ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} + p75_bs = {g: group_data for g in ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} + mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1,0,0)]*8 + + # Patch plt.subplots to return a grid of MagicMock axes (as np.array with dtype=object) + with patch.object(abm.plt, 'subplots') as mock_subplots: + fig_mock = MagicMock() + ax_mock = np.empty((6, 7), dtype=object) + for i in range(6): + for j in range(7): + ax_mock[i, j] = MagicMock() + mock_subplots.return_value = (fig_mock, ax_mock) + abm.plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs) + # Check that at least one ax's plot was called + assert any(ax_mock[i, j].plot.called for i in range(6) for j in range(7)) + assert fig_mock.suptitle.called + + def test__format_x_axis(self): + with patch('memilio.plot.plotAbmInfectionStates.plt') as mock_plt: + abm._format_x_axis(np.arange(10), '2021-03-01', 2) + assert mock_plt.gca.called + assert mock_plt.gcf.called + +if __name__ == '__main__': + unittest.main() \ No newline at end of file From 28c3187da683e884bc131415f114e306b34eaa8f Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Wed, 28 May 2025 00:25:02 +0200 Subject: [PATCH 06/38] Refactor setup.py: enhance PylintCommand documentation and update comments for clarity --- pycode/memilio-plot/setup.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pycode/memilio-plot/setup.py b/pycode/memilio-plot/setup.py index 9c6a85908a..b13449bd91 100644 --- a/pycode/memilio-plot/setup.py +++ b/pycode/memilio-plot/setup.py @@ -8,12 +8,13 @@ class PylintCommand(Command): - """Custom command to run pylint and get a report as html.""" + """ + Custom command to run pylint and get a report as html. + """ description = "Runs pylint and outputs the report as html." user_options = [] def initialize_options(self): - """ """ from pylint.reporters.json_reporter import JSONReporter from pylint.reporters.text import ParseableTextReporter, TextReporter from pylint_json2html import JsonExtendedReporter @@ -29,12 +30,10 @@ def initialize_options(self): } def finalize_options(self): - """ """ self.reporter, self.out_file = self.REPORTERS.get( self.out_format) # , self.REPORTERS.get("parseable")) def run(self): - """ """ os.makedirs("build_pylint", exist_ok=True) # Run pylint @@ -63,7 +62,7 @@ def run(self): test_suite='memilio.plot_test', install_requires=[ # smaller pandas versions contain a bug that sometimes prevents reading - # some excel files (e.g. population or mobility data) + # some excel files (e.g. population or twitter data) 'pandas>=1.2.2', 'matplotlib', # smaller numpy versions cause a security issue, 1.25 breaks testing with pyfakefs @@ -74,6 +73,7 @@ def run(self): 'pyxlsb', 'wget', 'folium', + 'scipy.ndimage', 'matplotlib', 'mapclassify', 'geopandas', From 2ba672e8459d087fdd5edeef7c859a91ce65f7b8 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Wed, 28 May 2025 00:26:06 +0200 Subject: [PATCH 07/38] Fix comment in setup.py: update example for Excel file types in install_requires --- pycode/memilio-plot/setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pycode/memilio-plot/setup.py b/pycode/memilio-plot/setup.py index b13449bd91..130996ac6a 100644 --- a/pycode/memilio-plot/setup.py +++ b/pycode/memilio-plot/setup.py @@ -62,7 +62,7 @@ def run(self): test_suite='memilio.plot_test', install_requires=[ # smaller pandas versions contain a bug that sometimes prevents reading - # some excel files (e.g. population or twitter data) + # some excel files (e.g. population or mobility data) 'pandas>=1.2.2', 'matplotlib', # smaller numpy versions cause a security issue, 1.25 breaks testing with pyfakefs From baa272b8ed10db30548423a2edf808731b182f42 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Wed, 28 May 2025 00:27:52 +0200 Subject: [PATCH 08/38] Remove plotAbmICUAndDeadComp.py: delete unused script for ICU and death comparison plots --- .../memilio/plot/plotAbmICUAndDeadComp.py | 266 ------------------ 1 file changed, 266 deletions(-) delete mode 100644 pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py diff --git a/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py b/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py deleted file mode 100644 index 510a88104e..0000000000 --- a/pycode/memilio-plot/memilio/plot/plotAbmICUAndDeadComp.py +++ /dev/null @@ -1,266 +0,0 @@ -# Python script to analyze bs runs -# input is a bs run folder with the following structure: -# bs_run_folder has a txt file for each bs run -# each txt file has a line for each time step -# each line has a column for each compartment as well as the timestep -# each column has the number of individuals in that compartment -# the first line of each txt file is the header - -import sys -import argparse -import os -import pandas as pd -import numpy as np -import matplotlib.pyplot as plt -import matplotlib -import matplotlib.colors as colors -import matplotlib.cm as cmx -import matplotlib.patches as mpatches -import matplotlib.lines as mlines -import h5py -from datetime import datetime -from matplotlib.dates import DateFormatter - -fontsize = 20 - - -def plot_dead(path): - # we will have a seperate plot the cumulative infected individuals, cumulative symptomatic individuals and cumulative dead individual - # we need to load the data - f_p50 = h5py.File( - path+"/infection_state_per_age_group/0/p50/Results.h5", 'r') - p50_bs = f_p50['0'] - - # do the same for 25 and 75 percentile - f_p25 = h5py.File( - path+"/infection_state_per_age_group/0/p25/Results.h5", 'r') - p25_bs = f_p25['0'] - - f_p75 = h5py.File( - path+"/infection_state_per_age_group/0/p75/Results.h5", 'r') - p75_bs = f_p75['0'] - - # do the same for 05 and 95 percentile - f_p05 = h5py.File( - path+"/infection_state_per_age_group/0/p05/Results.h5", 'r') - p05_bs = f_p05['0'] - - f_p95 = h5py.File( - path+"/infection_state_per_age_group/0/p95/Results.h5", 'r') - p95_bs = f_p95['0'] - - age_group_access = ['Group1', 'Group2', 'Group3', - 'Group4', 'Group5', 'Group6', 'Total'] - - # we need the real data json file cases_all_county_age - df_abb = pd.read_json( - path+"/../../../pydata/Germany/cases_all_county_age_ma1.json") - - # we just need the columns cases and date - # we need to offset the dates by 19 day - df_abb['Date'] = df_abb['Date'] + pd.DateOffset(days=18) - # we need just the dates bewteen 2021-03-01 and 2021-06-01 - df_abb = df_abb[(df_abb['Date'] >= '2021-03-01') & - (df_abb['Date'] <= '2021-06-01')] - # we just need the cases with id 3101 - df_abb = df_abb[df_abb['ID_County'] == 3101] - # df_abb['Deaths'] = np.round(df_abb[['Deaths']].to_numpy()) - - # we need the amount of dead persons for each age group: These are A00-A04, A05-A14, A15-A34, A35-A59, A60-A79, A80+ - age_groups = ['A00-A04', 'A05-A14', 'A15-A34', 'A35-A59', 'A60-A79', 'A80+'] - age_grous_string = ['Age 0-4', 'Age 5-14', 'Age 15-34', 'Age 35-59', 'Age 60-79', 'Age 80+'] - # we need to sum up the amount of dead persons for each age group - - # we want the deaths for the age groups - df_abb = df_abb[['Date', 'Deaths', 'Age_RKI']] - # we want a plot with 2 rows. Second row has a plot with each age group and the simulated and real dead persons - # First row has the cumulative dead persons - fig = plt.figure('Deaths') - fig.set_figwidth(20) - fig.set_figheight(9) - gs = fig.add_gridspec(2,6) - - # we need the cumulative dead persons - ax = fig.add_subplot(gs[0, :]) - df_total_dead = df_abb.groupby('Date').sum()[0:90] - y_real = df_total_dead['Deaths'].to_numpy() - # we need to substract the first value from the rest - y_real = y_real - y_real[0] - - y_sim = p50_bs['Total'][()][:, 7][::24][0:90] - y_sim = y_sim - y_sim[0] - - y_sim25 = p25_bs['Total'][()][:, 7][::24][0:90] - y_sim25 = y_sim25 - y_sim25[0] - - y_sim75 = p75_bs['Total'][()][:, 7][::24][0:90] - y_sim75 = y_sim75 - y_sim75[0] - - y_sim05 = p05_bs['Total'][()][:, 7][::24][0:90] - y_sim05 = y_sim05 - y_sim05[0] - - y_sim95 = p95_bs['Total'][()][:, 7][::24][0:90] - y_sim95 = y_sim95 - y_sim95[0] - - - - # we calculate the RMSE - rmse_dead = np.sqrt(((y_real- y_sim)**2).mean()) - # we need to plot the cumulative dead persons from the real world and from the simulation - - ax.plot(df_total_dead.index, y_sim, color='tab:blue',label='Simulated deaths') - ax.plot(df_total_dead.index, y_real, 'v',color='tab:red', linewidth=4, label='Extrapolated deaths from reported infection case data') - ax.fill_between(df_total_dead.index, y_sim75, y_sim25, - alpha=0.5, color='tab:blue', label='50% Confidence interval') - ax.fill_between(df_total_dead.index,y_sim95, y_sim05, - alpha=0.25, color='tab:blue', label='90% Confidence interval') - # ax.text(0.25, 0.8, 'RMSE: '+str(float("{:.2f}".format(rmse_dead))), horizontalalignment='center', - # verticalalignment='center', transform=plt.gca().transAxes, color='pink', fontsize=15) - ax.set_label('Number of individuals') - ax.set_title('Cumulative Deaths', fontsize=fontsize) - ax.set_ylabel('Number of individuals', fontsize=fontsize-8) - ax.legend(fontsize=fontsize-8) - - # now for each age group - for i, age_group in zip(range(6), age_group_access): - ax = fig.add_subplot(gs[1, i]) - # we need the amount of dead persons for each age group - df_abb_age_group = df_abb[df_abb['Age_RKI'] == age_groups[i]][0:90] - y_real = np.round(df_abb_age_group['Deaths'].to_numpy()) - # we need to plot the dead persons from the real world and from the simulation - ax.plot(df_abb_age_group['Date'], y_real-y_real[0], color='tab:red') - ax.plot(df_abb_age_group['Date'], p50_bs[age_group_access[i]][()][:, 7][::24][0:90]-p50_bs[age_group_access[i]][()][:, 7][::24][0], color='tab:blue') - ax.fill_between(df_abb_age_group['Date'], p75_bs[age_group_access[i]][()][:, 7][::24][0:90]-p75_bs[age_group_access[i]][()][:, 7][::24][0], p25_bs[age_group_access[i]][()][:, 7][::24][0:90]-p25_bs[age_group_access[i]][()][:, 7][::24][0], - alpha=0.5, color='tab:blue') - ax.set_title('Deaths, '+age_grous_string[i]) - ax.set_ybound(lower=0) - ax.set_xticks(df_abb_age_group['Date'][::50]) - ax.tick_params(axis='both', which='major', labelsize=fontsize-10) - ax.tick_params(axis='both', which='minor', labelsize=fontsize-10) - if i == 0: - ax.set_ylabel('Number of individuals',fontsize=fontsize-8) - ax.set_ybound(upper=1) - - plt.show() - -def plot_icu(path): - - df_abb = pd.read_json(path+"/../../../pydata/Germany/county_divi.json") - - perc_of_critical_in_icu_age = [0.55,0.55,0.55,0.56,0.54,0.46] - perc_of_critical_in_icu=0.55 - - age_group_access = ['Group1', 'Group2', 'Group3', - 'Group4', 'Group5', 'Group6', 'Total'] - - - # we just need the columns ICU_low and ICU_hig - df_abb = df_abb[['ID_County', 'ICU', 'Date']] - - df_abb = df_abb[df_abb['ID_County'] == 3101] - # we need just the dates bewteen 2021-03-01 and 2021-06-01 - df_abb = df_abb[(df_abb['Date'] >= '2021-03-01') & - (df_abb['Date'] <= '2021-06-01')] - - # we plot this against this the Amount of persons in the ICU from our model - f_p50 = h5py.File( - path+"/infection_state_per_age_group/0/p50/Results.h5", 'r') - total_50 = f_p50['0']['Total'][()][::24][0:90] - - total_50_age = f_p50['0'][age_group_access[0]][()] - for i in range(6): - total_50_age += f_p50['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] - total_50_age = total_50_age[::24][0:90] - - - # we plot this against this the Amount of persons in the ICU from our model - f_p75 = h5py.File( - path+"/infection_state_per_age_group/0/p75/Results.h5", 'r') - # total_75 = f_p75['0']['Total'][()][::24][0:90] - total_75_age = f_p75['0'][age_group_access[0]][()] - for i in range(6): - total_75_age += f_p75['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] - total_75_age = total_75_age[::24][0:90] - - # same with 25 percentile - f_p25 = h5py.File( - path+"/infection_state_per_age_group/0/p25/Results.h5", 'r') - # total_25 = f_p25['0']['Total'][()][::24][0:90] - total_25_age = f_p25['0'][age_group_access[0]][()] - for i in range(6): - total_25_age += f_p25['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] - total_25_age = total_25_age[::24][0:90] - - # same with 05 and 95 percentile - f_p05 = h5py.File( - path+"/infection_state_per_age_group/0/p05/Results.h5", 'r') - # total_05 = f_p05['0']['Total'][()][::24][0:90] - total_05_age = f_p05['0'][age_group_access[0]][()] - for i in range(6): - total_05_age += f_p05['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] - total_05_age = total_05_age[::24][0:90] - - f_p95 = h5py.File( - path+"/infection_state_per_age_group/0/p95/Results.h5", 'r') - # total_95 = f_p95['0']['Total'][()][::24][0:90] - total_95_age = f_p95['0'][age_group_access[0]][()] - for i in range(6): - total_95_age += f_p95['0'][age_group_access[i]][()]*perc_of_critical_in_icu_age[i] - total_95_age = total_95_age[::24][0:90] - - - ICU_Simulation_one_percentile = np.floor(total_50[:, 5]*perc_of_critical_in_icu) - ICU_Simulation = np.round(total_50_age[:, 5]) - ICU_Simulation75 = np.round(total_75_age[:, 5]) - ICU_Simulation25 = np.round(total_25_age[:, 5]) - ICU_Simulation05 = np.round(total_05_age[:, 5]) - ICU_Simulation95 = np.round(total_95_age[:, 5]) - ICU_Real = df_abb['ICU'][0:90] - - #smooth the data - # ICU_Real = gaussian_filter1d(ICU_Real, sigma=1, mode='nearest') - # ICU_Simulation = gaussian_filter1d(ICU_Simulation, sigma=1, mode='nearest') - - - - # we calculate the RMSE - rmse_ICU = np.sqrt(((ICU_Real - ICU_Simulation_one_percentile)**2).mean()) - - # plot the ICU beds and the ICU beds taken - fig, ax = plt.subplots(1, 1, constrained_layout=True) - fig.set_figwidth(12) - fig.set_figheight(9) - # we plot the ICU_low and the ICU_high - ax.plot(df_abb['Date'][0:90], ICU_Real,'x', color='tab:red', linewidth=10, label='Data') - ax.plot(df_abb['Date'][0:90], ICU_Simulation, color='tab:blue', label='Simulation') - # ax.plot(df_abb['Date'][0:90], ICU_Simulation_one_percentile, color='tab:green', label='Simulated ICU beds') - ax.fill_between(df_abb['Date'][0:90],ICU_Simulation75, ICU_Simulation25, - alpha=0.5, color='tab:blue', label='50% Confidence interval') - ax.fill_between(df_abb['Date'][0:90],ICU_Simulation05, ICU_Simulation95, - alpha=0.25, color='tab:blue', label='90% Confidence interval') - - - # we also write the rmse - # ax.text(0.25, 0.8, 'RMSE: '+str(float("{:.2f}".format(rmse_ICU))), horizontalalignment='center', - # verticalalignment='center', transform=plt.gca().transAxes, color='pink', fontsize=15) - ax.tick_params(axis='both', which='major', labelsize=fontsize-4) - ax.tick_params(axis='both', which='minor', labelsize=fontsize-4) - ax.set_ylabel('Occupied ICU beds', fontsize=fontsize) - ax.set_title('ICU beds', fontsize=fontsize+4) - ax.legend(fontsize=fontsize-4) - plt.show() - - - - -if __name__ == "__main__": - path = "" - - if (len(sys.argv) > 1): - n_runs = sys.argv[1] - else: - n_runs = len([entry for entry in os.listdir(path) - if os.path.isfile(os.path.join(path, entry))]) - - # plot_icu(path) - # plot_dead(path) From 20b651e8a99098a327fdf2550572951855ac933d Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Wed, 28 May 2025 00:34:25 +0200 Subject: [PATCH 09/38] formatting --- .../memilio/plot/plotAbmInfectionStates.py | 105 ++++++++++++------ .../test_plot_plotAbmInfectionStates.py | 40 ++++--- 2 files changed, 94 insertions(+), 51 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 91a6926c5e..00f2d749b3 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -37,19 +37,19 @@ # The used Loggers are: # struct LogInfectionStatePerAgeGroup : mio::LogAlways { # using Type = std::pair; -# /** +# /** # * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. # * @param[in] sim The simulation of the abm. # * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. # */ # static Type log(const mio::abm::Simulation& sim) # { -# +# # Eigen::VectorXd sum = Eigen::VectorXd::Zero( # Eigen::Index((size_t)mio::abm::InfectionState::Count * sim.get_world().parameters.get_num_groups())); # const auto curr_time = sim.get_time(); # const auto persons = sim.get_world().get_persons(); -# +# # // PRAGMA_OMP(parallel for) # for (auto i = size_t(0); i < persons.size(); ++i) { # auto& p = persons[i]; @@ -63,23 +63,23 @@ # return std::make_pair(curr_time, sum); # } # }; -# +# # struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { # using Type = std::pair; -# /** +# /** # * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. # * @param[in] sim The simulation of the abm. # * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. # */ # static Type log(const mio::abm::Simulation& sim) # { -# +# # Eigen::VectorXd sum = Eigen::VectorXd::Zero( # Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_world().parameters.get_num_groups())); # auto curr_time = sim.get_time(); # auto prev_time = sim.get_prev_time(); # const auto persons = sim.get_world().get_persons(); -# +# # // PRAGMA_OMP(parallel for) # for (auto i = size_t(0); i < persons.size(); ++i) { # auto& p = persons[i]; @@ -96,7 +96,7 @@ # return std::make_pair(curr_time, sum); # } # }; -# +# # The output of the loggers of several runs is stored in HDF5 files, with the memilio funciton mio::save_results in mio/io/result_io.h. @@ -112,6 +112,7 @@ def load_h5_results(base_path, percentile): data = {k: v[()] for k, v in f['0'].items()} return data + def plot_infections_loc_types_average( path_to_loc_types, start_date='2021-03-01', @@ -134,17 +135,22 @@ def plot_infections_loc_types_average( total_50 = p50['Total'] plt.figure('Infection_location_types') - plt.title('Infection per location type for the median run, rolling sum over 24 hours') + plt.title( + 'Infection per location type for the median run, rolling sum over 24 hours') color_plot = matplotlib.colormaps.get_cmap(colormap).colors # If you define further location types, you need to adjust this list - states_plot = [0, 1, 2, 3, 4, 5 , 6] - legend_plot = ['Home', 'School', 'Work', 'SocialEvent', 'BasicsShop', 'Hospital', 'ICU'] + states_plot = [0, 1, 2, 3, 4, 5, 6] + legend_plot = ['Home', 'School', 'Work', + 'SocialEvent', 'BasicsShop', 'Hospital', 'ICU'] for idx, i in enumerate(states_plot): - color = color_plot[i % len(color_plot)] if i < len(color_plot) else "black" + color = color_plot[i % len(color_plot)] if i < len( + color_plot) else "black" # Sum up every 24 hours, then smooth - indexer = pd.api.indexers.FixedForwardWindowIndexer(window_size=rolling_window) - y = pd.DataFrame(total_50[:, i]).rolling(window=indexer, min_periods=1).sum().to_numpy() + indexer = pd.api.indexers.FixedForwardWindowIndexer( + window_size=rolling_window) + y = pd.DataFrame(total_50[:, i]).rolling( + window=indexer, min_periods=1).sum().to_numpy() y = y[0::rolling_window].flatten() y = gaussian_filter1d(y, sigma=smooth_sigma, mode='nearest') plt.plot(time[0::rolling_window], y, color=color) @@ -155,6 +161,7 @@ def plot_infections_loc_types_average( plt.ylabel('Number of individuals') plt.show() + def plot_infection_states_results( path_to_infection_states, start_date='2021-03-01', @@ -177,7 +184,9 @@ def plot_infection_states_results( total_75 = p75['Total'] plot_infection_states_individual(time, p50, p25, p75, colormap) - plot_infection_states(time, total_50, total_25, total_75, start_date, colormap, xtick_step) + plot_infection_states(time, total_50, total_25, + total_75, start_date, colormap, xtick_step) + def plot_infection_states( x, y50, y25, y75, @@ -202,11 +211,14 @@ def plot_infection_states( for i in states_plot: plt.plot(x, y50[:, i], color=color_plot[i]) - plt.legend(legend_plot) # Needs to be done here, otherwise the percentage fill_between will not work correctly + # Needs to be done here, otherwise the percentage fill_between will not work correctly + plt.legend(legend_plot) for i in states_plot: - plt.fill_between(x, y50[:, i], y25[:, i], alpha=0.5, color=color_plot[i]) - plt.fill_between(x, y50[:, i], y75[:, i], alpha=0.5, color=color_plot[i]) + plt.fill_between(x, y50[:, i], y25[:, i], + alpha=0.5, color=color_plot[i]) + plt.fill_between(x, y50[:, i], y75[:, i], + alpha=0.5, color=color_plot[i]) plt.legend(legend_plot) _format_x_axis(x, start_date, xtick_step) @@ -214,6 +226,7 @@ def plot_infection_states( plt.ylabel('Number of individuals') plt.show() + def plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs, colormap='Set1'): """ Plots infection states for each age group. @@ -223,28 +236,37 @@ def plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs, colormap='Set1') @param[in] p75_bs 75th percentile values by group. @param[in] colormap Matplotlib colormap. """ - age_groups = ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total'] # Adjust as needed + age_groups = ['Group1', 'Group2', 'Group3', 'Group4', + 'Group5', 'Group6', 'Total'] # Adjust as needed color_plot = matplotlib.colormaps.get_cmap(colormap).colors - fig, ax = plt.subplots(6, len(age_groups), constrained_layout=True, figsize=(20, 9)) + fig, ax = plt.subplots( + 6, len(age_groups), constrained_layout=True, figsize=(20, 9)) for col_idx, group in enumerate(age_groups): y50 = p50_bs[group] y25 = p25_bs[group] y75 = p75_bs[group] # Infected no symptoms - _plot_state(ax[0, col_idx], x, y50[:, 1], y25[:, 1], y75[:, 1], color_plot[col_idx], '#Infected_no_symptoms, Age' + str(group)) + _plot_state(ax[0, col_idx], x, y50[:, 1], y25[:, 1], y75[:, 1], + color_plot[col_idx], '#Infected_no_symptoms, Age' + str(group)) # Infected symptoms - _plot_state(ax[1, col_idx], x, y50[:, 2], y25[:, 2], y75[:, 2], color_plot[col_idx], '#Infected_symptoms, Age' + str(group)) + _plot_state(ax[1, col_idx], x, y50[:, 2], y25[:, 2], y75[:, 2], + color_plot[col_idx], '#Infected_symptoms, Age' + str(group)) # Severe - _plot_state(ax[2, col_idx], x, y50[:, 4], y25[:, 4], y75[:, 4], color_plot[col_idx], '#Severe, Age' + str(group)) + _plot_state(ax[2, col_idx], x, y50[:, 4], y25[:, 4], y75[:, 4], + color_plot[col_idx], '#Severe, Age' + str(group)) # Critical - _plot_state(ax[3, col_idx], x, y50[:, 5], y25[:, 5], y75[:, 5], color_plot[col_idx], '#Critical, Age' + str(group)) + _plot_state(ax[3, col_idx], x, y50[:, 5], y25[:, 5], y75[:, 5], + color_plot[col_idx], '#Critical, Age' + str(group)) # Dead - _plot_state(ax[4, col_idx], x, y50[:, 7], y25[:, 7], y75[:, 7], color_plot[col_idx], '#Dead, Age' + str(group)) + _plot_state(ax[4, col_idx], x, y50[:, 7], y25[:, 7], + y75[:, 7], color_plot[col_idx], '#Dead, Age' + str(group)) # Recovered - _plot_state(ax[5, col_idx], x, y50[:, 6], y25[:, 6], y75[:, 6], color_plot[col_idx], '#Recovered, Age' + str(group)) - - fig.suptitle('Infection states per age group with 50% percentile', fontsize=16) + _plot_state(ax[5, col_idx], x, y50[:, 6], y25[:, 6], y75[:, 6], + color_plot[col_idx], '#Recovered, Age' + str(group)) + + fig.suptitle( + 'Infection states per age group with 50% percentile', fontsize=16) # We hide the Legend for the individual plots as it is too cluttered for ax_row in ax: @@ -253,6 +275,7 @@ def plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs, colormap='Set1') plt.show() + def _plot_state(ax, x, y50, y25, y75, color, title): """ Helper to plot a single state with fill_between. """ ax.set_xlabel('time (days)') @@ -263,6 +286,7 @@ def _plot_state(ax, x, y50, y25, y75, color, title): ax.set_title(title) ax.legend(['Simulation']) + def _format_x_axis(x, start_date, xtick_step): """ Helper to format x-axis as dates. """ start = datetime.strptime(start_date, '%Y-%m-%d') @@ -272,17 +296,24 @@ def _format_x_axis(x, start_date, xtick_step): plt.gca().set_xticklabels(xx_str[::xtick_step]) plt.gcf().autofmt_xdate() + def main(): """ Main function for CLI usage. """ - parser = argparse.ArgumentParser(description="Plot infection state and location type results.") - parser.add_argument("--path-to-infection-states", help="Path to infection states results") - parser.add_argument("--path-to-loc-types", help="Path to location types results") - parser.add_argument("--start-date", type=str, default='2021-03-01', help="Simulation start date (YYYY-MM-DD)") - parser.add_argument("--colormap", type=str, default='Set1', help="Matplotlib colormap") - parser.add_argument("--xtick-step", type=int, default=150, help="Step for x-axis ticks") + parser = argparse.ArgumentParser( + description="Plot infection state and location type results.") + parser.add_argument("--path-to-infection-states", + help="Path to infection states results") + parser.add_argument("--path-to-loc-types", + help="Path to location types results") + parser.add_argument("--start-date", type=str, default='2021-03-01', + help="Simulation start date (YYYY-MM-DD)") + parser.add_argument("--colormap", type=str, + default='Set1', help="Matplotlib colormap") + parser.add_argument("--xtick-step", type=int, + default=150, help="Step for x-axis ticks") args = parser.parse_args() - if args.path_to_infection_states: + if args.path_to_infection_states: plot_infection_states_results( args.path_to_infection_states, start_date=args.start_date, @@ -294,12 +325,12 @@ def main(): start_date=args.start_date, colormap=args.colormap, xtick_step=args.xtick_step) - + if not args.path_to_infection_states and not args.path_to_loc_types: print("Please provide a path to infection states or location types results.") sys.exit(1) plt.show() - + if __name__ == "__main__": main() diff --git a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py index 41d5a31da5..bb04cc54a9 100644 --- a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py @@ -25,6 +25,7 @@ import memilio.plot.plotAbmInfectionStates as abm + class TestPlotAbmInfectionStates(unittest.TestCase): @patch('memilio.plot.plotAbmInfectionStates.h5py.File') @@ -33,7 +34,8 @@ def test_load_h5_results(self, mock_h5file): mock_h5file().__enter__().get.return_value = {'0': mock_group} mock_h5file().__enter__().items.return_value = [('0', mock_group)] mock_h5file().__enter__().__getitem__.return_value = mock_group - mock_h5file().__enter__().items.return_value = [('Time', np.arange(10)), ('Total', np.ones((10, 8)))] + mock_h5file().__enter__().items.return_value = [ + ('Time', np.arange(10)), ('Total', np.ones((10, 8)))] with patch('memilio.plot.plotAbmInfectionStates.h5py.File', mock_h5file): result = abm.load_h5_results('dummy_path', 'p50') assert 'Time' in result @@ -46,9 +48,11 @@ def test_load_h5_results(self, mock_h5file): @patch('memilio.plot.plotAbmInfectionStates.gaussian_filter1d', side_effect=lambda x, sigma, mode: x) @patch('memilio.plot.plotAbmInfectionStates.pd.DataFrame') def test_plot_infections_loc_types_average(self, mock_df, mock_gauss, mock_matplotlib, mock_load): - mock_load.return_value = {'Time': np.arange(48), 'Total': np.ones((48, 7))} - mock_df.return_value.rolling.return_value.sum.return_value.to_numpy.return_value = np.ones((48, 1)) - mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1,0,0)]*7 + mock_load.return_value = { + 'Time': np.arange(48), 'Total': np.ones((48, 7))} + mock_df.return_value.rolling.return_value.sum.return_value.to_numpy.return_value = np.ones( + (48, 1)) + mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1, 0, 0)]*7 # Patch plt.gca().plot to a MagicMock with patch.object(abm.plt, 'gca') as mock_gca: @@ -64,9 +68,12 @@ def test_plot_infections_loc_types_average(self, mock_df, mock_gauss, mock_matpl @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states_individual') def test_plot_infection_states_results(self, mock_indiv, mock_states, mock_load): mock_load.side_effect = [ - {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones((10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))}, - {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones((10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))}, - {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones((10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))} + {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones( + (10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))}, + {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones( + (10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))}, + {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones( + (10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))} ] abm.plot_infection_states_results('dummy_path') assert mock_indiv.called @@ -78,7 +85,7 @@ def test_plot_infection_states(self, mock_matplotlib): y50 = np.ones((10, 8)) y25 = np.zeros((10, 8)) y75 = np.ones((10, 8))*2 - mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1,0,0)]*8 + mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1, 0, 0)]*8 # Patch plt.gca().plot and fill_between with patch.object(abm.plt, 'gca') as mock_gca: @@ -94,10 +101,13 @@ def test_plot_infection_states(self, mock_matplotlib): def test_plot_infection_states_individual(self, mock_matplotlib): x = np.arange(10) group_data = np.ones((10, 8)) - p50_bs = {g: group_data for g in ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} - p25_bs = {g: group_data for g in ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} - p75_bs = {g: group_data for g in ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} - mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1,0,0)]*8 + p50_bs = {g: group_data for g in [ + 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} + p25_bs = {g: group_data for g in [ + 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} + p75_bs = {g: group_data for g in [ + 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} + mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1, 0, 0)]*8 # Patch plt.subplots to return a grid of MagicMock axes (as np.array with dtype=object) with patch.object(abm.plt, 'subplots') as mock_subplots: @@ -109,7 +119,8 @@ def test_plot_infection_states_individual(self, mock_matplotlib): mock_subplots.return_value = (fig_mock, ax_mock) abm.plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs) # Check that at least one ax's plot was called - assert any(ax_mock[i, j].plot.called for i in range(6) for j in range(7)) + assert any(ax_mock[i, j].plot.called for i in range(6) + for j in range(7)) assert fig_mock.suptitle.called def test__format_x_axis(self): @@ -118,5 +129,6 @@ def test__format_x_axis(self): assert mock_plt.gca.called assert mock_plt.gcf.called + if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main() From 706aea060ca0519a05c6c62adbcaf384c793eae7 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Wed, 28 May 2025 20:47:21 +0200 Subject: [PATCH 10/38] Implement martins suggestions --- .../memilio/plot/plotAbmInfectionStates.py | 228 +++++++++++------- .../test_plot_plotAbmInfectionStates.py | 24 +- 2 files changed, 160 insertions(+), 92 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 00f2d749b3..0ff71347dd 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -99,6 +99,40 @@ # # The output of the loggers of several runs is stored in HDF5 files, with the memilio funciton mio::save_results in mio/io/result_io.h. +# Adjust these as needed. +# States and location type numbers need to match the infection states used in your simulation. +state_labels = { + 1: 'Exposed', + 2: 'I_Asymp', + 3: 'I_Symp', + 4: 'I_Severe', + 5: 'I_Critical', + 7: 'Dead' +} + +age_groups = ['Group1', 'Group2', 'Group3', 'Group4', + 'Group5', 'Group6', 'Total'] + +age_groups_dict = { + 'Group1': 'Ages 0-4', + 'Group2': 'Ages 5-14', + 'Group3': 'Ages 15-34', + 'Group4': 'Ages 35-59', + 'Group5': 'Ages 60-79', + 'Group6': 'Ages 80+', + 'Total': 'All Ages' +} + +location_type_labels = { + 0: 'Home', + 1: 'School', + 2: 'Work', + 3: 'SocialEvent', + 4: 'BasicsShop', + 5: 'Hospital', + 6: 'ICU' +} + def load_h5_results(base_path, percentile): """ Reads HDF5 results for a given group and percentile. @@ -138,12 +172,8 @@ def plot_infections_loc_types_average( plt.title( 'Infection per location type for the median run, rolling sum over 24 hours') color_plot = matplotlib.colormaps.get_cmap(colormap).colors - # If you define further location types, you need to adjust this list - states_plot = [0, 1, 2, 3, 4, 5, 6] - legend_plot = ['Home', 'School', 'Work', - 'SocialEvent', 'BasicsShop', 'Hospital', 'ICU'] - for idx, i in enumerate(states_plot): + for idx, i in enumerate(location_type_labels.keys()): color = color_plot[i % len(color_plot)] if i < len( color_plot) else "black" # Sum up every 24 hours, then smooth @@ -153,9 +183,9 @@ def plot_infections_loc_types_average( window=indexer, min_periods=1).sum().to_numpy() y = y[0::rolling_window].flatten() y = gaussian_filter1d(y, sigma=smooth_sigma, mode='nearest') - plt.plot(time[0::rolling_window], y, color=color) + plt.plot(time[0::rolling_window], y, color=color, linewidth=2.5) - plt.legend(legend_plot) + plt.legend(list(location_type_labels.values())) _format_x_axis(time, start_date, xtick_step) plt.xlabel('Date') plt.ylabel('Number of individuals') @@ -166,14 +196,10 @@ def plot_infection_states_results( path_to_infection_states, start_date='2021-03-01', colormap='Set1', - xtick_step=150): - """ Loads and plots infection state results. - - @param[in] base_path Path to results directory. - @param[in] start_date Start date as string. - @param[in] colormap Matplotlib colormap. - @param[in] xtick_step Step size for x-axis ticks. - """ + xtick_step=150, + show90=False +): + """ Loads and plots infection state results. """ # Load data p50 = load_h5_results(path_to_infection_states, "p50") p25 = load_h5_results(path_to_infection_states, "p25") @@ -182,109 +208,126 @@ def plot_infection_states_results( total_50 = p50['Total'] total_25 = p25['Total'] total_75 = p75['Total'] - - plot_infection_states_individual(time, p50, p25, p75, colormap) + p05 = p95 = None + total_05 = total_95 = None + if show90: + total_95 = load_h5_results(path_to_infection_states, "p95") + total_05 = load_h5_results(path_to_infection_states, "p05") + p95 = total_95['Total'] + p05 = total_05['Total'] + + plot_infection_states_individual( + time, p50, p25, p75, colormap, + p05_bs=total_05 if show90 else None, + p95_bs=total_95 if show90 else None, + show90=show90 + ) plot_infection_states(time, total_50, total_25, - total_75, start_date, colormap, xtick_step) + total_75, start_date, colormap, xtick_step, + y05=p05, y95=p95, show_90=show90) def plot_infection_states( x, y50, y25, y75, start_date='2021-03-01', colormap='Set1', - xtick_step=150): - """ Plots infection states with percentiles. + xtick_step=150, + y05=None, y95=None, show_90=False): + """ Plots infection states with percentiles and improved styling. """ + plt.figure('Infection_states') - @param[in] x Time array. - @param[in] y50 Median values. - @param[in] y25 25th percentile values. - @param[in] y75 75th percentile values. - @param[in] start_date Start date as string. - @param[in] colormap Matplotlib colormap. - @param[in] xtick_step Step size for x-axis ticks. - """ - plt.figure('Infection_states with 50% percentile') plt.title('Infection states with 50% percentile') + if show_90: + plt.title('Infection states with 50% and 90% percentiles') + color_plot = matplotlib.colormaps.get_cmap(colormap).colors - states_plot = [1, 2, 3, 4, 5, 7] - legend_plot = ['E', 'I_NSymp', 'I_Symp', 'I_Sev', 'I_Crit', 'Dead'] - for i in states_plot: - plt.plot(x, y50[:, i], color=color_plot[i]) - # Needs to be done here, otherwise the percentage fill_between will not work correctly - plt.legend(legend_plot) + states_plot = list(state_labels.keys()) for i in states_plot: + plt.plot(x, y50[:, i], color=color_plot[i], + linewidth=2.5, label=state_labels[i]) + # needs to be after the plot calls + plt.legend([state_labels[i] for i in states_plot]) + for i in states_plot: + plt.plot(x, y25[:, i], color=color_plot[i], + linestyle='dashdot', linewidth=1.2, alpha=0.7) + plt.plot(x, y75[:, i], color=color_plot[i], + linestyle='dashdot', linewidth=1.2, alpha=0.7) plt.fill_between(x, y50[:, i], y25[:, i], - alpha=0.5, color=color_plot[i]) + alpha=0.2, color=color_plot[i]) plt.fill_between(x, y50[:, i], y75[:, i], - alpha=0.5, color=color_plot[i]) + alpha=0.2, color=color_plot[i]) + # Optional: 90% percentile + if show_90 and y05 is not None and y95 is not None: + plt.plot(x, y05[:, i], color=color_plot[i], + linestyle='dashdot', linewidth=1.0, alpha=0.4) + plt.plot(x, y95[:, i], color=color_plot[i], + linestyle='dashdot', linewidth=1.0, alpha=0.4) + plt.fill_between(x, y05[:, i], y95[:, i], + # More transparent + alpha=0.25, color=color_plot[i]) - plt.legend(legend_plot) _format_x_axis(x, start_date, xtick_step) - plt.xlabel('Time') + plt.xlabel('Date') plt.ylabel('Number of individuals') plt.show() -def plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs, colormap='Set1'): - """ Plots infection states for each age group. +def plot_infection_states_individual( + x, p50_bs, p25_bs, p75_bs, colormap='Set1', + p05_bs=None, p95_bs=None, show90=False +): + """ Plots infection states for each age group, with optional 90% percentile. """ - @param[in] x Time array. - @param[in] p50_bs Median values by group. - @param[in] p25_bs 25th percentile values by group. - @param[in] p75_bs 75th percentile values by group. - @param[in] colormap Matplotlib colormap. - """ - age_groups = ['Group1', 'Group2', 'Group3', 'Group4', - 'Group5', 'Group6', 'Total'] # Adjust as needed color_plot = matplotlib.colormaps.get_cmap(colormap).colors + n_states = len(state_labels) fig, ax = plt.subplots( - 6, len(age_groups), constrained_layout=True, figsize=(20, 9)) + n_states, len(age_groups), constrained_layout=True, figsize=(20, 3 * n_states)) for col_idx, group in enumerate(age_groups): y50 = p50_bs[group] y25 = p25_bs[group] y75 = p75_bs[group] - # Infected no symptoms - _plot_state(ax[0, col_idx], x, y50[:, 1], y25[:, 1], y75[:, 1], - color_plot[col_idx], '#Infected_no_symptoms, Age' + str(group)) - # Infected symptoms - _plot_state(ax[1, col_idx], x, y50[:, 2], y25[:, 2], y75[:, 2], - color_plot[col_idx], '#Infected_symptoms, Age' + str(group)) - # Severe - _plot_state(ax[2, col_idx], x, y50[:, 4], y25[:, 4], y75[:, 4], - color_plot[col_idx], '#Severe, Age' + str(group)) - # Critical - _plot_state(ax[3, col_idx], x, y50[:, 5], y25[:, 5], y75[:, 5], - color_plot[col_idx], '#Critical, Age' + str(group)) - # Dead - _plot_state(ax[4, col_idx], x, y50[:, 7], y25[:, 7], - y75[:, 7], color_plot[col_idx], '#Dead, Age' + str(group)) - # Recovered - _plot_state(ax[5, col_idx], x, y50[:, 6], y25[:, 6], y75[:, 6], - color_plot[col_idx], '#Recovered, Age' + str(group)) - + y05 = p05_bs[group] if (show90 and p05_bs is not None) else None + y95 = p95_bs[group] if (show90 and p95_bs is not None) else None + for row_idx, (state_idx, label) in enumerate(state_labels.items()): + _plot_state( + ax[row_idx, col_idx], x, y50[:, state_idx], y25[:, + state_idx], y75[:, state_idx], + color_plot[col_idx], f'#{label}, {age_groups_dict[group]}', + y05=y05[:, state_idx] if y05 is not None else None, + y95=y95[:, state_idx] if y95 is not None else None, + show90=show90 + ) + # The legend should say: solid line = median, dashed line = 25% and 75% perc. and if show90 is True, dotted line = 5%, 25%, 75%, 95% perc. + perc_string = '25/75%' if not show90 else '5/25/75/95%' + ax[row_idx, col_idx].legend( + ['Median', f'{perc_string} perc.'], + loc='upper left', fontsize=8) + + string_short = ' and 90%' if show90 else '' fig.suptitle( - 'Infection states per age group with 50% percentile', fontsize=16) - - # We hide the Legend for the individual plots as it is too cluttered - for ax_row in ax: - for ax_col in ax_row: - ax_col.legend().set_visible(False) + 'Infection states per age group with 50' + string_short + ' percentile', + fontsize=16) plt.show() -def _plot_state(ax, x, y50, y25, y75, color, title): - """ Helper to plot a single state with fill_between. """ +def _plot_state(ax, x, y50, y25, y75, color, title, y05=None, y95=None, show90=False): + """ Helper to plot a single state with fill_between and optional 90% percentile. """ ax.set_xlabel('time (days)') ax.plot(x, y50, color=color, label='Median') ax.fill_between(x, y50, y25, alpha=0.5, color=color) ax.fill_between(x, y50, y75, alpha=0.5, color=color) + if show90 and y05 is not None and y95 is not None: + ax.plot(x, y05, color=color, linestyle='dotted', + linewidth=1.0, alpha=0.4) + ax.plot(x, y95, color=color, linestyle='dotted', + linewidth=1.0, alpha=0.4) + ax.fill_between(x, y05, y95, alpha=0.15, color=color) ax.tick_params(axis='y') ax.set_title(title) - ax.legend(['Simulation']) def _format_x_axis(x, start_date, xtick_step): @@ -311,20 +354,25 @@ def main(): default='Set1', help="Matplotlib colormap") parser.add_argument("--xtick-step", type=int, default=150, help="Step for x-axis ticks") + parser.add_argument("--90percentile", action="store_true", + help="If set, plot 90% percentile as well") args = parser.parse_args() - if args.path_to_infection_states: - plot_infection_states_results( - args.path_to_infection_states, - start_date=args.start_date, - colormap=args.colormap, - xtick_step=args.xtick_step) - if args.path_to_loc_types: - plot_infections_loc_types_average( - args.path_to_loc_types, - start_date=args.start_date, - colormap=args.colormap, - xtick_step=args.xtick_step) + path_to_infection_states = "/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_per_location_type_per_age_group/0" + path_to_loc_types = "/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_state_per_age_group/0" + + plot_infection_states_results( + path_to_loc_types, + start_date=args.start_date, + colormap=args.colormap, + xtick_step=args.xtick_step, + show90=True + ) + plot_infections_loc_types_average( + path_to_infection_states, + start_date=args.start_date, + colormap=args.colormap, + xtick_step=args.xtick_step) if not args.path_to_infection_states and not args.path_to_loc_types: print("Please provide a path to infection states or location types results.") diff --git a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py index bb04cc54a9..fbdba48535 100644 --- a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py @@ -85,13 +85,23 @@ def test_plot_infection_states(self, mock_matplotlib): y50 = np.ones((10, 8)) y25 = np.zeros((10, 8)) y75 = np.ones((10, 8))*2 + y05 = np.ones((10, 8))*-1 + y95 = np.ones((10, 8))*3 mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1, 0, 0)]*8 # Patch plt.gca().plot and fill_between with patch.object(abm.plt, 'gca') as mock_gca: mock_ax = MagicMock() mock_gca.return_value = mock_ax - abm.plot_infection_states(x, y50, y25, y75) + abm.plot_infection_states( + x, y50, y25, y75, + start_date='2021-03-01', + colormap='Set1', + xtick_step=2, + y05=y05, + y95=y95, + show_90=True + ) assert mock_ax.plot.called assert mock_ax.fill_between.called assert mock_ax.set_xticks.called @@ -107,6 +117,10 @@ def test_plot_infection_states_individual(self, mock_matplotlib): 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} p75_bs = {g: group_data for g in [ 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} + p05_bs = {g: group_data*-1 for g in [ + 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} + p95_bs = {g: group_data*3 for g in [ + 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1, 0, 0)]*8 # Patch plt.subplots to return a grid of MagicMock axes (as np.array with dtype=object) @@ -117,7 +131,13 @@ def test_plot_infection_states_individual(self, mock_matplotlib): for j in range(7): ax_mock[i, j] = MagicMock() mock_subplots.return_value = (fig_mock, ax_mock) - abm.plot_infection_states_individual(x, p50_bs, p25_bs, p75_bs) + abm.plot_infection_states_individual( + x, p50_bs, p25_bs, p75_bs, + colormap='Set1', + p05_bs=p05_bs, + p95_bs=p95_bs, + show90=True + ) # Check that at least one ax's plot was called assert any(ax_mock[i, j].plot.called for i in range(6) for j in range(7)) From ab2e0e8a572f7c4517f35031d3db4217ecfa1162 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Wed, 28 May 2025 20:48:23 +0200 Subject: [PATCH 11/38] Refactor plotAbmInfectionStates.py: update paths to use command line arguments for better flexibility --- .../memilio-plot/memilio/plot/plotAbmInfectionStates.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 0ff71347dd..cf2c0414d4 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -357,19 +357,15 @@ def main(): parser.add_argument("--90percentile", action="store_true", help="If set, plot 90% percentile as well") args = parser.parse_args() - - path_to_infection_states = "/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_per_location_type_per_age_group/0" - path_to_loc_types = "/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_state_per_age_group/0" - plot_infection_states_results( - path_to_loc_types, + args.path_to_infection_states, start_date=args.start_date, colormap=args.colormap, xtick_step=args.xtick_step, show90=True ) plot_infections_loc_types_average( - path_to_infection_states, + args.path_to_loc_types, start_date=args.start_date, colormap=args.colormap, xtick_step=args.xtick_step) From e4df6cb25a6c5e00c3e3e2da70b6df2bea220ae3 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Sat, 14 Jun 2025 12:02:28 +0200 Subject: [PATCH 12/38] Julias suggestions p1 --- .../memilio/plot/plotAbmInfectionStates.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index cf2c0414d4..a3d90ddd90 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -1,5 +1,5 @@ ############################################################################# -# Copyright (C) 2020-2024 MEmilio +# Copyright (C) 2020-2025 MEmilio # # Authors: Sascha Korf # @@ -30,17 +30,17 @@ from scipy.ndimage import gaussian_filter1d -# Module for plotting infection states and location types from ABM results. +# Module for plotting number of agents per infection state and number of infected agents per location type from ABM results. # This module provides functions to load and visualize infection states and -# location types from simulation results stored in HDF5 format and are output -# by the MEmilio agent-based model (ABM). +# location types from simulation results of the agent-based model (ABM) stored in HDF5 format. + # The used Loggers are: # struct LogInfectionStatePerAgeGroup : mio::LogAlways { # using Type = std::pair; # /** -# * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. +# * @brief Log the TimeSeries of the number of Person%s in an #InfectionState for every age group. # * @param[in] sim The simulation of the abm. -# * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. +# * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState for every age group. # */ # static Type log(const mio::abm::Simulation& sim) # { From c5cf7dce478579c93904c8a9e53262187dcb4b48 Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Sat, 14 Jun 2025 14:39:18 +0200 Subject: [PATCH 13/38] julia review p2 --- .../memilio/plot/plotAbmInfectionStates.py | 33 +++++++++---------- .../test_plot_plotAbmInfectionStates.py | 2 +- 2 files changed, 16 insertions(+), 19 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index a3d90ddd90..fbf5386bd6 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -53,12 +53,10 @@ # // PRAGMA_OMP(parallel for) # for (auto i = size_t(0); i < persons.size(); ++i) { # auto& p = persons[i]; -# if (p.get_should_be_logged()) { -# auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) + -# ((uint32_t)p.get_infection_state(curr_time)); -# // PRAGMA_OMP(atomic) -# sum[index] += 1; -# } +# auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) + +# ((uint32_t)p.get_infection_state(curr_time)); +# // PRAGMA_OMP(atomic) +# sum[index] += 1; # } # return std::make_pair(curr_time, sum); # } @@ -97,10 +95,9 @@ # } # }; # -# The output of the loggers of several runs is stored in HDF5 files, with the memilio funciton mio::save_results in mio/io/result_io.h. +# The output of the loggers of several runs is stored in HDF5 files using mio::save_results in mio/io/result_io.h. # Adjust these as needed. -# States and location type numbers need to match the infection states used in your simulation. state_labels = { 1: 'Exposed', 2: 'I_Asymp', @@ -139,7 +136,9 @@ def load_h5_results(base_path, percentile): @param[in] base_path Path to results directory. @param[in] percentile Subdirectory for percentile (e.g. 'p50'). - @return Dictionary with data arrays. + @return Dictionary with data arrays. Keys are dataset names from the HDF5 file + (e.g., 'Time', 'Total', age group names like 'Group1', 'Group2', etc.). + Values are numpy arrays containing the corresponding time series data. """ file_path = os.path.join(base_path, percentile, "Results.h5") with h5py.File(file_path, 'r') as f: @@ -170,7 +169,7 @@ def plot_infections_loc_types_average( plt.figure('Infection_location_types') plt.title( - 'Infection per location type for the median run, rolling sum over 24 hours') + 'Number of new infections per location type for the median run, rolling sum over 24 hours') color_plot = matplotlib.colormaps.get_cmap(colormap).colors for idx, i in enumerate(location_type_labels.keys()): @@ -216,7 +215,7 @@ def plot_infection_states_results( p95 = total_95['Total'] p05 = total_05['Total'] - plot_infection_states_individual( + plot_infection_states_by_age_group( time, p50, p25, p75, colormap, p05_bs=total_05 if show90 else None, p95_bs=total_95 if show90 else None, @@ -254,9 +253,7 @@ def plot_infection_states( linestyle='dashdot', linewidth=1.2, alpha=0.7) plt.plot(x, y75[:, i], color=color_plot[i], linestyle='dashdot', linewidth=1.2, alpha=0.7) - plt.fill_between(x, y50[:, i], y25[:, i], - alpha=0.2, color=color_plot[i]) - plt.fill_between(x, y50[:, i], y75[:, i], + plt.fill_between(x, y25[:, i], y75[:, i], alpha=0.2, color=color_plot[i]) # Optional: 90% percentile if show_90 and y05 is not None and y95 is not None: @@ -274,7 +271,7 @@ def plot_infection_states( plt.show() -def plot_infection_states_individual( +def plot_infection_states_by_age_group( x, p50_bs, p25_bs, p75_bs, colormap='Set1', p05_bs=None, p95_bs=None, show90=False ): @@ -318,8 +315,7 @@ def _plot_state(ax, x, y50, y25, y75, color, title, y05=None, y95=None, show90=F """ Helper to plot a single state with fill_between and optional 90% percentile. """ ax.set_xlabel('time (days)') ax.plot(x, y50, color=color, label='Median') - ax.fill_between(x, y50, y25, alpha=0.5, color=color) - ax.fill_between(x, y50, y75, alpha=0.5, color=color) + ax.fill_between(x, y25, y75, alpha=0.5, color=color) if show90 and y05 is not None and y95 is not None: ax.plot(x, y05, color=color, linestyle='dotted', linewidth=1.0, alpha=0.4) @@ -353,10 +349,11 @@ def main(): parser.add_argument("--colormap", type=str, default='Set1', help="Matplotlib colormap") parser.add_argument("--xtick-step", type=int, - default=150, help="Step for x-axis ticks") + default=150, help="Step for x-axis ticks (usually hours)") parser.add_argument("--90percentile", action="store_true", help="If set, plot 90% percentile as well") args = parser.parse_args() + plot_infection_states_results( args.path_to_infection_states, start_date=args.start_date, diff --git a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py index fbdba48535..d279f0a0a8 100644 --- a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py @@ -131,7 +131,7 @@ def test_plot_infection_states_individual(self, mock_matplotlib): for j in range(7): ax_mock[i, j] = MagicMock() mock_subplots.return_value = (fig_mock, ax_mock) - abm.plot_infection_states_individual( + abm.plot_infection_states_by_age_group( x, p50_bs, p25_bs, p75_bs, colormap='Set1', p05_bs=p05_bs, From f5dffd69ec4d4b42e15c18d7547b242e4dd20988 Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Sun, 15 Jun 2025 19:42:55 +0200 Subject: [PATCH 14/38] julia review pt3 --- .../memilio/plot/plotAbmInfectionStates.py | 33 +++- .../test_plot_plotAbmInfectionStates.py | 176 ++++++++++++++---- 2 files changed, 161 insertions(+), 48 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index fbf5386bd6..89ac9ee3f8 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -65,7 +65,7 @@ # struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { # using Type = std::pair; # /** -# * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. +# * @brief Log the TimeSeries of the number of newly infected Person%s for each Location Type and each age. # * @param[in] sim The simulation of the abm. # * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. # */ @@ -81,7 +81,6 @@ # // PRAGMA_OMP(parallel for) # for (auto i = size_t(0); i < persons.size(); ++i) { # auto& p = persons[i]; -# if (p.get_should_be_logged()) { # // PRAGMA_OMP(atomic) # if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) && # (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { @@ -89,7 +88,6 @@ # ((uint32_t)p.get_location().get_type()); # sum[index] += 1; # } -# } # } # return std::make_pair(curr_time, sum); # } @@ -153,7 +151,7 @@ def plot_infections_loc_types_average( smooth_sigma=1, rolling_window=24, xtick_step=150): - """ Plots rolling average infections per location type for the median run. + """ Plots rolling average new infections per location type for the median run. @param[in] base_path Path to results directory. @param[in] start_date Start date as string. @@ -179,7 +177,7 @@ def plot_infections_loc_types_average( indexer = pd.api.indexers.FixedForwardWindowIndexer( window_size=rolling_window) y = pd.DataFrame(total_50[:, i]).rolling( - window=indexer, min_periods=1).sum().to_numpy() + window=indexer, min_periods=1).mean().to_numpy() y = y[0::rolling_window].flatten() y = gaussian_filter1d(y, sigma=smooth_sigma, mode='nearest') plt.plot(time[0::rolling_window], y, color=color, linewidth=2.5) @@ -198,7 +196,15 @@ def plot_infection_states_results( xtick_step=150, show90=False ): - """ Loads and plots infection state results. """ + """ Loads and plots infection state results. + + @param[in] path_to_infection_states Path to results directory containing infection state data. + @param[in] start_date Start date as string (YYYY-MM-DD format). + @param[in] colormap Matplotlib colormap name. + @param[in] xtick_step Step size for x-axis ticks. + @param[in] show90 If True, plot 90% percentile (5% and 95%) in addition to 50% percentile. + """ + # Load data p50 = load_h5_results(path_to_infection_states, "p50") p25 = load_h5_results(path_to_infection_states, "p25") @@ -275,7 +281,17 @@ def plot_infection_states_by_age_group( x, p50_bs, p25_bs, p75_bs, colormap='Set1', p05_bs=None, p95_bs=None, show90=False ): - """ Plots infection states for each age group, with optional 90% percentile. """ + """ Plots infection states for each age group, with optional 90% percentile. + + @param[in] x Time array for x-axis. + @param[in] p50_bs Dictionary containing 50th percentile data for all age groups. + @param[in] p25_bs Dictionary containing 25th percentile data for all age groups. + @param[in] p75_bs Dictionary containing 75th percentile data for all age groups. + @param[in] colormap Matplotlib colormap name. + @param[in] p05_bs Dictionary containing 5th percentile data for all age groups (optional). + @param[in] p95_bs Dictionary containing 95th percentile data for all age groups (optional). + @param[in] show90 If True, plot 90% percentile bands in addition to 50% percentile. + """ color_plot = matplotlib.colormaps.get_cmap(colormap).colors n_states = len(state_labels) @@ -313,7 +329,6 @@ def plot_infection_states_by_age_group( def _plot_state(ax, x, y50, y25, y75, color, title, y05=None, y95=None, show90=False): """ Helper to plot a single state with fill_between and optional 90% percentile. """ - ax.set_xlabel('time (days)') ax.plot(x, y50, color=color, label='Median') ax.fill_between(x, y25, y75, alpha=0.5, color=color) if show90 and y05 is not None and y95 is not None: @@ -369,7 +384,7 @@ def main(): if not args.path_to_infection_states and not args.path_to_loc_types: print("Please provide a path to infection states or location types results.") - sys.exit(1) + plt.show() diff --git a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py index d279f0a0a8..2822f58a68 100644 --- a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py @@ -1,5 +1,5 @@ ############################################################################# -# Copyright (C) 2020-2024 MEmilio +# Copyright (C) 2020-2025 MEmilio # # Authors: Sascha Korf # @@ -34,8 +34,6 @@ def test_load_h5_results(self, mock_h5file): mock_h5file().__enter__().get.return_value = {'0': mock_group} mock_h5file().__enter__().items.return_value = [('0', mock_group)] mock_h5file().__enter__().__getitem__.return_value = mock_group - mock_h5file().__enter__().items.return_value = [ - ('Time', np.arange(10)), ('Total', np.ones((10, 8)))] with patch('memilio.plot.plotAbmInfectionStates.h5py.File', mock_h5file): result = abm.load_h5_results('dummy_path', 'p50') assert 'Time' in result @@ -67,17 +65,46 @@ def test_plot_infections_loc_types_average(self, mock_df, mock_gauss, mock_matpl @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states') @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states_individual') def test_plot_infection_states_results(self, mock_indiv, mock_states, mock_load): - mock_load.side_effect = [ - {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones( - (10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))}, - {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones( - (10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))}, - {'Time': np.arange(10), 'Total': np.ones((10, 8)), 'Group1': np.ones((10, 8)), 'Group2': np.ones((10, 8)), 'Group3': np.ones( - (10, 8)), 'Group4': np.ones((10, 8)), 'Group5': np.ones((10, 8)), 'Group6': np.ones((10, 8)), 'Total': np.ones((10, 8))} - ] + test_data = { + 'Time': np.arange(10), + 'Total': np.ones((10, 8)), + 'Group1': np.ones((10, 8)), + 'Group2': np.ones((10, 8)), + 'Group3': np.ones((10, 8)), + 'Group4': np.ones((10, 8)), + 'Group5': np.ones((10, 8)), + 'Group6': np.ones((10, 8)) + } + mock_load.side_effect = [test_data, test_data, test_data] + abm.plot_infection_states_results('dummy_path') - assert mock_indiv.called - assert mock_states.called + + # Verify functions are called with correct arguments + self.assertEqual(mock_load.call_count, 3, + "load_h5_results should be called 3 times (p25, p50, p75)") + + # Check that load_h5_results was called with correct percentiles + expected_calls = [ + unittest.mock.call('dummy_path', 'p25'), + unittest.mock.call('dummy_path', 'p50'), + unittest.mock.call('dummy_path', 'p75') + ] + mock_load.assert_has_calls(expected_calls, any_order=True) + + # Verify plotting functions were called with the loaded data + mock_states.assert_called_once() + mock_indiv.assert_called_once() + + # Check that plot_infection_states was called with correct data structure + states_call_args = mock_states.call_args + self.assertIsNotNone(states_call_args) + self.assertEqual(len(states_call_args[0]), 4) # x, y50, y25, y75 + + # Check that plot_infection_states_individual was called with correct data structure + indiv_call_args = mock_indiv.call_args + self.assertIsNotNone(indiv_call_args) + # x, p50_bs, p25_bs, p75_bs + self.assertEqual(len(indiv_call_args[0]), 4) @patch('memilio.plot.plotAbmInfectionStates.matplotlib') def test_plot_infection_states(self, mock_matplotlib): @@ -93,6 +120,7 @@ def test_plot_infection_states(self, mock_matplotlib): with patch.object(abm.plt, 'gca') as mock_gca: mock_ax = MagicMock() mock_gca.return_value = mock_ax + abm.plot_infection_states( x, y50, y25, y75, start_date='2021-03-01', @@ -102,35 +130,55 @@ def test_plot_infection_states(self, mock_matplotlib): y95=y95, show_90=True ) - assert mock_ax.plot.called - assert mock_ax.fill_between.called - assert mock_ax.set_xticks.called - assert mock_ax.set_xticklabels.called + + # Verify plot was called with correct data + self.assertTrue(mock_ax.plot.called, "Plot should be called") + plot_calls = mock_ax.plot.call_args_list + self.assertEqual(len(plot_calls), 8, + "Should plot 8 infection states") + + # Verify fill_between was called for confidence intervals + self.assertTrue(mock_ax.fill_between.called, + "fill_between should be called for confidence intervals") + fill_calls = mock_ax.fill_between.call_args_list + # Should have calls for both 50% and 90% confidence intervals if show_90=True + self.assertGreater(len(fill_calls), 0) + + # Verify axis formatting with correct parameters + mock_ax.set_xticks.assert_called_once() + xticks_call = mock_ax.set_xticks.call_args[0][0] + expected_ticks = np.arange(0, len(x), 2) # xtick_step=2 + np.testing.assert_array_equal(xticks_call, expected_ticks) + + # Verify xticklabels are set correctly + mock_ax.set_xticklabels.assert_called_once() + xticklabels_call = mock_ax.set_xticklabels.call_args[0][0] + self.assertEqual(len(xticklabels_call), len(expected_ticks)) + + # Verify that start_date is used in label formatting + if len(xticklabels_call) > 0: + # Labels should contain date strings when start_date is provided + self.assertTrue(all('2021' in str(label) + for label in xticklabels_call)) @patch('memilio.plot.plotAbmInfectionStates.matplotlib') def test_plot_infection_states_individual(self, mock_matplotlib): x = np.arange(10) group_data = np.ones((10, 8)) - p50_bs = {g: group_data for g in [ - 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} - p25_bs = {g: group_data for g in [ - 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} - p75_bs = {g: group_data for g in [ - 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} - p05_bs = {g: group_data*-1 for g in [ - 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} - p95_bs = {g: group_data*3 for g in [ - 'Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total']} + groups = ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total'] + p50_bs = {g: group_data for g in groups} + p25_bs = {g: group_data for g in groups} + p75_bs = {g: group_data for g in groups} + p05_bs = {g: group_data*-1 for g in groups} + p95_bs = {g: group_data*3 for g in groups} mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1, 0, 0)]*8 - # Patch plt.subplots to return a grid of MagicMock axes (as np.array with dtype=object) + # Patch plt.subplots to return a grid of MagicMock axes with patch.object(abm.plt, 'subplots') as mock_subplots: fig_mock = MagicMock() - ax_mock = np.empty((6, 7), dtype=object) - for i in range(6): - for j in range(7): - ax_mock[i, j] = MagicMock() + ax_mock = MagicMock() # Let the function determine the actual structure mock_subplots.return_value = (fig_mock, ax_mock) + abm.plot_infection_states_by_age_group( x, p50_bs, p25_bs, p75_bs, colormap='Set1', @@ -138,16 +186,66 @@ def test_plot_infection_states_individual(self, mock_matplotlib): p95_bs=p95_bs, show90=True ) - # Check that at least one ax's plot was called - assert any(ax_mock[i, j].plot.called for i in range(6) - for j in range(7)) - assert fig_mock.suptitle.called + + # Verify subplot was called (without assuming specific dimensions) + mock_subplots.assert_called_once() + subplot_call = mock_subplots.call_args + + # Verify that subplots was called with reasonable dimensions + if subplot_call and len(subplot_call[0]) >= 2: + rows, cols = subplot_call[0][:2] + self.assertGreater(rows, 0, "Should have at least one row") + self.assertGreater(cols, 0, "Should have at least one column") + + # The dimensions should be related to our data structure + total_subplots = rows * cols + expected_min_subplots = len(groups) * group_data.shape[1] + self.assertGreaterEqual(total_subplots, expected_min_subplots, + f"Should have at least {expected_min_subplots} subplots for {len(groups)} groups and {group_data.shape[1]} infection states") + + # Verify figure title is set + fig_mock.suptitle.assert_called_once() def test__format_x_axis(self): + test_x = np.arange(10) + test_start_date = '2021-03-01' + test_xtick_step = 2 + with patch('memilio.plot.plotAbmInfectionStates.plt') as mock_plt: - abm._format_x_axis(np.arange(10), '2021-03-01', 2) - assert mock_plt.gca.called - assert mock_plt.gcf.called + mock_ax = MagicMock() + mock_plt.gca.return_value = mock_ax + + abm._format_x_axis(test_x, test_start_date, test_xtick_step) + + # Verify that gca was called to get current axis + mock_plt.gca.assert_called_once() + + # Verify that gcf was called to get current figure + mock_plt.gcf.assert_called_once() + + # Verify axis formatting methods were called + self.assertTrue(mock_ax.set_xticks.called, + "set_xticks should be called") + self.assertTrue(mock_ax.set_xticklabels.called, + "set_xticklabels should be called") + + # Verify correct tick positions + xticks_call = mock_ax.set_xticks.call_args + if xticks_call and xticks_call[0]: + tick_positions = xticks_call[0][0] + expected_positions = np.arange(0, len(test_x), test_xtick_step) + np.testing.assert_array_equal( + tick_positions, expected_positions) + + # Verify tick labels are date strings + xticklabels_call = mock_ax.set_xticklabels.call_args + if xticklabels_call and xticklabels_call[0]: + tick_labels = xticklabels_call[0][0] + self.assertIsInstance(tick_labels, (list, np.ndarray)) + if len(tick_labels) > 0: + # Should contain date information when start_date is provided + self.assertTrue(any('2021' in str(label) + for label in tick_labels)) if __name__ == '__main__': From 75e1105dfd19c4801c0342cb0d7e02486f324e82 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Sun, 15 Jun 2025 19:48:27 +0200 Subject: [PATCH 15/38] format --- .../test_plot_plotAbmInfectionStates.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py index 2822f58a68..bf0b9b86b0 100644 --- a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py @@ -165,7 +165,8 @@ def test_plot_infection_states(self, mock_matplotlib): def test_plot_infection_states_individual(self, mock_matplotlib): x = np.arange(10) group_data = np.ones((10, 8)) - groups = ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Total'] + groups = ['Group1', 'Group2', 'Group3', + 'Group4', 'Group5', 'Group6', 'Total'] p50_bs = {g: group_data for g in groups} p25_bs = {g: group_data for g in groups} p75_bs = {g: group_data for g in groups} @@ -178,7 +179,7 @@ def test_plot_infection_states_individual(self, mock_matplotlib): fig_mock = MagicMock() ax_mock = MagicMock() # Let the function determine the actual structure mock_subplots.return_value = (fig_mock, ax_mock) - + abm.plot_infection_states_by_age_group( x, p50_bs, p25_bs, p75_bs, colormap='Set1', @@ -186,23 +187,23 @@ def test_plot_infection_states_individual(self, mock_matplotlib): p95_bs=p95_bs, show90=True ) - + # Verify subplot was called (without assuming specific dimensions) mock_subplots.assert_called_once() subplot_call = mock_subplots.call_args - + # Verify that subplots was called with reasonable dimensions if subplot_call and len(subplot_call[0]) >= 2: rows, cols = subplot_call[0][:2] self.assertGreater(rows, 0, "Should have at least one row") self.assertGreater(cols, 0, "Should have at least one column") - + # The dimensions should be related to our data structure total_subplots = rows * cols expected_min_subplots = len(groups) * group_data.shape[1] - self.assertGreaterEqual(total_subplots, expected_min_subplots, - f"Should have at least {expected_min_subplots} subplots for {len(groups)} groups and {group_data.shape[1]} infection states") - + self.assertGreaterEqual(total_subplots, expected_min_subplots, + f"Should have at least {expected_min_subplots} subplots for {len(groups)} groups and {group_data.shape[1]} infection states") + # Verify figure title is set fig_mock.suptitle.assert_called_once() From ddb9cb47b004d20a81b158c5f056c0ddfa0880be Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Sun, 15 Jun 2025 20:06:41 +0200 Subject: [PATCH 16/38] fix tests --- .../memilio/plot/plotAbmInfectionStates.py | 18 ++++++--- .../test_plot_plotAbmInfectionStates.py | 40 ++++++++++--------- 2 files changed, 34 insertions(+), 24 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 89ac9ee3f8..f936ed381e 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -151,7 +151,7 @@ def plot_infections_loc_types_average( smooth_sigma=1, rolling_window=24, xtick_step=150): - """ Plots rolling average new infections per location type for the median run. + """ Plots rolling sum of new infections per 24 hours location type for the median run. @param[in] base_path Path to results directory. @param[in] start_date Start date as string. @@ -177,7 +177,7 @@ def plot_infections_loc_types_average( indexer = pd.api.indexers.FixedForwardWindowIndexer( window_size=rolling_window) y = pd.DataFrame(total_50[:, i]).rolling( - window=indexer, min_periods=1).mean().to_numpy() + window=indexer, min_periods=1).sum().to_numpy() y = y[0::rolling_window].flatten() y = gaussian_filter1d(y, sigma=smooth_sigma, mode='nearest') plt.plot(time[0::rolling_window], y, color=color, linewidth=2.5) @@ -319,9 +319,17 @@ def plot_infection_states_by_age_group( ['Median', f'{perc_string} perc.'], loc='upper left', fontsize=8) + # Add y label for leftmost column + if col_idx == 0: + ax[row_idx, col_idx].set_ylabel('Number of individuals') + + # Add x label for bottom row + if row_idx == n_states - 1: + ax[row_idx, col_idx].set_xlabel('Time (days)') + string_short = ' and 90%' if show90 else '' fig.suptitle( - 'Infection states per age group with 50' + string_short + ' percentile', + 'Infection states per age group with 50%' + string_short + ' percentile', fontsize=16) plt.show() @@ -370,14 +378,14 @@ def main(): args = parser.parse_args() plot_infection_states_results( - args.path_to_infection_states, + "/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_state_per_age_group/0", start_date=args.start_date, colormap=args.colormap, xtick_step=args.xtick_step, show90=True ) plot_infections_loc_types_average( - args.path_to_loc_types, + "/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_per_location_type_per_age_group/0", start_date=args.start_date, colormap=args.colormap, xtick_step=args.xtick_step) diff --git a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py index bf0b9b86b0..2f616b7962 100644 --- a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py @@ -63,7 +63,7 @@ def test_plot_infections_loc_types_average(self, mock_df, mock_gauss, mock_matpl @patch('memilio.plot.plotAbmInfectionStates.load_h5_results') @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states') - @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states_individual') + @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states_by_age_group') def test_plot_infection_states_results(self, mock_indiv, mock_states, mock_load): test_data = { 'Time': np.arange(10), @@ -98,13 +98,13 @@ def test_plot_infection_states_results(self, mock_indiv, mock_states, mock_load) # Check that plot_infection_states was called with correct data structure states_call_args = mock_states.call_args self.assertIsNotNone(states_call_args) - self.assertEqual(len(states_call_args[0]), 4) # x, y50, y25, y75 + self.assertEqual(len(states_call_args[0]), 7) # x, y50, y25, y75 - # Check that plot_infection_states_individual was called with correct data structure + # Check that plot_infection_states_by_age_group was called with correct data structure indiv_call_args = mock_indiv.call_args self.assertIsNotNone(indiv_call_args) # x, p50_bs, p25_bs, p75_bs - self.assertEqual(len(indiv_call_args[0]), 4) + self.assertEqual(len(indiv_call_args[0]), 5) @patch('memilio.plot.plotAbmInfectionStates.matplotlib') def test_plot_infection_states(self, mock_matplotlib): @@ -134,8 +134,9 @@ def test_plot_infection_states(self, mock_matplotlib): # Verify plot was called with correct data self.assertTrue(mock_ax.plot.called, "Plot should be called") plot_calls = mock_ax.plot.call_args_list - self.assertEqual(len(plot_calls), 8, - "Should plot 8 infection states") + # From the actual function: 6 infection states * 5 plot calls each (median + 4 percentile lines) = 30 + self.assertEqual(len(plot_calls), 30, + "Should plot 30 lines total (6 states * 5 lines each)") # Verify fill_between was called for confidence intervals self.assertTrue(mock_ax.fill_between.called, @@ -162,7 +163,7 @@ def test_plot_infection_states(self, mock_matplotlib): for label in xticklabels_call)) @patch('memilio.plot.plotAbmInfectionStates.matplotlib') - def test_plot_infection_states_individual(self, mock_matplotlib): + def test_plot_infection_states_by_age_group(self, mock_matplotlib): x = np.arange(10) group_data = np.ones((10, 8)) groups = ['Group1', 'Group2', 'Group3', @@ -177,7 +178,11 @@ def test_plot_infection_states_individual(self, mock_matplotlib): # Patch plt.subplots to return a grid of MagicMock axes with patch.object(abm.plt, 'subplots') as mock_subplots: fig_mock = MagicMock() - ax_mock = MagicMock() # Let the function determine the actual structure + # From the actual function: n_states (6) rows, len(age_groups) (7) columns + ax_mock = np.empty((6, 7), dtype=object) + for i in range(6): + for j in range(7): + ax_mock[i, j] = MagicMock() mock_subplots.return_value = (fig_mock, ax_mock) abm.plot_infection_states_by_age_group( @@ -188,21 +193,17 @@ def test_plot_infection_states_individual(self, mock_matplotlib): show90=True ) - # Verify subplot was called (without assuming specific dimensions) + # Verify subplot was called with correct dimensions mock_subplots.assert_called_once() subplot_call = mock_subplots.call_args # Verify that subplots was called with reasonable dimensions if subplot_call and len(subplot_call[0]) >= 2: rows, cols = subplot_call[0][:2] - self.assertGreater(rows, 0, "Should have at least one row") - self.assertGreater(cols, 0, "Should have at least one column") - - # The dimensions should be related to our data structure - total_subplots = rows * cols - expected_min_subplots = len(groups) * group_data.shape[1] - self.assertGreaterEqual(total_subplots, expected_min_subplots, - f"Should have at least {expected_min_subplots} subplots for {len(groups)} groups and {group_data.shape[1]} infection states") + self.assertEqual( + rows, 6, "Should have 6 rows (number of infection states)") + self.assertEqual( + cols, 7, "Should have 7 columns (number of age groups)") # Verify figure title is set fig_mock.suptitle.assert_called_once() @@ -218,8 +219,9 @@ def test__format_x_axis(self): abm._format_x_axis(test_x, test_start_date, test_xtick_step) - # Verify that gca was called to get current axis - mock_plt.gca.assert_called_once() + # Verify that gca was called to get current axis (it's called twice in the function) + self.assertEqual(mock_plt.gca.call_count, 2, + "gca should be called twice") # Verify that gcf was called to get current figure mock_plt.gcf.assert_called_once() From c2433fcdff3af064a95886de2faa15d20eb697b0 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Sat, 21 Jun 2025 16:04:09 +0200 Subject: [PATCH 17/38] julia review pt5 --- .../memilio/plot/plotAbmInfectionStates.py | 21 ++++++-- .../test_plot_plotAbmInfectionStates.py | 50 +++++++++++++++++-- 2 files changed, 62 insertions(+), 9 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index f936ed381e..626c954ffa 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -67,7 +67,7 @@ # /** # * @brief Log the TimeSeries of the number of newly infected Person%s for each Location Type and each age. # * @param[in] sim The simulation of the abm. -# * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. +# * @return A pair of the TimePoint and the TimeSeries of newly infected Person%s for each Location Type and each age. # */ # static Type log(const mio::abm::Simulation& sim) # { @@ -238,7 +238,20 @@ def plot_infection_states( colormap='Set1', xtick_step=150, y05=None, y95=None, show_90=False): - """ Plots infection states with percentiles and improved styling. """ + """ Plots infection states with percentile bands. + + @param[in] x Time array for x-axis. + @param[in] y50 50th percentile data array. + @param[in] y25 25th percentile data array. + @param[in] y75 75th percentile data array. + @param[in] start_date Start date as string (YYYY-MM-DD format). + @param[in] colormap Matplotlib colormap name. + @param[in] xtick_step Step size for x-axis ticks. + @param[in] y05 5th percentile data array (optional). + @param[in] y95 95th percentile data array (optional). + @param[in] show_90 If True, plot 90% percentile bands in addition to 50% percentile. + """ + plt.figure('Infection_states') plt.title('Infection states with 50% percentile') @@ -378,14 +391,14 @@ def main(): args = parser.parse_args() plot_infection_states_results( - "/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_state_per_age_group/0", + path_to_infection_states=args.path_to_infection_states, start_date=args.start_date, colormap=args.colormap, xtick_step=args.xtick_step, show90=True ) plot_infections_loc_types_average( - "/Users/saschakorf/Nosynch/Arbeit/memilio/memilio/data/cluster_results/final_results/results_2024-09-20192904_best/infection_per_location_type_per_age_group/0", + path_to_loc_types=args.path_to_loc_types, start_date=args.start_date, colormap=args.colormap, xtick_step=args.xtick_step) diff --git a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py index 2f616b7962..dbf9d0a13a 100644 --- a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py @@ -52,15 +52,53 @@ def test_plot_infections_loc_types_average(self, mock_df, mock_gauss, mock_matpl (48, 1)) mock_matplotlib.colormaps.get_cmap.return_value.colors = [(1, 0, 0)]*7 - # Patch plt.gca().plot to a MagicMock - with patch.object(abm.plt, 'gca') as mock_gca: + # Patch plt methods + with patch.object(abm.plt, 'gca') as mock_gca, \ + patch.object(abm.plt, 'figure') as mock_figure, \ + patch.object(abm.plt, 'title') as mock_title, \ + patch.object(abm.plt, 'legend') as mock_legend, \ + patch.object(abm.plt, 'xlabel') as mock_xlabel, \ + patch.object(abm.plt, 'ylabel') as mock_ylabel, \ + patch.object(abm.plt, 'show') as mock_show: + mock_ax = MagicMock() mock_gca.return_value = mock_ax + abm.plot_infections_loc_types_average('dummy_path') + + # Test basic plotting functionality assert mock_ax.plot.called assert mock_ax.set_xticks.called assert mock_ax.set_xticklabels.called + # Test figure settings + self.assertEqual(mock_figure.call_count, 2, + "figure should be called twice") + # Verify first call is with the figure name + mock_figure.assert_any_call('Infection_location_types') + # Verify second call is without arguments (for autofmt_xdate) + mock_figure.assert_any_call() + mock_title.assert_called_once_with( + 'Number of new infections per location type for the median run, rolling sum over 24 hours') + mock_legend.assert_called_once() + mock_xlabel.assert_called_once_with('Date') + mock_ylabel.assert_called_once_with('Number of individuals') + mock_show.assert_called_once() + + # Verify legend was called with location type labels + legend_call_args = mock_legend.call_args + if legend_call_args and legend_call_args[0]: + legend_labels = legend_call_args[0][0] + # Should contain the location type labels from the function + expected_labels = list(abm.location_type_labels.values()) + self.assertEqual(legend_labels, expected_labels) + + # Verify that plot was called for each location type + plot_calls = mock_ax.plot.call_args_list + # Should plot 7 location types + self.assertEqual(len(plot_calls), 7, + "Should plot all 7 location types") + @patch('memilio.plot.plotAbmInfectionStates.load_h5_results') @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states') @patch('memilio.plot.plotAbmInfectionStates.plot_infection_states_by_age_group') @@ -98,12 +136,13 @@ def test_plot_infection_states_results(self, mock_indiv, mock_states, mock_load) # Check that plot_infection_states was called with correct data structure states_call_args = mock_states.call_args self.assertIsNotNone(states_call_args) - self.assertEqual(len(states_call_args[0]), 7) # x, y50, y25, y75 + # x, y50, y25, y75, y05, y95, show_90 + self.assertEqual(len(states_call_args[0]), 7) # Check that plot_infection_states_by_age_group was called with correct data structure indiv_call_args = mock_indiv.call_args self.assertIsNotNone(indiv_call_args) - # x, p50_bs, p25_bs, p75_bs + # x, p50_bs, p25_bs, p75_bs, p05_bs self.assertEqual(len(indiv_call_args[0]), 5) @patch('memilio.plot.plotAbmInfectionStates.matplotlib') @@ -143,7 +182,8 @@ def test_plot_infection_states(self, mock_matplotlib): "fill_between should be called for confidence intervals") fill_calls = mock_ax.fill_between.call_args_list # Should have calls for both 50% and 90% confidence intervals if show_90=True - self.assertGreater(len(fill_calls), 0) + self.assertEqual(len(fill_calls), 12, + "Should have 12 fill_between calls (6 states * 2 confidence intervals each)") # Verify axis formatting with correct parameters mock_ax.set_xticks.assert_called_once() From 606ba7f44ac6d031577a3ff34a64da61cdecac09 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Sat, 21 Jun 2025 16:09:48 +0200 Subject: [PATCH 18/38] julia review pt7 --- .../plot_test/test_plot_plotAbmInfectionStates.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py index dbf9d0a13a..4bc5dde5e1 100644 --- a/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot_test/test_plot_plotAbmInfectionStates.py @@ -196,11 +196,20 @@ def test_plot_infection_states(self, mock_matplotlib): xticklabels_call = mock_ax.set_xticklabels.call_args[0][0] self.assertEqual(len(xticklabels_call), len(expected_ticks)) + # Verify that the xticklabels contain the correct date formatting + expected_dates = ['2021-03-01', '2021-03-03', + '2021-03-05', '2021-03-07', '2021-03-09'] + for i, label in enumerate(xticklabels_call): + self.assertEqual(str(label), expected_dates[i], + f"Label at position {i} should be {expected_dates[i]}") + # Verify that start_date is used in label formatting if len(xticklabels_call) > 0: - # Labels should contain date strings when start_date is provided + # All labels should contain date strings when start_date is provided self.assertTrue(all('2021' in str(label) for label in xticklabels_call)) + # First label should match the start_date + self.assertEqual(str(xticklabels_call[0]), '2021-03-01') @patch('memilio.plot.plotAbmInfectionStates.matplotlib') def test_plot_infection_states_by_age_group(self, mock_matplotlib): @@ -233,7 +242,7 @@ def test_plot_infection_states_by_age_group(self, mock_matplotlib): show90=True ) - # Verify subplot was called with correct dimensions + # Verify that subplots was called to create a grid of axes mock_subplots.assert_called_once() subplot_call = mock_subplots.call_args From e316b6a2c1449284ca084568de8d3f1da35addc8 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Sat, 21 Jun 2025 19:12:13 +0200 Subject: [PATCH 19/38] fix: improve error message for missing file in fake filesystem --- .../surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py index 937e9a4efd..d560138582 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py @@ -505,7 +505,7 @@ def test_network_fit(self, mock_population, mock_baseline, mock_minimum): model.network_fit( path=self.path, model=md, modeltype='classic', training_parameter=training_parameter, filename="data_secir_groups.pickle") - error_message = "[Errno 2] No such file or directory in the fake filesystem: '" + \ + error_message = "[Errno 2] No such file or directory: '" + \ self.path + "data_secir_groups.pickle'" self.assertEqual(str(error.exception), error_message) From 8a002c78a7e8dcc796f17944280b879387e5c621 Mon Sep 17 00:00:00 2001 From: "Korf, Sascha Alexander" <51127093+xsaschako@users.noreply.github.com> Date: Mon, 23 Jun 2025 12:27:41 +0200 Subject: [PATCH 20/38] Update test_surrogatemodel_ode_secir_groups.py --- .../surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py index c0fca2d2be..2bc1bbcffa 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py @@ -505,7 +505,6 @@ def test_network_fit(self, mock_population, mock_baseline, mock_minimum): model.network_fit( path=self.path, model=md, modeltype='classic', training_parameter=training_parameter, filename="data_secir_groups.pickle") - # check error message error_message_part1 = "[Errno 2] No such file or directory" error_message_part2 = self.path From 392ec54f96a0cafb1d0a1dcef334ee3007d36eec Mon Sep 17 00:00:00 2001 From: "Korf, Sascha Alexander" <51127093+xsaschako@users.noreply.github.com> Date: Mon, 23 Jun 2025 12:28:08 +0200 Subject: [PATCH 21/38] Update test_surrogatemodel_ode_secir_groups.py --- .../surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py index 2bc1bbcffa..80f9c0f90d 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py @@ -514,7 +514,6 @@ def test_network_fit(self, mock_population, mock_baseline, mock_minimum): self.assertIn(error_message_part2, str(error.exception)) self.assertIn(error_message_part3, str(error.exception)) - # generate dataset with multiple output input_width = 5 label_width = 10 From 7ccdb45210b877e4da9587bb9147f5c199918815 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Mon, 9 Mar 2026 14:47:34 +0100 Subject: [PATCH 22/38] update --- cpp/examples/abm_history_object.cpp | 3 +- cpp/models/abm/common_abm_loggers.h | 61 +++++++++++++++++++++++++++++ cpp/models/abm/simulation.h | 13 ++++++ 3 files changed, 76 insertions(+), 1 deletion(-) diff --git a/cpp/examples/abm_history_object.cpp b/cpp/examples/abm_history_object.cpp index fb7e61b18f..15512bc8ab 100644 --- a/cpp/examples/abm_history_object.cpp +++ b/cpp/examples/abm_history_object.cpp @@ -25,6 +25,7 @@ #include "memilio/utils/abstract_parameter_distribution.h" #include "memilio/io/history.h" #include "memilio/utils/parameter_distributions.h" +#include "abm/common_abm_loggers.h" #include #include @@ -196,7 +197,7 @@ int main() } }; - mio::History history; + mio::History history; sim.advance(tmax, history); diff --git a/cpp/models/abm/common_abm_loggers.h b/cpp/models/abm/common_abm_loggers.h index 0d91c69c24..cbe7cb7620 100644 --- a/cpp/models/abm/common_abm_loggers.h +++ b/cpp/models/abm/common_abm_loggers.h @@ -189,6 +189,67 @@ struct LogInfectionState : mio::LogAlways { } }; +struct LogInfectionStatePerAgeGroup : mio::LogAlways { + using Type = std::pair; + /** + * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. + * @param[in] sim The simulation of the abm. + * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. + */ + static Type log(const mio::abm::Simulation<>& sim) + { + + Eigen::VectorXd sum = Eigen::VectorXd::Zero( + Eigen::Index((size_t)mio::abm::InfectionState::Count * sim.get_model().parameters.get_num_groups())); + const auto curr_time = sim.get_time(); + const auto persons = sim.get_model().get_persons(); + + + for (auto i = size_t(0); i < persons.size(); ++i) { + auto& p = persons[i]; + auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) + + ((uint32_t)p.get_infection_state(curr_time)); + // PRAGMA_OMP(atomic) + sum[index] += 1; + } + return std::make_pair(curr_time, sum); + } +}; + +struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { + using Type = std::pair; + /** + * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. + * @param[in] sim The simulation of the abm. + * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. + */ + static Type log(const mio::abm::Simulation<>& sim) + { + + Eigen::VectorXd sum = Eigen::VectorXd::Zero( + Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_model().parameters.get_num_groups())); + auto curr_time = sim.get_time(); + auto prev_time = sim.get_prev_time(); + const auto persons = sim.get_model().get_persons(); + + + for (auto i = size_t(0); i < persons.size(); ++i) { + auto& p = persons[i]; + + + if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) && + (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { + auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) + + ((uint32_t)p.get_location()); + sum[index] += 1; + + } + } + return std::make_pair(curr_time, sum); + } +}; + + /** * @brief This is like the DataWriterToMemory, but it only logs time series data. * @tparam Loggers The loggers that are used to log data. The loggers must return a touple with a TimePoint and a value. diff --git a/cpp/models/abm/simulation.h b/cpp/models/abm/simulation.h index 8e6ff52839..283fcd26bd 100644 --- a/cpp/models/abm/simulation.h +++ b/cpp/models/abm/simulation.h @@ -46,6 +46,7 @@ class Simulation Simulation(TimePoint t0, Model&& model) : m_model(std::move(model)) , m_t(t0) + , m_t_prev(t0-hours(1)) , m_dt(hours(1)) { } @@ -85,6 +86,15 @@ class Simulation return m_t; } + /** + * @brief Get the previous time of the Simulation. + */ + TimePoint get_prev_time() const + { + return m_t_prev; + } + + /** * @brief Get the Model that this Simulation evolves. */ @@ -103,11 +113,14 @@ class Simulation { auto dt = std::min(m_dt, tmax - m_t); m_model.evolve(m_t, dt); + m_t_prev = m_t; m_t += m_dt; + } Model m_model; ///< The Model to simulate. TimePoint m_t; ///< The current TimePoint of the Simulation. + TimePoint m_t_prev; ///< The previous TimePoint of the Simulation. TimeSpan m_dt; ///< The length of the time steps. }; From 2132e8ade625e2063814f83d763fb8affe4292fb Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:04:17 +0100 Subject: [PATCH 23/38] update example --- cpp/examples/abm_history_object.cpp | 175 +++++++++++----------------- 1 file changed, 71 insertions(+), 104 deletions(-) diff --git a/cpp/examples/abm_history_object.cpp b/cpp/examples/abm_history_object.cpp index 15512bc8ab..f577430502 100644 --- a/cpp/examples/abm_history_object.cpp +++ b/cpp/examples/abm_history_object.cpp @@ -17,83 +17,47 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#include "abm/common_abm_loggers.h" #include "abm/household.h" #include "abm/lockdown_rules.h" -#include "abm/simulation.h" -#include "abm/model.h" #include "abm/location_type.h" -#include "memilio/utils/abstract_parameter_distribution.h" +#include "abm/model.h" +#include "abm/simulation.h" #include "memilio/io/history.h" +#include "memilio/io/result_io.h" #include "memilio/utils/parameter_distributions.h" -#include "abm/common_abm_loggers.h" +#include "boost/filesystem.hpp" -#include -#include - -std::string convert_loc_id_to_string(std::tuple tuple_id) -{ - return std::to_string(static_cast(std::get<0>(tuple_id))) + "_" + - std::to_string(std::get<1>(tuple_id)); -} - -template -void write_log_to_file(const T& history) -{ - auto logg = history.get_log(); - // Write the results to a file. - auto loc_id = std::get<1>(logg); - auto time_points = std::get<0>(logg); - std::string input; - std::ofstream myfile("test_output.txt"); - myfile << "Locations as numbers:\n"; - for (auto&& id : loc_id[0]) { - myfile << convert_loc_id_to_string(id) << "\n"; - } - myfile << "Timepoints:\n"; - - for (auto&& t : time_points) { - input += std::to_string(t) + " "; - } - myfile << input << "\n"; - - myfile.close(); -} +#include int main() { mio::set_log_level(mio::LogLevel::warn); - // This is a minimal example with children and adults < 60y. - // We divided them into 4 different age groups, which are defined as follows: + + // Define 4 age groups. const size_t num_age_groups = 4; const auto age_group_0_to_4 = mio::AgeGroup(0); const auto age_group_5_to_14 = mio::AgeGroup(1); const auto age_group_15_to_34 = mio::AgeGroup(2); const auto age_group_35_to_59 = mio::AgeGroup(3); - // Create the model with 4 age groups. + // Create the model. auto model = mio::abm::Model(num_age_groups); - mio::ParameterDistributionLogNormal log_norm(4., 1.); - // Set same infection parameter for all age groups. For example, the incubation period is log normally distributed with parameters 4 and 1. model.parameters.get() = mio::ParameterDistributionLogNormal(4., 1.); - - // Set the age group the can go to school is AgeGroup(1) (i.e. 5-14) model.parameters.get()[age_group_5_to_14] = true; - // Set the age group the can go to work is AgeGroup(2) and AgeGroup(3) (i.e. 15-34 and 35-59) model.parameters.get().set_multiple({age_group_15_to_34, age_group_35_to_59}, true); - // There are 3 households for each household group. + // Set up households. int n_households = 3; - // For more than 1 family households we need families. These are parents and children and randoms (which are distributed like the data we have for these households). - auto child = mio::abm::HouseholdMember(num_age_groups); // A child is 50/50% 0-4 or 5-14. + auto child = mio::abm::HouseholdMember(num_age_groups); child.set_age_weight(age_group_0_to_4, 1); child.set_age_weight(age_group_5_to_14, 1); - auto parent = mio::abm::HouseholdMember(num_age_groups); // A parent is 50/50% 15-34 or 35-59. + auto parent = mio::abm::HouseholdMember(num_age_groups); parent.set_age_weight(age_group_15_to_34, 1); parent.set_age_weight(age_group_35_to_59, 1); - // Two-person household with one parent and one child. auto twoPersonHousehold_group = mio::abm::HouseholdGroup(); auto twoPersonHousehold_full = mio::abm::Household(); twoPersonHousehold_full.add_members(child, 1); @@ -101,7 +65,6 @@ int main() twoPersonHousehold_group.add_households(twoPersonHousehold_full, n_households); add_household_group_to_model(model, twoPersonHousehold_group); - // Three-person household with two parent and one child. auto threePersonHousehold_group = mio::abm::HouseholdGroup(); auto threePersonHousehold_full = mio::abm::Household(); threePersonHousehold_full.add_members(child, 1); @@ -109,59 +72,45 @@ int main() threePersonHousehold_group.add_households(threePersonHousehold_full, n_households); add_household_group_to_model(model, threePersonHousehold_group); - // Add one social event with 5 maximum contacts. - // Maximum contacs limit the number of people that a person can infect while being at this location. + // Add locations. auto event = model.add_location(mio::abm::LocationType::SocialEvent); model.get_location(event).get_infection_parameters().set(5); - // Add hospital and ICU with 5 maximum contacs. auto hospital = model.add_location(mio::abm::LocationType::Hospital); model.get_location(hospital).get_infection_parameters().set(5); auto icu = model.add_location(mio::abm::LocationType::ICU); model.get_location(icu).get_infection_parameters().set(5); - // Add one supermarket, maximum constacts are assumed to be 20. auto shop = model.add_location(mio::abm::LocationType::BasicsShop); model.get_location(shop).get_infection_parameters().set(20); - // At every school, the maximum contacts are 20. auto school = model.add_location(mio::abm::LocationType::School); model.get_location(school).get_infection_parameters().set(20); - // At every workplace, maximum contacts are 10. auto work = model.add_location(mio::abm::LocationType::Work); model.get_location(work).get_infection_parameters().set(10); - // People can get tested at work (and do this with 0.5 probability) from time point 0 to day 30. - auto validity_period = mio::abm::days(1); - auto probability = 0.5; - auto start_date = mio::abm::TimePoint(0); - auto end_date = mio::abm::TimePoint(0) + mio::abm::days(30); - auto test_type = mio::abm::TestType::Antigen; - auto test_parameters = model.parameters.get()[test_type]; - auto testing_criteria_work = mio::abm::TestingCriteria(); - auto testing_scheme_work = mio::abm::TestingScheme(testing_criteria_work, validity_period, start_date, end_date, - test_parameters, probability); - model.get_testing_strategy().add_scheme(mio::abm::LocationType::Work, testing_scheme_work); - - // Assign infection state to each person. - // The infection states are chosen randomly. - auto persons = model.get_persons(); - for (auto& person : persons) { - auto rng = mio::abm::PersonalRandomNumberGenerator(person); - mio::abm::InfectionState infection_state = - (mio::abm::InfectionState)(rand() % ((uint32_t)mio::abm::InfectionState::Count - 1)); + // Set up testing at work. + auto start_date = mio::abm::TimePoint(0); + auto end_date = mio::abm::TimePoint(0) + mio::abm::days(30); + auto test_type = mio::abm::TestType::Antigen; + auto test_parameters = model.parameters.get()[test_type]; + auto testing_scheme = mio::abm::TestingScheme(mio::abm::TestingCriteria(), mio::abm::days(1), start_date, + end_date, test_parameters, 0.5); + model.get_testing_strategy().add_scheme(mio::abm::LocationType::Work, testing_scheme); + + // Randomly assign infection states. + for (auto& person : model.get_persons()) { + auto rng = mio::abm::PersonalRandomNumberGenerator(person); + auto infection_state = (mio::abm::InfectionState)(rand() % ((uint32_t)mio::abm::InfectionState::Count - 1)); if (infection_state != mio::abm::InfectionState::Susceptible) person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(), model.parameters, start_date, infection_state)); } - // Assign locations to the people + // Assign locations to people. for (auto& person : model.get_persons()) { const auto pid = person.get_id(); - //assign shop and event model.assign_location(pid, event); model.assign_location(pid, shop); - //assign hospital and ICU model.assign_location(pid, hospital); model.assign_location(pid, icu); - //assign work/school to people depending on their age if (person.get_age() == age_group_5_to_14) { model.assign_location(pid, school); } @@ -170,36 +119,54 @@ int main() } } - // During the lockdown, social events are closed for 90% of people. - auto t_lockdown = mio::abm::TimePoint(0) + mio::abm::days(10); - mio::abm::close_social_events(t_lockdown, 0.9, model.parameters); - + // Run simulation. auto t0 = mio::abm::TimePoint(0); auto tmax = mio::abm::TimePoint(0) + mio::abm::days(30); auto sim = mio::abm::Simulation(t0, std::move(model)); - struct LogTimePoint : mio::LogAlways { - using Type = ScalarType; - static Type log(const mio::abm::Simulation<>& sim) - { - return sim.get_time().hours(); - } - }; - struct LogLocationIds : mio::LogOnce { - using Type = std::vector>; - static Type log(const mio::abm::Simulation<>& sim) - { - Type location_ids{}; - for (auto& location : sim.get_model().get_locations()) { - location_ids.push_back(std::make_tuple(location.get_type(), location.get_id().get())); - } - return location_ids; - } - }; - - mio::History history; + // Set up TimeSeriesWriter-based history loggers. + auto num_age = sim.get_model().parameters.get_num_groups(); + mio::History + historyInfectionStatePerAgeGroup{ + Eigen::Index((size_t)mio::abm::InfectionState::Count * num_age)}; + mio::History + historyInfectionPerLocationTypePerAgeGroup{ + Eigen::Index((size_t)mio::abm::LocationType::Count * num_age)}; + + sim.advance(tmax, historyInfectionStatePerAgeGroup, historyInfectionPerLocationTypePerAgeGroup); + + // Extract TimeSeries from histories. + auto infection_state_per_age_group = + std::get<0>(historyInfectionStatePerAgeGroup.get_log()); + auto infection_per_loc_type_per_age_group = + std::get<0>(historyInfectionPerLocationTypePerAgeGroup.get_log()); + + // Save results to HDF5 files. + auto num_groups = (int)num_age; + boost::filesystem::path result_dir("results"); + + if (!mio::create_directory(result_dir.string())) { + std::cerr << "Error creating result directory." << std::endl; + return 1; + } - sim.advance(tmax, history); + auto dir_infection_state = result_dir / "infection_state_per_age_group"; + if (!mio::create_directory(dir_infection_state.string())) { + std::cerr << "Error creating infection_state_per_age_group directory." << std::endl; + return 1; + } + if (!mio::save_result({infection_state_per_age_group}, {0}, num_groups, + (dir_infection_state / "results_run0.h5").string())) { + std::cerr << "Error saving infection_state_per_age_group results." << std::endl; + } - write_log_to_file(history); + auto dir_infection_per_loc = result_dir / "infection_per_location_type_per_age_group"; + if (!mio::create_directory(dir_infection_per_loc.string())) { + std::cerr << "Error creating infection_per_location_type_per_age_group directory." << std::endl; + return 1; + } + if (!mio::save_result({infection_per_loc_type_per_age_group}, {0}, num_groups, + (dir_infection_per_loc / "results_run0.h5").string())) { + std::cerr << "Error saving infection_per_location_type_per_age_group results." << std::endl; + } } From 0fdbb158b7d1e9329d0d9aedaea9a10e587389b9 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:17:51 +0100 Subject: [PATCH 24/38] update --- .../memilio/plot/plotAbmInfectionStates.py | 92 ++++--------------- 1 file changed, 20 insertions(+), 72 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 626c954ffa..ea9fe36c48 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -35,65 +35,9 @@ # location types from simulation results of the agent-based model (ABM) stored in HDF5 format. # The used Loggers are: -# struct LogInfectionStatePerAgeGroup : mio::LogAlways { -# using Type = std::pair; -# /** -# * @brief Log the TimeSeries of the number of Person%s in an #InfectionState for every age group. -# * @param[in] sim The simulation of the abm. -# * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState for every age group. -# */ -# static Type log(const mio::abm::Simulation& sim) -# { -# -# Eigen::VectorXd sum = Eigen::VectorXd::Zero( -# Eigen::Index((size_t)mio::abm::InfectionState::Count * sim.get_world().parameters.get_num_groups())); -# const auto curr_time = sim.get_time(); -# const auto persons = sim.get_world().get_persons(); -# -# // PRAGMA_OMP(parallel for) -# for (auto i = size_t(0); i < persons.size(); ++i) { -# auto& p = persons[i]; -# auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) + -# ((uint32_t)p.get_infection_state(curr_time)); -# // PRAGMA_OMP(atomic) -# sum[index] += 1; -# } -# return std::make_pair(curr_time, sum); -# } -# }; -# -# struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { -# using Type = std::pair; -# /** -# * @brief Log the TimeSeries of the number of newly infected Person%s for each Location Type and each age. -# * @param[in] sim The simulation of the abm. -# * @return A pair of the TimePoint and the TimeSeries of newly infected Person%s for each Location Type and each age. -# */ -# static Type log(const mio::abm::Simulation& sim) -# { -# -# Eigen::VectorXd sum = Eigen::VectorXd::Zero( -# Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_world().parameters.get_num_groups())); -# auto curr_time = sim.get_time(); -# auto prev_time = sim.get_prev_time(); -# const auto persons = sim.get_world().get_persons(); -# -# // PRAGMA_OMP(parallel for) -# for (auto i = size_t(0); i < persons.size(); ++i) { -# auto& p = persons[i]; -# // PRAGMA_OMP(atomic) -# if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) && -# (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { -# auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) + -# ((uint32_t)p.get_location().get_type()); -# sum[index] += 1; -# } -# } -# return std::make_pair(curr_time, sum); -# } -# }; -# -# The output of the loggers of several runs is stored in HDF5 files using mio::save_results in mio/io/result_io.h. +# LogInfectionStatePerAgeGroup +# LogInfectionPerLocationTypePerAgeGroup +# The output of the loggers of several runs is stored in HDF5 files using mio::save_results in mio/io/result_io.h, see abm_history_object.cpp. # Adjust these as needed. state_labels = { @@ -390,21 +334,25 @@ def main(): help="If set, plot 90% percentile as well") args = parser.parse_args() - plot_infection_states_results( - path_to_infection_states=args.path_to_infection_states, - start_date=args.start_date, - colormap=args.colormap, - xtick_step=args.xtick_step, - show90=True - ) - plot_infections_loc_types_average( - path_to_loc_types=args.path_to_loc_types, - start_date=args.start_date, - colormap=args.colormap, - xtick_step=args.xtick_step) - if not args.path_to_infection_states and not args.path_to_loc_types: print("Please provide a path to infection states or location types results.") + return + + if args.path_to_infection_states: + plot_infection_states_results( + path_to_infection_states=args.path_to_infection_states, + start_date=args.start_date, + colormap=args.colormap, + xtick_step=args.xtick_step, + show90=True + ) + + if args.path_to_loc_types: + plot_infections_loc_types_average( + path_to_loc_types=args.path_to_loc_types, + start_date=args.start_date, + colormap=args.colormap, + xtick_step=args.xtick_step) plt.show() From e22dd77003527ba9d0b69d45b030bade7a96ebd3 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:24:23 +0100 Subject: [PATCH 25/38] redo --- cpp/examples/abm_history_object.cpp | 174 ++++++++++++++++------------ 1 file changed, 103 insertions(+), 71 deletions(-) diff --git a/cpp/examples/abm_history_object.cpp b/cpp/examples/abm_history_object.cpp index f577430502..fb7e61b18f 100644 --- a/cpp/examples/abm_history_object.cpp +++ b/cpp/examples/abm_history_object.cpp @@ -17,47 +17,82 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#include "abm/common_abm_loggers.h" #include "abm/household.h" #include "abm/lockdown_rules.h" -#include "abm/location_type.h" -#include "abm/model.h" #include "abm/simulation.h" +#include "abm/model.h" +#include "abm/location_type.h" +#include "memilio/utils/abstract_parameter_distribution.h" #include "memilio/io/history.h" -#include "memilio/io/result_io.h" #include "memilio/utils/parameter_distributions.h" -#include "boost/filesystem.hpp" -#include +#include +#include + +std::string convert_loc_id_to_string(std::tuple tuple_id) +{ + return std::to_string(static_cast(std::get<0>(tuple_id))) + "_" + + std::to_string(std::get<1>(tuple_id)); +} + +template +void write_log_to_file(const T& history) +{ + auto logg = history.get_log(); + // Write the results to a file. + auto loc_id = std::get<1>(logg); + auto time_points = std::get<0>(logg); + std::string input; + std::ofstream myfile("test_output.txt"); + myfile << "Locations as numbers:\n"; + for (auto&& id : loc_id[0]) { + myfile << convert_loc_id_to_string(id) << "\n"; + } + myfile << "Timepoints:\n"; + + for (auto&& t : time_points) { + input += std::to_string(t) + " "; + } + myfile << input << "\n"; + + myfile.close(); +} int main() { mio::set_log_level(mio::LogLevel::warn); - - // Define 4 age groups. + // This is a minimal example with children and adults < 60y. + // We divided them into 4 different age groups, which are defined as follows: const size_t num_age_groups = 4; const auto age_group_0_to_4 = mio::AgeGroup(0); const auto age_group_5_to_14 = mio::AgeGroup(1); const auto age_group_15_to_34 = mio::AgeGroup(2); const auto age_group_35_to_59 = mio::AgeGroup(3); - // Create the model. + // Create the model with 4 age groups. auto model = mio::abm::Model(num_age_groups); + mio::ParameterDistributionLogNormal log_norm(4., 1.); + // Set same infection parameter for all age groups. For example, the incubation period is log normally distributed with parameters 4 and 1. model.parameters.get() = mio::ParameterDistributionLogNormal(4., 1.); + + // Set the age group the can go to school is AgeGroup(1) (i.e. 5-14) model.parameters.get()[age_group_5_to_14] = true; + // Set the age group the can go to work is AgeGroup(2) and AgeGroup(3) (i.e. 15-34 and 35-59) model.parameters.get().set_multiple({age_group_15_to_34, age_group_35_to_59}, true); - // Set up households. + // There are 3 households for each household group. int n_households = 3; - auto child = mio::abm::HouseholdMember(num_age_groups); + // For more than 1 family households we need families. These are parents and children and randoms (which are distributed like the data we have for these households). + auto child = mio::abm::HouseholdMember(num_age_groups); // A child is 50/50% 0-4 or 5-14. child.set_age_weight(age_group_0_to_4, 1); child.set_age_weight(age_group_5_to_14, 1); - auto parent = mio::abm::HouseholdMember(num_age_groups); + auto parent = mio::abm::HouseholdMember(num_age_groups); // A parent is 50/50% 15-34 or 35-59. parent.set_age_weight(age_group_15_to_34, 1); parent.set_age_weight(age_group_35_to_59, 1); + // Two-person household with one parent and one child. auto twoPersonHousehold_group = mio::abm::HouseholdGroup(); auto twoPersonHousehold_full = mio::abm::Household(); twoPersonHousehold_full.add_members(child, 1); @@ -65,6 +100,7 @@ int main() twoPersonHousehold_group.add_households(twoPersonHousehold_full, n_households); add_household_group_to_model(model, twoPersonHousehold_group); + // Three-person household with two parent and one child. auto threePersonHousehold_group = mio::abm::HouseholdGroup(); auto threePersonHousehold_full = mio::abm::Household(); threePersonHousehold_full.add_members(child, 1); @@ -72,45 +108,59 @@ int main() threePersonHousehold_group.add_households(threePersonHousehold_full, n_households); add_household_group_to_model(model, threePersonHousehold_group); - // Add locations. + // Add one social event with 5 maximum contacts. + // Maximum contacs limit the number of people that a person can infect while being at this location. auto event = model.add_location(mio::abm::LocationType::SocialEvent); model.get_location(event).get_infection_parameters().set(5); + // Add hospital and ICU with 5 maximum contacs. auto hospital = model.add_location(mio::abm::LocationType::Hospital); model.get_location(hospital).get_infection_parameters().set(5); auto icu = model.add_location(mio::abm::LocationType::ICU); model.get_location(icu).get_infection_parameters().set(5); + // Add one supermarket, maximum constacts are assumed to be 20. auto shop = model.add_location(mio::abm::LocationType::BasicsShop); model.get_location(shop).get_infection_parameters().set(20); + // At every school, the maximum contacts are 20. auto school = model.add_location(mio::abm::LocationType::School); model.get_location(school).get_infection_parameters().set(20); + // At every workplace, maximum contacts are 10. auto work = model.add_location(mio::abm::LocationType::Work); model.get_location(work).get_infection_parameters().set(10); - // Set up testing at work. - auto start_date = mio::abm::TimePoint(0); - auto end_date = mio::abm::TimePoint(0) + mio::abm::days(30); - auto test_type = mio::abm::TestType::Antigen; - auto test_parameters = model.parameters.get()[test_type]; - auto testing_scheme = mio::abm::TestingScheme(mio::abm::TestingCriteria(), mio::abm::days(1), start_date, - end_date, test_parameters, 0.5); - model.get_testing_strategy().add_scheme(mio::abm::LocationType::Work, testing_scheme); - - // Randomly assign infection states. - for (auto& person : model.get_persons()) { - auto rng = mio::abm::PersonalRandomNumberGenerator(person); - auto infection_state = (mio::abm::InfectionState)(rand() % ((uint32_t)mio::abm::InfectionState::Count - 1)); + // People can get tested at work (and do this with 0.5 probability) from time point 0 to day 30. + auto validity_period = mio::abm::days(1); + auto probability = 0.5; + auto start_date = mio::abm::TimePoint(0); + auto end_date = mio::abm::TimePoint(0) + mio::abm::days(30); + auto test_type = mio::abm::TestType::Antigen; + auto test_parameters = model.parameters.get()[test_type]; + auto testing_criteria_work = mio::abm::TestingCriteria(); + auto testing_scheme_work = mio::abm::TestingScheme(testing_criteria_work, validity_period, start_date, end_date, + test_parameters, probability); + model.get_testing_strategy().add_scheme(mio::abm::LocationType::Work, testing_scheme_work); + + // Assign infection state to each person. + // The infection states are chosen randomly. + auto persons = model.get_persons(); + for (auto& person : persons) { + auto rng = mio::abm::PersonalRandomNumberGenerator(person); + mio::abm::InfectionState infection_state = + (mio::abm::InfectionState)(rand() % ((uint32_t)mio::abm::InfectionState::Count - 1)); if (infection_state != mio::abm::InfectionState::Susceptible) person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(), model.parameters, start_date, infection_state)); } - // Assign locations to people. + // Assign locations to the people for (auto& person : model.get_persons()) { const auto pid = person.get_id(); + //assign shop and event model.assign_location(pid, event); model.assign_location(pid, shop); + //assign hospital and ICU model.assign_location(pid, hospital); model.assign_location(pid, icu); + //assign work/school to people depending on their age if (person.get_age() == age_group_5_to_14) { model.assign_location(pid, school); } @@ -119,54 +169,36 @@ int main() } } - // Run simulation. + // During the lockdown, social events are closed for 90% of people. + auto t_lockdown = mio::abm::TimePoint(0) + mio::abm::days(10); + mio::abm::close_social_events(t_lockdown, 0.9, model.parameters); + auto t0 = mio::abm::TimePoint(0); auto tmax = mio::abm::TimePoint(0) + mio::abm::days(30); auto sim = mio::abm::Simulation(t0, std::move(model)); - // Set up TimeSeriesWriter-based history loggers. - auto num_age = sim.get_model().parameters.get_num_groups(); - mio::History - historyInfectionStatePerAgeGroup{ - Eigen::Index((size_t)mio::abm::InfectionState::Count * num_age)}; - mio::History - historyInfectionPerLocationTypePerAgeGroup{ - Eigen::Index((size_t)mio::abm::LocationType::Count * num_age)}; - - sim.advance(tmax, historyInfectionStatePerAgeGroup, historyInfectionPerLocationTypePerAgeGroup); - - // Extract TimeSeries from histories. - auto infection_state_per_age_group = - std::get<0>(historyInfectionStatePerAgeGroup.get_log()); - auto infection_per_loc_type_per_age_group = - std::get<0>(historyInfectionPerLocationTypePerAgeGroup.get_log()); - - // Save results to HDF5 files. - auto num_groups = (int)num_age; - boost::filesystem::path result_dir("results"); - - if (!mio::create_directory(result_dir.string())) { - std::cerr << "Error creating result directory." << std::endl; - return 1; - } + struct LogTimePoint : mio::LogAlways { + using Type = ScalarType; + static Type log(const mio::abm::Simulation<>& sim) + { + return sim.get_time().hours(); + } + }; + struct LogLocationIds : mio::LogOnce { + using Type = std::vector>; + static Type log(const mio::abm::Simulation<>& sim) + { + Type location_ids{}; + for (auto& location : sim.get_model().get_locations()) { + location_ids.push_back(std::make_tuple(location.get_type(), location.get_id().get())); + } + return location_ids; + } + }; - auto dir_infection_state = result_dir / "infection_state_per_age_group"; - if (!mio::create_directory(dir_infection_state.string())) { - std::cerr << "Error creating infection_state_per_age_group directory." << std::endl; - return 1; - } - if (!mio::save_result({infection_state_per_age_group}, {0}, num_groups, - (dir_infection_state / "results_run0.h5").string())) { - std::cerr << "Error saving infection_state_per_age_group results." << std::endl; - } + mio::History history; - auto dir_infection_per_loc = result_dir / "infection_per_location_type_per_age_group"; - if (!mio::create_directory(dir_infection_per_loc.string())) { - std::cerr << "Error creating infection_per_location_type_per_age_group directory." << std::endl; - return 1; - } - if (!mio::save_result({infection_per_loc_type_per_age_group}, {0}, num_groups, - (dir_infection_per_loc / "results_run0.h5").string())) { - std::cerr << "Error saving infection_per_location_type_per_age_group results." << std::endl; - } + sim.advance(tmax, history); + + write_log_to_file(history); } From 7d544bf54cb58fc9c63f107eaba74218658aa177 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:48:41 +0100 Subject: [PATCH 26/38] update --- cpp/examples/abm_parameter_study.cpp | 75 +++++++++++++------ cpp/models/abm/result_simulation.h | 22 +++++- .../memilio/plot/plotAbmInfectionStates.py | 65 +++++++++++++--- 3 files changed, 127 insertions(+), 35 deletions(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index cfaba3ae94..079a6b3743 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -34,6 +34,7 @@ #include "memilio/utils/stl_util.h" #include +#include constexpr size_t num_age_groups = 4; @@ -173,9 +174,9 @@ int main() // Set start and end time for the simulation. auto t0 = mio::abm::TimePoint(0); - auto tmax = t0 + mio::abm::days(5); + auto tmax = t0 + mio::abm::days(10); // Set the number of simulations to run in the study - const size_t num_runs = 3; + const size_t num_runs = 10; // Create a parameter study. // Note that the study for the ABM currently does not make use of the arguments "parameters" or "dt", as we create @@ -187,42 +188,74 @@ int main() // study.get_rng().seed({12341234, 53456, 63451, 5232576, 84586, 52345}); const std::string result_dir = mio::path_join(mio::base_dir(), "example_results"); + // Remove old results to avoid HDF5 "name already exists" errors + std::filesystem::remove_all(result_dir); if (!mio::create_directory(result_dir)) { mio::log_error("Could not create result directory \"{}\".", result_dir); return 1; } + const std::string result_dir_new = mio::path_join(result_dir, "new"); + if (!mio::create_directory(result_dir_new)) { + mio::log_error("Could not create result directory \"{}\".", result_dir_new); + return 1; + } + + const std::string result_dir_detailed = mio::path_join(result_dir, "detailed"); + if (!mio::create_directory(result_dir_detailed)) { + mio::log_error("Could not create result directory \"{}\".", result_dir_detailed); + return 1; + } + auto ensemble_results = study.run( [](auto, auto t0_, auto, size_t) { return mio::abm::ResultSimulation(make_model(mio::thread_local_rng()), t0_); }, [result_dir](auto&& sim, auto&& run_idx) { - auto interpolated_result = mio::interpolate_simulation_result(sim.get_result()); + auto interpolated_result = mio::interpolate_simulation_result(sim.get_result()); + auto interpolated_result_detailed = mio::interpolate_simulation_result(sim.get_result_detailed()); + std::string outpath = mio::path_join(result_dir, "abm_minimal_run_" + std::to_string(run_idx) + ".txt"); std::ofstream outfile_run(outpath); sim.get_result().print_table(outfile_run, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4); + std::cout << "Results written to " << outpath << std::endl; - auto params = std::vector{}; - return std::vector{interpolated_result}; + return std::make_pair(std::vector{interpolated_result}, std::vector{interpolated_result_detailed}); }); if (ensemble_results.size() > 0) { - auto ensemble_results_p05 = ensemble_percentile(ensemble_results, 0.05); - auto ensemble_results_p25 = ensemble_percentile(ensemble_results, 0.25); - auto ensemble_results_p50 = ensemble_percentile(ensemble_results, 0.50); - auto ensemble_results_p75 = ensemble_percentile(ensemble_results, 0.75); - auto ensemble_results_p95 = ensemble_percentile(ensemble_results, 0.95); - - mio::unused(save_result(ensemble_results_p05, {0}, num_age_groups, - mio::path_join(result_dir, "Results_" + std::string("p05") + ".h5"))); - mio::unused(save_result(ensemble_results_p25, {0}, num_age_groups, - mio::path_join(result_dir, "Results_" + std::string("p25") + ".h5"))); - mio::unused(save_result(ensemble_results_p50, {0}, num_age_groups, - mio::path_join(result_dir, "Results_" + std::string("p50") + ".h5"))); - mio::unused(save_result(ensemble_results_p75, {0}, num_age_groups, - mio::path_join(result_dir, "Results_" + std::string("p75") + ".h5"))); - mio::unused(save_result(ensemble_results_p95, {0}, num_age_groups, - mio::path_join(result_dir, "Results_" + std::string("p95") + ".h5"))); + // Split the pairs into two separate ensemble vectors + std::vector>> ensemble_new, ensemble_detailed; + for (auto& [res_new, res_det] : ensemble_results) { + ensemble_new.push_back(std::move(res_new)); + ensemble_detailed.push_back(std::move(res_det)); + } + + // Percentiles for aggregated results + auto new_p05 = ensemble_percentile(ensemble_new, 0.05); + auto new_p25 = ensemble_percentile(ensemble_new, 0.25); + auto new_p50 = ensemble_percentile(ensemble_new, 0.50); + auto new_p75 = ensemble_percentile(ensemble_new, 0.75); + auto new_p95 = ensemble_percentile(ensemble_new, 0.95); + + mio::unused(save_result(new_p05, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p05.h5"))); + mio::unused(save_result(new_p25, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p25.h5"))); + mio::unused(save_result(new_p50, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p50.h5"))); + mio::unused(save_result(new_p75, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p75.h5"))); + mio::unused(save_result(new_p95, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p95.h5"))); + + // Percentiles for detailed results + auto det_p05 = ensemble_percentile(ensemble_detailed, 0.05); + auto det_p25 = ensemble_percentile(ensemble_detailed, 0.25); + auto det_p50 = ensemble_percentile(ensemble_detailed, 0.50); + auto det_p75 = ensemble_percentile(ensemble_detailed, 0.75); + auto det_p95 = ensemble_percentile(ensemble_detailed, 0.95); + + mio::unused(save_result(det_p05, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p05.h5"))); + mio::unused(save_result(det_p25, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p25.h5"))); + mio::unused(save_result(det_p50, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p50.h5"))); + mio::unused(save_result(det_p75, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p75.h5"))); + mio::unused(save_result(det_p95, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p95.h5"))); } mio::mpi::finalize(); diff --git a/cpp/models/abm/result_simulation.h b/cpp/models/abm/result_simulation.h index 2f1d02131a..d2f2128ea6 100644 --- a/cpp/models/abm/result_simulation.h +++ b/cpp/models/abm/result_simulation.h @@ -45,19 +45,33 @@ class ResultSimulation : public Simulation */ void advance(TimePoint tmax) { - Simulation::advance(tmax, history); + Simulation::advance(tmax, history, history_detailed); } /** * @brief Return the simulation result aggregated by infection states. */ - const mio::TimeSeries& get_result() const + mio::TimeSeries get_result() const { - return get<0>(history.get_log()); + return std::get<0>(history.get_log()); } + /** + * @brief Return the detailed simulation result aggregated by infection states. + */ + mio::TimeSeries get_result_detailed() const + { + return std::get<0>(history_detailed.get_log()); + } + + mio::History history{ - Eigen::Index(InfectionState::Count)}; ///< History used to create the result TimeSeries. + Eigen::Index(InfectionState::Count)}; + + mio::History history_detailed{ + Eigen::Index(InfectionState::Count) * this->get_model().parameters.get_num_groups()}; + + }; } // namespace abm diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index ea9fe36c48..1534fd73ea 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -82,7 +82,37 @@ def load_h5_results(base_path, percentile): (e.g., 'Time', 'Total', age group names like 'Group1', 'Group2', etc.). Values are numpy arrays containing the corresponding time series data. """ - file_path = os.path.join(base_path, percentile, "Results.h5") + # Try flat structure first (e.g., Results_p50.h5) + file_path_flat = os.path.join(base_path, f"Results_{percentile}.h5") + # Try subdirectory structure (e.g., p50/Results.h5) + file_path_nested = os.path.join(base_path, percentile, "Results.h5") + + # Determine which file exists + if os.path.exists(file_path_flat): + file_path = file_path_flat + elif os.path.exists(file_path_nested): + file_path = file_path_nested + else: + raise FileNotFoundError( + f"Could not find percentile results file. Tried:\n" + f" - {file_path_flat}\n" + f" - {file_path_nested}\n\n" + f"Expected directory structure (option 1):\n" + f" {base_path}/\n" + f" Results_p05.h5\n" + f" Results_p25.h5\n" + f" Results_p50.h5\n" + f" Results_p75.h5\n" + f" Results_p95.h5\n\n" + f"Or (option 2):\n" + f" {base_path}/\n" + f" p05/Results.h5\n" + f" p25/Results.h5\n" + f" p50/Results.h5\n" + f" p75/Results.h5\n" + f" p95/Results.h5\n" + ) + with h5py.File(file_path, 'r') as f: data = {k: v[()] for k, v in f['0'].items()} return data @@ -114,7 +144,17 @@ def plot_infections_loc_types_average( 'Number of new infections per location type for the median run, rolling sum over 24 hours') color_plot = matplotlib.colormaps.get_cmap(colormap).colors + # Check if data dimensions match expected location types + num_cols = total_50.shape[1] if len(total_50.shape) > 1 else 1 + num_location_types = len(location_type_labels) + + if num_cols < num_location_types: + print(f"Warning: Data has {num_cols} columns but {num_location_types} location types defined.") + print(f"Only plotting first {num_cols} location types.") + for idx, i in enumerate(location_type_labels.keys()): + if i >= num_cols: + break # Skip if we don't have data for this location type color = color_plot[i % len(color_plot)] if i < len( color_plot) else "black" # Sum up every 24 hours, then smooth @@ -124,13 +164,12 @@ def plot_infections_loc_types_average( window=indexer, min_periods=1).sum().to_numpy() y = y[0::rolling_window].flatten() y = gaussian_filter1d(y, sigma=smooth_sigma, mode='nearest') - plt.plot(time[0::rolling_window], y, color=color, linewidth=2.5) + plt.plot(time[0::rolling_window], y, color=color, linewidth=2.5, label=location_type_labels[i]) - plt.legend(list(location_type_labels.values())) + plt.legend() # Auto-generate legend from plot labels _format_x_axis(time, start_date, xtick_step) plt.xlabel('Date') plt.ylabel('Number of individuals') - plt.show() def plot_infection_states_results( @@ -231,7 +270,6 @@ def plot_infection_states( _format_x_axis(x, start_date, xtick_step) plt.xlabel('Date') plt.ylabel('Number of individuals') - plt.show() def plot_infection_states_by_age_group( @@ -250,22 +288,31 @@ def plot_infection_states_by_age_group( @param[in] show90 If True, plot 90% percentile bands in addition to 50% percentile. """ + # Dynamically detect available age groups from the data + available_groups = [key for key in p50_bs.keys() if key.startswith('Group') or key == 'Total'] + # Sort to ensure Group1, Group2, ..., Total order + available_groups = sorted(available_groups, key=lambda x: (x != 'Total', x)) + color_plot = matplotlib.colormaps.get_cmap(colormap).colors n_states = len(state_labels) fig, ax = plt.subplots( - n_states, len(age_groups), constrained_layout=True, figsize=(20, 3 * n_states)) + n_states, len(available_groups), constrained_layout=True, figsize=(20, 3 * n_states)) - for col_idx, group in enumerate(age_groups): + for col_idx, group in enumerate(available_groups): y50 = p50_bs[group] y25 = p25_bs[group] y75 = p75_bs[group] y05 = p05_bs[group] if (show90 and p05_bs is not None) else None y95 = p95_bs[group] if (show90 and p95_bs is not None) else None + + # Get group label, using default if not in predefined dict + group_label = age_groups_dict.get(group, group) + for row_idx, (state_idx, label) in enumerate(state_labels.items()): _plot_state( ax[row_idx, col_idx], x, y50[:, state_idx], y25[:, state_idx], y75[:, state_idx], - color_plot[col_idx], f'#{label}, {age_groups_dict[group]}', + color_plot[col_idx], f'#{label}, {group_label}', y05=y05[:, state_idx] if y05 is not None else None, y95=y95[:, state_idx] if y95 is not None else None, show90=show90 @@ -289,8 +336,6 @@ def plot_infection_states_by_age_group( 'Infection states per age group with 50%' + string_short + ' percentile', fontsize=16) - plt.show() - def _plot_state(ax, x, y50, y25, y75, color, title, y05=None, y95=None, show90=False): """ Helper to plot a single state with fill_between and optional 90% percentile. """ From 9724423f585bdce33673c809cb9b74b1363b3b62 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:49:11 +0100 Subject: [PATCH 27/38] fix comment formatting in plotAbmInfectionStates.py --- pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 1534fd73ea..3a4db0f75c 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -34,7 +34,7 @@ # This module provides functions to load and visualize infection states and # location types from simulation results of the agent-based model (ABM) stored in HDF5 format. -# The used Loggers are: +# The used Loggers are: # LogInfectionStatePerAgeGroup # LogInfectionPerLocationTypePerAgeGroup # The output of the loggers of several runs is stored in HDF5 files using mio::save_results in mio/io/result_io.h, see abm_history_object.cpp. From 22044c63823dda03f755184bcd5b4ca69eb2678e Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:54:48 +0100 Subject: [PATCH 28/38] update --- cpp/examples/abm_parameter_study.cpp | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 079a6b3743..92272f8302 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -195,9 +195,9 @@ int main() return 1; } - const std::string result_dir_new = mio::path_join(result_dir, "new"); - if (!mio::create_directory(result_dir_new)) { - mio::log_error("Could not create result directory \"{}\".", result_dir_new); + const std::string result_directory_standard = mio::path_join(result_dir, "standart_results"); + if (!mio::create_directory(result_directory_standard)) { + mio::log_error("Could not create result directory \"{}\".", result_directory_standard); return 1; } @@ -220,15 +220,16 @@ int main() sim.get_result().print_table(outfile_run, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4); std::cout << "Results written to " << outpath << std::endl; - return std::make_pair(std::vector{interpolated_result}, std::vector{interpolated_result_detailed}); + + return std::vector>{interpolated_result, interpolated_result_detailed}; }); if (ensemble_results.size() > 0) { - // Split the pairs into two separate ensemble vectors + std::vector>> ensemble_new, ensemble_detailed; - for (auto& [res_new, res_det] : ensemble_results) { - ensemble_new.push_back(std::move(res_new)); - ensemble_detailed.push_back(std::move(res_det)); + for (auto& run : ensemble_results) { + ensemble_new.push_back({std::move(run[0])}); + ensemble_detailed.push_back({std::move(run[1])}); } // Percentiles for aggregated results @@ -238,11 +239,11 @@ int main() auto new_p75 = ensemble_percentile(ensemble_new, 0.75); auto new_p95 = ensemble_percentile(ensemble_new, 0.95); - mio::unused(save_result(new_p05, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p05.h5"))); - mio::unused(save_result(new_p25, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p25.h5"))); - mio::unused(save_result(new_p50, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p50.h5"))); - mio::unused(save_result(new_p75, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p75.h5"))); - mio::unused(save_result(new_p95, {0}, num_age_groups, mio::path_join(result_dir_new, "Results_p95.h5"))); + mio::unused(save_result(new_p05, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p05.h5"))); + mio::unused(save_result(new_p25, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p25.h5"))); + mio::unused(save_result(new_p50, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p50.h5"))); + mio::unused(save_result(new_p75, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p75.h5"))); + mio::unused(save_result(new_p95, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p95.h5"))); // Percentiles for detailed results auto det_p05 = ensemble_percentile(ensemble_detailed, 0.05); From 6c6879f096fc7e42d8ede70d1dfb0af2a7fe8405 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Mon, 9 Mar 2026 17:08:15 +0100 Subject: [PATCH 29/38] update --- cpp/models/abm/common_abm_loggers.h | 2 +- cpp/models/abm/result_simulation.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/models/abm/common_abm_loggers.h b/cpp/models/abm/common_abm_loggers.h index cbe7cb7620..b0a40d5959 100644 --- a/cpp/models/abm/common_abm_loggers.h +++ b/cpp/models/abm/common_abm_loggers.h @@ -240,7 +240,7 @@ struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) && (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) + - ((uint32_t)p.get_location()); + ((uint32_t)p.get_location_type()); sum[index] += 1; } diff --git a/cpp/models/abm/result_simulation.h b/cpp/models/abm/result_simulation.h index d2f2128ea6..02aeffe069 100644 --- a/cpp/models/abm/result_simulation.h +++ b/cpp/models/abm/result_simulation.h @@ -69,7 +69,7 @@ class ResultSimulation : public Simulation Eigen::Index(InfectionState::Count)}; mio::History history_detailed{ - Eigen::Index(InfectionState::Count) * this->get_model().parameters.get_num_groups()}; + Eigen::Index(LocationType::Count) * this->get_model().parameters.get_num_groups()}; }; From 49d0587d83d1801fe1252951614b0f2e65e7663f Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Mon, 9 Mar 2026 17:15:28 +0100 Subject: [PATCH 30/38] update --- cpp/examples/abm_parameter_study.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 92272f8302..d6b44c4bb1 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -188,7 +188,6 @@ int main() // study.get_rng().seed({12341234, 53456, 63451, 5232576, 84586, 52345}); const std::string result_dir = mio::path_join(mio::base_dir(), "example_results"); - // Remove old results to avoid HDF5 "name already exists" errors std::filesystem::remove_all(result_dir); if (!mio::create_directory(result_dir)) { mio::log_error("Could not create result directory \"{}\".", result_dir); From 72fcfc59845003ebd81e76758caa00c4904b720a Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Fri, 31 Jul 2026 19:53:19 +0200 Subject: [PATCH 31/38] update --- cpp/examples/abm_parameter_study.cpp | 49 +++++++++---- cpp/models/abm/common_abm_loggers.h | 73 +++++++++++-------- cpp/models/abm/result_simulation.h | 18 ++--- cpp/models/abm/simulation.h | 11 +-- .../memilio/plot/plotAbmInfectionStates.py | 14 ++-- .../tests/test_plot_plotAbmInfectionStates.py | 37 +++++++--- 6 files changed, 126 insertions(+), 76 deletions(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index f9ea8c1720..61dd73d93e 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -173,9 +173,9 @@ int main() // Set start and end time for the simulation. auto t0 = mio::abm::TimePoint(0); - auto tmax = t0 + mio::abm::days(10); + auto tmax = t0 + mio::abm::days(30); // Set the number of simulations to run in the study - const size_t num_runs = 3; + const size_t num_runs = 11; // Set up an RNG mio::RandomNumberGenerator rng; @@ -200,6 +200,9 @@ int main() const auto result_dir_standard = mio::create_directories_or_exit(result_dir / "standard_results"); const auto result_dir_detailed = mio::create_directories_or_exit(result_dir / "detailed_results"); + // Collects the detailed result of each run, see the process_simulation_result lambda below. + std::vector>> ensemble_results_detailed; + // Run the study // The first lambda ("create_simulation" argument) sets up the simulation, the second ("process_simulation_result") // allows us to process each simulations result. Be mindful of the memory used for storing these results! @@ -212,7 +215,7 @@ int main() copy.reset_rng(ctr); return mio::abm::ResultSimulation(std::move(copy), t0_); }, - [&result_dir_standard, &result_dir_detailed](auto&& sim, auto&& run_idx) { + [&result_dir_standard, &result_dir_detailed, &ensemble_results_detailed](auto&& sim, auto&& run_idx) { auto interpolated_result = mio::interpolate_simulation_result(sim.get_result()); auto interpolated_result_detailed = mio::interpolate_simulation_result(sim.get_result_detailed()); @@ -228,27 +231,45 @@ int main() sim.get_result_detailed().print_table(outfile_run_detailed, {}, 7, 4); std::cout << "Detailed results written to " << outpath_detailed.string() << std::endl; + // The detailed result is collected separately, because ensemble_percentile requires all entries of an + // ensemble to have the same number of elements, which the two results do not have. + ensemble_results_detailed.push_back({interpolated_result_detailed}); + return std::vector{interpolated_result}; }); // The study collects all results on the root rank, so we only process the results there if (mio::mpi::is_root()) { - const auto write_percentile = [&](double p) { - auto ensemble_percentiles = ensemble_percentile(ensemble_results, p); - - std::ofstream out(result_dir / fmt::format("Results_p{:0<4.2}.txt", p)); + // The percentiles are written in the "Results_p05.h5" layout expected by + // pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py, so that the results of this example can be + // plotted directly. The percentile is given in whole percent to match that file name. + const auto write_percentile = [&](int p) { + auto ensemble_percentiles = ensemble_percentile(ensemble_results, p / 100.0); + auto ensemble_percentiles_detailed = ensemble_percentile(ensemble_results_detailed, p / 100.0); + + std::ofstream out(result_dir_standard / fmt::format("Results_p{:02d}.txt", p)); ensemble_percentiles.front().print_table(out, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4); - mio::unused(save_result(ensemble_percentiles, {0}, num_age_groups, - (result_dir / fmt::format("Results_p{:0<4.2}.h5", p)).string())); + // save_result splits each row into num_groups groups. Both results are a single vector per time point, + // so they are written as one group each, whose "Total" holds all entries of that vector. Passing the + // number of age groups here would instead split those entries across that many groups. + mio::unused(mio::save_result(ensemble_percentiles, {0}, 1, + (result_dir_standard / fmt::format("Results_p{:02d}.h5", p)).string())); + mio::unused(mio::save_result(ensemble_percentiles_detailed, {0}, 1, + (result_dir_detailed / fmt::format("Results_p{:02d}.h5", p)).string())); }; - write_percentile(0.05); - write_percentile(0.25); - write_percentile(0.50); - write_percentile(0.75); - write_percentile(0.95); + write_percentile(5); + write_percentile(25); + write_percentile(50); + write_percentile(75); + write_percentile(95); + + std::cout << "\nPercentiles written. Plot them with (in folder memilio-plot) &&:\n" + << " python -m memilio.plot.plotAbmInfectionStates \\\n" + << " --path-to-infection-states " << result_dir_standard.string() << " \\\n" + << " --path-to-loc-types " << result_dir_detailed.string() << std::endl; } mio::mpi::finalize(); diff --git a/cpp/models/abm/common_abm_loggers.h b/cpp/models/abm/common_abm_loggers.h index b0a40d5959..fb25e1d23d 100644 --- a/cpp/models/abm/common_abm_loggers.h +++ b/cpp/models/abm/common_abm_loggers.h @@ -37,10 +37,10 @@ namespace abm /** * @brief Struct to save specific mobility data of an agent. - * The data consists of: - * + * The data consists of the agent's id, the Location%s it moved between, the times it left and arrived, + * the TransportMode and ActivityType of the trip, and the agent's InfectionState. */ -struct mobility_data { +struct MobilityData { uint32_t agent_id; uint32_t from_id; uint32_t to_id; @@ -51,6 +51,11 @@ struct mobility_data { mio::abm::InfectionState infection_state; }; +/** + * @brief Deduce the ActivityType an agent is most likely pursuing from the LocationType it is at. + * @param[in] current_location The type of the Location the agent is currently at. + * @return The ActivityType associated with the given LocationType, or ActivityType::UnknownActivity if there is none. + */ constexpr mio::abm::ActivityType guess_activity_type(mio::abm::LocationType current_location) { switch (current_location) { @@ -153,7 +158,7 @@ struct LogDataForMobility : mio::LogAlways { static Type log(const mio::abm::Simulation<>& sim) { Type mobility_data{}; - for (Person p : sim.get_model().get_persons()) { + for (const Person& p : sim.get_model().get_persons()) { mobility_data.push_back( std::make_tuple(p.get_id(), p.get_location(), sim.get_time(), p.get_last_transport_mode(), guess_activity_type(p.get_location_type()), p.get_infection_state(sim.get_time()))); @@ -174,11 +179,9 @@ struct LogInfectionState : mio::LogAlways { */ static Type log(const mio::abm::Simulation<>& sim) { - Eigen::VectorX sum = Eigen::VectorX::Zero(Eigen::Index(mio::abm::InfectionState::Count)); auto curr_time = sim.get_time(); - PRAGMA_OMP(for) for (auto& location : sim.get_model().get_locations()) { for (uint32_t inf_state = 0; inf_state < (int)mio::abm::InfectionState::Count; inf_state++) { sum[inf_state] += sim.get_model().get_subpopulation(location.get_id(), curr_time, @@ -189,60 +192,68 @@ struct LogInfectionState : mio::LogAlways { } }; +/** +* @brief Logger to log the TimeSeries of the number of Person%s in an #InfectionState per AgeGroup. +* This is a finer grained variant of LogInfectionState. It is not used by ResultSimulation, but provided here so that +* a Simulation can be given an age resolved History where that resolution is needed. +*/ struct LogInfectionStatePerAgeGroup : mio::LogAlways { using Type = std::pair; - /** - * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. + /** + * @brief Log the TimeSeries of the number of Person%s in an #InfectionState per AgeGroup. * @param[in] sim The simulation of the abm. - * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. + * @return A pair of the TimePoint and a vector counting the Person%s per AgeGroup and #InfectionState, + * indexed by `age_group * InfectionState::Count + infection_state`. */ static Type log(const mio::abm::Simulation<>& sim) { - Eigen::VectorXd sum = Eigen::VectorXd::Zero( Eigen::Index((size_t)mio::abm::InfectionState::Count * sim.get_model().parameters.get_num_groups())); const auto curr_time = sim.get_time(); - const auto persons = sim.get_model().get_persons(); - - for (auto i = size_t(0); i < persons.size(); ++i) { - auto& p = persons[i]; - auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) + - ((uint32_t)p.get_infection_state(curr_time)); - // PRAGMA_OMP(atomic) - sum[index] += 1; + for (const Person& p : sim.get_model().get_persons()) { + auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) + + ((uint32_t)p.get_infection_state(curr_time)); + sum[index] += 1; } return std::make_pair(curr_time, sum); } }; +/** +* @brief Logger to log the TimeSeries of new #Infection%s per LocationType and AgeGroup. +*/ struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { using Type = std::pair; - /** - * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. + /** + * @brief Log the TimeSeries of new #Infection%s per LocationType and AgeGroup. + * A Person is counted if it became #InfectionState::Exposed since the previous time step. It is attributed to the + * LocationType it is at when the transition is observed. At the first time step there is no previous time step to + * compare against, so no Person is counted. This avoids miscounting the initially infected Person%s of the + * Model as new #Infection%s. * @param[in] sim The simulation of the abm. - * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState. + * @return A pair of the TimePoint and a vector counting the newly exposed Person%s per AgeGroup and LocationType, + * indexed by `age_group * LocationType::Count + location_type`. */ static Type log(const mio::abm::Simulation<>& sim) { - Eigen::VectorXd sum = Eigen::VectorXd::Zero( Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_model().parameters.get_num_groups())); - auto curr_time = sim.get_time(); - auto prev_time = sim.get_prev_time(); - const auto persons = sim.get_model().get_persons(); - - - for (auto i = size_t(0); i < persons.size(); ++i) { - auto& p = persons[i]; + auto curr_time = sim.get_time(); + auto prev_time = sim.get_prev_time(); + // Before the first time step is evolved there is no previous state to compare against. Persons that are + // infected at initialization must not be counted as new Infections. + if (prev_time >= curr_time) { + return std::make_pair(curr_time, sum); + } + for (const Person& p : sim.get_model().get_persons()) { if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) && (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) + - ((uint32_t)p.get_location_type()); + ((uint32_t)p.get_location_type()); sum[index] += 1; - } } return std::make_pair(curr_time, sum); diff --git a/cpp/models/abm/result_simulation.h b/cpp/models/abm/result_simulation.h index 02aeffe069..c0073d562a 100644 --- a/cpp/models/abm/result_simulation.h +++ b/cpp/models/abm/result_simulation.h @@ -49,29 +49,27 @@ class ResultSimulation : public Simulation } /** - * @brief Return the simulation result aggregated by infection states. + * @brief Return the simulation result aggregated by #InfectionState. */ - mio::TimeSeries get_result() const + const mio::TimeSeries& get_result() const { return std::get<0>(history.get_log()); } - /** - * @brief Return the detailed simulation result aggregated by infection states. + /** + * @brief Return the detailed simulation result, i.e. the new #Infection%s aggregated by LocationType and AgeGroup. */ - mio::TimeSeries get_result_detailed() const + const mio::TimeSeries& get_result_detailed() const { return std::get<0>(history_detailed.get_log()); } - mio::History history{ - Eigen::Index(InfectionState::Count)}; + Eigen::Index(InfectionState::Count)}; ///< History used to create the result TimeSeries. mio::History history_detailed{ - Eigen::Index(LocationType::Count) * this->get_model().parameters.get_num_groups()}; - - + Eigen::Index(LocationType::Count) * + this->get_model().parameters.get_num_groups()}; ///< History used to create the detailed result TimeSeries. }; } // namespace abm diff --git a/cpp/models/abm/simulation.h b/cpp/models/abm/simulation.h index 283fcd26bd..da66e734eb 100644 --- a/cpp/models/abm/simulation.h +++ b/cpp/models/abm/simulation.h @@ -46,7 +46,9 @@ class Simulation Simulation(TimePoint t0, Model&& model) : m_model(std::move(model)) , m_t(t0) - , m_t_prev(t0-hours(1)) + // Before the first time step is evolved there is no previous time. m_t_prev is set to m_t so that loggers can + // detect this case via get_prev_time() >= get_time() and skip transition-based counting. + , m_t_prev(t0) , m_dt(hours(1)) { } @@ -87,14 +89,14 @@ class Simulation } /** - * @brief Get the previous time of the Simulation. + * @brief Get the time of the Simulation before the last evolve step. + * Equal to get_time() as long as no time step has been evolved yet. */ TimePoint get_prev_time() const { return m_t_prev; } - /** * @brief Get the Model that this Simulation evolves. */ @@ -115,12 +117,11 @@ class Simulation m_model.evolve(m_t, dt); m_t_prev = m_t; m_t += m_dt; - } Model m_model; ///< The Model to simulate. TimePoint m_t; ///< The current TimePoint of the Simulation. - TimePoint m_t_prev; ///< The previous TimePoint of the Simulation. + TimePoint m_t_prev; ///< The TimePoint of the Simulation before the last evolve step. TimeSpan m_dt; ///< The length of the time steps. }; diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 3a4db0f75c..4ce5b3a2cc 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -35,8 +35,8 @@ # location types from simulation results of the agent-based model (ABM) stored in HDF5 format. # The used Loggers are: -# LogInfectionStatePerAgeGroup -# LogInfectionPerLocationTypePerAgeGroup +# LogInfectionState +# LogInfectionPerLocationTypePerAgeGroup # The output of the loggers of several runs is stored in HDF5 files using mio::save_results in mio/io/result_io.h, see abm_history_object.cpp. # Adjust these as needed. @@ -127,7 +127,7 @@ def plot_infections_loc_types_average( xtick_step=150): """ Plots rolling sum of new infections per 24 hours location type for the median run. - @param[in] base_path Path to results directory. + @param[in] path_to_loc_types Path to results directory. @param[in] start_date Start date as string. @param[in] colormap Matplotlib colormap. @param[in] smooth_sigma Sigma for Gaussian smoothing. @@ -291,7 +291,7 @@ def plot_infection_states_by_age_group( # Dynamically detect available age groups from the data available_groups = [key for key in p50_bs.keys() if key.startswith('Group') or key == 'Total'] # Sort to ensure Group1, Group2, ..., Total order - available_groups = sorted(available_groups, key=lambda x: (x != 'Total', x)) + available_groups = sorted(available_groups, key=lambda x: (x == 'Total', x)) color_plot = matplotlib.colormaps.get_cmap(colormap).colors n_states = len(state_labels) @@ -375,8 +375,8 @@ def main(): default='Set1', help="Matplotlib colormap") parser.add_argument("--xtick-step", type=int, default=150, help="Step for x-axis ticks (usually hours)") - parser.add_argument("--90percentile", action="store_true", - help="If set, plot 90% percentile as well") + parser.add_argument("--show-90percentile", action="store_true", + help="If set, plot 90%% percentile as well") args = parser.parse_args() if not args.path_to_infection_states and not args.path_to_loc_types: @@ -389,7 +389,7 @@ def main(): start_date=args.start_date, colormap=args.colormap, xtick_step=args.xtick_step, - show90=True + show90=args.show_90percentile ) if args.path_to_loc_types: diff --git a/pycode/memilio-plot/tests/test_plot_plotAbmInfectionStates.py b/pycode/memilio-plot/tests/test_plot_plotAbmInfectionStates.py index 4bc5dde5e1..8e137de955 100644 --- a/pycode/memilio-plot/tests/test_plot_plotAbmInfectionStates.py +++ b/pycode/memilio-plot/tests/test_plot_plotAbmInfectionStates.py @@ -18,6 +18,7 @@ # limitations under the License. ############################################################################# +import os import unittest from unittest.mock import patch, MagicMock import numpy as np @@ -30,17 +31,37 @@ class TestPlotAbmInfectionStates(unittest.TestCase): @patch('memilio.plot.plotAbmInfectionStates.h5py.File') def test_load_h5_results(self, mock_h5file): + # load_h5_results checks that the results file exists before opening it, so os.path.exists has to be + # patched as well. Only the flat layout (Results_p50.h5) is reported as existing here. mock_group = {'Time': np.arange(10), 'Total': np.ones((10, 8))} mock_h5file().__enter__().get.return_value = {'0': mock_group} mock_h5file().__enter__().items.return_value = [('0', mock_group)] mock_h5file().__enter__().__getitem__.return_value = mock_group - with patch('memilio.plot.plotAbmInfectionStates.h5py.File', mock_h5file): + with patch('memilio.plot.plotAbmInfectionStates.os.path.exists', + side_effect=lambda p: p.endswith('Results_p50.h5')): result = abm.load_h5_results('dummy_path', 'p50') assert 'Time' in result assert 'Total' in result np.testing.assert_array_equal(result['Time'], np.arange(10)) np.testing.assert_array_equal(result['Total'], np.ones((10, 8))) + @patch('memilio.plot.plotAbmInfectionStates.h5py.File') + def test_load_h5_results_nested_layout(self, mock_h5file): + # The nested layout (p50/Results.h5) is supported as a fallback. + mock_group = {'Time': np.arange(10), 'Total': np.ones((10, 8))} + mock_h5file().__enter__().items.return_value = [('0', mock_group)] + mock_h5file().__enter__().__getitem__.return_value = mock_group + with patch('memilio.plot.plotAbmInfectionStates.os.path.exists', + side_effect=lambda p: p.endswith(os.path.join('p50', 'Results.h5'))): + result = abm.load_h5_results('dummy_path', 'p50') + np.testing.assert_array_equal(result['Time'], np.arange(10)) + + def test_load_h5_results_missing_file(self): + # Neither layout exists, so a FileNotFoundError explaining both is raised. + with patch('memilio.plot.plotAbmInfectionStates.os.path.exists', return_value=False): + with self.assertRaises(FileNotFoundError): + abm.load_h5_results('dummy_path', 'p50') + @patch('memilio.plot.plotAbmInfectionStates.load_h5_results') @patch('memilio.plot.plotAbmInfectionStates.matplotlib') @patch('memilio.plot.plotAbmInfectionStates.gaussian_filter1d', side_effect=lambda x, sigma, mode: x) @@ -71,19 +92,17 @@ def test_plot_infections_loc_types_average(self, mock_df, mock_gauss, mock_matpl assert mock_ax.set_xticks.called assert mock_ax.set_xticklabels.called - # Test figure settings - self.assertEqual(mock_figure.call_count, 2, - "figure should be called twice") - # Verify first call is with the figure name - mock_figure.assert_any_call('Infection_location_types') - # Verify second call is without arguments (for autofmt_xdate) - mock_figure.assert_any_call() + # Test figure settings. The figure is created once, with its name. Showing the figure is left to the + # caller, so plot_infections_loc_types_average must not call plt.show itself. + self.assertEqual(mock_figure.call_count, 1, + "figure should be called once") + mock_figure.assert_called_once_with('Infection_location_types') mock_title.assert_called_once_with( 'Number of new infections per location type for the median run, rolling sum over 24 hours') mock_legend.assert_called_once() mock_xlabel.assert_called_once_with('Date') mock_ylabel.assert_called_once_with('Number of individuals') - mock_show.assert_called_once() + mock_show.assert_not_called() # Verify legend was called with location type labels legend_call_args = mock_legend.call_args From d40ead4859e5eb259dc364c1641c962848b0a791 Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Fri, 31 Jul 2026 19:58:13 +0200 Subject: [PATCH 32/38] delete comment --- cpp/models/abm/simulation.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/cpp/models/abm/simulation.h b/cpp/models/abm/simulation.h index da66e734eb..03ea23e2ce 100644 --- a/cpp/models/abm/simulation.h +++ b/cpp/models/abm/simulation.h @@ -46,8 +46,6 @@ class Simulation Simulation(TimePoint t0, Model&& model) : m_model(std::move(model)) , m_t(t0) - // Before the first time step is evolved there is no previous time. m_t_prev is set to m_t so that loggers can - // detect this case via get_prev_time() >= get_time() and skip transition-based counting. , m_t_prev(t0) , m_dt(hours(1)) { From 19b306e9194556496ec55771082cbc09823ca4c2 Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Fri, 31 Jul 2026 20:01:57 +0200 Subject: [PATCH 33/38] format --- .../memilio/plot/plotAbmInfectionStates.py | 28 +++++++++++-------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py index 4ce5b3a2cc..b948851ffa 100644 --- a/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py +++ b/pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py @@ -86,7 +86,7 @@ def load_h5_results(base_path, percentile): file_path_flat = os.path.join(base_path, f"Results_{percentile}.h5") # Try subdirectory structure (e.g., p50/Results.h5) file_path_nested = os.path.join(base_path, percentile, "Results.h5") - + # Determine which file exists if os.path.exists(file_path_flat): file_path = file_path_flat @@ -112,7 +112,7 @@ def load_h5_results(base_path, percentile): f" p75/Results.h5\n" f" p95/Results.h5\n" ) - + with h5py.File(file_path, 'r') as f: data = {k: v[()] for k, v in f['0'].items()} return data @@ -147,11 +147,12 @@ def plot_infections_loc_types_average( # Check if data dimensions match expected location types num_cols = total_50.shape[1] if len(total_50.shape) > 1 else 1 num_location_types = len(location_type_labels) - + if num_cols < num_location_types: - print(f"Warning: Data has {num_cols} columns but {num_location_types} location types defined.") + print( + f"Warning: Data has {num_cols} columns but {num_location_types} location types defined.") print(f"Only plotting first {num_cols} location types.") - + for idx, i in enumerate(location_type_labels.keys()): if i >= num_cols: break # Skip if we don't have data for this location type @@ -164,7 +165,8 @@ def plot_infections_loc_types_average( window=indexer, min_periods=1).sum().to_numpy() y = y[0::rolling_window].flatten() y = gaussian_filter1d(y, sigma=smooth_sigma, mode='nearest') - plt.plot(time[0::rolling_window], y, color=color, linewidth=2.5, label=location_type_labels[i]) + plt.plot(time[0::rolling_window], y, color=color, + linewidth=2.5, label=location_type_labels[i]) plt.legend() # Auto-generate legend from plot labels _format_x_axis(time, start_date, xtick_step) @@ -289,10 +291,12 @@ def plot_infection_states_by_age_group( """ # Dynamically detect available age groups from the data - available_groups = [key for key in p50_bs.keys() if key.startswith('Group') or key == 'Total'] + available_groups = [key for key in p50_bs.keys( + ) if key.startswith('Group') or key == 'Total'] # Sort to ensure Group1, Group2, ..., Total order - available_groups = sorted(available_groups, key=lambda x: (x == 'Total', x)) - + available_groups = sorted( + available_groups, key=lambda x: (x == 'Total', x)) + color_plot = matplotlib.colormaps.get_cmap(colormap).colors n_states = len(state_labels) fig, ax = plt.subplots( @@ -304,10 +308,10 @@ def plot_infection_states_by_age_group( y75 = p75_bs[group] y05 = p05_bs[group] if (show90 and p05_bs is not None) else None y95 = p95_bs[group] if (show90 and p95_bs is not None) else None - + # Get group label, using default if not in predefined dict group_label = age_groups_dict.get(group, group) - + for row_idx, (state_idx, label) in enumerate(state_labels.items()): _plot_state( ax[row_idx, col_idx], x, y50[:, state_idx], y25[:, @@ -391,7 +395,7 @@ def main(): xtick_step=args.xtick_step, show90=args.show_90percentile ) - + if args.path_to_loc_types: plot_infections_loc_types_average( path_to_loc_types=args.path_to_loc_types, From faf0f4048fcfd7dfb71a501b670321b5191128f0 Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Fri, 31 Jul 2026 20:38:36 +0200 Subject: [PATCH 34/38] fix test --- cpp/tests/test_abm_mobility_rules.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/tests/test_abm_mobility_rules.cpp b/cpp/tests/test_abm_mobility_rules.cpp index f5bd9fcc3a..469bcba68f 100644 --- a/cpp/tests/test_abm_mobility_rules.cpp +++ b/cpp/tests/test_abm_mobility_rules.cpp @@ -495,8 +495,9 @@ TEST_F(TestMobilityRules, hospital) EXPECT_EQ(mio::abm::go_to_hospital(rng_inf, p_inf, t, dt, mio::abm::Parameters(num_age_groups)), mio::abm::LocationType::Hospital); + // The infection has to start at t, otherwise it may already have progressed to InfectedSevere by then. auto p_car = - make_test_person(this->get_rng(), home, age_group_15_to_34, mio::abm::InfectionState::InfectedSymptoms); + make_test_person(this->get_rng(), home, age_group_15_to_34, mio::abm::InfectionState::InfectedSymptoms, t); auto rng_car = mio::abm::PersonalRandomNumberGenerator(this->get_rng(), p_car); // Ensure person has infection symptoms still stay at home EXPECT_EQ(mio::abm::go_to_hospital(rng_car, p_car, t, dt, mio::abm::Parameters(num_age_groups)), From 650985ffcbe562ebcd7cb75ce9136f43bee133bb Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Thu, 13 Aug 2026 13:57:32 +0200 Subject: [PATCH 35/38] update --- cpp/tests/test_abm_simulation.cpp | 82 +++++++++++++++++++++++++++++++ 1 file changed, 82 insertions(+) diff --git a/cpp/tests/test_abm_simulation.cpp b/cpp/tests/test_abm_simulation.cpp index 8fd25b66aa..fd1fa9d80e 100644 --- a/cpp/tests/test_abm_simulation.cpp +++ b/cpp/tests/test_abm_simulation.cpp @@ -171,4 +171,86 @@ TEST(TestSimulation, ResultSimulation) } EXPECT_EQ(sim.get_result().get_value(0)[(Eigen::Index)mio::abm::InfectionState::Susceptible], 1.0); EXPECT_EQ(sim.get_result().get_value(N - 1)[(Eigen::Index)mio::abm::InfectionState::Susceptible], 1.0); + + // the detailed result has the same time points, but is resolved by LocationType and AgeGroup + ASSERT_EQ(sim.get_result_detailed().get_num_time_points(), N); + EXPECT_THAT(sim.get_result_detailed().get_times(), ElementsAreLinspace(t0.days(), tmax.days(), N)); + EXPECT_EQ(sim.get_result_detailed().get_num_elements(), + Eigen::Index((size_t)mio::abm::LocationType::Count * num_age_groups)); + // nobody gets infected in this setup, so no new Infection is ever counted + for (const auto& tp : sim.get_result_detailed()) { + EXPECT_EQ(tp.sum(), 0.0); + } +} + +TEST(TestSimulation, resultSimulationDetailedCountsNewInfections) +{ + // a Person that becomes Exposed during the Simulation must show up in the detailed result at its LocationType and + // AgeGroup, while a Person that is already Exposed at t0 must not be counted. + auto model = mio::abm::Model(num_age_groups); + auto home = model.add_location(mio::abm::LocationType::Home); + const auto t0 = mio::abm::TimePoint(0); + + // an Exposed Person at initialization, and a second one that only becomes Exposed after the first time step + auto initially_exposed = add_test_person(model, home, age_group_15_to_34, mio::abm::InfectionState::Exposed, t0); + auto newly_exposed = add_test_person(model, home, age_group_35_to_59, mio::abm::InfectionState::Susceptible, t0); + model.assign_location(initially_exposed, home); + model.assign_location(newly_exposed, home); + + // give the second Person an Infection that starts in the Exposed state one hour in, i.e. after the first step + const auto t_infection = t0 + mio::abm::hours(1); + auto& person = model.get_person(newly_exposed); + auto rng_person = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person); + person.add_new_infection(mio::abm::Infection(rng_person, mio::abm::VirusVariant::Wildtype, person.get_age(), + model.parameters, t_infection, + mio::abm::InfectionState::Exposed)); + + auto sim = mio::abm::ResultSimulation(std::move(model), t0); + sim.advance(t0 + mio::abm::hours(2)); + + const auto& detailed = sim.get_result_detailed(); + const auto index = Eigen::Index((size_t)mio::abm::LocationType::Count * (size_t)age_group_35_to_59.get() + + (size_t)mio::abm::LocationType::Home); + + // the initially Exposed Person is not counted, since there is no previous time step at t0 + EXPECT_EQ(detailed.get_value(0).sum(), 0.0); + // the new Infection of the second Person is counted exactly once, at Home and in its AgeGroup + double total = 0.0; + for (Eigen::Index i = 1; i < detailed.get_num_time_points(); ++i) { + total += detailed.get_value(i).sum(); + } + EXPECT_EQ(total, 1.0); + EXPECT_EQ(detailed.get_value(1)[index], 1.0); +} + +TEST(TestSimulation, logInfectionStatePerAgeGroup) +{ + // LogInfectionStatePerAgeGroup counts every Person by AgeGroup and InfectionState + auto model = mio::abm::Model(num_age_groups); + auto home = model.add_location(mio::abm::LocationType::Home); + + auto p1 = add_test_person(model, home, age_group_5_to_14, mio::abm::InfectionState::Susceptible); + auto p2 = add_test_person(model, home, age_group_5_to_14, mio::abm::InfectionState::Susceptible); + auto p3 = add_test_person(model, home, age_group_60_to_79, mio::abm::InfectionState::Dead); + model.assign_location(p1, home); + model.assign_location(p2, home); + model.assign_location(p3, home); + + auto sim = mio::abm::Simulation(mio::abm::TimePoint(0), std::move(model)); + mio::History history{ + Eigen::Index((size_t)mio::abm::InfectionState::Count * num_age_groups)}; + sim.advance(mio::abm::TimePoint(0) + mio::abm::hours(1), history); + + const auto& log = std::get<0>(history.get_log()); + ASSERT_EQ(log.get_num_time_points(), 2); + EXPECT_EQ(log.get_num_elements(), Eigen::Index((size_t)mio::abm::InfectionState::Count * num_age_groups)); + + const auto idx = [](mio::AgeGroup age, mio::abm::InfectionState state) { + return Eigen::Index((size_t)mio::abm::InfectionState::Count * (size_t)age.get() + (size_t)state); + }; + for (Eigen::Index i = 0; i < log.get_num_time_points(); ++i) { + EXPECT_EQ(log.get_value(i).sum(), 3.0); + EXPECT_EQ(log.get_value(i)[idx(age_group_5_to_14, mio::abm::InfectionState::Susceptible)], 2.0); + EXPECT_EQ(log.get_value(i)[idx(age_group_60_to_79, mio::abm::InfectionState::Dead)], 1.0); + } } From d8102a23fee4981944e6b308d17f8168ef52573b Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Thu, 8 Oct 2026 18:28:45 +0200 Subject: [PATCH 36/38] 1210 Address review: move study histories to example, fix new infection logger - Revert ResultSimulation and Simulation to main; the parameter study example now uses a small StudySimulation wrapper holding its histories. - LogInfectionPerLocationTypePerAgeGroup keeps its previous log time and location types itself and attributes new infections to the location where the interaction happened, not to where the person moved afterwards. - Gather the detailed results across MPI ranks by returning them from the processing lambda. - Handle the IOResult of save_result instead of discarding it. - LogInfectionState counts in a single pass over all persons (#1406). --- cpp/examples/abm_parameter_study.cpp | 81 ++++++++++++++++++++------ cpp/models/abm/common_abm_loggers.h | 57 ++++++++++++------- cpp/models/abm/result_simulation.h | 18 +----- cpp/models/abm/simulation.h | 12 ---- cpp/tests/test_abm_simulation.cpp | 85 ++++++++++++++++------------ 5 files changed, 152 insertions(+), 101 deletions(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 61dd73d93e..52015b5409 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -17,10 +17,11 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#include "abm/result_simulation.h" +#include "abm/common_abm_loggers.h" #include "abm/household.h" #include "abm/lockdown_rules.h" #include "abm/model.h" +#include "abm/simulation.h" #include "abm/time.h" #include "memilio/compartments/parameter_studies.h" @@ -165,6 +166,44 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) return model; } +/** + * @brief Simulation for this study, holding the Histories needed to evaluate a run. + * ParameterStudy only calls advance(tmax), so the Simulation has to own its Histories. Besides the result aggregated by + * #InfectionState, it logs the new Infections per LocationType and AgeGroup. + */ +class StudySimulation : public mio::abm::Simulation<> +{ +public: + StudySimulation(mio::abm::Model&& model, mio::abm::TimePoint t0) + : Simulation(t0, std::move(model)) + , history_detailed{Eigen::Index((size_t)mio::abm::LocationType::Count * + get_model().parameters.get_num_groups())} + { + } + + void advance(mio::abm::TimePoint tmax) + { + Simulation::advance(tmax, history, history_detailed); + } + + /// @brief The result aggregated by #InfectionState. + const mio::TimeSeries& get_result() const + { + return std::get<0>(history.get_log()); + } + + /// @brief The new Infections per LocationType and AgeGroup. + const mio::TimeSeries& get_result_detailed() const + { + return std::get<0>(history_detailed.get_log()); + } + +private: + mio::History history{ + Eigen::Index(mio::abm::InfectionState::Count)}; + mio::History history_detailed; +}; + int main() { mio::mpi::init(); @@ -200,9 +239,6 @@ int main() const auto result_dir_standard = mio::create_directories_or_exit(result_dir / "standard_results"); const auto result_dir_detailed = mio::create_directories_or_exit(result_dir / "detailed_results"); - // Collects the detailed result of each run, see the process_simulation_result lambda below. - std::vector>> ensemble_results_detailed; - // Run the study // The first lambda ("create_simulation" argument) sets up the simulation, the second ("process_simulation_result") // allows us to process each simulations result. Be mindful of the memory used for storing these results! @@ -213,9 +249,9 @@ int main() // per person to 2^16 = 65536 const auto ctr = mio::Counter(static_cast(run_idx) << 16); copy.reset_rng(ctr); - return mio::abm::ResultSimulation(std::move(copy), t0_); + return StudySimulation(std::move(copy), t0_); }, - [&result_dir_standard, &result_dir_detailed, &ensemble_results_detailed](auto&& sim, auto&& run_idx) { + [&result_dir_standard, &result_dir_detailed](auto&& sim, auto&& run_idx) { auto interpolated_result = mio::interpolate_simulation_result(sim.get_result()); auto interpolated_result_detailed = mio::interpolate_simulation_result(sim.get_result_detailed()); @@ -231,20 +267,25 @@ int main() sim.get_result_detailed().print_table(outfile_run_detailed, {}, 7, 4); std::cout << "Detailed results written to " << outpath_detailed.string() << std::endl; - // The detailed result is collected separately, because ensemble_percentile requires all entries of an - // ensemble to have the same number of elements, which the two results do not have. - ensemble_results_detailed.push_back({interpolated_result_detailed}); - - return std::vector{interpolated_result}; + // Both results are returned, so that the study gathers them from all MPI ranks on the root rank. + return std::vector{interpolated_result, interpolated_result_detailed}; }); // The study collects all results on the root rank, so we only process the results there if (mio::mpi::is_root()) { + // Split the gathered results into one ensemble per result, because ensemble_percentile requires all entries + // of an ensemble to have the same number of elements, which the two results do not have. + std::vector>> ensemble_results_standard, ensemble_results_detailed; + for (auto& run_result : ensemble_results) { + ensemble_results_standard.push_back({std::move(run_result[0])}); + ensemble_results_detailed.push_back({std::move(run_result[1])}); + } + // The percentiles are written in the "Results_p05.h5" layout expected by // pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py, so that the results of this example can be // plotted directly. The percentile is given in whole percent to match that file name. const auto write_percentile = [&](int p) { - auto ensemble_percentiles = ensemble_percentile(ensemble_results, p / 100.0); + auto ensemble_percentiles = ensemble_percentile(ensemble_results_standard, p / 100.0); auto ensemble_percentiles_detailed = ensemble_percentile(ensemble_results_detailed, p / 100.0); std::ofstream out(result_dir_standard / fmt::format("Results_p{:02d}.txt", p)); @@ -254,10 +295,18 @@ int main() // save_result splits each row into num_groups groups. Both results are a single vector per time point, // so they are written as one group each, whose "Total" holds all entries of that vector. Passing the // number of age groups here would instead split those entries across that many groups. - mio::unused(mio::save_result(ensemble_percentiles, {0}, 1, - (result_dir_standard / fmt::format("Results_p{:02d}.h5", p)).string())); - mio::unused(mio::save_result(ensemble_percentiles_detailed, {0}, 1, - (result_dir_detailed / fmt::format("Results_p{:02d}.h5", p)).string())); + const auto save_status = mio::save_result( + ensemble_percentiles, {0}, 1, (result_dir_standard / fmt::format("Results_p{:02d}.h5", p)).string()); + if (!save_status) { + std::cout << "Error writing result: " << save_status.error().formatted_message() << std::endl; + } + const auto save_status_detailed = + mio::save_result(ensemble_percentiles_detailed, {0}, 1, + (result_dir_detailed / fmt::format("Results_p{:02d}.h5", p)).string()); + if (!save_status_detailed) { + std::cout << "Error writing detailed result: " << save_status_detailed.error().formatted_message() + << std::endl; + } }; write_percentile(5); diff --git a/cpp/models/abm/common_abm_loggers.h b/cpp/models/abm/common_abm_loggers.h index fb25e1d23d..618933abbd 100644 --- a/cpp/models/abm/common_abm_loggers.h +++ b/cpp/models/abm/common_abm_loggers.h @@ -30,6 +30,8 @@ #include "abm/mobility_data.h" #include "memilio/utils/mioomp.h" +#include + namespace mio { namespace abm @@ -181,11 +183,14 @@ struct LogInfectionState : mio::LogAlways { { Eigen::VectorX sum = Eigen::VectorX::Zero(Eigen::Index(mio::abm::InfectionState::Count)); - auto curr_time = sim.get_time(); - for (auto& location : sim.get_model().get_locations()) { - for (uint32_t inf_state = 0; inf_state < (int)mio::abm::InfectionState::Count; inf_state++) { - sum[inf_state] += sim.get_model().get_subpopulation(location.get_id(), curr_time, - mio::abm::InfectionState(inf_state)); + const auto curr_time = sim.get_time(); + const auto& model = sim.get_model(); + // A single pass over all Person%s, instead of calling get_subpopulation for every Location and + // #InfectionState, which would iterate over all Person%s each time. Person%s currently located in another + // Model (e.g. in a graph ABM) are not counted, same as in get_subpopulation. + for (const Person& p : model.get_persons()) { + if (p.get_location_model_id() == model.get_id()) { + sum[(Eigen::Index)p.get_infection_state(curr_time)] += 1; } } return std::make_pair(curr_time, sum); @@ -222,42 +227,52 @@ struct LogInfectionStatePerAgeGroup : mio::LogAlways { /** * @brief Logger to log the TimeSeries of new #Infection%s per LocationType and AgeGroup. +* Unlike the other loggers, this one keeps a state: the time and the LocationType of each Person at its previous log. +* Hence, it has to log after every time step, which is the case when used in a History passed to Simulation::advance. */ struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { using Type = std::pair; /** * @brief Log the TimeSeries of new #Infection%s per LocationType and AgeGroup. - * A Person is counted if it became #InfectionState::Exposed since the previous time step. It is attributed to the - * LocationType it is at when the transition is observed. At the first time step there is no previous time step to - * compare against, so no Person is counted. This avoids miscounting the initially infected Person%s of the - * Model as new #Infection%s. + * A Person is counted if it became #InfectionState::Exposed since the previous log. A Person is infected during the + * interaction at the beginning of a time step and may move to another Location afterwards. Hence, it is attributed + * to the LocationType it was at during the previous log, not to its current one. + * At the first log there is no previous log to compare against, so no Person is counted. This avoids miscounting + * the initially infected Person%s of the Model as new #Infection%s. * @param[in] sim The simulation of the abm. * @return A pair of the TimePoint and a vector counting the newly exposed Person%s per AgeGroup and LocationType, * indexed by `age_group * LocationType::Count + location_type`. */ - static Type log(const mio::abm::Simulation<>& sim) + Type log(const mio::abm::Simulation<>& sim) { Eigen::VectorXd sum = Eigen::VectorXd::Zero( Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_model().parameters.get_num_groups())); - auto curr_time = sim.get_time(); - auto prev_time = sim.get_prev_time(); + const auto curr_time = sim.get_time(); + const auto persons = sim.get_model().get_persons(); + const size_t num_persons = persons.size(); + // Only Persons that were already present at the previous log can be compared against it. + const size_t num_previous = m_prev_location_types.size(); + m_prev_location_types.resize(num_persons); - // Before the first time step is evolved there is no previous state to compare against. Persons that are - // infected at initialization must not be counted as new Infections. - if (prev_time >= curr_time) { - return std::make_pair(curr_time, sum); - } - - for (const Person& p : sim.get_model().get_persons()) { - if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) && + size_t i = 0; + for (const Person& p : persons) { + if (i < num_previous && (p.get_infection_state(m_prev_time) != mio::abm::InfectionState::Exposed) && (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) + - ((uint32_t)p.get_location_type()); + ((uint32_t)m_prev_location_types[i]); sum[index] += 1; } + // Remember where the Person is now, which is where the interaction of the next time step happens. + m_prev_location_types[i] = p.get_location_type(); + ++i; } + m_prev_time = curr_time; return std::make_pair(curr_time, sum); } + +private: + mio::abm::TimePoint m_prev_time{0}; ///< Time of the previous log. + std::vector m_prev_location_types; ///< LocationType of each Person at the previous log. }; diff --git a/cpp/models/abm/result_simulation.h b/cpp/models/abm/result_simulation.h index c0073d562a..2f1d02131a 100644 --- a/cpp/models/abm/result_simulation.h +++ b/cpp/models/abm/result_simulation.h @@ -45,31 +45,19 @@ class ResultSimulation : public Simulation */ void advance(TimePoint tmax) { - Simulation::advance(tmax, history, history_detailed); + Simulation::advance(tmax, history); } /** - * @brief Return the simulation result aggregated by #InfectionState. + * @brief Return the simulation result aggregated by infection states. */ const mio::TimeSeries& get_result() const { - return std::get<0>(history.get_log()); - } - - /** - * @brief Return the detailed simulation result, i.e. the new #Infection%s aggregated by LocationType and AgeGroup. - */ - const mio::TimeSeries& get_result_detailed() const - { - return std::get<0>(history_detailed.get_log()); + return get<0>(history.get_log()); } mio::History history{ Eigen::Index(InfectionState::Count)}; ///< History used to create the result TimeSeries. - - mio::History history_detailed{ - Eigen::Index(LocationType::Count) * - this->get_model().parameters.get_num_groups()}; ///< History used to create the detailed result TimeSeries. }; } // namespace abm diff --git a/cpp/models/abm/simulation.h b/cpp/models/abm/simulation.h index 03ea23e2ce..8e6ff52839 100644 --- a/cpp/models/abm/simulation.h +++ b/cpp/models/abm/simulation.h @@ -46,7 +46,6 @@ class Simulation Simulation(TimePoint t0, Model&& model) : m_model(std::move(model)) , m_t(t0) - , m_t_prev(t0) , m_dt(hours(1)) { } @@ -86,15 +85,6 @@ class Simulation return m_t; } - /** - * @brief Get the time of the Simulation before the last evolve step. - * Equal to get_time() as long as no time step has been evolved yet. - */ - TimePoint get_prev_time() const - { - return m_t_prev; - } - /** * @brief Get the Model that this Simulation evolves. */ @@ -113,13 +103,11 @@ class Simulation { auto dt = std::min(m_dt, tmax - m_t); m_model.evolve(m_t, dt); - m_t_prev = m_t; m_t += m_dt; } Model m_model; ///< The Model to simulate. TimePoint m_t; ///< The current TimePoint of the Simulation. - TimePoint m_t_prev; ///< The TimePoint of the Simulation before the last evolve step. TimeSpan m_dt; ///< The length of the time steps. }; diff --git a/cpp/tests/test_abm_simulation.cpp b/cpp/tests/test_abm_simulation.cpp index fd1fa9d80e..d350dcbcf6 100644 --- a/cpp/tests/test_abm_simulation.cpp +++ b/cpp/tests/test_abm_simulation.cpp @@ -171,56 +171,67 @@ TEST(TestSimulation, ResultSimulation) } EXPECT_EQ(sim.get_result().get_value(0)[(Eigen::Index)mio::abm::InfectionState::Susceptible], 1.0); EXPECT_EQ(sim.get_result().get_value(N - 1)[(Eigen::Index)mio::abm::InfectionState::Susceptible], 1.0); - - // the detailed result has the same time points, but is resolved by LocationType and AgeGroup - ASSERT_EQ(sim.get_result_detailed().get_num_time_points(), N); - EXPECT_THAT(sim.get_result_detailed().get_times(), ElementsAreLinspace(t0.days(), tmax.days(), N)); - EXPECT_EQ(sim.get_result_detailed().get_num_elements(), - Eigen::Index((size_t)mio::abm::LocationType::Count * num_age_groups)); - // nobody gets infected in this setup, so no new Infection is ever counted - for (const auto& tp : sim.get_result_detailed()) { - EXPECT_EQ(tp.sum(), 0.0); - } } -TEST(TestSimulation, resultSimulationDetailedCountsNewInfections) +TEST(TestSimulation, logInfectionPerLocationTypePerAgeGroup) { - // a Person that becomes Exposed during the Simulation must show up in the detailed result at its LocationType and + // a Person that becomes Exposed must be counted at the LocationType it was at when it got infected and in its // AgeGroup, while a Person that is already Exposed at t0 must not be counted. auto model = mio::abm::Model(num_age_groups); auto home = model.add_location(mio::abm::LocationType::Home); + auto work = model.add_location(mio::abm::LocationType::Work); const auto t0 = mio::abm::TimePoint(0); + // the Persons only move when we tell them to + model.use_mobility_rules(false); - // an Exposed Person at initialization, and a second one that only becomes Exposed after the first time step + // an Exposed Person at initialization, one that gets infected at Home and one that gets infected at Work auto initially_exposed = add_test_person(model, home, age_group_15_to_34, mio::abm::InfectionState::Exposed, t0); - auto newly_exposed = add_test_person(model, home, age_group_35_to_59, mio::abm::InfectionState::Susceptible, t0); - model.assign_location(initially_exposed, home); - model.assign_location(newly_exposed, home); - - // give the second Person an Infection that starts in the Exposed state one hour in, i.e. after the first step - const auto t_infection = t0 + mio::abm::hours(1); - auto& person = model.get_person(newly_exposed); - auto rng_person = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person); - person.add_new_infection(mio::abm::Infection(rng_person, mio::abm::VirusVariant::Wildtype, person.get_age(), - model.parameters, t_infection, - mio::abm::InfectionState::Exposed)); - - auto sim = mio::abm::ResultSimulation(std::move(model), t0); + auto stays_home = add_test_person(model, home, age_group_35_to_59, mio::abm::InfectionState::Susceptible, t0); + auto moves_home = add_test_person(model, work, age_group_60_to_79, mio::abm::InfectionState::Susceptible, t0); + for (auto id : {initially_exposed, stays_home, moves_home}) { + model.assign_location(id, home); + model.assign_location(id, work); + } + + // both new Infections start in the Exposed state during the first time step + const auto t_infection = t0 + mio::abm::minutes(30); + for (auto id : {stays_home, moves_home}) { + auto& person = model.get_person(id); + auto rng_person = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person); + person.add_new_infection(mio::abm::Infection(rng_person, mio::abm::VirusVariant::Wildtype, person.get_age(), + model.parameters, t_infection, + mio::abm::InfectionState::Exposed)); + } + + auto sim = mio::abm::Simulation(t0, std::move(model)); + mio::History history{ + Eigen::Index((size_t)mio::abm::LocationType::Count * num_age_groups)}; + + // log the initial state, evolve one time step, then move the Person infected at Work back Home before logging. + // This is what happens when a Person gets infected during the interaction of a time step and moves on in the + // mobility of the same time step. + history.log(sim); + sim.advance(t0 + mio::abm::hours(1)); + sim.get_model().change_location(moves_home, home); + history.log(sim); + // a further time step without new Infections sim.advance(t0 + mio::abm::hours(2)); + history.log(sim); - const auto& detailed = sim.get_result_detailed(); - const auto index = Eigen::Index((size_t)mio::abm::LocationType::Count * (size_t)age_group_35_to_59.get() + - (size_t)mio::abm::LocationType::Home); + const auto& detailed = std::get<0>(history.get_log()); + ASSERT_EQ(detailed.get_num_time_points(), 3); + EXPECT_EQ(detailed.get_num_elements(), Eigen::Index((size_t)mio::abm::LocationType::Count * num_age_groups)); + const auto idx = [](mio::AgeGroup age, mio::abm::LocationType type) { + return Eigen::Index((size_t)mio::abm::LocationType::Count * (size_t)age.get() + (size_t)type); + }; - // the initially Exposed Person is not counted, since there is no previous time step at t0 + // the initially Exposed Person is not counted, since there is no previous log at t0 EXPECT_EQ(detailed.get_value(0).sum(), 0.0); - // the new Infection of the second Person is counted exactly once, at Home and in its AgeGroup - double total = 0.0; - for (Eigen::Index i = 1; i < detailed.get_num_time_points(); ++i) { - total += detailed.get_value(i).sum(); - } - EXPECT_EQ(total, 1.0); - EXPECT_EQ(detailed.get_value(1)[index], 1.0); + // both new Infections are counted exactly once, where they happened, even though one Person has moved on since + EXPECT_EQ(detailed.get_value(1).sum(), 2.0); + EXPECT_EQ(detailed.get_value(1)[idx(age_group_35_to_59, mio::abm::LocationType::Home)], 1.0); + EXPECT_EQ(detailed.get_value(1)[idx(age_group_60_to_79, mio::abm::LocationType::Work)], 1.0); + EXPECT_EQ(detailed.get_value(2).sum(), 0.0); } TEST(TestSimulation, logInfectionStatePerAgeGroup) From 64998ea525aeef3cce01e9aa5e451b12a2f9c908 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Thu, 8 Oct 2026 18:31:02 +0200 Subject: [PATCH 37/38] update comment --- cpp/models/abm/common_abm_loggers.h | 4 ---- 1 file changed, 4 deletions(-) diff --git a/cpp/models/abm/common_abm_loggers.h b/cpp/models/abm/common_abm_loggers.h index 618933abbd..86c0cc820e 100644 --- a/cpp/models/abm/common_abm_loggers.h +++ b/cpp/models/abm/common_abm_loggers.h @@ -185,9 +185,6 @@ struct LogInfectionState : mio::LogAlways { Eigen::VectorX::Zero(Eigen::Index(mio::abm::InfectionState::Count)); const auto curr_time = sim.get_time(); const auto& model = sim.get_model(); - // A single pass over all Person%s, instead of calling get_subpopulation for every Location and - // #InfectionState, which would iterate over all Person%s each time. Person%s currently located in another - // Model (e.g. in a graph ABM) are not counted, same as in get_subpopulation. for (const Person& p : model.get_persons()) { if (p.get_location_model_id() == model.get_id()) { sum[(Eigen::Index)p.get_infection_state(curr_time)] += 1; @@ -275,7 +272,6 @@ struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { std::vector m_prev_location_types; ///< LocationType of each Person at the previous log. }; - /** * @brief This is like the DataWriterToMemory, but it only logs time series data. * @tparam Loggers The loggers that are used to log data. The loggers must return a touple with a TimePoint and a value. From bab6f5084ff6ddca72558affc52a1d60eda22801 Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Wed, 14 Oct 2026 10:43:59 +0200 Subject: [PATCH 38/38] 1210 Track new infections per person instead of per previous time - LogInfectionPerLocationTypePerAgeGroup remembers the InfectionState, LocationType and PersonId of each Person at the previous log, instead of the previous time. This removes the dependence on a start time of 0. - Persons without a previous log, at the first log or added since the previous log, are not counted but tracked from then on. An entry is only used if it belongs to the same PersonId. - Test a Simulation not starting at 0 and a Person added during the Simulation. --- cpp/models/abm/common_abm_loggers.h | 48 +++++++++++++++++------------ cpp/tests/test_abm_simulation.cpp | 15 +++++++-- 2 files changed, 40 insertions(+), 23 deletions(-) diff --git a/cpp/models/abm/common_abm_loggers.h b/cpp/models/abm/common_abm_loggers.h index 86c0cc820e..a97c8a681f 100644 --- a/cpp/models/abm/common_abm_loggers.h +++ b/cpp/models/abm/common_abm_loggers.h @@ -224,18 +224,20 @@ struct LogInfectionStatePerAgeGroup : mio::LogAlways { /** * @brief Logger to log the TimeSeries of new #Infection%s per LocationType and AgeGroup. -* Unlike the other loggers, this one keeps a state: the time and the LocationType of each Person at its previous log. -* Hence, it has to log after every time step, which is the case when used in a History passed to Simulation::advance. +* Unlike the other loggers, this one keeps a state: the #InfectionState and the LocationType of each Person at its +* previous log. Hence, it has to log after every time step, which is the case when used in a History passed to +* Simulation::advance. As History stores its own Logger instances, use a new History for every Simulation. */ struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { using Type = std::pair; /** * @brief Log the TimeSeries of new #Infection%s per LocationType and AgeGroup. - * A Person is counted if it became #InfectionState::Exposed since the previous log. A Person is infected during the - * interaction at the beginning of a time step and may move to another Location afterwards. Hence, it is attributed - * to the LocationType it was at during the previous log, not to its current one. - * At the first log there is no previous log to compare against, so no Person is counted. This avoids miscounting - * the initially infected Person%s of the Model as new #Infection%s. + * A Person is counted if it was not #InfectionState::Exposed at the previous log, but is now. A Person is infected + * during the interaction at the beginning of a time step and may move to another Location afterwards. Hence, it is + * attributed to the LocationType it was at during the previous log, not to its current one. + * A Person without a previous log is not counted, but only remembered for the next log. This is the case for all + * Person%s at the first log, which avoids miscounting the initially infected Person%s of the Model as new + * #Infection%s, and for Person%s that were added to the Model since the previous log. * @param[in] sim The simulation of the abm. * @return A pair of the TimePoint and a vector counting the newly exposed Person%s per AgeGroup and LocationType, * indexed by `age_group * LocationType::Count + location_type`. @@ -244,32 +246,38 @@ struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways { { Eigen::VectorXd sum = Eigen::VectorXd::Zero( Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_model().parameters.get_num_groups())); - const auto curr_time = sim.get_time(); - const auto persons = sim.get_model().get_persons(); - const size_t num_persons = persons.size(); - // Only Persons that were already present at the previous log can be compared against it. - const size_t num_previous = m_prev_location_types.size(); - m_prev_location_types.resize(num_persons); + const auto curr_time = sim.get_time(); + const auto persons = sim.get_model().get_persons(); + // Entries beyond the previous number of Persons are new and have no previous log. + const size_t num_previous = m_previous.size(); + m_previous.resize(persons.size()); size_t i = 0; for (const Person& p : persons) { - if (i < num_previous && (p.get_infection_state(m_prev_time) != mio::abm::InfectionState::Exposed) && - (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) { + const auto curr_state = p.get_infection_state(curr_time); + // The entry only belongs to this Person if the Persons were not reordered since the previous log. + const bool has_previous_log = i < num_previous && m_previous[i].id == p.get_id(); + if (has_previous_log && m_previous[i].infection_state != mio::abm::InfectionState::Exposed && + curr_state == mio::abm::InfectionState::Exposed) { auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) + - ((uint32_t)m_prev_location_types[i]); + ((uint32_t)m_previous[i].location_type); sum[index] += 1; } // Remember where the Person is now, which is where the interaction of the next time step happens. - m_prev_location_types[i] = p.get_location_type(); + m_previous[i] = {p.get_id(), p.get_location_type(), curr_state}; ++i; } - m_prev_time = curr_time; return std::make_pair(curr_time, sum); } private: - mio::abm::TimePoint m_prev_time{0}; ///< Time of the previous log. - std::vector m_prev_location_types; ///< LocationType of each Person at the previous log. + /// @brief State of a Person at the previous log. + struct PreviousLog { + PersonId id; + mio::abm::LocationType location_type; + mio::abm::InfectionState infection_state; + }; + std::vector m_previous; ///< State of each Person at the previous log, in the order of the Persons. }; /** diff --git a/cpp/tests/test_abm_simulation.cpp b/cpp/tests/test_abm_simulation.cpp index d350dcbcf6..9d7e6b8fe3 100644 --- a/cpp/tests/test_abm_simulation.cpp +++ b/cpp/tests/test_abm_simulation.cpp @@ -180,7 +180,8 @@ TEST(TestSimulation, logInfectionPerLocationTypePerAgeGroup) auto model = mio::abm::Model(num_age_groups); auto home = model.add_location(mio::abm::LocationType::Home); auto work = model.add_location(mio::abm::LocationType::Work); - const auto t0 = mio::abm::TimePoint(0); + // the Simulation does not have to start at 0 + const auto t0 = mio::abm::TimePoint(0) + mio::abm::days(10); // the Persons only move when we tell them to model.use_mobility_rules(false); @@ -214,12 +215,19 @@ TEST(TestSimulation, logInfectionPerLocationTypePerAgeGroup) sim.advance(t0 + mio::abm::hours(1)); sim.get_model().change_location(moves_home, home); history.log(sim); - // a further time step without new Infections + // a further time step without new Infections, but a Person that is added to the Model already Exposed. It has no + // previous log, so it is not counted as a new Infection, but tracked from now on. sim.advance(t0 + mio::abm::hours(2)); + auto late_arrival = add_test_person(sim.get_model(), home, age_group_15_to_34, mio::abm::InfectionState::Exposed, + t0 + mio::abm::hours(2)); + sim.get_model().assign_location(late_arrival, home); + sim.get_model().assign_location(late_arrival, work); + history.log(sim); + sim.advance(t0 + mio::abm::hours(3)); history.log(sim); const auto& detailed = std::get<0>(history.get_log()); - ASSERT_EQ(detailed.get_num_time_points(), 3); + ASSERT_EQ(detailed.get_num_time_points(), 4); EXPECT_EQ(detailed.get_num_elements(), Eigen::Index((size_t)mio::abm::LocationType::Count * num_age_groups)); const auto idx = [](mio::AgeGroup age, mio::abm::LocationType type) { return Eigen::Index((size_t)mio::abm::LocationType::Count * (size_t)age.get() + (size_t)type); @@ -232,6 +240,7 @@ TEST(TestSimulation, logInfectionPerLocationTypePerAgeGroup) EXPECT_EQ(detailed.get_value(1)[idx(age_group_35_to_59, mio::abm::LocationType::Home)], 1.0); EXPECT_EQ(detailed.get_value(1)[idx(age_group_60_to_79, mio::abm::LocationType::Work)], 1.0); EXPECT_EQ(detailed.get_value(2).sum(), 0.0); + EXPECT_EQ(detailed.get_value(3).sum(), 0.0); } TEST(TestSimulation, logInfectionStatePerAgeGroup)