Commit bbd0edb8 authored by Pablo Iturrieta Rebolledo's avatar Pablo Iturrieta Rebolledo
Browse files

Readded n_test. Added result figures to repo. Modified forecast figures and...

Readded n_test. Added result figures to repo. Modified forecast figures and finalized K_function gen
parent 9d57b109
......@@ -5,6 +5,9 @@
<orderEntry type="jdk" jdkName="Python 3.8 (venv)" jdkType="Python SDK" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
<component name="PyDocumentationSettings">
<option name="renderExternalDocumentation" value="true" />
</component>
<component name="TestRunnerService">
<option name="PROJECT_TEST_RUNNER" value="pytest" />
</component>
......
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="PySciProjectComponent">
<option name="PY_SCI_VIEW" value="true" />
<option name="PY_SCI_VIEW_SUGGESTED" value="true" />
</component>
</project>
\ No newline at end of file
......@@ -151,7 +151,7 @@ def run(use_saved=False, catalog_name='bollettino'):
else:
Result = process_forecasts_batch(models[yr], catalogs[yr], save_obj=result_path)
Results[yr] = Result
plot_results(Result, years=yr, savepath=paths.get_csep_figpath(test, yr, type))
plot_results(Result, years=yr, savepath=paths.get_csep_figpath(test, yr, catalog_name))
plot_merged_batchs(merge_results_batchs(Results[years[0]], Results[years[1]]),
savepath=paths.get_csep_figpath(test, "%i-%i" % tuple(years), catalog_name))
......
......@@ -53,7 +53,8 @@ def plot_results(Results, years, savepath=False):
'figsize':(7,6),
'xlabel': 'Log-likelihood (Spatial)',
'linewidth': 0.5,
'capsize': 0}
'capsize': 0,
'percentile': 95}
ax = plots.plot_poisson_consistency_test(Results, plot_args=plot_args, one_sided_lower=True)
......@@ -135,10 +136,10 @@ def run(use_saved=False, catalog_name='bollettino'):
else:
Result = process_forecasts_batch(models[yr], catalogs[yr], save_obj=result_path)
Results[yr] = Result
plot_results(Result, years=yr, savepath=paths.get_csep_figpath(test, yr, catalog_name))
plot_results(Result, years=yr, savepath=paths.get_csep_figpath(test, yr, catalog_name+'2'))
plot_merged_batchs(merge_results_batchs(Results[years[0]], Results[years[1]]),
savepath=paths.get_csep_figpath(test, "%i-%i" % tuple(years), catalog_name))
# plot_merged_batchs(merge_results_batchs(Results[years[0]], Results[years[1]]),
# savepath=paths.get_csep_figpath(test, "%i-%i22" % tuple(years), catalog_name))
print('Scores 5yr')
get_score(Results[5])
print('Scores 10yr')
......@@ -146,8 +147,8 @@ def run(use_saved=False, catalog_name='bollettino'):
if __name__ == '__main__':
use_saved = False
catalogs = ['bollettino', 'emrcmt', 'horus']
use_saved = True
catalogs = ['bollettino', 'emrcmt', 'horus'][:1]
for cat in catalogs:
run(use_saved, cat)
......
......@@ -51,24 +51,27 @@ 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),
'figsize':(8,6),
'xlabel': '',
'ylabel': 'Information Gain',
'ylim': (-5, 5),
'linewidth': 1.5,
'capsize': 0,
'markersize': 6,
'ylims': (-4,4)}
'percentile': 95}
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',
legend_elements = [Line2D([0], [0], marker='o', lw=0, label='Sign. $(p<0.05)$',
markerfacecolor='black', markeredgecolor='black', markersize=5),
Line2D([0], [0], marker='o', lw=0, label='W-test not significant',
Line2D([0], [0], marker='o', lw=0, label='Not sign. $(p\geq0.05)$',
markerfacecolor='white', markeredgecolor='black', markersize=5)
]
ax.legend(handles=legend_elements, loc='lower right')
ax.legend(handles=legend_elements, title="W-test", loc='upper right',
handletextpad=0.5, fontsize=8, title_fontsize=10, framealpha=0.4)
if folder:
plt.savefig(os.path.join(folder, 'TW_%i_ref%s.png' % (years, name_ref)), dpi=300)
......@@ -89,7 +92,6 @@ def plot_all_results(Results, years, p=0.05, order=False, savepath=False):
data_t = T_value[arg_ind, :][:, arg_ind] ## Flip rows/cols if ordered by value
data_w = W_quantile[arg_ind, :][:, arg_ind]
print(T_quantile)
data_tq = T_quantile[arg_ind, :][:, arg_ind]
fig, ax = plt.subplots(1, 1, figsize=(12,10))
......
......@@ -3,7 +3,7 @@ import pickle
import os
import cartopy
import paths
import get_models
import get_models, get_catalogs
import numpy as np
# forecasts = get_models.five_years()
......
......@@ -229,42 +229,39 @@ def K_r_spatstat(points, model, normpower=2,plot=False, r=None):
return K
def plot_results(Results, alpha=0.05):
def plot_results(Results, cat='bollettino', alpha=0.05):
n = 0
L_diff = dict.fromkeys(Results.keys())
pcf_diff = dict.fromkeys(Results.keys())
for key, value in Results.items():
# n+=1
# if n == 2:
# break
K_sims = value['K_sims'][2]
K_sims = value['K_sims']
K_cat = value['K_cat']
r_cat = value['rk_cat']
r_sims = value['Rk_sims'][0]
# K_avg = np.nanmean(K_sims, axis=0)
# K_up = np.nanquantile(K_sims, 1-alpha/2, axis=0)
# K_down = np.nanquantile(K_sims, alpha/2, axis=0)
# print(K_avg)
# # sns.set_theme(font='sans-serif')
K_avg = np.nanmean(K_sims, axis=0)
K_up = np.nanquantile(K_sims, 1-alpha/2, axis=0)
K_down = np.nanquantile(K_sims, alpha/2, axis=0)
sns.set_style("darkgrid", {"axes.facecolor": ".9", 'font.family':'Ubuntu'})
# for i in K_sims:
# sns.lineplot(r_sims, i, lw=0.1, color='black')
# sns.lineplot(r_cat, K_cat, color='r', label='Observed catalog')
# sns.lineplot(r_sims, K_avg, label='Sim. average')
#
# # sns.lin
# plt.fill_between(r_sims, K_down, K_up, color='gray', alpha=0.2, label=r'Sim. envelope ($\alpha=%.2f$)' % (1-2*alpha))
# plt.title("Model - %s" % key, fontsize=16)
# plt.xlabel(r"$r~~\mathrm{[km]}$")
# plt.ylabel(r"$\hat{K}(r)$")
# plt.legend(loc=2)
#
# plt.savefig(paths.get_ripleyk_figpath('K', 10, key))
# plt.show()
for i in K_sims:
sns.lineplot(x=r_sims, y=i, lw=0.05, color='black')
sns.lineplot(x=r_cat, y=K_cat, color='r', label='Observed catalog')
sns.lineplot(x=r_sims, y=K_avg, label='Sim. average')
plt.fill_between(r_sims, K_down, K_up, color='gray', alpha=0.2, label=r'Sim. envelope ($\alpha=%.2f$)' % (1-2*alpha))
plt.title("Model: %s - Catalog: %s " % (key, cat), fontsize=16)
plt.xlabel(r"$r~~\mathrm{[km]}$")
plt.ylabel(r"$\hat{K}(r)$")
plt.legend(loc=2)
plt.savefig(paths.get_kripley_figpath('K', 10, key + '_' + cat )))
plt.show()
L_sims = value['L_sims']
L_cat = value['L_cat']
......@@ -286,21 +283,21 @@ def plot_results(Results, alpha=0.05):
plt.xlabel(r"$r~~\mathrm{[km]}$")
plt.ylabel(r"$L(r) = \sqrt{\frac{\hat{K}(r)}{\pi}}$")
plt.legend(loc=4)
plt.savefig(paths.get_kripley_figpath('L', 10, key))
plt.savefig(paths.get_kripley_figpath('L', 10, key + '_' + cat ))
plt.show()
for i in L_sims:
sns.lineplot(x=r_sims, y=i-r_sims, lw=0.05, color='blue', ls='--')
sns.lineplot(x=r_cat, y=L_cat-r_sims, color='r', label='Observed catalog')
sns.lineplot(x=r_sims, y=L_avg-r_sims, label='Sim. average')
plt.fill_between(r_sims, L_down-r_sims, L_up-r_sims, color='gray', alpha=0.4, label=r'Sim. envelope ($\alpha=%.2f$)' % (1-alpha))
plt.title("Model - %s" % key, fontsize=16)
plt.xlabel(r"$r~~\mathrm{[km]}$")
plt.ylabel(r"$L(r)-r = \sqrt{\frac{\hat{K}(r)}{\pi}}-r$")
plt.legend(loc=4)
plt.savefig(paths.get_kripley_figpath('Lr', 10, key + '_' + cat))
plt.show()
# for i in L_sims:
# sns.lineplot(x=r_sims, y=i-r_sims, lw=0.05, color='blue', ls='--')
# sns.lineplot(x=r_cat, y=L_cat-r_sims, color='r', label='Observed catalog')
# sns.lineplot(x=r_sims, y=L_avg-r_sims, label='Sim. average')
# plt.fill_between(r_sims, L_down-r_sims, L_up-r_sims, color='gray', alpha=0.4, label=r'Sim. envelope ($\alpha=%.2f$)' % (1-alpha))
# plt.title("Model - %s" % key, fontsize=16)
# plt.xlabel(r"$r~~\mathrm{[km]}$")
# plt.ylabel(r"$L(r)-r = \sqrt{\frac{\hat{K}(r)}{\pi}}-r$")
# plt.legend(loc=4)
# plt.savefig(paths.get_kripley_figpath('Lr', 10, key))
# plt.show()
#
pcf_sims = value['PCF_sims']
pcf_cat = value['pcf_cat']
r_cat = value['rpcf_cat']
......@@ -314,19 +311,19 @@ def plot_results(Results, alpha=0.05):
(pcf_cat<pcf_down)*(pcf_cat-pcf_down)
pcf_over[np.isnan(pcf_over)] = 0.
pcf_diff[key] = pcf_over
# sns.lineplot(x=pcfr_sims, y=pcf_cat, color='r', label='Observed catalog')
# sns.lineplot(x=pcfr_sims, y=pcf_avg, label='Sim. average')
# plt.fill_between(pcfr_sims, pcf_down, pcf_up, color='gray', alpha=0.2, label=r'Sim. envelope ($\alpha=%.2f$)' % (1-alpha))
# plt.title("Model - %s" % key, fontsize=16)
# plt.ylim([-20, 20])
# plt.xlabel(r"$r~~\mathrm{[km]}$")
# plt.ylabel(r"$g(r) = \frac{1}{2\pi r}\,\frac{dK(r)}{dr}$")
# plt.legend(loc=4)
# plt.savefig(paths.get_ripleyk_figpath('pcf', 10, key))
# plt.show()
# env_diff['r'] = r_sims
# b = plot_combined(L_diff, pcf_diff, r_sims, pcfr_sims)
sns.lineplot(x=pcfr_sims, y=pcf_cat, color='r', label='Observed catalog')
sns.lineplot(x=pcfr_sims, y=pcf_avg, label='Sim. average')
plt.fill_between(pcfr_sims, pcf_down, pcf_up, color='gray', alpha=0.2, label=r'Sim. envelope ($\alpha=%.2f$)' % (1-alpha))
plt.title("Model - %s" % key, fontsize=16)
plt.ylim([-20, 20])
plt.xlabel(r"$r~~\mathrm{[km]}$")
plt.ylabel(r"$g(r) = \frac{1}{2\pi r}\,\frac{dK(r)}{dr}$")
plt.legend(loc=4)
plt.savefig(paths.get_kripley_figpath('pcf', 10, key + '_' + cat))
plt.show()
b = plot_combined(L_diff, pcf_diff, r_sims, pcfr_sims)
# return b
plt.show()
......@@ -368,7 +365,7 @@ def plot_combined(A,B, x, xx, order='inc'):
extent=[xmin, xmax, ymin, ymax],
vmin=20, vmax=80, zorder=n)
b = ax.fill(y, x, facecolor='none', edgecolor='none', zorder=n)
b = ax.fill(np.append(y, index), np.append(x, x[-1]), facecolor='none', edgecolor='none', zorder=n)
for i in b:
im.set_clip_path(i)
n += 1
......@@ -435,7 +432,9 @@ def plot_combined(A,B, x, xx, order='inc'):
ax.set_ylabel("$r\,[\mathrm{km}]$")
plt.title("Ripley's L and Pair Correlation Functions")
plt.tight_layout()
plt.savefig(paths.get_kripley_figpath('Total', 10, key + '_' + cat))
plt.show()
return b
def run():
......@@ -455,7 +454,7 @@ def run():
print(model.name)
points = np.stack((catalog.get_longitudes(), catalog.get_latitudes())).T
r = np.linspace(0,800,200)
r = np.linspace(0, 800, 200)
K = K_r_spatstat(points, model, r=r)
rk_cat = K[0]
K_cat = K[2]
......@@ -516,9 +515,9 @@ def run():
if __name__ == "__main__":
Results = run()
# for cat in ['bollettino', 'horus', 'emrcmt'][:1]:
# with open(paths.get_kripley_result_path('K_%s' % cat, 10), 'rb') as file_:
# Results = pickle.load(file_)
# Results = run()
for cat in ['bollettino', 'horus', 'emrcmt']:
with open(paths.get_kripley_result_path('K_%s' % cat, 10), 'rb') as file_:
Results = pickle.load(file_)
# a = plot_results(Results)
a = plot_results(Results, cat)
results/figures/csep/T_10_emrcmt.png

241 KB | W: | H:

results/figures/csep/T_10_emrcmt.png

255 KB | W: | H:

results/figures/csep/T_10_emrcmt.png
results/figures/csep/T_10_emrcmt.png
results/figures/csep/T_10_emrcmt.png
results/figures/csep/T_10_emrcmt.png
  • 2-up
  • Swipe
  • Onion skin
results/figures/csep/T_10_horus.png

240 KB | W: | H:

results/figures/csep/T_10_horus.png

256 KB | W: | H:

results/figures/csep/T_10_horus.png
results/figures/csep/T_10_horus.png
results/figures/csep/T_10_horus.png
results/figures/csep/T_10_horus.png
  • 2-up
  • Swipe
  • Onion skin
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