From 94aa79b3bc65a02231e76a5f3055adc5b124a56f Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 16 Jul 2019 17:39:31 +0200
Subject: [PATCH] [PAPER1] modify the strength criterion for the hypercube (use
the 0.98 quantile and its evolution in 10 years)
---
.../scm_models_data/abstract_study.py | 2 +-
.../2-main_choice_starting_years.py | 85 -------------------
.../abstract_hypercube_visualizer.py | 8 +-
.../altitude_hypercube_visualizer.py | 49 ++++++++---
.../study_visualization/study_visualizer.py | 6 +-
.../comparisons_visualization.py | 8 +-
.../paper1_steps/__init__.py | 0
.../__init__.py | 0
.../main1_good_stationary_gev_fit.py} | 0
.../main2_choice_to_not_use_starting_years.py | 36 ++++++++
...main3_non_stationary_strength_evolution.py | 40 +++++++++
experiment/paper1_steps/utils.py | 29 +++++++
.../abstract_gev_change_point_test.py | 43 ++++++++--
.../abstract_univariate_test.py | 2 +-
.../trend_analysis/univariate_test/utils.py | 2 +-
.../temporal_linear_margin_model.py | 2 +-
.../extreme_models/result_from_fit.py | 6 +-
.../margin_fits/gev/gev_params.py | 18 ++++
18 files changed, 219 insertions(+), 117 deletions(-)
delete mode 100644 experiment/meteo_france_data/scm_models_data/paper1_steps/2-main_choice_starting_years.py
rename experiment/{meteo_france_data/scm_models_data => }/paper1_steps/__init__.py (100%)
create mode 100644 experiment/paper1_steps/hard extreme evolution - annual maxima/__init__.py
rename experiment/{meteo_france_data/scm_models_data/paper1_steps/1-main_good_stationary_gev_fit.py => paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py} (100%)
create mode 100644 experiment/paper1_steps/hard extreme evolution - annual maxima/main2_choice_to_not_use_starting_years.py
create mode 100644 experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py
create mode 100644 experiment/paper1_steps/utils.py
diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py
index df09b2ed..4232785c 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -299,7 +299,7 @@ class AbstractStudy(object):
x, y = list(row)
massif_name = row.name
value = massif_name_to_value[massif_name]
- str_value = str(round(value, 1)) if isinstance(value, str) else str(value)
+ str_value = str(value)
ax.text(x, y, str_value, horizontalalignment='center', verticalalignment='center', fontsize=7)
if scaled:
diff --git a/experiment/meteo_france_data/scm_models_data/paper1_steps/2-main_choice_starting_years.py b/experiment/meteo_france_data/scm_models_data/paper1_steps/2-main_choice_starting_years.py
deleted file mode 100644
index 51e242af..00000000
--- a/experiment/meteo_france_data/scm_models_data/paper1_steps/2-main_choice_starting_years.py
+++ /dev/null
@@ -1,85 +0,0 @@
-import time
-
-from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusRecentSwe
-from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_hypercube_visualizer import \
- AltitudeHypercubeVisualizer
-from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_hypercube_visualizer_extended import \
- AltitudeHypercubeVisualizerBisExtended, QuantityHypercubeWithoutTrendExtended, \
- AltitudeHypercubeVisualizerWithoutTrendExtended, QuantityHypercubeWithoutTrend
-from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
- Altitude_Hypercube_Year_Visualizer
-from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.main_files.main_fast_hypercube_one_altitudes import \
- get_fast_parameters
-from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.main_files.main_full_hypercube import \
- get_full_parameters
-from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.quantity_altitude_visualizer import \
- QuantityAltitudeHypercubeVisualizer
-from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.utils_hypercube import \
- load_altitude_visualizer, load_quantity_visualizer
-from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
- ALL_ALTITUDES, SCM_STUDIES
-from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import GevLocationChangePointTest, \
- GevScaleChangePointTest
-
-
-def get_fast_altitude_visualizer(altitude_hypercube_class):
- altitudes, last_starting_year, nb_data_reduced_for_speed, only_first_one, save_to_file, trend_test_class = get_fast_parameters()
- study_classes = [CrocusRecentSwe]
- visualizer = load_altitude_visualizer(altitude_hypercube_class, altitudes, last_starting_year,
- nb_data_reduced_for_speed, only_first_one, save_to_file, study_classes,
- trend_test_class)
- return visualizer
-
-
-def main_fast_old_spatial_repartition():
- # Simply the main graph
- get_fast_altitude_visualizer(Altitude_Hypercube_Year_Visualizer).visualize_massif_trend_test_one_altitude()
-
-
-def get_full_altitude_visualizer(altitude_hypercube_class, exact_starting_year=None, altitude=900):
- altitudes, first_starting_year, last_starting_year, nb_data_reduced_for_speed, only_first_one, save_to_file, trend_test_class = get_full_parameters(
- altitude=altitude)
- if exact_starting_year is not None:
- first_starting_year, last_starting_year = None, None
- study_classes = [CrocusRecentSwe]
- trend_test_class = GevScaleChangePointTest
- visualizer = load_altitude_visualizer(altitude_hypercube_class, altitudes, last_starting_year,
- nb_data_reduced_for_speed, only_first_one, save_to_file, study_classes,
- trend_test_class, first_starting_year=first_starting_year,
- exact_starting_year=exact_starting_year)
- return visualizer
-
-
-FULL_ALTITUDES = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
-HALF_ALTITUDES = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][::2]
-
-
-def main_fast_spatial_repartition():
- for altitude in FULL_ALTITUDES[-1:]:
- vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
- exact_starting_year=1958)
- vizualiser.save_to_file = False
- vizualiser.visualize_massif_trend_test_one_altitude()
-
-
-def main_full_spatial_repartition():
- for altitude in FULL_ALTITUDES[:]:
- # Compute for the most likely starting year
- # vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude)
- # vizualiser.visualize_massif_trend_test_one_altitude()
- # Compute the trend for a linear trend
- vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
- exact_starting_year=1958)
- vizualiser.visualize_massif_trend_test_one_altitude()
-
-
-def main_run():
- main_full_spatial_repartition()
- # main_fast_spatial_repartition()
-
-
-if __name__ == '__main__':
- start = time.time()
- main_run()
- duration = time.time() - start
- print('Full run took {}s'.format(round(duration, 1)))
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
index abffff98..0b0d2689 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
@@ -20,11 +20,13 @@ class AbstractHypercubeVisualizer(object):
def __init__(self, tuple_to_study_visualizer: Dict[Tuple, StudyVisualizer],
trend_test_class,
nb_data_reduced_for_speed=False,
+ reduce_strength_array=False,
save_to_file=False,
first_starting_year=None,
last_starting_year=None,
exact_starting_year=None,
verbose=True):
+ self.reduce_strength_array = reduce_strength_array
self.verbose = verbose
self.save_to_file = save_to_file
self.trend_test_class = trend_test_class
@@ -87,13 +89,17 @@ class AbstractHypercubeVisualizer(object):
)
@cached_property
- def df_hypercube_trend_strength(self) -> pd.DataFrame:
+ def df_hypercube_trend_slope_relative_strength(self) -> pd.DataFrame:
return self._df_hypercube_trend_meta(idx=1)
@cached_property
def df_hypercube_trend_nllh(self) -> pd.DataFrame:
return self._df_hypercube_trend_meta(idx=2)
+ @cached_property
+ def df_hypercube_trend_constant_quantile(self) -> pd.DataFrame:
+ return self._df_hypercube_trend_meta(idx=3)
+
# Some properties
@property
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
index 55b1c08d..c4079ad2 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
@@ -9,7 +9,9 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
SCM_STUDY_NAME_TO_COLOR, SCM_STUDY_NAME_TO_ABBREVIATION, SCM_STUDY_CLASS_TO_ABBREVIATION, SCM_STUDIES_NAMES
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
StudyVisualizer
+from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import AbstractGevChangePointTest
from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
+from extreme_estimator.margin_fits.gev.gev_params import GevParams
ALTITUDES_XLABEL = 'altitudes'
@@ -71,11 +73,13 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
assert not isinstance(s_trend_type_percentage.index, pd.MultiIndex)
s_trend_type_percentage *= 100
series = [s_trend_type_percentage]
- # # Reduce df_strength to a serie s_trend_strength
- # df_strength = self.df_hypercube_trend_strength[df_bool]
- # s_trend_strength = reduction_function(df_strength)
- # # Group result
- # series = [s_trend_type_percentage, s_trend_strength]
+ if self.reduce_strength_array:
+ # Reduce df_strength to a serie s_trend_strength
+ df_strength = self.df_hypercube_trend_slope_relative_strength[df_bool]
+ s_trend_strength = reduction_function(df_strength)
+ df_constant = self.df_hypercube_trend_constant_quantile[df_bool]
+ s_trend_constant = reduction_function(df_constant)
+ series.extend([s_trend_strength, s_trend_constant])
return series
def subtitle_to_reduction_function(self, reduction_function, level=None, add_detailed_plot=False, subtitle=None):
@@ -280,6 +284,8 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
massif_to_color = {}
add_text = self.nb_rows > 1
massif_to_year = {}
+ massif_to_strength = {}
+ massif_to_constant = {}
poster_trend_types = [AbstractUnivariateTest.SIGNIFICATIVE_POSITIVE_TREND,
AbstractUnivariateTest.SIGNIFICATIVE_NEGATIVE_TREND,
AbstractUnivariateTest.NEGATIVE_TREND,
@@ -292,15 +298,36 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
massif_to_color_for_trend_type = {k: color for k, v in dict(serie).items() if not np.isnan(v)}
massif_to_color.update(massif_to_color_for_trend_type)
if add_text:
- massif_to_year_for_trend_type = {k: int(v) for k, v in
- self.trend_type_to_series(reduction_function, isin_parameters)[
- display_trend_type][1].items()
- if k in massif_to_color_for_trend_type}
- massif_to_year.update(massif_to_year_for_trend_type)
+ if self.reduce_strength_array:
+ massif_to_value_for_trend_type = [{k: v for k, v in
+ self.trend_type_to_series(reduction_function,
+ isin_parameters)[
+ display_trend_type][i].items()
+ if k in massif_to_color_for_trend_type} for i in [1, 2]]
+ massif_to_strength.update(massif_to_value_for_trend_type[0])
+ massif_to_constant.update(massif_to_value_for_trend_type[1])
+ else:
+ massif_to_value_for_trend_type = {k: int(v) for k, v in
+ self.trend_type_to_series(reduction_function,
+ isin_parameters)[
+ display_trend_type][1].items()
+ if k in massif_to_color_for_trend_type}
+ massif_to_year.update(massif_to_value_for_trend_type)
+ # Compute massif_to_value
+ if self.reduce_strength_array:
+ massif_name_to_value = {m: "{} {}{} / {} year(s)".format(
+ int(massif_to_constant[m]),
+ "+" if massif_to_strength[m] > 0 else "",
+ round(massif_to_strength[m] * massif_to_constant[m], 1),
+ AbstractGevChangePointTest.nb_years_for_quantile_evolution)
+ for m in massif_to_strength}
+ else:
+ massif_name_to_value = massif_to_year
self.study.visualize_study(None, massif_name_to_color=massif_to_color, show=False,
show_label=False, scaled=True, add_text=add_text,
- massif_name_to_value=massif_to_year)
+ massif_name_to_value=massif_name_to_value)
+ print(subtitle)
title = self.set_trend_test_reparition_title(subtitle, set=False)
# row_title = self.get_title_plot(xlabel='massifs', ax_idx=i)
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
index 74dbf0a1..4e696e1a 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
@@ -463,11 +463,11 @@ class StudyVisualizer(VisualizationParameters):
sample_one_massif_from_each_region)
for massif_name, gev_change_point_test_results in massif_name_to_gev_change_point_test_results.items():
trend_test_res, best_idxs = gev_change_point_test_results
- trend_test_res = [(a, b, c) if i in best_idxs else (np.nan, np.nan, c)
- for i, (a, b, c, *_) in enumerate(trend_test_res)]
+ trend_test_res = [(a, b, c, d) if i in best_idxs else (np.nan, np.nan, c, np.nan)
+ for i, (a, b, c, d, *_) in enumerate(trend_test_res)]
massif_name_to_trend_res[massif_name] = list(zip(*trend_test_res))
nb_res = len(list(massif_name_to_trend_res.values())[0])
- assert nb_res == 3
+ assert nb_res == 4
all_massif_name_to_res = [{k: v[idx_res] for k, v in massif_name_to_trend_res.items()}
for idx_res in range(nb_res)]
return [pd.DataFrame(massif_name_to_res, index=self.starting_years_for_change_point_test).transpose()
diff --git a/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py b/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py
index 89bc4334..5ccecce6 100644
--- a/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py
+++ b/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py
@@ -251,10 +251,10 @@ class ComparisonsVisualization(VisualizationParameters):
label += "\n {} starting in {}".format(display_trend_type, most_likely_year)
ordered_dict[TREND_TYPE_CNAME] = display_trend_type
ordered_dict['most likely year'] = most_likely_year
- # Display the nllh against the starting year
+ # Display the deviance against the starting year
step = 1
- ax2.plot(starting_years[::step], [t[3] for t in trend_test_res][::step], color=plot_color, marker='o')
- ax2.plot(starting_years[::step], [t[4] for t in trend_test_res][::step], color=plot_color, marker='x')
+ ax2.plot(starting_years[::step], [t[4] for t in trend_test_res][::step], color=plot_color, marker='o')
+ ax2.plot(starting_years[::step], [t[5] for t in trend_test_res][::step], color=plot_color, marker='x')
# Plot maxima
ax.grid()
ax.plot(years, maxima, label=label, color=plot_color)
@@ -270,7 +270,7 @@ class ComparisonsVisualization(VisualizationParameters):
res = safe_run_r_estimator(function=r('gev.fit'), xdat=ro.FloatVector(data),
use_start=True)
res = ResultFromIsmev(res, {})
- gev_params = res.stationary_gev_params
+ gev_params = res.constant_gev_params
lim = 1.5 * max(data)
x = np.linspace(0, lim, 1000)
diff --git a/experiment/meteo_france_data/scm_models_data/paper1_steps/__init__.py b/experiment/paper1_steps/__init__.py
similarity index 100%
rename from experiment/meteo_france_data/scm_models_data/paper1_steps/__init__.py
rename to experiment/paper1_steps/__init__.py
diff --git a/experiment/paper1_steps/hard extreme evolution - annual maxima/__init__.py b/experiment/paper1_steps/hard extreme evolution - annual maxima/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/experiment/meteo_france_data/scm_models_data/paper1_steps/1-main_good_stationary_gev_fit.py b/experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py
similarity index 100%
rename from experiment/meteo_france_data/scm_models_data/paper1_steps/1-main_good_stationary_gev_fit.py
rename to experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py
diff --git a/experiment/paper1_steps/hard extreme evolution - annual maxima/main2_choice_to_not_use_starting_years.py b/experiment/paper1_steps/hard extreme evolution - annual maxima/main2_choice_to_not_use_starting_years.py
new file mode 100644
index 00000000..c43a29e7
--- /dev/null
+++ b/experiment/paper1_steps/hard extreme evolution - annual maxima/main2_choice_to_not_use_starting_years.py
@@ -0,0 +1,36 @@
+import time
+
+from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
+ Altitude_Hypercube_Year_Visualizer
+from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+
+
+def main_fast_spatial_repartition():
+ for altitude in FULL_ALTITUDES[-1:]:
+ vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
+ exact_starting_year=1958)
+ vizualiser.save_to_file = False
+ vizualiser.visualize_massif_trend_test_one_altitude()
+
+
+def main_full_spatial_repartition():
+ for altitude in FULL_ALTITUDES[:]:
+ # Compute for the most likely starting year
+ # vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude)
+ # vizualiser.visualize_massif_trend_test_one_altitude()
+ # Compute the trend for a linear trend
+ vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
+ exact_starting_year=1958)
+ vizualiser.visualize_massif_trend_test_one_altitude()
+
+
+def main_run():
+ # main_full_spatial_repartition()
+ main_fast_spatial_repartition()
+
+
+if __name__ == '__main__':
+ start = time.time()
+ main_run()
+ duration = time.time() - start
+ print('Full run took {}s'.format(round(duration, 1)))
diff --git a/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py b/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py
new file mode 100644
index 00000000..eccbb692
--- /dev/null
+++ b/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py
@@ -0,0 +1,40 @@
+import time
+
+from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
+ Altitude_Hypercube_Year_Visualizer
+
+"""
+Visualize the 0.99 quantile initial value and its evolution
+"""
+from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+
+
+def main_fast_spatial_risk_evolution():
+ for altitude in FULL_ALTITUDES[-1:]:
+ vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
+ exact_starting_year=1958, reduce_strength_array=True)
+ vizualiser.save_to_file = False
+ vizualiser.visualize_massif_trend_test_one_altitude()
+
+
+def main_full_spatial_repartition():
+ for altitude in FULL_ALTITUDES[:]:
+ # Compute for the most likely starting year
+ # vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude)
+ # vizualiser.visualize_massif_trend_test_one_altitude()
+ # Compute the trend for a linear trend
+ vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
+ exact_starting_year=1958)
+ vizualiser.visualize_massif_trend_test_one_altitude()
+
+
+def main_run():
+ # main_full_spatial_repartition()
+ main_fast_spatial_risk_evolution()
+
+
+if __name__ == '__main__':
+ start = time.time()
+ main_run()
+ duration = time.time() - start
+ print('Full run took {}s'.format(round(duration, 1)))
diff --git a/experiment/paper1_steps/utils.py b/experiment/paper1_steps/utils.py
new file mode 100644
index 00000000..699d710f
--- /dev/null
+++ b/experiment/paper1_steps/utils.py
@@ -0,0 +1,29 @@
+import time
+
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusRecentSwe
+from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
+ Altitude_Hypercube_Year_Visualizer
+from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.main_files.main_full_hypercube import \
+ get_full_parameters
+from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.utils_hypercube import \
+ load_altitude_visualizer
+from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import GevScaleChangePointTest, \
+ GevLocationChangePointTest
+
+FULL_ALTITUDES = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
+
+
+def get_full_altitude_visualizer(altitude_hypercube_class, exact_starting_year=None, altitude=900,
+ reduce_strength_array=False):
+ altitudes, first_starting_year, last_starting_year, nb_data_reduced_for_speed, only_first_one, save_to_file, trend_test_class = get_full_parameters(
+ altitude=altitude)
+ if exact_starting_year is not None:
+ first_starting_year, last_starting_year = None, None
+ study_classes = [CrocusRecentSwe]
+ trend_test_class = GevLocationChangePointTest
+ visualizer = load_altitude_visualizer(altitude_hypercube_class, altitudes, last_starting_year,
+ nb_data_reduced_for_speed, only_first_one, save_to_file, study_classes,
+ trend_test_class, first_starting_year=first_starting_year,
+ exact_starting_year=exact_starting_year)
+ visualizer.reduce_strength_array = reduce_strength_array
+ return visualizer
diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
index 93a49a12..f7acb27e 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
@@ -21,6 +21,9 @@ from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_tempor
class AbstractGevChangePointTest(AbstractUnivariateTest):
RRunTimeError_TREND = 'R RunTimeError trend'
+ # I should use the quantile from the Eurocode for the buildings
+ quantile_for_strength = 0.98
+ nb_years_for_quantile_evolution = 10
def __init__(self, years, maxima, starting_year, non_stationary_model_class, gev_param_name):
super().__init__(years, maxima, starting_year)
@@ -28,8 +31,12 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
df = pd.DataFrame({AbstractCoordinates.COORDINATE_T: years})
df_maxima_gev = pd.DataFrame(maxima, index=df.index)
observations = AbstractSpatioTemporalObservations(df_maxima_gev=df_maxima_gev)
- self.coordinates = AbstractTemporalCoordinates.from_df(df, transformation_class=CenteredScaledNormalization) # type: AbstractTemporalCoordinates
+ self.coordinates = AbstractTemporalCoordinates.from_df(df,
+ transformation_class=CenteredScaledNormalization) # type: AbstractTemporalCoordinates
self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
+ # For the moment, we chose:
+ # -the 0.99 quantile (even if for building maybe we should use the 1/50 so 0.98 quantile)
+ # -to see the evolution between two successive years
try:
# Fit stationary model
@@ -60,11 +67,8 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
return self.non_stationary_deviance - self.stationary_deviance
@property
- def non_stationary_nllh(self):
- if self.crashed:
- return np.nan
- else:
- return self.non_stationary_estimator.result_from_fit.nllh
+ def non_stationary_constant_gev_params(self) -> GevParams:
+ return self.non_stationary_estimator.result_from_fit.constant_gev_params
@property
def stationary_deviance(self):
@@ -80,6 +84,13 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
else:
return self.non_stationary_estimator.result_from_fit.deviance
+ @property
+ def non_stationary_nllh(self):
+ if self.crashed:
+ return np.nan
+ else:
+ return self.non_stationary_estimator.result_from_fit.nllh
+
@property
def is_significant(self) -> bool:
return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=1)
@@ -111,11 +122,23 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
return self.get_coef(self.non_stationary_estimator, AbstractCoordinates.COORDINATE_T)
@property
- def test_trend_strength(self):
+ def test_trend_slope_strength(self):
if self.crashed:
return 0.0
else:
- return self.percentage_of_change_per_year
+ slope = self._slope_strength()
+ slope *= self.nb_years_for_quantile_evolution * self.coordinates.transformed_distance_between_two_successive_years[0]
+ return slope
+
+ def _slope_strength(self):
+ raise NotImplementedError
+
+ @property
+ def test_trend_constant_quantile(self):
+ if self.crashed:
+ return 0.0
+ else:
+ return self.non_stationary_constant_gev_params.quantile(p=self.quantile_for_strength)
@property
def percentage_of_change_per_year(self):
@@ -135,6 +158,10 @@ class GevLocationChangePointTest(AbstractGevChangePointTest):
super().__init__(years, maxima, starting_year,
NonStationaryLocationStationModel, GevParams.LOC)
+ def _slope_strength(self):
+ return self.non_stationary_constant_gev_params.quantile_strength_evolution_ratio(p=self.quantile_for_strength,
+ mu1=self.non_stationary_linear_coef)
+
class GevScaleChangePointTest(AbstractGevChangePointTest):
diff --git a/experiment/trend_analysis/univariate_test/abstract_univariate_test.py b/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
index 52863f98..7381e71c 100644
--- a/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
@@ -105,7 +105,7 @@ class AbstractUnivariateTest(object):
return len(self.years)
@property
- def test_trend_strength(self):
+ def test_trend_slope_strength(self):
return 0.0
@property
diff --git a/experiment/trend_analysis/univariate_test/utils.py b/experiment/trend_analysis/univariate_test/utils.py
index 32703989..475a8855 100644
--- a/experiment/trend_analysis/univariate_test/utils.py
+++ b/experiment/trend_analysis/univariate_test/utils.py
@@ -9,7 +9,7 @@ from utils import NB_CORES
def compute_gev_change_point_test_result(smooth_maxima, starting_year, trend_test_class, years):
trend_test = trend_test_class(years, smooth_maxima, starting_year) # type: AbstractGevChangePointTest
assert isinstance(trend_test, AbstractGevChangePointTest)
- return trend_test.test_trend_type, trend_test.test_trend_strength, trend_test.non_stationary_nllh, trend_test.non_stationary_deviance, trend_test.stationary_deviance
+ return trend_test.test_trend_type, trend_test.test_trend_slope_strength, trend_test.non_stationary_nllh, trend_test.test_trend_constant_quantile, trend_test.non_stationary_deviance, trend_test.stationary_deviance
def compute_gev_change_point_test_results(multiprocessing, maxima, starting_years, trend_test_class,
diff --git a/extreme_estimator/extreme_models/margin_model/temporal_linear_margin_model.py b/extreme_estimator/extreme_models/margin_model/temporal_linear_margin_model.py
index 9f9a8057..6c7b064c 100644
--- a/extreme_estimator/extreme_models/margin_model/temporal_linear_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/temporal_linear_margin_model.py
@@ -18,7 +18,7 @@ class TemporalLinearMarginModel(LinearMarginModel):
def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,
df_coordinates_temp: pd.DataFrame) -> ResultFromFit:
- # Gev Fit
+ # Gev Fit with isMev package
assert data.shape[1] == len(df_coordinates_temp.values)
res = safe_run_r_estimator(function=r('gev.fit'), use_start=self.use_start_value,
xdat=ro.FloatVector(data[0]), y=df_coordinates_temp.values, mul=self.mul,
diff --git a/extreme_estimator/extreme_models/result_from_fit.py b/extreme_estimator/extreme_models/result_from_fit.py
index 006d5575..78d9d83a 100644
--- a/extreme_estimator/extreme_models/result_from_fit.py
+++ b/extreme_estimator/extreme_models/result_from_fit.py
@@ -48,6 +48,10 @@ class ResultFromFit(object):
def convergence(self) -> str:
raise NotImplementedError
+ @property
+ def constant_gev_params(self) -> GevParams:
+ raise NotImplementedError
+
class ResultFromIsmev(ResultFromFit):
@@ -76,7 +80,7 @@ class ResultFromIsmev(ResultFromFit):
return coef_dict
@property
- def stationary_gev_params(self) -> GevParams:
+ def constant_gev_params(self) -> GevParams:
params = {k.split('Coeff1')[0]: v for k, v in self.margin_coef_dict.items()
if 'Coeff1' in k and 'temp' not in k}
return GevParams.from_dict(params)
diff --git a/extreme_estimator/margin_fits/gev/gev_params.py b/extreme_estimator/margin_fits/gev/gev_params.py
index 2094bbeb..90f6f946 100644
--- a/extreme_estimator/margin_fits/gev/gev_params.py
+++ b/extreme_estimator/margin_fits/gev/gev_params.py
@@ -45,6 +45,24 @@ class GevParams(ExtremeParams):
def __str__(self):
return self.to_dict().__str__()
+ def quantile_strength_evolution_ratio(self, p=0.99, mu1=0.0, sigma1=0.0):
+ """
+ Compute the relative evolution of some quantile with respect to time.
+ (when mu1 and sigma1 can be modified with time)
+
+ :param p: level of the quantile
+ :param mu1: temporal slope of the location parameter
+ :param sigma1: temporal slope of the scale parameter
+ :return: A string summarizing the evolution ratio
+ """
+ initial_quantile = self.quantile(p)
+ quantity_increased = mu1
+ if sigma1 != 0:
+ quantity_increased += (sigma1 / self.shape) * (1 - (- np.float_power(np.log(p), -self.shape)))
+ return quantity_increased / initial_quantile
+
+ # Compute some indicators (such as the mean and the variance)
+
def g(self, k):
# Compute the g_k parameters as defined in wikipedia
# https://fr.wikipedia.org/wiki/Loi_d%27extremum_g%C3%A9n%C3%A9ralis%C3%A9e
--
GitLab