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

updated brier score and multi score figure

parent 9af0f80b
......@@ -21,9 +21,9 @@ plt.rcParams['legend.fontsize'] = 14
plt.rcParams['figure.titlesize'] = 18
def pgs_i(Forecasts, ref_forecast, Catalog, bin='spatial_magnitude'):
def pgs_i(Forecasts, ref_forecast, Catalog, bin='spatial_magnitude', alpha=0.05):
N = Catalog.get_number_of_events()
k = len(Forecasts)
if bin == 'spatial_magnitude':
catalog_counts = Catalog.spatial_magnitude_counts()
......@@ -38,25 +38,64 @@ def pgs_i(Forecasts, ref_forecast, Catalog, bin='spatial_magnitude'):
probabilities_forecasts.append(prob_i)
idx_eq = Catalog.get_spatial_idx
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)
s_idx = Catalog.get_spatial_idx()
m_idx = Catalog.get_mag_idx()
s_idx, m_idx = np.nonzero(catalog_counts)
print(payout[1].shape)
results = []
print('Model', 'Net Return', 'Net Return in EQ bins')
for model, i in enumerate(payout):
net = np.sum(i)
s_2 = 1/(N-1) *(np.sum(np.power(payout[ref_forecast][s_idx, m_idx] - i[s_idx, m_idx], 2))) - \
1 / (N**2 - N) * np.power(np.sum(payout[ref_forecast][s_idx, m_idx] - i[s_idx, m_idx]), 2)
# s_2 = 1/(N-1) *(np.sum(np.power(payout[ref_forecast] - i, 2))) - \
# 1 / (N**2 - N) * np.power(np.sum(payout[ref_forecast] - i), 2)
std = np.sqrt(s_2)
df = N - 1
t_critical = scipy.stats.t.ppf(1 - (alpha / 2), df)
pgs_lower = net_avg - (t_critical * std / np.sqrt(N))
pgs_upper = net_avg + (t_critical * std / np.sqrt(N))
median = Forecasts[ref_forecast].sum() - Forecasts[model].sum()
diff = payout[ref_forecast][s_idx, m_idx] - i[s_idx, m_idx]
d = diff - median
d = np.compress(np.not_equal(d, 0), d, axis=-1)
count = len(d)
r = scipy.stats.rankdata(abs(d))
r_plus = np.sum((d > 0) * r, axis=0)
r_minus = np.sum((d < 0) * r, axis=0)
t = min(r_plus, r_minus)
mn = count * (count + 1.) * 0.25
se = count * (count + 1.) * (2. * count + 1.)
replist, repnum = scipy.stats.find_repeats(r)
if repnum.size != 0:
se -= 0.5 * (repnum * (repnum * repnum - 1)).sum()
se = np.sqrt(se / 24)
z = (t - mn) / se
prob = 2. * scipy.stats.distributions.norm.sf(abs(z))
net_avg = np.sum(i )
result_i = csep.models.EvaluationResult(
name='Parimutuel Gambling Score',
test_distribution=net,
test_distribution=net_avg,
quantile=(pgs_lower, pgs_upper),
observed_statistic = prob,
sim_name=(Forecasts[model].name, Forecasts[ref_forecast].name))
results.append(result_i)
print(result_i.sim_name,net)
return results
......@@ -212,9 +251,12 @@ def plot_all_results(Results, years, range_ = (-10, 10), order=False, savepath=F
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
PGS_q = np.array([[(PGS_i.quantile[0] > 0 and PGS_i.quantile[1] > 0 ) or
(PGS_i.quantile[0] < 0 and PGS_i.quantile[1] < 0) for PGS_i in Ref_model] for Ref_model in Results]).T
PGS_wq = np.array([[PGS_i.observed_statistic for PGS_i in Ref_model] for Ref_model in Results]).T
print(PGS_q)
score = np.sum(PGS_value, axis=1)/PGS_value.shape[0]
score = np.sum(PGS_value > 0, axis=1)
if order:
arg_ind = np.flip(np.argsort(score))
......@@ -226,19 +268,23 @@ def plot_all_results(Results, years, range_ = (-10, 10), order=False, savepath=F
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=ax, cbar_kws={'pad': 0.01, 'shrink': 0.7, 'label':' Net Return',
'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')
for n, i in enumerate(PGS_q):
for m, j in enumerate(i):
if j > 0 and PGS_wq[n,m] < 0.05:
ax.scatter(n + 0.5, m + 0.5, marker='o', s=5, color='black')
# ax.scatter(n + 0.5, m + 0.5, marker='o', s=75, facecolor='None', edgecolor='black')
# for n, i in enumerate(PGS_wq):
# for m, j in enumerate(i):
# if j < 0.05:
# ax.scatter(n + 0.5, m + 0.5, marker='o', s=5, 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)
ax.set_title('Paired PGS')
plt.tight_layout()
if savepath:
......@@ -248,7 +294,7 @@ def plot_all_results(Results, years, range_ = (-10, 10), order=False, savepath=F
def run_cross(catalog_name='bollettino'):
test = 'PGS'
years = [5, 10]
years = [10]
models = {5 : get_models.five_years(),
10 : get_models.ten_years()}
n_models = len(models[5])
......@@ -272,7 +318,7 @@ def run_cross(catalog_name='bollettino'):
Results.append(Result_ij)
# plot_pgs(Result, yr, typ, savepath=paths.get_csep_fig_path(test, yr, 'ref%s' % typ))
plot_all_results(Results, yr, range_=(-10, 10), savepath=paths.get_csep_figpath('pgs', yr, catalog_name), order=True)
plot_all_results(Results, yr, range_=(-10, 10), savepath=paths.get_csep_figpath('pgs', yr, catalog_name, 'svg'), order=True)
return Results
# for typ in types:
# plot_merged_batchs(Results[(10, typ)], Results[(5, typ)], typ,
......@@ -281,7 +327,7 @@ def run_cross(catalog_name='bollettino'):
if __name__ == '__main__':
catalogs = ['bollettino', 'emrcmt', 'horus']
for cat in catalogs:
for cat in catalogs[:1]:
# run(cat)
run_cross(cat)
a = run_cross(cat)
......@@ -97,34 +97,39 @@ 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, 'shrink': 0.7,
ax=ax, cbar_kws={'pad': 0.01, 'shrink': 0.7, 'label':'Information Gain',
'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')
for n, i in enumerate(data_tq):
for m, j in enumerate(i):
if j > 0:
ax.scatter(n + 0.5, m + 0.5, marker='o', s=75, facecolor='None', edgecolor='black')
for n, i in enumerate(data_w):
for m, j in enumerate(i):
if j < p:
if j > 0 and data_w[n,m] < p :
# ax.scatter(n + 0.5, m + 0.5, marker='o', s=75, facecolor='None', edgecolor='black')
ax.scatter(n + 0.5, m + 0.5, marker='o', s=5, color='black')
# for n, i in enumerate(data_w):
# for m, j in enumerate(i):
# if j < p:
# ax.scatter(n + 0.5, m + 0.5, marker='o', s=5, 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)
ax.set_title('Information Gain $T-$test')
# ax.set_title('Information Gain $T-$test')
legend_elements = [Line2D([0], [0], marker='o', lw=0, label='$W$-significant',
markerfacecolor="black", markeredgecolor='black', markersize=4),
Line2D([0], [0], marker='o', lw=0, label='$T$-significant',
markerfacecolor="None", markeredgecolor='black', markersize=10)
# legend_elements = [Line2D([0], [0], marker='o', lw=0, label='$W$-significant',
# markerfacecolor="black", markeredgecolor='black', markersize=4),
# Line2D([0], [0], marker='o', lw=0, label='$T$-significant',
# markerfacecolor="None", markeredgecolor='black', markersize=10)
# ]
legend_elements = [Line2D([0], [0], marker='o', lw=0, label='$T$- and $W$-significant',
markerfacecolor="black", markeredgecolor='black', markersize=4)
]
fig.legend(handles=legend_elements, loc='lower right',
bbox_to_anchor=(0.79, 0.08, 0.2, 0.2), handletextpad=0)
bbox_to_anchor=(0.79, 0.0, 0.2, 0.2), handletextpad=0)
plt.tight_layout()
plt.savefig(os.path.join(savepath), dpi=300, transparent=True)
plt.show()
......@@ -165,7 +170,7 @@ def run(use_saved=False, catalog_name='bollettino'):
# 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_figpath(test, yr, type=catalog_name + typ * 'order'))
savepath=paths.get_csep_figpath(test, yr, type=catalog_name + typ * 'order', format='svg'))
if __name__ =='__main__':
......
......@@ -167,7 +167,7 @@ def get_csep_figpath(test, year, type=None, format='png'):
type = ''
return join(csep_figs, "%s_%s%s.%s" % (test, str(year), type, format))
def get_kripley_figpath(analysis, year, model):
def get_kripley_figpath(analysis, year, model, format='png'):
"""
Get the figure path of a particular test
:param test: Name of the test
......@@ -175,7 +175,7 @@ def get_kripley_figpath(analysis, year, model):
:param type: Particular type of test
:return:
"""
return join(kripley_figs, "%s_%s_%s.png" % (analysis, str(year), model))
return join(kripley_figs, "%s_%s_%s.%s" % (analysis, str(year), model, format))
def get_sentropy_figpath(analysis, year, model):
......
......@@ -23,6 +23,8 @@ import seaborn as sns
import time
import pandas as pd
##3 CHECK ORDER OF SIMULATIONS
## CHECK LENGTH OF VECTOR
def get_simulations(model, catalog, lons, lats):
......@@ -247,21 +249,21 @@ def plot_results(Results, cat='bollettino', alpha=0.05):
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(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()
# sns.set_style("darkgrid", {"axes.facecolor": ".9", 'font.family':'Ubuntu'})
# 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']
......@@ -274,38 +276,44 @@ def plot_results(Results, cat='bollettino', alpha=0.05):
(L_cat<L_down)*(L_cat-L_down)
L_over[np.isnan(L_over)] = 0.
L_diff[key] = L_over
for i in L_sims:
sns.lineplot(x=r_sims, y=i, lw=0.05, color='blue', ls='--')
sns.lineplot(x=r_cat, y=L_cat, color='r', label='Observed catalog')
sns.lineplot(x=r_sims, y=L_avg, label='Sim. average')
plt.fill_between(r_sims, L_down, L_up, 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) = \sqrt{\frac{\hat{K}(r)}{\pi}}$")
plt.legend(loc=4)
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, lw=0.05, color='black', ls='--')
# sns.lineplot(x=r_cat, y=L_cat, color='r', label='Observed catalog')
# sns.lineplot(x=r_sims, y=L_avg, label='Sim. average')
# plt.fill_between(r_sims, L_down, L_up, 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) = \sqrt{\frac{\hat{K}(r)}{\pi}}$")
# plt.legend(loc=4)
# 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()
pcf_sims = value['PCF_sims']
pcf_cat = value['pcf_cat']
pcf_sims = [np.diff(K_i)/np.diff(r_sims)/2/np.pi/r_sims[:-1] for K_i in K_sims] #todo
# pcf_cat = value['pcf_cat'] #todo
pcf_cat = np.diff(K_cat) / np.diff(r_cat) / 2 / np.pi / r_cat[1:] #todo
# pcf_cat = np.hstack((pcf_cat[0], pcf_cat)) #todo
r_cat = value['rpcf_cat']
pcfr_sims = value['Rpcf_sims'][0]
pcfr_sims = value['Rpcf_sims'][0] #todo
pcfr_sims = r_sims[:-1] #todo
pcf_avg = np.mean(pcf_sims, axis=0)
pcf_up = np.quantile(pcf_sims, 1-alpha/2, axis=0)
pcf_down = np.quantile(pcf_sims, alpha/2, axis=0)
print(len(pcf_cat))
pcf_over = (pcf_cat>pcf_up)*(pcf_cat-pcf_up) +\
(pcf_down<pcf_cat)*(pcf_cat<pcf_up)*0 +\
(pcf_cat<pcf_down)*(pcf_cat-pcf_down)
......@@ -315,7 +323,7 @@ def plot_results(Results, cat='bollettino', alpha=0.05):
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.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)
......@@ -323,7 +331,7 @@ def plot_results(Results, cat='bollettino', alpha=0.05):
plt.show()
b = plot_combined(L_diff, pcf_diff, r_sims, pcfr_sims)
b = plot_combined(L_diff, pcf_diff, r_sims, pcfr_sims, order=None)
# return b
plt.show()
......@@ -332,7 +340,7 @@ def plot_combined(A,B, x, xx, order='inc'):
sns.set_style("dark", {"axes.facecolor": ".9", 'font.family': 'Ubuntu'})
from matplotlib import cm
figsize = (10,6)
figsize = (10,5)
fig, ax = plt.subplots(figsize=figsize)
xlims = []
range_ = [np.min(np.array([i for i in A.values()]).ravel()),
......@@ -343,10 +351,13 @@ def plot_combined(A,B, x, xx, order='inc'):
if order is None:
iter_array = A.keys()
print('HI')
elif order == 'inc':
model_order = np.argsort([np.sum(i != 0) for i in A.values()])
iter_array = [list(A.keys())[i] for i in model_order]
range_ = [-np.max(np.array([i for i in A.values()])), np.max(np.array([i for i in A.values()]))]
print(range_)
for index, key in enumerate(iter_array):
vals = A[key]
vals[np.isnan(vals)] = 0
......@@ -359,24 +370,23 @@ def plot_combined(A,B, x, xx, order='inc'):
y = index -0.5 - new_vals
img_data = np.arange(0, 100, 0.01)
img_data = np.flip(np.linspace(range_[0], range_[1], 100))
img_data = img_data.reshape(img_data.size, 1).T
im = ax.imshow(img_data, aspect='auto', origin='lower', cmap=plt.cm.coolwarm_r,
extent=[xmin, xmax, ymin, ymax],
vmin=20, vmax=80, zorder=n)
im = ax.imshow(img_data, aspect='auto', origin='lower', cmap=plt.cm.coolwarm,
extent=[xmin, xmax, ymin, ymax], zorder=n, vmin=-250, vmax=250)
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
ax.fill_betweenx(x, index - 0.5 - new_vals, [index ]*len(vals),
ax.fill_betweenx(x, index - 0.5 - new_vals, [index]*len(vals),
facecolor='none', edgecolor='gray', alpha=0.8,
where = (vals!=0), interpolate=True, zorder=n)
plt.axvline(index, c='black',zorder=n, lw=0.7)
plt.grid()
plt.axvline(index, c='black',zorder=n, lw=0.9)
plt.grid()
cba = plt.colorbar(im, shrink=0.5, fraction=0.03, pad=0.03)
cba.ax.set_title('$L(r)$', pad=15)
# plot and clip the imag
......@@ -385,31 +395,31 @@ def plot_combined(A,B, x, xx, order='inc'):
n = 0
range_ = [np.sqrt(np.min(np.array([i for i in B.values()])).ravel()),
np.sqrt(np.max(np.array([i for i in B.values()]).ravel()))]
range_ = [0,
np.sqrt(np.abs(np.max(np.array([i for i in B.values()]).ravel())))]
for index, key in enumerate(iter_array):
vals = B[key]
vals[np.isnan(vals)] = 0
vals = np.sqrt(vals)*np.sign(vals)
vals = np.sqrt(np.abs(vals))
vals[np.isnan(vals)] = 0
vals[-1] = 0
xamp = 1.5
xamp = 1.2
new_vals = (vals) / (range_[1]) * xamp - 0.5
xmin, xmax = index - xamp, index
ymin, ymax = min(xx), max(xx)
# plt.axvline(index, c='white', zorder=n + 1)
y = index - 0.5 - new_vals
img_data = np.arange(0, 100, 0.01)
img_data = np.linspace( range_[1], range_[0], 100)
img_data = np.flip(np.linspace(range_[0], range_[1], num=100))
img_data = img_data.reshape(img_data.size, 1).T
# plot and clip the imag
im = ax.imshow(img_data, aspect='auto', origin='lower', cmap=plt.cm.Greens,
extent=[xmin, xmax, ymin, ymax],
vmin=0, vmax=100, alpha=0.9, zorder=n)
extent=[xmin, xmax, ymin, ymax], alpha=0.9, zorder=n, vmax=25)
b = ax.fill(y, xx, facecolor='none', edgecolor='none', zorder=n)
for i in b:
im.set_clip_path(i)
......@@ -418,11 +428,12 @@ def plot_combined(A,B, x, xx, order='inc'):
facecolor='none', edgecolor='gray', alpha=0.8,
where=(vals != 0), interpolate=True, zorder=n)
cba = fig.colorbar(im, shrink=0.5, fraction=0.03, pad=0.03)
cba.ax.set_title('$\sqrt{g(r)}$', pad=15)
ax.set_xlim([-1, 19])
ax.set_ylim([-10, 550])
ax.set_ylim([-10, 800])
# plt.gca().invert_yaxis()
......@@ -430,9 +441,10 @@ def plot_combined(A,B, x, xx, order='inc'):
ax.set_xticks(numpy.arange(len(A)))
# ax.xaxis.tick_top()
ax.set_ylabel("$r\,[\mathrm{km}]$")
plt.title("Ripley's L and Pair Correlation Functions")
# plt.title("Ripley's L and Pair Correlation Functions")
plt.tight_layout()
plt.savefig(paths.get_kripley_figpath('Total', 10, key + '_' + cat))
print(paths.get_kripley_figpath('Total', 10, '_' + cat))
plt.savefig(paths.get_kripley_figpath('Total', 10, '_' + cat, format='png'), dpi=300)
plt.show()
return b
......@@ -515,9 +527,40 @@ def run():
if __name__ == "__main__":
Results = run()
for cat in ['bollettino', 'horus', 'emrcmt']:
# 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_)
#
a = plot_results(Results, cat)
## ########## Test for plotting simulations ######
# models = get_models.ten_years()
# model = models[13]
# lons = model.get_longitudes()
# lats = model.get_latitudes()
# catalog = get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino())
# sim_catalogs = [get_simulations(model, catalog, lons, lats)[0] for i in range(1)]
# ax = model.plot(plot_args={'clim': (-6, -2)})
# for i in sim_catalogs:
# ax.plot(i[:, 0], i[:, 1], '.', color='white', markersize=5)
# plt.show()
models = get_models.ten_years()
model = models[11]
lons = model.get_longitudes()
lats = model.get_latitudes()
catalog = get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino())
points = np.stack((catalog.get_longitudes(), catalog.get_latitudes())).T
r = np.linspace(0, 800, 200)
K = K_r_spatstat(points, model, r=r)
rk_cat = K[0]
K_cat = K[2]
PCF = pcf_spatstat(K)
plt.plot(PCF[2])
plt.plot(np.diff(K[2])/np.diff(K[0])/2/np.pi/K[0][:-1])
plt.show()
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
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