Commit 82b27067 authored by Pablo Iturrieta Rebolledo's avatar Pablo Iturrieta Rebolledo
Browse files

Modified T-test to include New Model object definition, and modified Results...

Modified T-test to include New Model object definition,  and modified Results Data structure to [ [ (Model_i, ref) for Model_i in Models ] for ref in Models ]
parent ae66fce9
import seaborn as sns
from csep.core import forecasts, catalogs
from csep.core import poisson_evaluations as poisson
from csep.utils import datasets, time_utils, plots, readers
import matplotlib.ticker as tick
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from csep.core import poisson_evaluations as poisson
from csep.utils import plots
import numpy as np
import copy
import scipy.stats
import pickle
import os
from codes import get_catalogs, get_models, filepaths
......@@ -23,14 +15,14 @@ def process_forecasts_batch(Model_batch, catalog, save_obj=False, load_obj=False
if not load_obj:
Results_T = []
Results_W = []
for f_i, f_j in itertools.product(Model_batch, Model_batch):
T_ij = poisson.paired_t_test(f_j, f_i, catalog)
Results_T.append(T_ij)
W_ij = poisson.w_test(f_j, f_i, catalog)
Results_W.append(W_ij)
for f_i in Model_batch:
T_ref_i = []
W_ref_i = []
for f_j in Model_batch:
T_ref_i.append(poisson.paired_t_test(f_j, f_i, catalog))
W_ref_i.append(poisson.w_test(f_j, f_i, catalog))
Results_T.append(T_ref_i)
Results_W.append(W_ref_i)
if save_obj:
with open(save_obj, 'wb') as obj:
pickle.dump((Results_T, Results_W), obj)
......@@ -42,22 +34,13 @@ def process_forecasts_batch(Model_batch, catalog, save_obj=False, load_obj=False
raise Exception('Results object not found')
return Results_T, Results_W
def merge_results_batchs(batch1, batch2):
L_merged = []
for j in batch1:
for i in batch2:
if i.sim_name == j.sim_name:
Forecast_5 = copy.deepcopy(i)
Forecast_10 = copy.deepcopy(j)
Forecast_5.sim_name = Forecast_5.sim_name + ' - 5yr'
Forecast_10.sim_name = '10yr'
L_merged.append(Forecast_5)
L_merged.append(Forecast_10)
return L_merged
def plot_single_results(Results, Forecasts, years, ref_model=0,folder=False):
plot_args = {'title': 'Paired T-test & W-test\n %i years experiment\nBase Model: %s' % (years, Forecasts[ref_model][0].name) ,
def plot_single_results(Results, years, ref_model=0,folder=False):
name_ref = Results[0][ref_model][0].sim_name[1]
plot_args = {'title': 'Paired T-test & W-test\n %i years experiment\nBase Model: %s' %
(years, name_ref) ,
'figsize':(6,8),
'xlabel': '',
'ylabel': 'Information Gain',
......@@ -67,8 +50,8 @@ def plot_single_results(Results, Forecasts, years, ref_model=0,folder=False):
'ylims': (-3,3)}
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))]
T_Results_per_model = Results[0][ref_model]
W_Results_per_model = Results[1][ref_model]
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),
......@@ -77,89 +60,50 @@ def plot_single_results(Results, Forecasts, years, ref_model=0,folder=False):
]
ax.legend(handles=legend_elements, loc='lower right')
if folder:
plt.savefig(os.path.join(folder, 'TW_%i_ref%s.png' % (years, Forecasts[ref_model].name)), dpi=300)
plt.savefig(os.path.join(folder, 'TW_%i_ref%s.png' % (years, name_ref)), dpi=300)
plt.show()
def plot_all_results(Results, Forecasts, years, order=False, savepath=False):
names = []
T = []
W = []
for t_i in Results[0]:
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).reshape((len(Forecasts),len(Forecasts))).T
score = np.sum(T, axis=1)/T.shape[0]
for w_i in Results[1]:
if w_i.quantile < 0.05:
W.append(1)
else:
W.append(0)
W = np.array(W).reshape((len(Forecasts),len(Forecasts))).T
def plot_all_results(Results, years, w_alpha=0.05, order=False, savepath=False):
names = [i[0].sim_name[1] for i in Results[0]]
T_value = np.array([[T_i.observed_statistic for T_i in Ref_model] for Ref_model in Results[0]]).T
W_value = np.array([[W_i.quantile for W_i in Ref_model] for Ref_model in Results[1]]).T
score = np.sum(T_value, axis=1)/T_value.shape[0]
if order:
arg_ind = np.flip(np.argsort(score))
data_t = T[arg_ind,:]
data_t = data_t[:, arg_ind]
data_w = W[arg_ind,:]
data_w = data_w[:, arg_ind]
else:
arg_ind = np.arange(len(score))
data_t = T_value[arg_ind, :][:, arg_ind] ## Flip rows/cols if ordered by value
data_w = W_value[arg_ind, :][:, arg_ind]
fig2, ax2 = plt.subplots(1, 1, figsize=(12,10))
fig, ax = plt.subplots(1, 1, figsize=(12,10))
cmap = sns.diverging_palette(220, 20, as_cmap=True)
sns.heatmap(data_t, vmin=-3, vmax=3, center= 0, cmap= cmap,
ax=ax2, cbar_kws={'pad': 0.01, 'label': 'Information Gain/EQ\nModel $i$ respect to Model $j$', 'shrink': 0.6,
ax=ax, cbar_kws={'pad': 0.01, 'label': 'Information Gain/EQ\nModel $i$ respect to Model $j$', 'shrink': 0.6,
'anchor': (0., 0.)})
ax2.set_yticklabels([names[i] for i in arg_ind], rotation='horizontal')
ax2.set_xticklabels([names[i] for i in arg_ind],rotation='vertical')
ax.set_yticklabels([names[i] for i in arg_ind], rotation='horizontal')
ax.set_xticklabels([names[i] for i in arg_ind],rotation='vertical')
for n, i in enumerate(data_w):
for m, j in enumerate(i):
if j ==1:
ax2.scatter(n + 0.5, m + 0.5,marker='o', s=2, color='black')
if j < w_alpha:
ax.scatter(n + 0.5, m + 0.5,marker='o', s=2, color='black')
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)
ax2.set_title('H-clustering of Paired T-test\n(%i years experiment)' % years, fontsize=16)
ax.set_title('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), dpi=300)
plt.show()
else:
fig, ax = plt.subplots(1, 1, figsize=(12, 10))
cmap = sns.diverging_palette(220, 20, as_cmap=True)
sns.heatmap(T, vmin=-4, vmax=2, center=0, cmap=cmap, ax=ax,
cbar_kws={'pad': 0.01, 'label': 'Information Gain/EQ\n$i$ respect to $j$', 'shrink': 0.6,
'anchor': (0., 0.)})
np.where(T == 0)
for n, i in enumerate(W):
for m, j in enumerate(i):
if j ==1:
ax.scatter(n + 0.5, m + 0.5,marker='o', s=2, color='black')
ax.set_xticklabels(names, rotation='vertical')
ax.set_yticklabels(names, rotation='horizontal')
ax.set_title('Heatmap Paired T-test \n(%i years)' % years, fontsize=16)
plt.tight_layout()
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.tight_layout()
plt.savefig(os.path.join(savepath), dpi=300)
plt.show()
......@@ -187,13 +131,13 @@ def run(use_saved=False):
for yr in years:
for model_i in range(len(models[yr])):
plot_single_results(Results[yr], models[yr], years=yr, ref_model=model_i,
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], models[yr], years=yr, order=typ,
plot_all_results(Results[yr], years=yr, order=typ,
savepath=filepaths.get_csep_fig_path(test, yr, type=typ*'order'))
if __name__ =='__main__':
use_saved=False
use_saved=True
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