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

Added K ripley inhom

parent 3e3a395c
......@@ -5,7 +5,7 @@
<paths name="pciturri@rs1:22">
<serverdata>
<mappings>
<mapping deploy="/tmp/pycharm_project_983" local="$PROJECT_DIR$" />
<mapping deploy="/home/pciturri/pycharmprojects/italy-experiment-testing" local="$PROJECT_DIR$" web="/" />
</mappings>
</serverdata>
</paths>
......
......@@ -2,7 +2,7 @@
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="jdk" jdkName="Remote Python 2.7.17 (sftp://pciturri@rs1:22/usr/bin/python)" jdkType="Python SDK" />
<orderEntry type="jdk" jdkName="Python 3.7 (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="Remote Python 2.7.17 (sftp://pciturri@rs1:22/usr/bin/python)" project-jdk-type="Python SDK" />
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.7 (venv)" project-jdk-type="Python SDK" />
</project>
\ No newline at end of file
......@@ -167,10 +167,10 @@ def run(use_saved=False):
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_fig_path(test, yr, type))
plot_results(Result, years=yr, savepath=paths.get_csep_figpath(test, yr, type))
plot_merged_batchs(merge_results_batchs(Results[years[1]], Results[years[0]]),
savepath=paths.get_csep_fig_path(test, "%i-%i" % tuple(years), types[0]))
savepath=paths.get_csep_figpath(test, "%i-%i" % tuple(years), types[0]))
print('Scores 5yr')
get_score(Results[5])
print('Scores 10yr')
......
......@@ -143,10 +143,10 @@ def run(use_saved=False):
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_fig_path(test, yr, type))
plot_results(Result, years=yr, savepath=paths.get_csep_figpath(test, yr, type))
plot_merged_batchs(merge_results_batchs(Results[years[0]], Results[years[1]]),
savepath=paths.get_csep_fig_path(test, "%i-%i" % tuple(years), types[0]))
savepath=paths.get_csep_figpath(test, "%i-%i" % tuple(years), types[0]))
print('Scores 5yr')
get_score(Results[5], catalogs[5])
print('Scores 10yr')
......
......@@ -130,11 +130,11 @@ def run(use_saved=False):
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_fig_path(test, yr, type))
plot_results(Result, years=yr, savepath=paths.get_csep_figpath(test, yr, type))
plot_merged_batchs(merge_results_batchs(Results[years[0]], Results[years[1]]),
savepath=paths.get_csep_fig_path(test, "%i-%i" % tuple(years), types[0]))
savepath=paths.get_csep_figpath(test, "%i-%i" % tuple(years), types[0]))
print('Scores 5yr')
get_score(Results[5])
print('Scores 10yr')
......
......@@ -142,11 +142,11 @@ def run(use_saved=False):
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_fig_path(test, yr, type))
plot_results(Result, years=yr, savepath=paths.get_csep_figpath(test, yr, type))
plot_merged_batchs(merge_results_batchs(Results[years[0]], Results[years[1]]),
savepath=paths.get_csep_fig_path(test, "%i-%i" % tuple(years), types[0]))
savepath=paths.get_csep_figpath(test, "%i-%i" % tuple(years), types[0]))
print('Scores 5yr')
get_score(Results[5])
print('Scores 10yr')
......
......@@ -190,11 +190,11 @@ def run(use_saved=False):
if typ == int:
typ = models[yr][typ][0].name
plot_pgs(Result, yr, typ, savepath=paths.get_csep_fig_path(test, yr, 'ref%s' % typ))
plot_pgs(Result, yr, typ, savepath=paths.get_csep_figpath(test, yr, 'ref%s' % typ))
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))
savepath=paths.get_csep_figpath(test, '%i-%i' % (5, 10), typ))
def plot_all_results(Results, years, range_ = (-10, 10), order=False, savepath=False):
......@@ -253,7 +253,7 @@ def run_cross(yr=10,use_saved=False):
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)
plot_all_results(Results, 10, range_=(-10, 10), savepath=paths.get_csep_figpath('pgs', 10, 'all'), order=True)
return Results
# for typ in types:
# plot_merged_batchs(Results[(10, typ)], Results[(5, typ)], typ,
......
......@@ -126,10 +126,10 @@ def run(use_saved=False):
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_fig_path(test, yr, type))
plot_results(Result, years=yr, savepath=paths.get_csep_figpath(test, yr, type))
plot_merged_batchs(merge_results_batchs(Results[years[0]], Results[years[1]]),
savepath=paths.get_csep_fig_path(test, "%i-%i" % tuple(years), types[0]))
savepath=paths.get_csep_figpath(test, "%i-%i" % tuple(years), types[0]))
print('Scores 5yr')
get_score(Results[5])
print('Scores 10yr')
......
......@@ -138,7 +138,7 @@ def run(use_saved=False):
# 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'))
savepath=paths.get_csep_figpath(test, yr, type=typ * 'order'))
if __name__ =='__main__':
use_saved = True
......
......@@ -4,48 +4,7 @@ import os
import cartopy
import paths
import get_models
import numpy as np
plot_properties = {'grid_labels': True,
'borders': True,
'feature_lw': 0.5,
'basemap': 'stock_img',
'cmap': 'rainbow',
'clim' : [-8, 0],
'projection': cartopy.crs.Mercator()}
# forecasts = get_models.five_years()
# figures_dir = paths.forecast_figures[5]
#
# for model_i in forecasts:
# plot_properties['filename'] = os.path.join(figures_dir, model_i.name)
# plot_properties['title'] = model_i.name + ' ' + model_i.author
# print(plot_properties['title'])
# model_i.plot(plot_args=plot_properties)
#
# forecasts = get_models.ten_years()
# figures_dir = paths.forecast_figures[10]
# for model_i in forecasts:
# plot_properties['filename'] = os.path.join(figures_dir, model_i.name)
# plot_properties['title'] = model_i.name + ' ' + model_i.author
# print(plot_properties['title'])
# model_i.plot(plot_args=plot_properties)
### Contours
plot_properties = {'grid_labels': True,
'borders': True,
'feature_lw': 0.5,
'basemap': 'stock_img',
'cmap': 'rainbow',
'clim' : [-8, 0],
'projection': cartopy.crs.Mercator()}
forecasts = get_models.ten_years()
figures_dir = paths.forecast_figures[10]
for model_i in forecasts:
plot_properties['filename'] = os.path.join(figures_dir, 'fig_bins/%s_8-0-15' % model_i.name)
plot_properties['title'] = model_i.name + ' ' + model_i.author
print(plot_properties['title'])
model_i.plot(plot_args=plot_properties)
......@@ -96,6 +96,7 @@ def forecast_xmlreader(filename_):
forecast.m_bin_width = m_bin_width
forecast.cell_dim = cell_dim
forecast.depth = depth
forecast.years = filename.split('.')[-3][:-2]
os.remove(filename_dat)
return forecast
......@@ -131,6 +132,7 @@ def ten_years():
if __name__ == '__main__':
process_models()
# path = paths.five_years_src
# model = forecast_xmlreader(path[0])
......@@ -132,7 +132,7 @@ def plot_results(years, labels, ll_type=None, n_eqk = 1,
### Information Gain
T = np.array([[T_i.observed_statistic for T_i in Ref_model] for
Ref_model in paths.get_csep_result('T', 5)[0]]).T
Ref_model in paths.get_csep_result('T', years)[0]]).T
if ref == 'avg':
t_score = np.sum(T, axis=1) / T.shape[0]
else:
......@@ -224,17 +224,17 @@ if __name__ == '__main__':
labels = ['Information Gain (t-test)', 'Net Return (PGS)',
'$\log(L_N)$ quantile', 'Poisson Consistency $(p=%.2f)$' % p]
plot_results(5, l_cl='CL', ll_type='quantile', labels=labels, ref='avg', lowcuts=[False, 0.1, False],
savepath=paths.get_csep_fig_path('Total', 5, 'avg_CLq'))
savepath=paths.get_csep_figpath('Total', 5, 'avg_CLq'))
plot_results(10, l_cl='CL', ll_type='quantile', labels=labels, ref='avg', lowcuts=[False, 0.1, False],
savepath=paths.get_csep_fig_path('Total', 10, 'avg_CLq'))
savepath=paths.get_csep_figpath('Total', 10, 'avg_CLq'))
labels = ['Information Gain (t-test)', 'Net Return (PGS)',
'$\log(L_N)$ - $\log(L_N)_{ref}$', 'Poisson Consistency $(p=%.2f)$' % p]
plot_results(5, ref=11, l_cl='CL', labels=labels, lowcuts=[False, -1, -15],
savepath=paths.get_csep_fig_path('Total', 5, 'ref11_CL'))
savepath=paths.get_csep_figpath('Total', 5, 'ref11_CL'))
plot_results(10, ref=11, l_cl='CL', labels=labels, lowcuts=[False, -2, -20],
savepath=paths.get_csep_fig_path('Total', 10, 'ref11_CL'))
savepath=paths.get_csep_figpath('Total', 10, 'ref11_CL'))
......@@ -242,16 +242,16 @@ if __name__ == '__main__':
labels = ['Information Gain (t-test)', 'Net Return (PGS)',
'$\log(L)$ quantile', 'Poisson Consistency $(p=%.2f)$' % p]
plot_results(5, l_cl='L', ll_type='quantile', labels=labels, ref='avg', lowcuts=[False, 0.1, False],
savepath=paths.get_csep_fig_path('Total', 5, 'avg_Lq'))
savepath=paths.get_csep_figpath('Total', 5, 'avg_Lq'))
plot_results(10, l_cl='L', ll_type='quantile', labels=labels, ref='avg', lowcuts=[False, -1, False],
savepath=paths.get_csep_fig_path('Total', 10, 'avg_Lq'))
savepath=paths.get_csep_figpath('Total', 10, 'avg_Lq'))
labels = ['Information Gain (t-test)', 'Net Return (PGS)',
'$\log(L)$ - $\log(L)_{ref}$', 'Poisson Consistency $(p=%.2f)$' % p]
plot_results(5, ref=11, l_cl='L', labels=labels, lowcuts=[False, -1, -15],
savepath=paths.get_csep_fig_path('Total', 5, 'ref11_L'))
savepath=paths.get_csep_figpath('Total', 5, 'ref11_L'))
plot_results(10, ref=11, l_cl='L', labels=labels, lowcuts=[False, -1, -15],
savepath=paths.get_csep_fig_path('Total', 10, 'ref11_L'))
savepath=paths.get_csep_figpath('Total', 10, 'ref11_L'))
......@@ -2,40 +2,44 @@ import csep
import os
import itertools
import pickle
from os.path import join, dirname
# Data years
years = [5, 10]
# Get absolute paths of Package
codes = os.path.dirname(os.path.realpath(__file__))
main = os.path.dirname(codes)
codes = dirname(os.path.realpath(__file__))
main = dirname(codes)
# Data paths
data = os.path.join(main, 'data')
data = join(main, 'data')
# Catalogs
catalogs = os.path.join(data, 'catalogs')
horus = os.path.join(catalogs, 'horus.txt')
emrcmt = os.path.join(catalogs, 'emrcmt.csv')
bollettino = os.path.join(catalogs, 'bollettino.txt')
taroni2018 = os.path.join(catalogs, 'taroni2018.txt')
catalogs = join(data, 'catalogs')
horus = join(catalogs, 'horus.txt')
emrcmt = join(catalogs, 'emrcmt.csv')
bollettino = join(catalogs, 'bollettino.txt')
taroni2018 = join(catalogs, 'taroni2018.txt')
# Region
regions = join(data, "region")
italy_region = join(regions, "polygon.shp")
# 5 year models
dir_models_5yr = os.path.join(data, 'five-year-models-10.01/source/forecasts/')
five_years_src = [os.path.join(dir_models_5yr, i) for i in
dir_models_5yr = join(data, 'five-year-models-10.01/source/forecasts/')
five_years_src = [join(dir_models_5yr, i) for i in
sorted(os.listdir(dir_models_5yr)) if i.endswith('xml')]
# 10 year models
dir_models_10yr = os.path.join(data, 'ten-year-models-10.01/source/forecasts/')
ten_years_src = [os.path.join(dir_models_10yr, i) for i in
dir_models_10yr = join(data, 'ten-year-models-10.01/source/forecasts/')
ten_years_src = [join(dir_models_10yr, i) for i in
sorted(os.listdir(dir_models_10yr)) if i.endswith('xml')]
# Processed models
models_five_years = os.path.join(data, 'five-year-models-10.01/models_5yr.obj')
models_ten_years = os.path.join(data, 'ten-year-models-10.01/models_10yr.obj')
models_five_years = join(data, 'five-year-models-10.01/models_5yr.obj')
models_ten_years = join(data, 'ten-year-models-10.01/models_10yr.obj')
# Results folder
results = os.path.join(main, 'results')
results = join(main, 'results')
# Path to store the results
......@@ -48,33 +52,38 @@ results = os.path.join(main, 'results')
$testtype ...
"""
store = os.path.join(results, 'stored_results')
stored_csep = os.path.join(store, 'csep_tests')
store = join(results, 'stored_results')
stored_csep = join(store, 'csep_tests')
spatent_logs = join(store, 'spatent')
store_years = {i: os.path.join(stored_csep, '%iyr') % i for i in years}
store_years = {i: join(stored_csep, '%iyr') % i for i in years}
spatent_years = {i: join(spatent_logs, '%iyr') % i for i in years}
tests = ['N', 'M', 'S', 'L', 'CL', 'T', 'PGS']
for i, j in store_years.items():
os.makedirs(j, exist_ok=True)
for i, j in spatent_years.items():
os.makedirs(j, exist_ok=True)
# Figures
figs = os.path.join(results, 'figures')
figs = join(results, 'figures')
## Forecast figures
forecast_figures = {5: os.path.join(data, 'five-year-models-10.01/figures'),
10: os.path.join(data, 'ten-year-models-10.01/figures')}
forecast_figures = {5: join(data, 'five-year-models-10.01/figures'),
10: join(data, 'ten-year-models-10.01/figures')}
for dir_ in forecast_figures.values():
os.makedirs(dir_, exist_ok=True)
## CSEP figures
csep_figs = os.path.join(figs, 'csep')
csep_figs = join(figs, 'csep')
os.makedirs(csep_figs, exist_ok=True)
## Spatial_analysis figures
spatial_figs = os.path.join(figs, 'spatial')
os.makedirs(spatial_figs, exist_ok=True)
spatial_figs = join(figs, 'spatial')
spat_entropy_figs = join(spatial_figs, 'spat_entropy')
ripley_k = join(spatial_figs, 'ripley_k')
os.makedirs(spat_entropy_figs, exist_ok=True)
os.makedirs(ripley_k, exist_ok=True)
# Create dirs
# os.makedirs(stored_csep, exist_ok=True)
......@@ -92,7 +101,18 @@ def get_csep_result_path(test, year, type=False):
type = '_' + str(type)
else:
type = ''
return os.path.join(store_years[year], test + type)
return join(store_years[year], test + type)
def get_spatent_result_path(method, year):
"""
Get the file path of a result object
:param test: Name of the test
:param year: Year duration of the test
:param type: Particular type of test
:return:
"""
return join(spatent_years[int(year)], 'log_%s_%s.txt' % (method, str(year)))
def get_csep_result(test, year, type=False):
"""
......@@ -104,7 +124,18 @@ def get_csep_result(test, year, type=False):
"""
return pickle.load(open(get_csep_result_path(test, year, type), 'rb'))
def get_csep_fig_path(test, year, type=None):
def get_csep_result(test, year, type=False):
"""
Get the saved CSEPEvaluationResult object of a test
:param test: Name of the test
:param year: Year duration of the test
:param type: Particular type of test
:return:
"""
return pickle.load(open(get_csep_result_path(test, year, type), 'rb'))
def get_csep_figpath(test, year, type=None):
"""
Get the figure path of a particular test
:param test: Name of the test
......@@ -116,9 +147,9 @@ def get_csep_fig_path(test, year, type=None):
type = '_' + str(type)
else:
type = ''
return os.path.join(csep_figs, "%s_%s%s.png" % (test, str(year), type))
return join(csep_figs, "%s_%s%s.png" % (test, str(year), type))
def get_spatial_fig_path(analysis, year, model):
def get_ripleyk_figpath(analysis, year, model):
"""
Get the figure path of a particular test
:param test: Name of the test
......@@ -126,5 +157,15 @@ def get_spatial_fig_path(analysis, year, model):
:param type: Particular type of test
:return:
"""
return os.path.join(spatial_figs, "%s_%s_%s.png" % (analysis, str(year), model))
return join(ripley_k, "%s_%s_%s.png" % (analysis, str(year), model))
def get_sentropy_figpath(analysis, year, model):
"""
Get the figure path of a particular test
:param test: Name of the test
:param year: Year duration of the test
:param type: Particular type of test
:return:
"""
return join(spat_entropy_figs, "%s_%s_%s.png" % (analysis, str(year), model))
import os
import cartopy
import paths
import get_models
plot_properties = {'grid_labels': True,
'borders': True,
'feature_lw': 0.5,
'basemap': 'stock_img',
'cmap': 'rainbow',
'clim' : [-8, 0],
'projection': cartopy.crs.Mercator()}
# forecasts = get_models.five_years()
# figures_dir = paths.forecast_figures[5]
#
# for model_i in forecasts:
# plot_properties['filename'] = os.path.join(figures_dir, model_i.name)
# plot_properties['title'] = model_i.name + ' ' + model_i.author
# print(plot_properties['title'])
# model_i.plot(plot_args=plot_properties)
#
# forecasts = get_models.ten_years()
# figures_dir = paths.forecast_figures[10]
# for model_i in forecasts:
# plot_properties['filename'] = os.path.join(figures_dir, model_i.name)
# plot_properties['title'] = model_i.name + ' ' + model_i.author
# print(plot_properties['title'])
# model_i.plot(plot_args=plot_properties)
### Contours
plot_properties = {'grid_labels': True,
'borders': True,
'feature_lw': 0.5,
'basemap': 'stock_img',
'cmap': 'rainbow',
'clim' : [-8, 0],
'projection': cartopy.crs.Mercator()}
forecasts = get_models.ten_years()
figures_dir = paths.forecast_figures[10]
for model_i in forecasts:
plot_properties['filename'] = os.path.join(figures_dir, 'fig_bins/%s_8-0-15' % model_i.name)
plot_properties['title'] = model_i.name + ' ' + model_i.author
print(plot_properties['title'])
model_i.plot(plot_args=plot_properties)
import numpy
import resource
usage = resource.getrusage(resource.RLIMIT_NICE).ru_maxrss
print(usage)
\ No newline at end of file
......@@ -13,106 +13,111 @@ from multiprocessing import Pool
import os
import time
def get_Hxy_2d(A, B, bins, range_):
hist2d = np.histogram2d(A.ravel(), B.ravel(), bins, [range_,range_],normed=True)
def get_Hxy(A, B, bins=100, range_='global', log=True):
if log:
a = A.ravel()[~np.isnan(A.ravel())]
b = B.ravel()[~np.isnan(B.ravel())]
a[a==0] = 1e-17
b[b == 0] = 1e-17
a = np.log10(a)
b = np.log10(b)
else:
a = A.ravel()[~np.isnan(A.ravel())]
b = B.ravel()[~np.isnan(B.ravel())]
if range_ == 'global':
range_ = [[np.min([np.min(a), np.min(b)]),
np.max([np.max(a), np.max(b)])],
[np.min([np.min(a), np.min(b)]),
np.max([np.max(a), np.max(b)])]]
elif range_ == 'local':
range_ = [[np.min(a), np.max(a)],
[np.min(b), np.max(b)]]
elif range_ == 'q_local':
range_ = [[np.quantile(a,0.02), np.quantile(a, 0.98)],
[np.quantile(b,0.02), np.quantile(b, 0.98)]]
elif range_ == 'q_global':
new_array = np.hstack((a,b))
range_ = [[np.quantile(new_array,0.05), np.quantile(new_array, 0.95)],
[np.quantile(new_array,0.05), np.quantile(new_array, 0.95)]]
hist2d = np.histogram2d(a, b, bins, range_, density=True)
dx = np.unique(np.diff(hist2d[1]))[0]*np.unique(np.diff(hist2d[2]))[0]
hist2d_dx = hist2d[0]*dx
Hxy = 0
for i in hist2d[0]:
for j in i[~np.isnan(i)]:
if j or j!=0.0:
Hxy -= j*np.log(j)
return Hxy
def get_Hx(A, bins, range_):
hist = np.histogram(A.ravel(),bins,range_)
Hx = 0
for i in hist[0][~np.isnan(hist)]:
if i!=0:
Hx-= i*np.log(i)
for i in hist2d_dx:
for j in i:
if j != 0.0:
Hxy -= j*np.log10(j)
return Hxy
def get_Hx(A, bins=30, range_='global', log=True):
HX = A.shape[0]*A.shape[1]
def entropy2d(a):
b = a.ravel()[~np.isnan(a.ravel())]
if b.size:
b = b[np.where(b.ravel()!=0.)]
return reduce((lambda x,y: x-y*np.log(y)),b)
if log:
a = A.ravel()[~np.isnan(A.ravel())]
a[a==0] = 1e-17
a = np.log10(a)
else:
return 0
def entropy1d(a):
a = A.ravel()[~np.isnan(A.ravel())]
if range_ == 'local' or range_ == 'global':
range_ = [np.min(a), np.max(a)]
elif range_ == 'q_local' or range_ == 'q_global':
range_ = [np.quantile(a, 0.02), np.quantile(a, 0.98)]
hist = np.histogram(a, bins,range_, density=True)
hist_dx = hist[0]*np.unique(np.diff(hist[1]))[0]
Hx = 0
for i in a:
for i in hist_dx:
if i!=0:
Hx += -i*np.log2(i)
else:
pass
Hx-= i*np.log10(i)
return Hx
def get_NSMI_2D(A, B, bins, range_):
m = np.size(A, 0)
n = np.size(A, 1)
p = np.size(A, 2)