diff --git a/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/study_visualizer_for_shape_repartition.py b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/study_visualizer_for_shape_repartition.py
index 0a25e17896e9a627e57846f3185fe0e7fa8d18f1..bc0f04da20a9dcdd3ae3c06943d99fa8c4847ffb 100644
--- a/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/study_visualizer_for_shape_repartition.py
+++ b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/study_visualizer_for_shape_repartition.py
@@ -16,7 +16,7 @@ class StudyVisualizerForShape(StudyVisualizerForNonStationaryTrends):
@cached_property
def massif_name_to_unconstrained_shape_parameter(self):
return {m: t.unconstrained_estimator_gev_params.shape
- for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+ for m, t in self.massif_name_to_trend_test_that_minimized_aic.items()}
@cached_property
def massif_name_to_change_value(self):
diff --git a/experiment/paper_past_snow_loads/paper_main_utils.py b/experiment/paper_past_snow_loads/paper_main_utils.py
index 17de30ad67be6bc8dc68bd128df747e0912ce2b1..8ba05340dfb5c582d86881a8c5342ddef7d8373a 100644
--- a/experiment/paper_past_snow_loads/paper_main_utils.py
+++ b/experiment/paper_past_snow_loads/paper_main_utils.py
@@ -1,4 +1,5 @@
from collections import OrderedDict
+from enum import Enum
from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
@@ -6,14 +7,20 @@ from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear
TemporalMarginFitMethod
-def load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty, study_class, uncertainty_methods,
+def load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncertainty, study_class,
+ uncertainty_methods,
study_visualizer_class=StudyVisualizerForNonStationaryTrends,
save_to_file=True):
fit_method = TemporalMarginFitMethod.extremes_fevd_mle
+ select_only_acceptable_shape_parameter = True
+ print('Fit method: {}, Select shape parameter: {}'.format(fit_method, select_only_acceptable_shape_parameter))
altitude_to_visualizer = OrderedDict()
for altitude in altitudes:
altitude_to_visualizer[altitude] = study_visualizer_class(
study=study_class(altitude=altitude), multiprocessing=True, save_to_file=save_to_file,
uncertainty_massif_names=massif_names, uncertainty_methods=uncertainty_methods,
- non_stationary_contexts=non_stationary_uncertainty, fit_method=fit_method)
+ model_subsets_for_uncertainty=model_subsets_for_uncertainty, fit_method=fit_method,
+ select_only_acceptable_shape_parameter=select_only_acceptable_shape_parameter)
return altitude_to_visualizer
+
+
diff --git a/experiment/paper_past_snow_loads/paper_utils.py b/experiment/paper_past_snow_loads/paper_utils.py
index bfbeae1939fbc7424f6e60003ca5f3b77384dd6b..66dd821a91924de8628bd8a922a9d0b7e38ac1f9 100644
--- a/experiment/paper_past_snow_loads/paper_utils.py
+++ b/experiment/paper_past_snow_loads/paper_utils.py
@@ -1,8 +1,19 @@
+from enum import Enum
+
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoadEurocode, \
CrocusSnowLoad3Days
+from root_utils import get_display_name_from_object_type
paper_altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700]
paper_study_classes = [CrocusSnowLoadTotal, CrocusSnowLoadEurocode, CrocusSnowLoad3Days]
# dpi_paper1_figure = 700
dpi_paper1_figure = None
+class ModelSubsetForUncertainty(Enum):
+ stationary_gumbel = 0
+ stationary_gumbel_and_gev = 1
+ non_stationary_gumbel = 2
+ non_stationary_gumbel_and_gev = 3
+
+
+
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index ac139f9fef9a855e4ac4c633129b47aa7e1dcfbf..2f81491c503d4606ca54e78d8ee5c038fd2cf24c 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -5,7 +5,7 @@ import matplotlib as mpl
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoadEurocode, \
CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days
from experiment.paper_past_snow_loads.paper_main_utils import load_altitude_to_visualizer
-from experiment.paper_past_snow_loads.paper_utils import paper_study_classes, paper_altitudes
+from experiment.paper_past_snow_loads.paper_utils import paper_study_classes, paper_altitudes, ModelSubsetForUncertainty
from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import \
plot_uncertainty_massifs
from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_histogram import \
@@ -30,12 +30,12 @@ def minor_result(altitude):
def compute_minimized_aic(visualizer):
- _ = visualizer.massif_name_to_minimized_aic_non_stationary_trend_test
+ _ = visualizer.massif_name_to_trend_test_that_minimized_aic
return True
def intermediate_result(altitudes, massif_names=None,
- non_stationary_uncertainty=None, uncertainty_methods=None,
+ model_subsets_for_uncertainty=None, uncertainty_methods=None,
study_class=CrocusSnowLoadTotal,
multiprocessing=False):
"""
@@ -49,7 +49,7 @@ def intermediate_result(altitudes, massif_names=None,
:return:
"""
# Load altitude to visualizer
- altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty,
+ altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncertainty,
study_class, uncertainty_methods)
# Load variable object efficiently
for v in altitude_to_visualizer.values():
@@ -83,23 +83,25 @@ def major_result():
uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
massif_names = None
- study_classes = paper_study_classes[:1]
+ study_classes = paper_study_classes[:2]
+ model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
+ ModelSubsetForUncertainty.stationary_gumbel_and_gev,
+ ModelSubsetForUncertainty.non_stationary_gumbel,
+ ModelSubsetForUncertainty.non_stationary_gumbel_and_gev]
+ # model_subsets_for_uncertainty = None
# study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
for study_class in study_classes:
- if study_class == CrocusSnowLoadEurocode:
- non_stationary_uncertainty = [False]
- else:
- non_stationary_uncertainty = [False, True][:]
- intermediate_result(paper_altitudes, massif_names, non_stationary_uncertainty, uncertainty_methods, study_class)
+ intermediate_result(paper_altitudes, massif_names, model_subsets_for_uncertainty,
+ uncertainty_methods, study_class)
if __name__ == '__main__':
- # major_result()
- intermediate_result(altitudes=[900, 1200], massif_names=['Vercors'],
- uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
- ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
- non_stationary_uncertainty=[False, True][1:],
- multiprocessing=True)
+ major_result()
+ # intermediate_result(altitudes=[900, 1200], massif_names=['Vercors'],
+ # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+ # ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+ # non_stationary_uncertainty=[False, True][1:],
+ # multiprocessing=True)
# intermediate_result(altitudes=[900, 1200], massif_names=['Maurienne'],
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
# ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
index 9fbb0793d34728bab6678347c967f9cb50e555b9..3f274fdb691f9dd23811a688f808cc8fa8336907 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
@@ -7,7 +7,7 @@ from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_A
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
SCM_STUDY_CLASS_TO_ABBREVIATION
-from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure
+from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure, ModelSubsetForUncertainty
from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
@@ -55,21 +55,26 @@ def plot_single_uncertainty_massif(altitude_to_visualizer: Dict[int, StudyVisual
massif_name):
visualizer = list(altitude_to_visualizer.values())[0]
- for non_stationary_context in visualizer.non_stationary_contexts:
+ model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
+ ModelSubsetForUncertainty.non_stationary_gumbel_and_gev]
+ print('Subsets for uncertainty curves:{}'.format(model_subsets_for_uncertainty))
+ for model_subset_for_uncertainty in model_subsets_for_uncertainty:
ax = create_adjusted_axes(1, 1)
- plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, non_stationary_context,
+ plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, model_subset_for_uncertainty,
altitude_to_visualizer)
# Save plot
massif_names_str = massif_name
- model_names_str = 'NonStationarity={}'.format(non_stationary_context)
+ model_names_str = get_display_name_from_object_type(model_subset_for_uncertainty)
visualizer.plot_name = model_names_str + '_' + massif_names_str
visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
plt.close()
-def get_label_name(non_stationary_context, ci_method_name, add_method_suffix):
- model_symbol = 'N' if non_stationary_context else '0'
- parameter = ', 2017' if non_stationary_context else ''
+def get_label_name(model_subset_for_uncertainty, ci_method_name, add_method_suffix):
+ model_symbol = 'N' if model_subset_for_uncertainty is not ModelSubsetForUncertainty.stationary_gumbel else '0'
+ parameter = ', 2017' if model_subset_for_uncertainty not in [ModelSubsetForUncertainty.stationary_gumbel,
+ ModelSubsetForUncertainty.stationary_gumbel_and_gev] \
+ else ''
model_name = ' $ z_p(\\boldsymbol{\\widehat{\\theta}_{\\mathcal{M}'
# model_name += '_' + model_symbol
model_name += '}}'
@@ -81,7 +86,7 @@ def get_label_name(non_stationary_context, ci_method_name, add_method_suffix):
return model_name
-def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, non_stationary_context,
+def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, model_subset_for_uncertainty,
altitude_to_visualizer: Dict[
int, StudyVisualizerForNonStationaryTrends]):
""" Generic function that might be used by many other more global functions"""
@@ -104,12 +109,12 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
# Plot uncertainties
color = ci_method_to_color[uncertainty_method]
valid_altitudes = plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color,
- massif_name, non_stationary_context, uncertainty_method,
+ massif_name, model_subset_for_uncertainty, uncertainty_method,
nb_uncertainty_methods)
# Plot some data for the non valid altitudes
- # Plot bars of TDRL only in the non stationary case
- if j == 0 and non_stationary_context:
+ # Plot bars of TDRL only in the general non stationary case
+ if model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel_and_gev:
plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, legend_size, legend_size)
ax.legend(loc=2, prop={'size': legend_size})
@@ -174,11 +179,11 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, leg
def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
- non_stationary_context, uncertainty_method, nb_uncertainty_methods):
+ model_subset_for_uncertainty, uncertainty_method, nb_uncertainty_methods):
# Compute ordered_return_level_uncertaines for a given massif_name, uncertainty methods, and non stationary context
ordered_return_level_uncertainties = []
for visualizer in altitude_to_visualizer.values():
- u = visualizer.triplet_to_eurocode_uncertainty[(uncertainty_method, non_stationary_context, massif_name)]
+ u = visualizer.triplet_to_eurocode_uncertainty[(uncertainty_method, model_subset_for_uncertainty, massif_name)]
ordered_return_level_uncertainties.append(u)
# Display
mean = [r.mean_estimate for r in ordered_return_level_uncertainties]
@@ -189,7 +194,7 @@ def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitud
ordered_return_level_uncertainties = list(np.array(ordered_return_level_uncertainties)[not_nan_index])
ci_method_name = str(uncertainty_method).split('.')[1].replace('_', ' ')
add_method_suffix = nb_uncertainty_methods > 1 or 'mle' not in ci_method_name
- label_name = get_label_name(non_stationary_context, ci_method_name, add_method_suffix=add_method_suffix)
+ label_name = get_label_name(model_subset_for_uncertainty, ci_method_name, add_method_suffix=add_method_suffix)
ax.plot(valid_altitudes, mean, linestyle='--', marker='o', color=color,
label=label_name)
lower_bound = [r.confidence_interval[0] for r in ordered_return_level_uncertainties]
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
index 68e535ebb8d4856e78c65f783b38874dd8b7e84a..dd5ff29783a363a5767e66cfcc56859ffecd3a6f 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
@@ -8,6 +8,7 @@ from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends
StudyVisualizerForNonStationaryTrends
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color, \
ci_method_to_label, ConfidenceIntervalMethodFromExtremes
+from root_utils import get_display_name_from_object_type
def plot_uncertainty_histogram(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
@@ -16,15 +17,15 @@ def plot_uncertainty_histogram(altitude_to_visualizer: Dict[int, StudyVisualizer
"""
altitude_to_visualizer = {a: v for a, v in altitude_to_visualizer.items() if a in EUROCODE_ALTITUDES}
visualizer = list(altitude_to_visualizer.values())[0]
- for non_stationary_context in visualizer.non_stationary_contexts:
- plot_histogram(altitude_to_visualizer, non_stationary_context)
+ for model_subset_for_uncertainty in visualizer.model_subsets_for_uncertainty:
+ plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty)
-def plot_histogram(altitude_to_visualizer, non_stationary_context):
+def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
"""
Plot a single graph for potentially several confidence interval method
:param altitude_to_visualizer:
- :param non_stationary_context:
+ :param model_subset_for_uncertainty:
:return:
"""
visualizers = list(altitude_to_visualizer.values())
@@ -39,7 +40,7 @@ def plot_histogram(altitude_to_visualizer, non_stationary_context):
width = 100
else:
width = 200
- plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width=width)
+ plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width=width)
fontsize_label = 15
legend_size = 15
ax.set_xticks(altitudes)
@@ -51,14 +52,14 @@ def plot_histogram(altitude_to_visualizer, non_stationary_context):
ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
ax.set_ylim([0, 100])
ax.set_yticks([10 * i for i in range(11)])
- visualizer.plot_name = 'Percentages of exceedance with non_stationary={}'.format(non_stationary_context)
+ visualizer.plot_name = 'Percentages of exceedance with {}'.format(get_display_name_from_object_type(model_subset_for_uncertainty))
# visualizer.show = True
visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
ax.clear()
-def plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width):
- three_percentages_of_excess = [v.three_percentages_of_excess(ci_method, non_stationary_context) for v in
+def plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width):
+ three_percentages_of_excess = [v.three_percentages_of_excess(ci_method, model_subset_for_uncertainty) for v in
visualizers]
epsilon = 0.5
three_percentages_of_excess = [(a, b, c) if a == b else (max(epsilon, a), b, c) for (a, b, c) in three_percentages_of_excess]
diff --git a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
index 4db5de7c26a993d31c5724268ffb8e0d617e2a8d..7406dfffc37401b7574f56dcca46731f1a7a810f 100644
--- a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -1,9 +1,8 @@
from collections import OrderedDict
-
-import matplotlib.pyplot as plt
from multiprocessing.pool import Pool
-from typing import Dict, List
+from typing import Dict, List, Tuple
+import matplotlib.pyplot as plt
import numpy as np
from cached_property import cached_property
@@ -15,20 +14,18 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
StudyVisualizer
from experiment.paper_past_snow_loads.check_mcmc_convergence_for_return_levels.gelman_convergence_test import \
compute_gelman_convergence_value
+from experiment.paper_past_snow_loads.paper_utils import ModelSubsetForUncertainty
from experiment.trend_analysis.abstract_score import MeanScore
from experiment.trend_analysis.univariate_test.extreme_trend_test.abstract_gev_trend_test import AbstractGevTrendTest
-from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gev_trend_test_one_parameter import \
- GevLocationTrendTest, GevScaleTrendTest
from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter import \
GumbelLocationTrendTest, GevStationaryVersusGumbel, GumbelScaleTrendTest, GumbelVersusGumbel
from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_three_parameters.gev_trend_test_three_parameters import \
GevLocationAndScaleTrendTestAgainstGumbel
from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters import \
- GevLocationAndScaleTrendTest, GevLocationAgainstGumbel, GevScaleAgainstGumbel
+ GevLocationAgainstGumbel, GevScaleAgainstGumbel
from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gumbel_test_two_parameters import \
GumbelLocationAndScaleTrendTest
-from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
- GumbelTemporalModel
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import GumbelTemporalModel
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
ConfidenceIntervalMethodFromExtremes
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
@@ -44,27 +41,30 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True,
score_class=MeanScore,
uncertainty_methods=None,
- non_stationary_contexts=None,
+ model_subsets_for_uncertainty=None,
uncertainty_massif_names=None,
effective_temporal_covariate=2017,
relative_change_trend_plot=True,
non_stationary_trend_test_to_marker=None,
- fit_method=None):
+ fit_method=None,
+ select_only_acceptable_shape_parameter=False):
super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis,
normalization_under_one_observations, score_class)
# Add some attributes
+ self.select_only_acceptable_shape_parameter = select_only_acceptable_shape_parameter
self.fit_method = fit_method
self.non_stationary_trend_test_to_marker = non_stationary_trend_test_to_marker
self.relative_change_trend_plot = relative_change_trend_plot
self.effective_temporal_covariate = effective_temporal_covariate
- self.non_stationary_contexts = non_stationary_contexts
+ self.model_subsets_for_uncertainty = model_subsets_for_uncertainty
self.uncertainty_methods = uncertainty_methods
self.uncertainty_massif_names = uncertainty_massif_names
# Assign some default arguments
- if self.non_stationary_contexts is None:
- self.non_stationary_contexts = [False, True][:1]
+ if self.model_subsets_for_uncertainty is None:
+ self.model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
+ ModelSubsetForUncertainty.non_stationary_gumbel_and_gev][:]
if self.uncertainty_methods is None:
self.uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
@@ -73,10 +73,11 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
if self.non_stationary_trend_test_to_marker is None:
# Assign default argument for the non stationary trends
self.non_stationary_trend_test = [GumbelVersusGumbel,
- GumbelLocationTrendTest, GumbelScaleTrendTest, GumbelLocationAndScaleTrendTest,
+ GumbelLocationTrendTest, GumbelScaleTrendTest,
+ GumbelLocationAndScaleTrendTest,
GevStationaryVersusGumbel,
- GevLocationAgainstGumbel, GevScaleAgainstGumbel, GevLocationAndScaleTrendTestAgainstGumbel
- ]
+ GevLocationAgainstGumbel, GevScaleAgainstGumbel,
+ GevLocationAndScaleTrendTestAgainstGumbel]
self.non_stationary_trend_test_to_marker = {t: t.marker for t in self.non_stationary_trend_test}
else:
self.non_stationary_trend_test = list(self.non_stationary_trend_test_to_marker.keys())
@@ -112,20 +113,50 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
d[m] = (years[mask], maxima[mask])
return d
+ @property
+ def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
+ return self.massif_name_to_trend_test_tuple[0]
+
+ @property
+ def massif_name_to_stationary_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
+ return self.massif_name_to_trend_test_tuple[1]
+
+ @property
+ def massif_name_to_gumbel_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
+ return self.massif_name_to_trend_test_tuple[2]
+
@cached_property
- def massif_name_to_minimized_aic_non_stationary_trend_test(self) -> Dict[str, AbstractGevTrendTest]:
+ def massif_name_to_trend_test_tuple(self) -> Tuple[
+ Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest]]:
starting_year = None
massif_name_to_trend_test_that_minimized_aic = {}
+ massif_name_to_stationary_trend_test_that_minimized_aic = {}
+ massif_name_to_gumbel_trend_test_that_minimized_aic = {}
for massif_name, (x, y) in self.massif_name_to_non_null_years_and_maxima.items():
quantile_level = self.massif_name_to_eurocode_quantile_level_in_practice[massif_name]
- non_stationary_trend_test = [
+ all_trend_test = [
t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level,
fit_method=self.fit_method)
for t in self.non_stationary_trend_test] # type: List[AbstractGevTrendTest]
- # Extract the model with minimized AIC
- trend_test_that_minimized_aic = sorted(non_stationary_trend_test, key=lambda t: t.aic)[0]
+ # Exclude GEV models whose shape parameter is not in the support of the prior distribution for GMLE
+ if self.select_only_acceptable_shape_parameter:
+ acceptable_shape_parameter = lambda s: -0.5 <= s <= 0.5 # physically acceptable prior
+ all_trend_test = [t for t in all_trend_test
+ if acceptable_shape_parameter(t.unconstrained_estimator_gev_params.shape)]
+ sorted_trend_test = sorted(all_trend_test, key=lambda t: t.aic)
+ # Extract the stationary or non-stationary model that minimized AIC
+ trend_test_that_minimized_aic = sorted_trend_test[0]
massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
- return massif_name_to_trend_test_that_minimized_aic
+ # Extract the stationary model that minimized AIC
+ stationary_trend_test_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
+ [GumbelVersusGumbel, GevStationaryVersusGumbel]][0]
+ massif_name_to_stationary_trend_test_that_minimized_aic[massif_name] = stationary_trend_test_that_minimized_aic
+ # Extract the Gumbel model that minimized AIC
+ gumbel_trend_test_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
+ [GumbelVersusGumbel, GumbelLocationTrendTest, GumbelScaleTrendTest,
+ GumbelLocationAndScaleTrendTest]][0]
+ massif_name_to_gumbel_trend_test_that_minimized_aic[massif_name] = gumbel_trend_test_that_minimized_aic
+ return massif_name_to_trend_test_that_minimized_aic, massif_name_to_stationary_trend_test_that_minimized_aic, massif_name_to_gumbel_trend_test_that_minimized_aic
# Part 1 - Trends
@@ -140,7 +171,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
max_abs_change = 1e-10
return max_abs_change
- def plot_trends(self, max_abs_tdrl=None, add_colorbar=True):
+ def plot_trends(self, max_abs_tdrl=None, add_colorbar=True):
if max_abs_tdrl is not None:
self.global_max_abs_change = max_abs_tdrl
ax = self.study.visualize_study(massif_name_to_value=self.massif_name_to_change_value,
@@ -208,12 +239,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
@cached_property
def massif_name_to_tdrl_value(self):
return {m: t.time_derivative_of_return_level for m, t in
- self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+ self.massif_name_to_trend_test_that_minimized_aic.items()}
@cached_property
def massif_name_to_relative_change_value(self):
return {m: t.relative_change_in_return_level(initial_year=self.initial_year, final_year=self.final_year)
- for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+ for m, t in self.massif_name_to_trend_test_that_minimized_aic.items()}
@property
def initial_year(self):
@@ -242,7 +273,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
@cached_property
def massif_name_to_marker_style(self):
d = {}
- for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items():
+ for m, t in self.massif_name_to_trend_test_that_minimized_aic.items():
d[m] = {'marker': self.non_stationary_trend_test_to_marker[type(t)],
'color': 'k',
'markersize': 7,
@@ -251,22 +282,24 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
# Part 2 - Uncertainty return level plot
- def massif_name_and_non_stationary_context_to_model_class(self, massif_name, non_stationary_context):
- if not non_stationary_context:
+ def massif_name_and_model_subset_to_model_class(self, massif_name, model_subset_for_uncertainty):
+ if model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gumbel:
return GumbelTemporalModel
+ elif model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gumbel_and_gev:
+ return self.massif_name_to_stationary_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
+ elif model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel:
+ return self.massif_name_to_gumbel_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
+ elif model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel_and_gev:
+ return self.massif_name_to_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
else:
- return self.massif_name_to_minimized_aic_non_stationary_trend_test[massif_name].unconstrained_model_class
-
- @property
- def nb_contexts(self):
- return len(self.non_stationary_contexts)
+ raise ValueError(model_subset_for_uncertainty)
- def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_context) \
+ def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, model_subset_for_uncertainty) \
-> Dict[str, EurocodeConfidenceIntervalFromExtremes]:
# Compute for the uncertainty massif names
arguments = [
[self.massif_name_to_non_null_years_and_maxima[m],
- self.massif_name_and_non_stationary_context_to_model_class(m, non_stationary_context),
+ self.massif_name_and_model_subset_to_model_class(m, model_subset_for_uncertainty),
ci_method, self.effective_temporal_covariate,
self.massif_name_to_eurocode_quantile_level_in_practice[m]
] for m in self.uncertainty_massif_names]
@@ -287,13 +320,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
@cached_property
def triplet_to_eurocode_uncertainty(self):
- # -> Dict[(str, bool, str), EurocodeConfidenceIntervalFromExtremes]
d = {}
for ci_method in self.uncertainty_methods:
- for non_stationary_uncertainty in self.non_stationary_contexts:
+ for model_subset_for_uncertainty in self.model_subsets_for_uncertainty:
for massif_name, eurocode_uncertainty in self.all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(
- ci_method, non_stationary_uncertainty).items():
- d[(ci_method, non_stationary_uncertainty, massif_name)] = eurocode_uncertainty
+ ci_method, model_subset_for_uncertainty).items():
+ d[(ci_method, model_subset_for_uncertainty, massif_name)] = eurocode_uncertainty
return d
def model_name_to_uncertainty_method_to_ratio_above_eurocode(self):
@@ -319,18 +351,13 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
return {m: r().valeur_caracteristique(altitude=self.study.altitude)
for m, r in massif_name_to_eurocode_region.items() if m in self.uncertainty_massif_names}
- def three_percentages_of_excess(self, ci_method, non_stationary_context):
+ def three_percentages_of_excess(self, ci_method, model_subset_for_uncertainty):
eurocode_and_uncertainties = [(self.massif_name_to_eurocode_values[massif_name],
self.triplet_to_eurocode_uncertainty[
- (ci_method, non_stationary_context, massif_name)])
+ (ci_method, model_subset_for_uncertainty, massif_name)])
for massif_name in self.uncertainty_massif_names]
a = np.array([(uncertainty.confidence_interval[0] > eurocode,
uncertainty.mean_estimate > eurocode,
uncertainty.confidence_interval[1] > eurocode)
for eurocode, uncertainty in eurocode_and_uncertainties])
return 100 * np.mean(a, axis=0)
-
-
-
-
-