Commit 4a21fd4a authored by Pablo Iturrieta Rebolledo's avatar Pablo Iturrieta Rebolledo
Browse files

edited T-test figures to include legend

created get_results_figs() code to process all csep results batch
parent 106e65a2
......@@ -132,7 +132,7 @@ def run(used_saved=False):
if __name__ == '__main__':
used_saved = True
used_saved = False
run(used_saved)
......
......@@ -71,6 +71,12 @@ def plot_single_results(Results, Forecasts, years, ref_model=0,savepath=False):
T_Results_per_model = Results[0][ref_model*len(Forecasts):(ref_model+1)*(len(Forecasts))]
W_Results_per_model = Results[1][ref_model * len(Forecasts):(ref_model + 1) * (len(Forecasts))]
ax = plots.plot_comparison_test(T_Results_per_model, W_Results_per_model, plot_args=plot_args)
legend_elements = [Line2D([0], [0], marker='o', lw=0,label='W-test significant',
markerfacecolor='black', markeredgecolor='black', markersize=5),
Line2D([0], [0], marker='o', lw=0, label='W-test not significant',
markerfacecolor='white', markeredgecolor='black', markersize=5)
]
ax.legend(handles=legend_elements, loc='lower right')
if savepath:
plt.savefig(os.path.join(savepath, 'TW_%i_ref%s.png' % (years, Forecasts[ref_model][0].name)), dpi=300)
......@@ -121,6 +127,12 @@ def plot_all_results(Results, Forecasts, years, order=False, p=0.05, savepath=Fa
plt.text(len(score) + 2.5 , i+0.7, '%.3f' % j, fontsize=14)
ax2.set_title('H-clustering of Paired T-test\n(%i years experiment)' % years, fontsize=16)
legend_elements = [Line2D([0], [0], marker='o', lw=0, label='w-significant',
markerfacecolor='black', markeredgecolor='black', markersize=1)]
fig2.legend(handles=legend_elements, loc='lower right',
bbox_to_anchor=(0.78, 0.08, 0.2, 0.2))
plt.tight_layout()
plt.savefig(os.path.join(savepath,'T_%syr_order.png' % years), dpi=300)
plt.show()
......@@ -145,7 +157,10 @@ def plot_all_results(Results, Forecasts, years, order=False, p=0.05, savepath=Fa
plt.text(len(score) + 2.5, -0.2, '$G_t$', fontsize=20)
for i, j in enumerate(score):
plt.text(len(score) + 2.5 , i+0.7, '%.3f' % j, fontsize=14)
legend_elements = [Line2D([0], [0], marker='o', lw=0, label='w-significant',
markerfacecolor='black', markeredgecolor='black', markersize=1)]
fig.legend(handles=legend_elements, loc='lower right',
bbox_to_anchor=(0.78, 0.08, 0.2, 0.2))
plt.savefig(os.path.join(savepath, 'T_%syr.png' % years), dpi=300)
plt.show()
......@@ -153,11 +168,11 @@ if __name__ =='__main__':
Models5 = get_models.get_models_5yr()
cat5 = get_catalogs.filter_5yr_08_2009(get_catalogs.taroni2008())
# Models10 = get_models.get_models_10yr()
# cat10 = get_catalogs.filter_10yr_01_2010(get_catalogs.buolletino())
Models10 = get_models.get_models_10yr()
cat10 = get_catalogs.filter_10yr_01_2010(get_catalogs.buolletino())
Results5 = process_forecasts_batch(Models5, cat5, save_obj=filepaths.test_path['TW_5'])
# Results10 = process_forecasts_batch(None, None, load_obj=filepaths.test_path['L_10'])
Results5 = process_forecasts_batch(Models5, cat5, load_obj=filepaths.test_path['TW_5'])
Results10 = process_forecasts_batch(Models10, cat10, load_obj=filepaths.test_path['TW_10'])
tw_figs_path = os.path.join(filepaths.csep_figs, 'tw_tests')
......@@ -165,318 +180,11 @@ if __name__ =='__main__':
for i in range(len(Models5)):
plot_single_results(Results5, Models5,ref_model=i, years=5, savepath=tw_figs_path)
a = plot_all_results(Results5, Models5, years=5, savepath=filepaths.csep_figs)
a = plot_all_results(Results5, Models5, years=5, order=True, savepath=filepaths.csep_figs)
for i in range(len(Models10)):
plot_single_results(Results10, Models5,ref_model=i, years=10, savepath=tw_figs_path)
a = plot_all_results(Results10, Models10, years=10, savepath=filepaths.csep_figs)
a = plot_all_results(Results10, Models10, years=10, order=True, savepath=filepaths.csep_figs)
# plot_results(Results10, years=10, savepath=filepaths.csep_figs)
# plot_merged_batchs(merge_results_batchs(Results10, Results5), savepath=filepaths.csep_figs)
# get_score(Results5, cat5)
# get_IG(Results5)
#
#
#
# process_cat = False
# process_results = False
# plot_single = False
# plot_multiple = True
#
# Dir = '/media/pablo/Pandapan/gfz/PhD/codes/'
# data_dir = os.path.join(Dir,'data')
# results_dir = os.path.join(Dir,'results/figs/results_ercmt')
#
# model_dir = os.path.join(data_dir, '10yr')
#
# # cat_dir = os.path.join(data_dir, 'horus.txt')
# # cat_obj = os.path.join(data_dir, 'cat_horus.obj.txt')
# cat_dir = os.path.join(data_dir, 'emrcmt.csv')
# cat_obj = os.path.join(data_dir, 'cat_ercmt.obj')
#
# # cat_dir = os.path.join(data_dir, 'horus.txt')
# # cat_obj = os.path.join(data_dir, 'cat_horus.obj.txt')
# res_obj = os.path.join(data_dir, 't_test.obj')
#
# Forecasts_5 = []
# Forecasts_10 = []
# T_10 = []
# W_10 =[]
# T_5 = []
# W_5 = []
#
# if process_cat:
# start_date = time_utils.strptime_to_utc_epoch('2010-01-01 00:00:00.00')
# end_date = time_utils.strptime_to_utc_epoch('2020-01-01 00:00:00.00')
# cat_data = readers.ingv_emrcmt(cat_dir)
# cat = csep.core.catalogs.CSEPCatalog()
# cat.catalog = cat_data
# cat.update_catalog_stats()
#
# cat.filter(f'origin_time >= {start_date}')
# cat.filter(f'origin_time < {end_date}')
# cat.filter('magnitude >= 4.95')
#
# with open(cat_obj, 'wb') as obj:
# pickle.dump(cat, obj)
# else:
# with open(cat_obj, 'rb') as obj:
# cat = pickle.load(obj)
#
# if process_results:
#
# ### 10 yr
# Forecasts = []
# for i, model_i in enumerate([i for i in os.listdir(model_dir) if i.endswith('obj')]):
# with open(os.path.join(model_dir, model_i), 'rb') as obj:
# Forecast, Metadata = pickle.load(obj)
# Forecast.name = Forecast.name.split(' - ')[0]
# Forecasts_10.append(Forecast)
# if i == 0:
# cat.filter_spatial(region=Forecast.region)
#
# for f_i, f_j in itertools.product(Forecasts_10, Forecasts_10):
#
# T_ij = poisson.paired_t_test(f_j, f_i, cat)
# T_10.append(T_ij)
# W_ij = poisson.w_test(f_j, f_i, cat)
# W_10.append(W_ij)
# ### 5 yr
# start_date = time_utils.strptime_to_utc_epoch('2010-01-01 00:00:00.00')
# end_date = time_utils.strptime_to_utc_epoch('2015-01-01 00:00:00.00')
# cat.filter(f'origin_time >= {start_date}')
# cat.filter(f'origin_time < {end_date}')
#
# ### 5 yr
# Forecasts_5 = []
# for i, model_i in enumerate([i for i in os.listdir(model_dir) if i.endswith('obj')]):
# with open(os.path.join(model_dir, model_i), 'rb') as obj:
# Forecast, Metadata = pickle.load(obj)
# Forecast.name = Forecast.name.split(' - ')[0]
# Forecasts_5.append(Forecast)
# if i == 0:
# cat.filter_spatial(region=Forecast.region)
#
# for f_i, f_j in itertools.product(Forecasts_5, Forecasts_5):
#
# T_ij = poisson.paired_t_test(f_j, f_i, cat)
# T_5.append(T_ij)
# W_ij = poisson.w_test(f_j, f_i, cat)
# W_5.append(W_ij)
# with open(res_obj, 'wb') as obj:
# pickle.dump((T_5, T_10, W_5, W_10, Forecasts_5, Forecasts_10), obj)
#
# else:
# with open(res_obj, 'rb') as obj:
# T_5, T_10, W_5, W_10, Forecasts_5, Forecasts_10 = pickle.load(obj)
#
#
#
# if plot_single:
#
#
# ### Select base model
# model_id = 0
#
# ### 10 years
# ax = plots.plot_comparison_test(T_10[model_id*len(Forecasts_10):(model_id+1)*(len(Forecasts_10))],
# plot_args={'title': 'Paired T-test & W-test\n 10 years experiment\nBase Model: %s' % Forecasts_10[model_id].name,
# 'xlabel': '',
# 'ylabel': 'Information Gain/EQ\n(respect to base model)'} )
#
# for n, ind in enumerate(np.arange(model_id*len(Forecasts_10),(model_id+1)*(len(Forecasts_10)))):
# w = '%.3f' % W_10[ind].quantile
# t = T_10[ind].test_distribution[1]
# if np.isinf(t) or np.isinf(-t) or np.isnan(t) or t==0:
# t = 0
# w = ''
#
# if np.isinf(t) or np.isinf(-t) or np.isnan(t) or t == 0:
# t = 0
# w = ''
#
# if w != '' and W_5[ind].quantile < 0.05:
# sig = '$\checkmark$'
# else:
# sig = '$\circ$'
#
# plt.text(n, t + 0.25, sig, fontsize=18, horizontalalignment='center', verticalalignment='center')
#
# ax.yaxis.grid()
# xTickPos, _ = plt.xticks()
# ax.yaxis.set_major_locator(tick.MaxNLocator(integer=True))
# ax.set_ylim([ax.get_ylim()[0]+3, ax.get_ylim()[1]+0.5])
# ax.bar(xTickPos, np.array([9999] * len(xTickPos)), bottom=-2000,
# width=(xTickPos[1] - xTickPos[0]), color=['gray', 'w'], alpha=0.2)
# ax.set_xlim([ax.get_xlim()[0] + 0.5, ax.get_xlim()[1] - 0.5])
# plt.savefig(os.path.join(results_dir,'t_%s_10yr.png' % Forecasts_10[model_id].name), dpi=300)
# plt.show()
#
# #
# # ### 5 years
# # ax = plots.plot_comparison_test(T_5[model_id*len(Forecasts_5):(model_id+1)*(len(Forecasts_5))],
# # plot_args={'title': 'Paired T-test & W-test\n 5 years experiment\nBase Model: %s' % Forecasts_5[model_id].name,
# # 'xlabel': '',
# # 'ylabel': 'Information Gain/EQ\n(respect to base model)'} )
# #
# # for n, ind in enumerate(np.arange(model_id*len(Forecasts_5),(model_id+1)*(len(Forecasts_5)))):
# #
# # w = '%.3f' % W_5[ind].quantile
# # t = T_5[ind].test_distribution[1]
# #
# # if np.isinf(t) or np.isinf(-t) or np.isnan(t) or t==0:
# # t = 0
# # w = ''
# #
# # if w != '' and W_5[ind].quantile < 0.05:
# # sig = '$\checkmark$'
# # else:
# # sig = '$\circ$'
# #
# # plt.text(n,t+0.25, sig, fontsize=18, horizontalalignment='center',verticalalignment='center')
# #
# #
# # ax.yaxis.grid()
# #
# # xTickPos, _ = plt.xticks()
# # ax.yaxis.set_major_locator(tick.MaxNLocator(integer=True))
# # ax.set_ylim([ax.get_ylim()[0]+5, ax.get_ylim()[1]+0.5])
# # ax.bar(xTickPos, np.array([9999] * len(xTickPos)), bottom=-2000,
# # width=(xTickPos[1] - xTickPos[0]), color=['gray', 'w'], alpha=0.2)
# # ax.set_xlim([ax.get_xlim()[0] + 0.5, ax.get_xlim()[1] - 0.5])
# # plt.savefig(os.path.join(results_dir,'t_%s_5yr.png' % Forecasts_10[model_id].name), dpi=300)
# # plt.show()
#
# if plot_multiple:
#
#
# T = []
# names = []
# for t_i in T_10:
# T.append(t_i.observed_statistic)
# if t_i.sim_name[0] == t_i.sim_name[1]:
# names.append(t_i.sim_name[0])
# T = np.array(T)
# T = T.reshape((len(Forecasts_10),len(Forecasts_10))).T
#
#
# T[T==np.inf] = np.nan
# T[T==-np.inf] = np.nan
# mask = np.isnan(T)
#
# plt.close('all')
#
# fig, ax = plt.subplots(1, 1, figsize=(12,10))
# # cbar_ax = fig.add_axes([.93, .15, .03, .5])
#
# g = sns.heatmap(T, vmin=-4, vmax=2, center= 0, mask=mask,cmap= 'RdGy_r', ax=ax,
# cbar_kws={'pad':0.01, 'label':'Information Gain/EQ\n$i$ respect to $j$', 'shrink':0.6, 'anchor':(0.,0.)})
# ind_zero = np.where(T==0)
# ax.scatter(ind_zero[0] + 0.5, ind_zero[1] + 0.5, s=1)
#
# ax.set_xticklabels(names, rotation='vertical')
# ax.set_yticklabels(names,rotation='horizontal')
# g.set_facecolor('xkcd:yellow')
# ax.set_title('Heatmap Paired T-test \n(10 years)', fontsize=16)
# plt.tight_layout()
# # plt.savefig(os.path.join(results_dir,'T_heatmap_10yr.png'), dpi=300)
# plt.show()
#
# score = np.nansum(T, axis=1)/T.shape[0]
# with open(os.path.join(results_dir, 'g_t.txt'), 'w') as file_:
# for i in zip(names, score):
# file_.write('%s\t%.3f\n' % i)
# T = np.delete(T, 16, 0)
# T = np.delete(T, 16, 1)
# names.pop(16)
# T = np.delete(T, 13, 0)
# T = np.delete(T, 13, 1)
# names.pop(13)
# g = sns.clustermap(T, method='median')
# data = g.data2d.to_numpy()
# idx = [str(i) for i in np.array(g.data2d.index)]
# labels = [names[int(i)] for i in idx]
#
# fig2, ax2 = plt.subplots(1, 1, figsize=(15,10))
# d = sns.heatmap(data, vmin=-1.5, vmax=1.5, center= 0, cmap= 'RdGy_r',
# ax=ax2, cbar_kws={'pad': 0.01, 'label': 'Information Gain/EQ\nModel $i$ respect to Model $j$', 'shrink': 0.6,
# 'anchor': (0., 0.)})
# ax2.set_yticklabels(labels, rotation='horizontal')
# ax2.set_xticklabels(labels,rotation='vertical')
# ind_zero = np.where(data==0)
#
# ax2.scatter(ind_zero[0] + 1, ind_zero[1] + 0.5, s=2)
#
#
# score =np.sum(data,axis=1)/data.shape[0]
# plt.text(len(score) + 2.5, -0.2, '$G_t$', fontsize=20)
# for i, j in enumerate(score):
# plt.text(len(score) + 2.5 , i+0.7, '%.3f' % j, fontsize=14)
#
# ax2.set_title('H-clustering of Paired T-test\n(10 years experiment)', fontsize=16)
# plt.tight_layout()
# plt.savefig(os.path.join(results_dir,'T_hc_10yr_score.png'), dpi=300)
# plt.show()
#
# ##### 5 years
#
# #
# # T = []
# #
# # for t_i in T_5:
# # T.append(t_i.observed_statistic)
# #
# # T = np.array(T)
# # T = T.reshape((len(Forecasts_5),len(Forecasts_5))).T
# #
# #
# # T[T==np.inf] = np.nan
# # T[T==-np.inf] = np.nan
# # mask = np.isnan(T)
# #
# # plt.close('all')
# #
# # fig, ax = plt.subplots(1, 1, figsize=(12,10))
# #
# # g = sns.heatmap(T, vmin=-4, vmax=2, center= 0, mask=mask,cmap= 'RdGy_r', ax=ax,
# # cbar_kws={'pad':0.01, 'label':'Information Gain/EQ\nModel $i$ respect to Model $j$', 'shrink':0.6, 'anchor':(0.,0.)})
# # ind_zero = np.where(T==0)
# # ax.scatter(ind_zero[0] + 0.5, ind_zero[1] + 0.5, s=1)
# # ax.set_ylabel('Model $i$')
# # ax.set_xlabel('Model $j$')
# # g.set_facecolor('xkcd:yellow')
# # ax.set_title('Heatmap $i-j$ Paired T-test \n(5 years experiment)', fontsize=16)
# # plt.tight_layout()
# # plt.savefig(os.path.join(results_dir,'T_heatmap_5yr.png'), dpi=300)
# # plt.show()
# #
# #
# # T = np.delete(T, 13, 0)
# # T = np.delete(T, 13, 1)
# # T = np.delete(T, 8, 0)
# # T = np.delete(T, 8, 1)
# # T = np.delete(T, 7, 0)
# # T = np.delete(T, 7, 1)
# # T = np.delete(T, 6, 0)
# # T = np.delete(T, 6, 1)
# #
# # # T[T<0] = 0.
# # g = sns.clustermap(T, method='median')
# # data = g.data2d.to_numpy()
# # idx = [str(i) for i in np.array(g.data2d.index)]
# #
# #
# # fig2, ax2 = plt.subplots(1, 1, figsize=(12,10))
# # d = sns.heatmap(data, vmin=-1.5, vmax=1.5, center= 0, cmap= 'RdGy_r',
# # ax=ax2, cbar_kws={'pad': 0.01, 'label': 'Information Gain/EQ\nModel $i$ respect to Model $j$', 'shrink': 0.6,
# # 'anchor': (0., 0.)})
# # ax2.set_ylabel('Model $i$')
# # ax2.set_xlabel('Model $j$')
# # ax2.set_yticklabels(idx, rotation='horizontal')
# # ax2.set_xticklabels(idx)
# # ind_zero = np.where(data==0)
# #
# # ax2.scatter(ind_zero[0] + 0.5, ind_zero[1] + 0.5, s=2)
# # ax2.set_title('H-clustering of $i-j$ Paired T-test\n(5 years)', fontsize=16)
# # plt.tight_layout()
# # plt.savefig(os.path.join(results_dir,'T_hc_5yr.png'), dpi=300)
# # plt.show()
# #
#
from codes import get_catalogs, get_models, filepaths
from csep_tests import n_test, m_test, s_test, l_test, cl_test, pgs
import pickle
import itertools
def get_results(years):
results = {key : pickle.load(open(item, 'rb')) for key,item in filepaths.test_path.items()}
print(filepaths.test_path)
get_results(5)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment