[PAPER 1] major refactor for drawing uncertainty

experiment/eurocode_data/main_eurocode_drawing.py

@@ -8,7 +8,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
 from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
     ConfidenceIntervalMethodFromExtremes
 from experiment.paper_past_snow_loads.result_trends_and_return_levels.eurocode_visualizer import \
-    plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict, get_model_name
+    plot_uncertainty_massifs, get_model_name
 from experiment.eurocode_data.massif_name_to_departement import MASSIF_NAMES_ALPS
 from experiment.eurocode_data.utils import EUROCODE_ALTITUDES
 from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days, CrocusSweTotal
@@ -73,7 +73,7 @@ def plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, tem
               model_name_to_massif_name_to_ordered_return_level.keys()}
         massif_name_to_model_name_to_ordered_return_level_uncertainties[massif_name] = d2
     # Plot graph
-    plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(
+    plot_uncertainty_massifs(
         massif_name_to_model_name_to_ordered_return_level_uncertainties, nb_massif_names=len(massif_names),
         nb_model_names=len(model_class_and_last_year))
     if show:

experiment/meteo_france_data/scm_models_data/abstract_study.py

@@ -225,6 +225,7 @@ class AbstractStudy(object):
 
     @property
     def study_massif_names(self) -> List[str]:
+        # Massif names that are present in the current study (i.e. for the current altitude)
         return self.altitude_to_massif_names[self.altitude]
 
     @property
@@ -490,6 +491,7 @@ class AbstractStudy(object):
         for massif_name in self.massif_name_to_altitudes.keys():
             for altitude in self.massif_name_to_altitudes[massif_name]:
                 altitude_to_massif_names[altitude].append(massif_name)
+        # massif_names are ordered in the same way as all_massif_names
         return altitude_to_massif_names
 
     """ Visualization methods """

experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py

@@ -3,6 +3,8 @@ from typing import Dict, List, Tuple
 import matplotlib.pyplot as plt
 import numpy as np
 
+from experiment.paper_past_snow_loads.result_trends_and_return_levels.study_visualizer_for_non_stationary_trends import \
+    StudyVisualizerForNonStationaryTrends
 from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
     EurocodeConfidenceIntervalFromExtremes
 from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
@@ -10,21 +12,6 @@ from experiment.meteo_france_data.scm_models_data.visualization.utils import cre
 from root_utils import get_display_name_from_object_type
 
 
-def get_label_name(model_name, ci_method_name: str):
-    is_non_stationary = model_name == 'NonStationary'
-    model_symbol = 'N' if is_non_stationary else '0'
-    parameter = ', 2017' if is_non_stationary else ''
-    model_name = ' $ \widehat{z_p}(\\boldsymbol{\\theta_{\mathcal{M}_'
-    model_name += model_symbol
-    model_name += '}}'
-    model_name += parameter
-    model_name += ')_{ \\textrm{' + ci_method_name.upper().split(' ')[1] + '}} $ '
-    return model_name
-
-
-def get_model_name(model_class):
-    return get_display_name_from_object_type(model_class).split('Stationary')[0] + 'Stationary'
-
 def massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, massif_names,
                                                       uncertainty_methods, temporal_covariate,
                                                       non_stationary_model):
@@ -33,7 +20,7 @@ def massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, ma
     for altitude, visualizer in altitude_to_visualizer.items():
         print('Processing altitude = {} '.format(altitude))
         for ci_method in uncertainty_methods:
-            d = visualizer.massif_name_to_altitude_and_eurocode_level_uncertainty_for_minimized_aic_model_class(
+            d = visualizer.massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(
                 massif_names, ci_method,
                 temporal_covariate, non_stationary_model)
             # Append the altitude one by one
@@ -45,78 +32,88 @@ def massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, ma
     return massif_name_to_ordered_eurocode_level_uncertainty
 
 
-def plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(altitude_to_visualizer,
-                                                                         massif_names,
-                                                                         non_stationary_models_for_uncertainty,
-                                                                         uncertainty_methods):
-    """
-    Rows correspond to massif names
-    Columns correspond to stationary/non stationary model name for a given date
-    Uncertainty result_trends_and_return_levels correpsond to the different plot on the graph
+def plot_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+    """ Plot several uncertainty plots
     :return:
     """
-    # Compute the dictionary of interest
-    # Plot uncertainties
-    model_name_to_massif_name_to_ordered_return_level = {}
-    for non_stationary_model in non_stationary_models_for_uncertainty:
-        d = massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, massif_names,
-                                                              uncertainty_methods,
-                                                              temporal_covariate=2017,
-                                                              non_stationary_model=non_stationary_model)
-        model_name_to_massif_name_to_ordered_return_level[non_stationary_model] = d
-
-    # Transform the dictionary into the desired format
-    d = {}
-    for massif_name in massif_names:
-        d2 = {model_name: model_name_to_massif_name_to_ordered_return_level[model_name][massif_name] for model_name
-              in
-              model_name_to_massif_name_to_ordered_return_level.keys()}
-        d[massif_name] = d2
-
+    visualizer = list(altitude_to_visualizer.values())[0]
+    # Subdivide massif names in group of 5
+    m = 4
+    uncertainty_massif_names = visualizer.uncertainty_massif_names
+    n = (len(uncertainty_massif_names) // m) + 1
+    for i in list(range(n))[:]:
+        massif_names = uncertainty_massif_names[m * i: m * (i + 1)]
+        print(massif_names)
+        plot_subgroup_uncertainty_massifs(altitude_to_visualizer, massif_names)
+
+
+def plot_subgroup_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
+                                      massif_names):
+    """Create a plot with a maximum of 4 massif names
+    We will save the plot at this level
+    """
+    visualizer = list(altitude_to_visualizer.values())[0]
     nb_massif_names = len(massif_names)
-    nb_model_names = len(non_stationary_models_for_uncertainty)
-    axes = create_adjusted_axes(nb_massif_names, nb_model_names)
+    assert nb_massif_names <= 5
+    axes = create_adjusted_axes(nb_massif_names, visualizer.nb_contexts)
     if nb_massif_names == 1:
         axes = [axes]
-    for ax, (massif_name, model_name_to_uncertainty_level) in zip(axes, d.items()):
-        plot_model_name_to_uncertainty_method_to_ordered_dict(model_name_to_uncertainty_level,
-                                                              massif_name, ax)
+    for ax, massif_name in zip(axes, massif_names):
+        plot_single_uncertainty_massif(altitude_to_visualizer,
+                                       massif_name, ax)
 
     # Save plot
-    visualizer = list(altitude_to_visualizer.values())[0]
     massif_names_str = '_'.join(massif_names)
-    model_names_str = 'NonStationarity=' + '_'.join([str(e) for e in non_stationary_models_for_uncertainty])
+    model_names_str = 'NonStationarity=' + '_'.join([str(e) for e in visualizer.non_stationary_contexts])
     visualizer.plot_name = model_names_str + '_' + massif_names_str
     visualizer.show_or_save_to_file(no_title=True)
 
-    # plt.suptitle('50-year return levels of extreme snow loads in France for several confiance interval methods.')
 
-def plot_model_name_to_uncertainty_method_to_ordered_dict(d, massif_name, axes):
-    if len(d) == 1:
+def plot_single_uncertainty_massif(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
+                                   massif_name, axes):
+    visualizer = list(altitude_to_visualizer.values())[0]
+    if visualizer.nb_contexts == 1:
         axes = [axes]
-    for ax, (model_name, uncertainty_method_to_ordered_dict) in zip(axes, d.items()):
-        plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name,
-                                                         uncertainty_method_to_ordered_dict)
+    for ax, non_stationary_context in zip(axes, visualizer.non_stationary_contexts):
+        plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, non_stationary_context,
+                                                                  altitude_to_visualizer)
 
 
-def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name,
-                                                     label_to_ordered_return_level_uncertainties:
-                                                     Dict[str, List[
-                                                         EurocodeConfidenceIntervalFromExtremes]]):
+def get_label_name(model_name, ci_method_name: str):
+    is_non_stationary = model_name == 'NonStationary'
+    model_symbol = 'N' if is_non_stationary else '0'
+    parameter = ', 2017' if is_non_stationary else ''
+    model_name = ' $ \widehat{z_p}(\\boldsymbol{\\theta_{\mathcal{M}_'
+    model_name += model_symbol
+    model_name += '}}'
+    model_name += parameter
+    model_name += ')_{ \\textrm{' + ci_method_name.upper().split(' ')[1] + '}} $ '
+    return model_name
+
+
+def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, non_stationary_context,
+                                                              altitude_to_visualizer: Dict[
+                                                                  int, StudyVisualizerForNonStationaryTrends]):
     """ Generic function that might be used by many other more global functions"""
+    altitudes = list(altitude_to_visualizer.keys())
+    visualizer = list(altitude_to_visualizer.values())[0]
     colors = ['tab:green', 'tab:olive']
     alpha = 0.2
     # Display the EUROCODE return level
     eurocode_region = massif_name_to_eurocode_region[massif_name]()
 
     # Display the return level from model class
-    for j, (color, (label, l)) in enumerate(zip(colors, label_to_ordered_return_level_uncertainties.items())):
-        l = list(zip(*l))
-        altitudes = l[0]
-        ordered_return_level_uncertaines = l[1]  # type: List[EurocodeConfidenceIntervalFromExtremes]
+    for j, (color, uncertainty_method) in enumerate(zip(colors, visualizer.uncertainty_methods)):
         # Plot eurocode standards only for the first loop
         if j == 0:
-            eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax, altitudes=altitudes)
+            eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax,
+                                                                                                   altitudes=altitudes)
+        # Compute ordered_return_level_uncertaines for a given massif_name, uncertainty methods, and non stationary context
+        ordered_return_level_uncertaines = []
+        for visualizer in altitude_to_visualizer.values():
+            u = visualizer.triplet_to_eurocode_uncertainty[(uncertainty_method, non_stationary_context, massif_name)]
+            ordered_return_level_uncertaines.append(u)
+        # Display
         mean = [r.mean_estimate for r in ordered_return_level_uncertaines]
         # Filter and keep only non nan values
         not_nan_index = [not np.isnan(m) for m in mean]
@@ -124,8 +121,9 @@ def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name
         altitudes = list(np.array(altitudes)[not_nan_index])
         ordered_return_level_uncertaines = list(np.array(ordered_return_level_uncertaines)[not_nan_index])
 
-        ci_method_name = str(label).split('.')[1].replace('_', ' ')
-        ax.plot(altitudes, mean, linestyle='--', marker='o', color=color, label=get_label_name(model_name, ci_method_name))
+        ci_method_name = str(uncertainty_method).split('.')[1].replace('_', ' ')
+        ax.plot(altitudes, mean, linestyle='--', marker='o', color=color,
+                label=get_label_name(non_stationary_context, ci_method_name))
         lower_bound = [r.confidence_interval[0] for r in ordered_return_level_uncertaines]
         upper_bound = [r.confidence_interval[1] for r in ordered_return_level_uncertaines]
         ax.fill_between(altitudes, lower_bound, upper_bound, color=color, alpha=alpha)
@@ -133,14 +131,15 @@ def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name
     ax.set_ylim([0.0, 16])
     massif_name_str = massif_name.replace('_', ' ')
     eurocode_region_str = get_display_name_from_object_type(type(eurocode_region))
-    is_non_stationary_model = model_name if isinstance(model_name, bool) else 'Non' in model_name
+    is_non_stationary_model = non_stationary_context if isinstance(non_stationary_context,
+                                                                   bool) else 'Non' in non_stationary_context
     if is_non_stationary_model:
-        model_name = 'non-stationary'
+        non_stationary_context = 'non-stationary'
     else:
-        model_name = 'stationary'
-    title = '{} ({} Eurocodes area) with a {} model'.format(massif_name_str, eurocode_region_str, model_name)
+        non_stationary_context = 'stationary'
+    title = '{} ({} Eurocodes area) with a {} model'.format(massif_name_str, eurocode_region_str,
+                                                            non_stationary_context)
     ax.set_title(title)
     ax.set_ylabel('50-year return level of SL (kN $m^-2$)')
     ax.set_xlabel('Altitude (m)')
     ax.grid()
-

experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -7,7 +7,7 @@ import matplotlib.pyplot as plt
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
     StudyVisualizer
 from experiment.paper_past_snow_loads.result_trends_and_return_levels.eurocode_visualizer import \
-    plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict
+    plot_uncertainty_massifs
 from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
 from experiment.paper_past_snow_loads.result_trends_and_return_levels.study_visualizer_for_non_stationary_trends import \
     StudyVisualizerForNonStationaryTrends
@@ -19,33 +19,57 @@ mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
 
-def draw_snow_load_map(altitude):
+def minor_result(altitude):
+    """Plot trends for a single altitude to be fast"""
     visualizer = StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude), multiprocessing=True)
     visualizer.plot_trends()
 
 
-def main_results():
-    altitudes = [[1500, 1800]][0]
-    uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
-                           ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
-    massif_names = ['Chartreuse']
-    non_stationary_models_for_uncertainty = [False, True][:1]
+def intermediate_result(altitudes, massif_names=None,
+                        non_stationary_uncertainty=None, uncertainty_methods=None,
+                        study_class=CrocusSnowLoadTotal):
+    """
+    Plot all the trends for all altitudes
+    And enable to plot uncertainty plot for some specific massif_names, uncertainty methods to be fast
+    :param altitudes:
+    :param massif_names:
+    :param non_stationary_uncertainty:
+    :param uncertainty_methods:
+    :param study_class:
+    :return:
+    """
     # Load altitude to visualizer
     altitude_to_visualizer = OrderedDict()
     for altitude in altitudes:
         altitude_to_visualizer[altitude] = StudyVisualizerForNonStationaryTrends(
-            study=CrocusSnowLoadTotal(altitude=altitude), multiprocessing=True, save_to_file=True)
+            study=study_class(altitude=altitude), multiprocessing=True, save_to_file=True,
+            uncertainty_massif_names=massif_names, uncertainty_methods=uncertainty_methods,
+            non_stationary_contexts=non_stationary_uncertainty)
     # Plot trends
     max_abs_tdrl = max([visualizer.max_abs_tdrl for visualizer in altitude_to_visualizer.values()])
     for visualizer in altitude_to_visualizer.values():
         visualizer.plot_trends(max_abs_tdrl)
     # Plot graph
-    plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(altitude_to_visualizer,
-                                                                         massif_names,
-                                                                         non_stationary_models_for_uncertainty,
-                                                                         uncertainty_methods)
+    plot_uncertainty_massifs(altitude_to_visualizer)
+    return altitude_to_visualizer
+
+
+def major_result():
+    altitudes = [[1500, 1800]][0]
+    uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+                           ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
+    massif_names = ['Chartreuse']
+    non_stationary_models_for_uncertainty = [False, True][:1]
+    #
+    # altitudes
+    # study_class
 
 
 if __name__ == '__main__':
-    # draw_snow_load_map(altitude=1800)
-    main_results()
+    # minor_result(altitude=1800)
+    intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
+                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
+                        non_stationary_uncertainty=[False])
+    intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
+                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
+                        non_stationary_uncertainty=[False])

experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py

@@ -15,6 +15,8 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter impo
     GevLocationTrendTest
 from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
 from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
+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 \
     compute_eurocode_confidence_interval, EurocodeConfidenceIntervalFromExtremes
 from root_utils import NB_CORES
@@ -26,11 +28,29 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                  vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False,
                  temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False,
                  complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True,
-                 score_class=MeanScore):
+                 score_class=MeanScore,
+                 uncertainty_methods=None,
+                 non_stationary_contexts=None,
+                 uncertainty_massif_names=None,
+                 effective_temporal_covariate=2017):
         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.effective_temporal_covariate = effective_temporal_covariate
+        self.non_stationary_contexts = non_stationary_contexts
+        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.uncertainty_methods is None:
+            self.uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+                                        ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
+        if self.uncertainty_massif_names is None:
+            self.uncertainty_massif_names = self.study.study_massif_names
+        # Assign default argument for the non stationary trends
         self.non_stationary_trend_test = [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest]
         self.non_stationary_trend_test_to_marker = dict(zip(self.non_stationary_trend_test, ["s", "^", "D"]))
 
@@ -103,7 +123,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                     'fillstyle': 'full' if t.is_significant else 'none'}
         return d
 
-    # Part 1 - Uncertainty return level plot
+    # Part 2 - Uncertainty return level plot
 
     def massif_name_to_model_class(self, massif_name, non_stationary_model):
         if not non_stationary_model:
@@ -111,23 +131,37 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         else:
             return self.massif_name_to_minimized_aic_non_stationary_trend_test[massif_name].unconstrained_model_class
 
-    def massif_name_to_uncertainty(self):
-        pass
+    @property
+    def nb_contexts(self):
+        return len(self.non_stationary_contexts)
+
+    @property
+    def nb_uncertainty_method(self):
+        return len(self.uncertainty_methods)
 
-    def massif_name_to_altitude_and_eurocode_level_uncertainty_for_minimized_aic_model_class(self, massif_names,
-                                                                                             ci_method,
-                                                                                             effective_temporal_covariate,
-                                                                                             non_stationary_model) \
+    def massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_model) \
             -> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]:
         arguments = [
             [self.massif_name_to_non_null_years_and_maxima[m],
              self.massif_name_to_model_class(m, non_stationary_model),
-             ci_method, effective_temporal_covariate] for m in massif_names]
+             ci_method, self.effective_temporal_covariate] for m in self.uncertainty_massif_names]
         if self.multiprocessing:
             with Pool(NB_CORES) as p:
                 res = p.starmap(compute_eurocode_confidence_interval, arguments)
         else:
             res = [compute_eurocode_confidence_interval(*argument) for argument in arguments]
-        altitudes_and_res = [(self.study.altitude, r) for r in res]
-        massif_name_to_eurocode_return_level_uncertainty = OrderedDict(zip(massif_names, altitudes_and_res))
+        massif_name_to_eurocode_return_level_uncertainty = OrderedDict(zip(self.uncertainty_massif_names, res))
         return massif_name_to_eurocode_return_level_uncertainty
+
+    @cached_property
+    def triplet_to_eurocode_uncertainty(self):
+        d = {}
+        for ci_method in self.uncertainty_methods:
+            for non_stationary_uncertainty in self.non_stationary_contexts:
+                for uncertainty_massif_name, eurocode_uncertainty in self.massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(
+                        ci_method, non_stationary_uncertainty).items():
+                    d[(ci_method, non_stationary_uncertainty, uncertainty_massif_name)] = eurocode_uncertainty
+        return d
+
+    def model_name_to_uncertainty_method_to_ratio_above_eurocode(self):
+        assert self.uncertainty_massif_names == self.study.study_massif_names