Commit 51608a9f authored by Pablo Iturrieta Rebolledo's avatar Pablo Iturrieta Rebolledo
Browse files

Shortened name of forecast to HRSS

Modified CL_test to include new fonts
Modified t_test toinclude font.
Created function to plot PGS as t-style heatmap matrix
parent 91c4b83c
......@@ -49,12 +49,25 @@ def merge_results_batchs(batch1, batch2):
return CL_merged
def plot_results(Results, years, savepath=False):
plot_args = {'title': 'CL-test (%s years)' % years,
'figsize':(7,6),
plot_args = {'title': '$L_N-$test (%s years)' % years,
'figsize':(6,7),
'xlabel': 'Log-likelihood',
'title_fontsize' : 14,
'linewidth': 0.5,
'capsize': 0}
'capsize': 0,
'hbars':True}
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = 'Ubuntu'
plt.rcParams['font.monospace'] = 'Ubuntu Mono'
plt.rcParams['font.size'] = 10
plt.rcParams['axes.labelsize'] = 10
plt.rcParams['axes.labelweight'] = 'bold'
plt.rcParams['xtick.labelsize'] = 8
plt.rcParams['ytick.labelsize'] = 8
plt.rcParams['legend.fontsize'] = 10
plt.rcParams['figure.titlesize'] = 12
ax = plots.plot_poisson_consistency_test(Results, plot_args=plot_args, one_sided_lower=True)
for i,j in enumerate(Results):
......@@ -62,24 +75,37 @@ def plot_results(Results, years, savepath=False):
ax.xaxis.set_major_locator(tick.MaxNLocator(integer=True))
legend_elements = [Line2D([0], [0], marker='o', lw=0,label='Sim. expected value',
color='k', markersize=3),
color='k', markersize=2),
Line2D([0], [0], color='k', lw=1, label='Sim. 95% conf.'),
Line2D([0], [1], color='k', ls='o',lw=0, markersize=4, label='Obs. Likelihood')]
Line2D([0], [1], color='green', marker='o',lw=0, markersize=4, label='Obs. Likelihood (passes)'),
Line2D([0], [1], color='red', marker='o',lw=0, markersize=4, label='Obs. Likelihood (fails)')]
ax.legend(handles=legend_elements, loc='lower left')
ax.legend(handles=legend_elements, loc='lower left', fontsize=9)
if savepath:
plt.savefig(os.path.join(savepath), dpi=300)
plt.show()
def plot_merged_batchs(batch, savepath):
plot_args = {'title': 'CL-test (5-10 years)',
plot_args = {'title': '$L_N-$test',
'figsize':(6,7),
'xlabel': 'Log-likelihood',
'title_fontsize' : 14,
'linewidth': 0.5,
'capsize': 0,
'hbars':False}
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = 'Ubuntu'
plt.rcParams['font.monospace'] = 'Ubuntu Mono'
plt.rcParams['font.size'] = 10
plt.rcParams['axes.labelsize'] = 10
plt.rcParams['axes.labelweight'] = 'bold'
plt.rcParams['xtick.labelsize'] = 8
plt.rcParams['ytick.labelsize'] = 8
plt.rcParams['legend.fontsize'] = 10
plt.rcParams['figure.titlesize'] = 12
ax = plots.plot_poisson_consistency_test(batch, plot_args=plot_args, one_sided_lower=True)
for i,j in enumerate(batch):
plt.plot(np.mean(j.test_distribution),len(batch)-i-1, 'ko', markersize=2)
......@@ -87,12 +113,12 @@ def plot_merged_batchs(batch, savepath):
yTickPos = yTickPos[:-1]
ax.barh(yTickPos[::2] + 0.5, 1000*np.array([45] * len(yTickPos[::2])), left=-1000,
height=2 * (yTickPos[1] - yTickPos[0]), color=['w', 'gray'], alpha=0.3)
legend_elements = [Line2D([0], [0], marker='o', lw=0,label='Sim. expected value',
color='k', markersize=3),
legend_elements = [Line2D([0], [0], marker='o', lw=0,label='Sim. expected value',
color='k', markersize=2),
Line2D([0], [0], color='k', lw=1, label='Sim. 95% conf.'),
Line2D([0], [1], color='k', ls='o',lw=0, markersize=4, label='Obs. Likelihood')]
ax.legend(handles=legend_elements, loc='lower left')
Line2D([0], [1], color='green', marker='o',lw=0, markersize=4, label='Obs. Likelihood (passes)'),
Line2D([0], [1], color='red', marker='o',lw=0, markersize=4, label='Obs. Likelihood (fails)')]
ax.legend(handles=legend_elements, loc='lower left', fontsize=9)
if savepath:
plt.savefig(os.path.join(savepath), dpi=300)
......
......@@ -7,6 +7,46 @@ import numpy as np
import scipy.stats
from codes import get_catalogs, get_models, paths
import itertools
import seaborn as sns
def pgs_i(Forecasts, ref_forecast, Catalog, bin='spatial_magnitude'):
k = len(Forecasts)
if bin == 'spatial_magnitude':
catalog_counts = Catalog.spatial_magnitude_counts()
elif bin == 'spatial':
catalog_counts = Catalog.spatial_counts().reshape(-1,1)
probabilities_forecasts = []
for f_i in Forecasts:
data = np.abs(f_i.data)
if bin == 'spatial':
data = np.sum(data, axis=1).reshape(-1, 1)
prob_i = scipy.stats.poisson.pmf(catalog_counts, data)
probabilities_forecasts.append(prob_i)
ref = probabilities_forecasts[ref_forecast]
payout = np.zeros((k, catalog_counts.shape[0], catalog_counts.shape[1]))
for i in range(payout.shape[0]):
payout[i] = -1 + 2*probabilities_forecasts[i] / (probabilities_forecasts[i]+ref)
results = []
print('Model', 'Net Return', 'Net Return in EQ bins')
for model, i in enumerate(payout):
net = np.sum(i)
result_i = csep.models.EvaluationResult(
name='Parimutuel Gambling Score',
test_distribution=net,
sim_name=(Forecasts[model].name, Forecasts[ref_forecast].name))
results.append(result_i)
print(result_i.sim_name,net)
return results
def pgs(Forecasts, Catalog, ref='avg', bin='spatial_magnitude', save_obj=False, load_obj=False):
......@@ -156,5 +196,78 @@ def run(use_saved=False):
plot_merged_batchs(Results[(10, typ)], Results[(5, typ)], typ,
savepath=paths.get_csep_fig_path(test, '%i-%i' % (5, 10), typ))
def plot_all_results(Results, years, range_ = (-10, 10), order=False, savepath=False):
names = [i[0].sim_name[1] for i in Results]
PGS_value = np.array([[PGS_i.test_distribution for PGS_i in Ref_model] for Ref_model in Results]).T
score = np.sum(PGS_value, axis=1)/PGS_value.shape[0]
if order:
arg_ind = np.flip(np.argsort(score))
else:
arg_ind = np.arange(len(score))
data_t = PGS_value[arg_ind, :][:, arg_ind] ## Flip rows/cols if ordered by value
fig, ax = plt.subplots(1, 1, figsize=(12,10))
cmap = sns.diverging_palette(275, 150, s=80, l=55, as_cmap=True)
sns.heatmap(data_t, vmin=range_[0], vmax=range_[1], center= 0, cmap= cmap,
ax=ax, cbar_kws={'pad': 0.01, 'shrink': 0.7,
'anchor': (0., 0.),'ticks':[-10,-5,0,5,10]})
ax.set_yticklabels([names[i] for i in arg_ind], rotation='horizontal')
ax.set_xticklabels([names[i] for i in arg_ind], rotation='vertical')
# plt.text(len(score) + 2.5, -0.2, '$G_t$', fontsize=20)
# for i, j in enumerate([score[i] for i in arg_ind]):
# plt.text(len(score) + 2.5 , i+0.7, '%.3f' % j, fontsize=14)
ax.set_title('Paired PGS')
plt.tight_layout()
if savepath:
plt.savefig(os.path.join(savepath), dpi=300, transparent=True)
plt.show()
def run_cross(yr=10,use_saved=False):
test = 'PGS'
models = {5 : get_models.five_years(),
10 : get_models.ten_years()}
n_models = len(models[5])
catalogs = {5 : get_catalogs.filter_5yr_08_2009(get_catalogs.taroni2008()),
10 : get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino())}
Results = []
for i in range(n_models):
Result_ij = pgs_i(models[yr], i, catalogs[yr], bin='spatial_magnitude')
Results.append(Result_ij)
# plot_pgs(Result, yr, typ, savepath=paths.get_csep_fig_path(test, yr, 'ref%s' % typ))
plot_all_results(Results, 10, range_=(-10, 10), savepath=paths.get_csep_fig_path('pgs', 10, 'all'), order=True)
return Results
# for typ in types:
# plot_merged_batchs(Results[(10, typ)], Results[(5, typ)], typ,
# savepath=paths.get_csep_fig_path(test, '%i-%i' % (5, 10), typ))
if __name__ == '__main__':
run()
a = run_cross(yr=10)
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = 'Ubuntu'
plt.rcParams['font.monospace'] = 'Ubuntu Mono'
plt.rcParams['font.size'] = 14
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['axes.labelweight'] = 'bold'
plt.rcParams['xtick.labelsize'] = 14
plt.rcParams['ytick.labelsize'] = 14
plt.rcParams['legend.fontsize'] = 14
plt.rcParams['figure.titlesize'] = 18
......@@ -85,8 +85,9 @@ def plot_all_results(Results, years, p=0.05, order=False, savepath=False):
cmap = sns.diverging_palette(220, 20, as_cmap=True)
sns.heatmap(data_t, vmin=-3, vmax=3, center= 0, cmap= cmap,
ax=ax, cbar_kws={'pad': 0.01, 'label': 'Information Gain/EQ\nModel $i$ respect to Model $j$', 'shrink': 0.6,
'anchor': (0., 0.)})
ax=ax, cbar_kws={'pad': 0.01, 'shrink': 0.7,
'anchor': (0., 0.)}),
ax.set_yticklabels([names[i] for i in arg_ind], rotation='horizontal')
ax.set_xticklabels([names[i] for i in arg_ind],rotation='vertical')
......@@ -94,23 +95,19 @@ def plot_all_results(Results, years, p=0.05, order=False, savepath=False):
for m, j in enumerate(i):
if j < p:
ax.scatter(n + 0.5, m + 0.5,marker='o', s=2, color='black')
# for n, i in enumerate(data_tq):
# for m, j in enumerate(i):
# if j > 0:
# ax.scatter(n + 0.5, m + 0.75,marker='o', s=4, color='green')
plt.text(len(score) + 2.5, -0.2, '$G_t$', fontsize=20)
for i, j in enumerate([score[i] for i in arg_ind]):
plt.text(len(score) + 2.5 , i+0.7, '%.3f' % j, fontsize=14)
# plt.text(len(score) + 2.5, -0.2, '$G_t$', fontsize=20)
# for i, j in enumerate([score[i] for i in arg_ind]):
# plt.text(len(score) + 2.5 , i+0.7, '%.3f' % j, fontsize=14)
ax.set_title('Paired T-test\n(%i years experiment)' % years, fontsize=16)
ax.set_title('Information Gain $T-$test')
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))
bbox_to_anchor=(0.79, 0.08, 0.2, 0.2), handletextpad=0)
plt.tight_layout()
plt.savefig(os.path.join(savepath), dpi=300)
plt.savefig(os.path.join(savepath), dpi=300, transparent=True)
plt.show()
def run(use_saved=False):
......@@ -135,15 +132,25 @@ def run(use_saved=False):
tw_figs_path = os.path.join(paths.csep_figs, 'tw_single')
os.makedirs(tw_figs_path, exist_ok=True)
for yr in years:
for model_i in range(len(models[yr])):
plot_single_results(Results[yr], years=yr, ref_model=model_i,
folder=tw_figs_path)
# for yr in years:
# for model_i in range(len(models[yr])):
# plot_single_results(Results[yr], years=yr, ref_model=model_i,
# folder=tw_figs_path)
for yr, typ in itertools.product(years, types):
plot_all_results(Results[yr], years=yr, order=typ,
savepath=paths.get_csep_fig_path(test, yr, type=typ * 'order'))
if __name__ =='__main__':
use_saved = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = 'Ubuntu'
plt.rcParams['font.monospace'] = 'Ubuntu Mono'
plt.rcParams['font.size'] = 14
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['axes.labelweight'] = 'bold'
plt.rcParams['xtick.labelsize'] = 14
plt.rcParams['ytick.labelsize'] = 14
plt.rcParams['legend.fontsize'] = 14
plt.rcParams['figure.titlesize'] = 18
run(use_saved)
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