[paper 1] refactor study visualizer using "model subset for uncertainty" name...

[paper 1] refactor study visualizer using "model subset for uncertainty" name instead of unclear names.

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):

experiment/paper_past_snow_loads/paper_main_utils.py

 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
+
+

experiment/paper_past_snow_loads/paper_utils.py

+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
+
+
+

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:],

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]

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]

experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py

 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)
-
-
-
-
-