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

added inhomogeneous pointprocess

parent d986d221
......@@ -2,7 +2,7 @@
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="jdk" jdkName="Python 3.8 (venv)" jdkType="Python SDK" />
<orderEntry type="inheritedJdk" />
<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.8 (venv)" 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
......@@ -41,8 +41,8 @@ def merge_results_batchs(batch1, 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'
Forecast_5.sim_name = Forecast_5.sim_name + ' - 10yr'
Forecast_10.sim_name = '5yr'
CL_merged.append(Forecast_5)
CL_merged.append(Forecast_10)
......@@ -157,7 +157,7 @@ def run(use_saved=False):
if __name__ == '__main__':
use_saved = False
use_saved = True
run(use_saved)
......
......@@ -47,7 +47,7 @@ def plot_single_results(Results, years, ref_model=0,folder=False):
'linewidth': 1.5,
'capsize': 0,
'markersize': 6,
'ylims': (-3,3)}
'ylims': (-4,4)}
T_Results_per_model = Results[0][ref_model]
......@@ -94,10 +94,10 @@ def plot_all_results(Results, years, p=0.05, order=False, savepath=False):
for m, j in enumerate(i):
if j < p:
ax.scatter(n + 0.5, m + 0.5,marker='o', s=2, color='black')
for n, i in enumerate(data_tq):
for m, j in enumerate(i):
if j > 0:
ax.scatter(n + 0.5, m + 0.75,marker='o', s=4, color='green')
# for n, i in enumerate(data_tq):
# for m, j in enumerate(i):
# if j > 0:
# ax.scatter(n + 0.5, m + 0.75,marker='o', s=4, color='green')
plt.text(len(score) + 2.5, -0.2, '$G_t$', fontsize=20)
for i, j in enumerate([score[i] for i in arg_ind]):
......@@ -135,10 +135,10 @@ def run(use_saved=False):
tw_figs_path = os.path.join(paths.csep_figs, 'tw_single')
os.makedirs(tw_figs_path, exist_ok=True)
# for yr in years:
# for model_i in range(len(models[yr])):
# plot_single_results(Results[yr], years=yr, ref_model=model_i,
# folder=tw_figs_path)
for yr in years:
for model_i in range(len(models[yr])):
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], years=yr, order=typ,
savepath=paths.get_csep_fig_path(test, yr, type=typ * 'order'))
......
......@@ -211,7 +211,7 @@ def plot_results(years, labels, ll_type=None, n_eqk = 1,
Axes = plot_legends(Axes, colors, labels)
Axes.set_title('Test results - %i years (%s)' % (years, ref_name) , pad=15, fontsize=14)
if savepath:
plt.savefig(savepath, dpi=300)
plt.savefig(savepath, dpi=300, transparent=True)
plt.show()
......
......@@ -71,6 +71,11 @@ for dir_ in forecast_figures.values():
csep_figs = os.path.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)
# Create dirs
# os.makedirs(stored_csep, exist_ok=True)
# os.makedirs(csep_figs, exist_ok=True)
......@@ -112,3 +117,14 @@ def get_csep_fig_path(test, year, type=None):
else:
type = ''
return os.path.join(csep_figs, "%s_%s%s.png" % (test, str(year), type))
def get_spatial_fig_path(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 os.path.join(spatial_figs, "%s_%s_%s.png" % (analysis, str(year), model))
......@@ -37,13 +37,14 @@ def K_r(Events, rates, Dists, r):
for i in range(n):
for j in range(n):
if j != i:
count = np.sum(Dists[i, :] < r)
mask = (Dists[i, :] < r) & (Dists[i, :] > 0)
count = np.sum(mask)
K_r += count / rates[i] / rates[j]
return K_r/899900
#################### CATALOG
plt.close('all')
cat = get_catalogs.filter_10yr_01_2010(get_catalogs.emrcmt())
cat = get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino())
region = cat.region
counts= cat.spatial_counts()
......@@ -59,8 +60,19 @@ Events = np.array([[i, j] for i,j in zip(cat.get_longitudes(), cat.get_latitudes
## Forecast
models = get_models.ten_years()
R_s = np.arange(0, 400, 20)
for model_i in range(6):
n_r = 25
r_max = 1000
range_models = range(19)
savefig = True
n_sim = 50
R_s = (np.logspace(0, 3, num =n_r)-1).astype(int)
T_x = []
T = []
for model_i in range_models:
print('Analyzing model %s' % models[model_i].name)
plt.figure(figsize=(6,4))
plt.title('$K-$function\n %s' % models[model_i].name)
forecast_data = models[model_i].spatial_counts()
......@@ -80,16 +92,14 @@ for model_i in range(6):
scale = n_obs / n_fore
expected_forecast_count = int(n_obs)
num_events_to_simulate = expected_forecast_count
print(num_events_to_simulate)
sampling_weights = numpy.cumsum(forecast_data.ravel()) / numpy.sum(forecast_data)
n_sim = 50
K_sim = []
for i in range(n_sim):
sim_fore = numpy.zeros(sampling_weights.shape)
sim = csep.core.poisson_evaluations._simulate_catalog(num_events_to_simulate, sampling_weights, sim_fore)
ind_sim = np.where(sim != 0)[0]
ind_sim = np.where(sim != 0)[0] #todo fix 2 events per cell
events_sim = np.array([models[model_i].get_longitudes()[ind_sim], models[model_i].get_latitudes()[ind_sim]]).T
dists_sim = dists(events_sim)
rates_sim = forecast_data[ind_sim]
......@@ -98,7 +108,9 @@ for model_i in range(6):
K.append(K_r(events_sim, rates_sim, dists_sim, i))
K_sim.append(K)
K_sim = np.array(K_sim)
K_avg = np.median(K_sim, axis=0)
K_std = np.std(K_sim)
......@@ -112,21 +124,36 @@ for model_i in range(6):
# color='gray', alpha=0.2)
L = np.divide(np.sqrt(K_sim), K_std)
L_avg = np.mean(L, axis=0)
L_up = np.quantile(L, 0.9, axis=0)
L_down = np.quantile(L, 0.1, axis=0)
L_std = np.std(L, axis=0)
L_up = np.quantile(L, 0.95, axis=0)
# L_up = L_avg + L_std
L_down = np.quantile(L, 0.05, axis=0)
# L_down = L_avg - L_std
L_cat = np.divide(np.sqrt(K_cat), K_std)
plt.plot(R_s, L_cat, color='red')
plt.plot(R_s, L_avg)
plt.plot(R_s, L_avg, '--', color='blue')
plt.fill_between(R_s, L_down, L_up,
color='gray', alpha=0.2)
color='gray', alpha=0.3)
plt.xlabel('$r$ [km]')
plt.ylabel('$L(r)=\frac{K(r)}{\pi}$')
plt.ylabel(r'$L(r)= \frac{K(r)}{\pi}$')
plt.legend([r'Observed $L(r)$', r'Simulated average $\overline{L(r)}$', '$0.95$ confidence intervals'], loc=4, fontsize=8)
plt.tight_layout()
plt.grid()
if savefig:
plt.savefig(paths.get_spatial_fig_path('K', 10, models[model_i].name), dpi=300)
plt.show()
# plt.figure()
# plt.plot(R_s, T_x_i)
# plt.show()
# cmap = plt.get_cmap('jet')
# fig = plt.figure()
# # a=plt.imshow(np.flipud(A), cmap=cmap)
......
import matplotlib.pyplot as plt
import numpy as np
import os
import entropy_calc as ec
from codes import get_models
from codes import paths
import numpy
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator
import numpy as np
from codes import get_catalogs
import csep
import itertools
import cartopy
#################### CATALOG
plot_properties = {'grid_labels': True,
'borders': True,
'feature_lw': 0.5,
'basemap': 'ESRI_terrain',
'cmap': 'rainbow',
'alpha_exp': 0.8,
'clim' : [-8, 0],
'projection': cartopy.crs.Mercator()}
cat_props = {'markersize': 3,
'markercolor': 'red',
'alpha': 0.6}
sim_props = {'markersize': 3,
'markercolor': 'blue',
'alpha': 0.6}
cat = get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino())
region = cat.region
counts= cat.spatial_counts()
A = np.zeros(region.idx_map.shape)
A[ A == 0] = np.nan
for i in range(A.shape[0]):
for j in range(A.shape[1]):
if not np.isnan(region.idx_map[i,j]):
A[i, j] = counts[int(region.idx_map[i,j])]
cat_ind = cat.get_spatial_idx()
Events = np.array([[i, j] for i,j in zip(cat.get_longitudes(), cat.get_latitudes())])
models = get_models.ten_years()
range_models = range(1,2)
savefig = True
n_sim = 100
T_x = []
T = []
for model_i in range_models:
forecast_data = models[model_i].spatial_counts()
rates_eqk = forecast_data[cat_ind]
plot_properties['title'] = models[model_i].name
forecast_data = models[model_i].spatial_counts()
observed_data = cat.spatial_counts()
n_obs = numpy.sum(observed_data)
n_fore = numpy.sum(forecast_data)
scale = n_obs / n_fore
expected_forecast_count = int(n_obs)
num_events_to_simulate = expected_forecast_count
sampling_weights = numpy.cumsum(forecast_data.ravel()) / numpy.sum(forecast_data)
#
Ax = models[model_i].plot(plot_args=plot_properties)
cat_props['title'] = models[model_i].name
cat_props['filename'] = paths.figs + '/sims/' + models[model_i].name+'_CAT.png'
Ax1 = cat.plot(ax=Ax, plot_args=cat_props)
for i in range(n_sim):
sim_fore = numpy.zeros(sampling_weights.shape)
sim = csep.core.poisson_evaluations._simulate_catalog(num_events_to_simulate, sampling_weights, sim_fore)
ind_sim = np.where(sim != 0)[0] #todo fix 2 events per cell
events_sim = np.array([models[model_i].get_longitudes()[ind_sim], models[model_i].get_latitudes()[ind_sim]]).T
rates_sim = forecast_data[ind_sim]
tuple = [ (1, 1, i[1], i[0], 5,10) for i in events_sim]
sim_cat = csep.catalogs.CSEPCatalog(data=tuple)
Ax = models[model_i].plot(plot_args=plot_properties)
sim_props['title'] = models[model_i].name
sim_props['filename'] = paths.figs + '/sims/' + models[model_i].name+'_%i.png'%i
Ax1 = cat.plot(ax=Ax, plot_args=cat_props)
sim_cat.region=region
sim_cat.plot(ax=Ax,plot_args=sim_props, show=True)
\documentclass[9pt,ignorenonframetext,mathserif,aspectratio=169]{beamer}
\usetheme{CambridgeUS}
\usepackage[3D]{movie15}
%\documentclass[9pt,ignorenonframetext,mathserif,handout]{beamer}
%\usetheme[height=7mm]{Madrid}
%\usepackage{pgfpages}
%\pgfpagesuselayout{2 on 1}[letterpaper,border shrink=12mm]
%\pgfpageslogicalpageoptions{1}{border code=\pgfusepath{stroke}}
%\pgfpageslogicalpageoptions{2}{border code=\pgfusepath{stroke}}
%%%%%%%%%%%%%%
%% Beamer options %%
%%%%%%%%%%%%%%
\setbeamertemplate{items}[ball]
\setbeamertemplate{blocks}[rounded][shadow=true]
\setbeamertemplate{navigation symbols}{}
\usepackage[makeroom]{cancel}
%\usepackage[T1]{fontenc}
%\usepackage{fourier}
% Uncomment for Article/Notes mode - comments are included
%\documentclass{article}
%\usepackage{beamerarticle}
\usepackage{relsize,multirow}
\usepackage{natbib}
%\usepackage{multimedia}
%\usepackage{algpseudocode}
% User-def commands
% Calculus symbols
\newcommand{\pd}[2]{\frac{\partial #1}{ \partial #2}} % First partial derivative command
\newcommand{\td}[2]{\frac{\mathrm{d} #1}{ \mathrm{d} #2}} % First derivative
\newcommand{\md}[2]{\frac{\mathrm{D} #1}{ \mathrm{D} #2}} % Material derivative
\newcommand{\pdd}[2]{\frac{\partial^2 #1}{ \partial #2 ^2}} % Second partial derivative command
\newcommand{\pddc}[3]{\frac{\partial^2 #1}{ \partial #2 \partial #3}} % Second partial derivative command
\newcommand{\bs}[1]{\boldsymbol{#1}}
\newcommand{\mbf}[1]{\mathbf{#1}}
\newcommand{\mbfh}[1]{\hat{\mathbf{#1}}}
\newcommand{\etal}{\textit{et al}. }
% Continuum mechanics & FEM symbols
\def\sca #1{\mbox{\rm{#1}}{}}
\def\mat #1{\mbox{\boldmath $\mathsf #1$}}
\def\vec #1{\mbox{\boldmath $#1$}{}}
\def\ten #1{\mbox{\boldmath $#1$}{}}
\def\Ass {\overset{\hspace*{0.4cm} n_{\scas{el}}}
{\underset{\scaf{c},\scaf{d}=1}{\msf{A}}}}
\def\ltr #1{\mbox{\sf{#1}}}
\def\bltr #1{\mbox{\sffamily{\bfseries{{#1}}}}}
% math operators and symbols
\DeclareMathOperator{\trace}{trace}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\symm}{symm}
\DeclareMathOperator{\sk}{skew}
\DeclareMathOperator{\grad}{grad}
\DeclareMathOperator{\Grad}{Grad}
\DeclareMathOperator{\dive}{div}
\DeclareMathOperator{\Dive}{Div}
\DeclareMathOperator{\curl}{curl}
\DeclareMathOperator{\Curl}{Curl}
%\DeclareMathOperator{\Assembly}{\mbox{\mathbb{A}} }
\DeclareMathOperator*{\A}{ \mathlarger{\mathlarger{\boldsymbol{\mathbb{A}}}} }
\newcommand{\vectwo}[2]{\begin{bmatrix} #1 \\ #2 \end{bmatrix}}
\newcommand{\mattwo}[4]{\begin{bmatrix} #1 & #2\\ #3 & #4 \end{bmatrix}}
\newcommand{\vecthree}[3]{\begin{bmatrix} #1 \\ #2 \\ #3 \end{bmatrix}}
\newcommand{\arginf}{\mathop{\mathrm{arg\,inf}}}
\newcommand{\argsup}{\mathop{\mathrm{arg\,sup}}}
\def\d#1{\mbox{$\,\mathrm{d}#1$}}
\def\order#1{\mbox{$\mathcal{O}(#1)$}}
\def\python{\mbox{\tt Python}}
\def\R{\mbox{\(\mathbb{R}\)}}
\def\dx{\mbox{\(\,\mathrm{d}x\)}}
\def\TODO{\mbox{\tt TODO}}
\usepackage[applemac]{inputenc} % Allows for accents!! Useful for spanish documents
%%% Theorem declaration %%%
\newtheorem{thm}{\bf{Teorema}}
\newtheorem{lema}{\bf{Lema}}
%\newtheorem{example}{\bf{Ejemplo}}
\newtheorem{defn}{\bf{Definici\'{o}n}}
\newtheorem{prop}{\bf{Proposici\'{o}n}}
%\newtheorem{corollary}{\bf{Corolario}}
%\newtheorem{proof}{\bf{Demostraci\'{o}n}}
\newtheorem{soln}{\bf{Soluci\'{o}n}}
\newtheorem{rmk}{\bf{Nota}}
\newtheorem{objgen}{\bf{Objetivo}}
% Different font families
%\usepackage{helvet}
%\usepackage{palatino}
%\usepackage{xmpmulti}
%\usepackage{multimedia}
%\usepackage{movie15}
\usepackage{moreverb}
\usepackage{amsmath, amsfonts, amsthm,alltt}
\usepackage{multirow, hyperref}
\usepackage[absolute,showboxes,verbose,overlay]{textpos}
% Paquete para escribir codigo
\usepackage{listings}
\lstset{basicstyle=\footnotesize\ttfamily,breaklines=true}
\usepackage{alltt}
% page numbering
%\setbeamertemplate{footline}{\begin{center} \small \insertframenumber \end{center}}
%\setbeamertemplate{headerline}{\begin{center} \small \insertframenumber \end{center}}
%\setbeamertemplate{footline}{\hfill\small\insertframenumber/\inserttotalframenumber}1
%\usepackage{beamerarticle} // not installed
%\usepackage{beamerthemesplit} // Activate for custom appearance
\makeatother
\setbeamertemplate{footline}
{\hfill\insertframenumber/ \inserttotalframenumber\hspace*{1ex}}
\date{\tiny 12 of November, 2020}
%%%%%%%%
% Author data %
%%%%%%%%
\vspace{2cm}
\title{Testing the first Italy Experiment:\\ Results and analysis of the spatial distribution of long-term
forecasts}
\author{Pablo Iturrieta\inst{1,2}, William Savran\inst{3}, Asim Khawaja\inst{1,2}, W. Marzocchi\inst{4}, M.
Taroni\inst{5}, G. Falcone \inst{5} \& Danijel Schorlemmer\inst{1}\\ [6pt]}
\institute[GFZ]{\inst{1} GFZ German Research Centre for Geosciences \and \inst{2} University of Potsdam, Germany
\and \inst{3} University of Southern California, Southern California Earthquake Center
\and \inst{4} University of Naples \and \inst{5} Istituto Nazionale di Geofisica e Vulcanologia }
%%%%%%%%
% Document %
%%%%%%%%
\begin{document}
\begin{frame}[t]
\mbox{}%
\begin{minipage}{0.25\textwidth}
\centering
\includegraphics[width=0.6\textwidth]{../../figures/old/logo_rise.jpeg}
\end{minipage}%
\hfill%
\begin{minipage}{0.25\textwidth}
\centering
\includegraphics[width=0.8\textwidth]{../../figures/old/logo_csep.png}
\end{minipage}%
\hfill%
\begin{minipage}{0.25\textwidth}
\centering
\includegraphics[width=0.65\textwidth]{../../figures/old/logo_gfz.png}
\end{minipage}%
\begin{minipage}{0.25\textwidth}
\centering
\includegraphics[width=0.5\textwidth]{../../figures/old/logo_potsdam.png}
\end{minipage}%
\mbox{}%
\vspace{0.cm}
\maketitle
\end{frame}
\begin{frame}[t]
\frametitle{The 10-year Italy Experiment}
\begin{columns}[t]
\column{.45\textwidth}
\textbf{\large Prospective testing }
\begin{itemize}
\item Experiment definition by CSEP\footnotemark as a consensus with researchers, modelers and catalog makers.
\begin{itemize}
\item Temporal windows: 1-day, 3-month, 5-year, 10-year.
\item Discretization: $0.1^\circ \times 0.1^\circ$ and $\Delta m_w=0.1$ spatial and magnitude binning.
\item $M_{target}\geq 4.95$ for 10-year forecast
\item Evaluation metrics
\item Learning data
\item Observation data
\end{itemize}
\end{itemize}
\column{.55\textwidth}
\begin{figure}
\vspace*{-1.5cm}
\includegraphics[width=0.8\textwidth]{../../figures/old/events_epsg3857.png}
\end{figure}
\end{columns}
\footnotetext[1]{Schorlemmer et al., 2010}
\footnotetext[2]{Taroni et al., 2018}
\end{frame}
\begin{frame}[t]
\frametitle{The 10-year Italy Experiment}
\begin{columns}[t]
\column{.4\textwidth}
\textbf{\large The forecasts}
\begin{itemize}
\item Rate-based forecast
\item Rate variability based on different underlying assumptions. E.g:
\begin{itemize}
\item Time dependence
\item Historical/Instrumental seismicity
\item Physics based
\item GR parameters variability
\item etc.
\end{itemize}
\item Strong variability in the models' spatial distribution
\end{itemize}
\column{.6\textwidth}
\only<1>{
\begin{figure}
\vspace*{-1cm}
\includegraphics[width=1.1\textwidth]{../../figures/old/map_TripleS-hybrid_10yr.png}
\end{figure}
}
\only<2>{
\begin{figure}
\vspace*{-1cm}
\includegraphics[width=1.1\textwidth]{../../figures/old/map_ALM_IT_10yr.png}
\end{figure}
}
\only<3>{
\begin{figure}
\vspace*{-1cm}
\includegraphics[width=1.1\textwidth]{../../figures/old/map_MPS04_10yr.png}
\end{figure}
}
\end{columns}
\end{frame}
\begin{frame}[t]
\frametitle{The 10-year Italy Experiment}
\begin{columns}[t]
\column{.4\textwidth}
\textbf{\large The results }
\begin{itemize}
\item The $\mathcal{L}_N$-test (Conditional joint Log-likelihood)
\begin{itemize}
\item Based in the Log-likelihood of a catalog $\Omega$ under a forecast $\Lambda$:\\
$\mathcal{L}_{\Omega|\Lambda}=\sum\limits_{i=1}^n\omega_i\ln\lambda_i-\lambda_i-\ln\omega_i!$\\
where $\omega_i$ is the event count and $\lambda_i$ the expected value in the bin $i$