Commit 64b18876 authored by Pablo Iturrieta Rebolledo's avatar Pablo Iturrieta Rebolledo
Browse files

Created multi plot for visualizing all model's results

parent ba725e2f
......@@ -2,7 +2,7 @@
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="jdk" jdkName="Python 3.7 (venv)" jdkType="Python SDK" />
<orderEntry type="jdk" jdkName="Python 3.8 (venv)" jdkType="Python SDK" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
<component name="TestRunnerService">
......
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.7 (venv)" project-jdk-type="Python SDK" />
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8 (venv)" project-jdk-type="Python SDK" />
</project>
\ No newline at end of file
......@@ -19,7 +19,7 @@ from rpy2.robjects import globalenv
import cartopy
import seaborn as sns
import time
import pandas as pd
def get_simulations(model, catalog, lons, lats):
......@@ -229,6 +229,8 @@ def K_r_spatstat(points, model, normpower=2,plot=False, r=None):
def plot_results(Results, 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:
......@@ -261,13 +263,17 @@ def plot_results(Results, alpha=0.05):
# plt.savefig(paths.get_ripleyk_figpath('K', 10, key))
# plt.show()
# L_sims = value['L_sims']
# L_cat = value['L_cat']
# L_avg = np.mean(L_sims, axis=0)
# L_up = np.quantile(L_sims, 1-alpha/2, axis=0)
# L_down = np.quantile(L_sims, alpha/2, axis=0)
#
#
L_sims = value['L_sims']
L_cat = value['L_cat']
L_avg = np.mean(L_sims, axis=0)
L_up = np.quantile(L_sims, 1-alpha/2, axis=0)
L_down = np.quantile(L_sims, alpha/2, axis=0)
L_over = (L_cat>L_up)*(L_cat-L_up) +\
(L_down<L_cat)*(L_cat<L_up)*0 +\
(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')
......@@ -279,25 +285,155 @@ def plot_results(Results, alpha=0.05):
# plt.legend(loc=4)
# plt.savefig(paths.get_ripleyk_figpath('L', 10, key))
# 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_ripleyk_figpath('Lr', 10, key))
# plt.show()
#
pcf_sims = value['PCF_sims']
pcf_cat = value['pcf_cat']
r_cat = value['rpcf_cat']
r_sims = value['Rpcf_sims'][0]
pcfr_sims = value['Rpcf_sims'][0]
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)
sns.lineplot(x=r_cat, y=pcf_cat, color='r', label='Observed catalog')
sns.lineplot(x=r_sims, y=pcf_avg, label='Sim. average')
plt.fill_between(r_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()
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)
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)
return b
plt.show()
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)
fig, ax = plt.subplots(figsize=figsize)
xlims = []
range_ = [np.min(np.array([i for i in A.values()]).ravel()),
np.max(np.array([i for i in A.values()]).ravel())]
n = 0
if order is None:
iter_array = A.keys()
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]
for index, key in enumerate(iter_array):
vals = A[key]
vals[np.isnan(vals)] = 0
xamp = 1.5
new_vals = (vals)/(range_[1])*xamp - 0.5
# ax.fill([index]*len(new_vals), new_vals)
xmin, xmax = index - xamp, index + xamp
ymin, ymax = min(x), max(x)
y = index -0.5 - new_vals
img_data = np.arange(0, 100, 0.01)
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)
b = ax.fill(y, x, 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),
facecolor='none', edgecolor='gray', alpha=0.8,
where = (vals!=0), interpolate=True, zorder=n)
plt.axvline(index, c='white',zorder=n)
# plot and clip the imag
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()))]
for index, key in enumerate(iter_array):
vals = B[key]
vals[np.isnan(vals)] = 0
vals = np.sqrt(vals)*np.sign(vals)
vals[np.isnan(vals)] = 0
vals[-1] = 0
xamp = 1.5
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 = 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)
b = ax.fill(y, xx, facecolor='none', edgecolor='none', zorder=n)
for i in b:
im.set_clip_path(i)
n += 1
a = ax.fill_betweenx(xx, y, [index] * len(vals),
facecolor='none', edgecolor='gray', alpha=0.8,
where=(vals != 0), interpolate=True, zorder=n)
ax.set_xlim([-1, 19])
ax.set_ylim([-10, 550])
# plt.gca().invert_yaxis()
ax.set_xticklabels([key for key in iter_array], rotation=45, ha='right')
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.tight_layout()
plt.show()
return b
def run_parallel():
......@@ -443,9 +579,9 @@ def run():
if __name__ == "__main__":
Results = run_parallel()
# Results = run_parallel()
#
with open(paths.get_kripley_result_path('KRipley_bolletino', 10), 'rb') as file_:
Results = pickle.load(file_)
plot_results(Results)
a = plot_results(Results)
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