[PAPER 1] fix ticks for colorbar. improve labels and cases where there is some missing labels

experiment/eurocode_data/utils.py

 
+# Eurocode quantile str
+EUROCODE_RETURN_LEVEL_STR = '50-year return level of SL (kN $m^-2$)'
 # Eurocode quantile correspond to a 50 year return period
 EUROCODE_QUANTILE = 0.98
 # Altitudes (between low and mid altitudes) < 2000m and should be > 200m

experiment/paper_past_snow_loads/data/main_eurocode_plot.py

@@ -3,10 +3,13 @@ import numpy as np
 
 from experiment.eurocode_data.eurocode_region import C2, E, C1
 from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
+from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR
 from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from root_utils import get_display_name_from_object_type
 
 
+
+
 def main_eurocode_norms(ax=None):
     if ax is None:
         ax = plt.gca()
@@ -20,7 +23,7 @@ def main_eurocode_norms(ax=None):
                 ax.legend()
                 ax.xaxis.set_ticks([250 * i for i in range(1, 9)])
                 ax.tick_params(axis='both', which='major', labelsize=13)
-                ax.set_ylabel('50-year return level of SL (kN $m^-2$)')
+                ax.set_ylabel(EUROCODE_RETURN_LEVEL_STR)
                 ax.set_xlabel('Altitude (m)')
                 ax.set_ylim([0.0, 11.0])
                 ax.grid()

experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py

@@ -3,8 +3,11 @@ from typing import Dict, List, Tuple
 import matplotlib.pyplot as plt
 import numpy as np
 
+from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR
 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.abstract_extract_eurocode_return_level import \
+    AbstractExtractEurocodeReturnLevel
 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
@@ -12,31 +15,11 @@ from experiment.meteo_france_data.scm_models_data.visualization.utils import cre
 from root_utils import get_display_name_from_object_type
 
 
-def massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, massif_names,
-                                                      uncertainty_methods, temporal_covariate,
-                                                      non_stationary_model):
-    massif_name_to_ordered_eurocode_level_uncertainty = {
-        massif_name: {ci_method: [] for ci_method in uncertainty_methods} for massif_name in massif_names}
-    for altitude, visualizer in altitude_to_visualizer.items():
-        print('Processing altitude = {} '.format(altitude))
-        for ci_method in uncertainty_methods:
-            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
-            for massif_name, return_level_uncertainty in d.items():
-                print(massif_name, return_level_uncertainty[0], return_level_uncertainty[1].confidence_interval,
-                      return_level_uncertainty[1].mean_estimate)
-                massif_name_to_ordered_eurocode_level_uncertainty[massif_name][ci_method].append(
-                    return_level_uncertainty)
-    return massif_name_to_ordered_eurocode_level_uncertainty
-
-
 def plot_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
     """ Plot several uncertainty plots
     :return:
     """
-    visualizer = list(altitude_to_visualizer.values())[0]
+    visualizer = list(altitude_to_visualizer.values())[-1]
     # Subdivide massif names in group of 3
     m = 1
     uncertainty_massif_names = visualizer.uncertainty_massif_names
@@ -80,10 +63,9 @@ def plot_single_uncertainty_massif(altitude_to_visualizer: Dict[int, StudyVisual
                                                                   altitude_to_visualizer)
 
 
-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 ''
+def get_label_name(non_stationary_context, ci_method_name):
+    model_symbol = 'N' if non_stationary_context else '0'
+    parameter = ', 2017' if non_stationary_context else ''
     model_name = ' $ \widehat{z_p}(\\boldsymbol{\\theta_{\mathcal{M}_'
     model_name += model_symbol
     model_name += '}}'
@@ -105,26 +87,18 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
 
     # Display the return level from model class
     for j, (color, uncertainty_method) in enumerate(zip(colors, visualizer.uncertainty_methods)):
-        # Plot eurocode standards only for the first loop
         if j == 0:
             # Plot eurocode norm
             eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax,
                                                                                                    altitudes=altitudes)
-            # Plot bars of TDRL only in the non stationary case
-            if non_stationary_context:
-                tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in altitude_to_visualizer.values()]
-                # Plot bars
-                colors = [v.massif_name_to_tdrl_color[massif_name] for v in altitude_to_visualizer.values()]
-                ax.bar(altitudes, tdrl_values, width=150, color=colors)
-                # Plot markers
-                markers_kwargs = [v.massif_name_to_marker_style(markersize=4)[massif_name]
-                                  for v, tdrl_value in zip(altitude_to_visualizer.values(), tdrl_values)]
-                for altitude, marker_kwargs, value in zip(altitudes, markers_kwargs, tdrl_values):
-                    ax.plot([altitude], [value / 2], **marker_kwargs)
 
         # Plot uncertainties
-        plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
-                                              non_stationary_context, uncertainty_method)
+        valid_altitudes = plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color,
+                                                                massif_name, non_stationary_context, uncertainty_method)
+
+        # Plot bars of TDRL only in the non stationary case
+        if j == 0 and non_stationary_context:
+            plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes)
 
     ax.legend(loc=2)
     ax.set_ylim([-1, 16])
@@ -140,11 +114,27 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
                                                             non_stationary_context)
     ax.set_title(title)
     ax.set_xticks(altitudes)
-    ax.set_ylabel('50-year return level of SL (kN $m^-2$)')
+    ax.set_ylabel(EUROCODE_RETURN_LEVEL_STR)
     ax.set_xlabel('Altitude (m)')
     ax.grid()
 
 
+def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes):
+    visualizers = [v for a, v in altitude_to_visualizer.items() if a in valid_altitudes and massif_name in v.uncertainty_massif_names]
+    if len(visualizers) > 0:
+        tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in visualizers]
+        # Plot bars
+        colors = [v.massif_name_to_tdrl_color[massif_name] for v in visualizers]
+        ax.bar(valid_altitudes, tdrl_values, width=150, color=colors, label=visualizers[0].label_tdrl_bar,
+               edgecolor='black', hatch='//')
+        # Plot markers
+        markers_kwargs = [v.massif_name_to_marker_style[massif_name] for v in visualizers]
+        for altitude, marker_kwargs, value in zip(valid_altitudes, markers_kwargs, tdrl_values):
+            # ax.plot([altitude], [value / 2], **marker_kwargs)
+            # Better to plot all the markers on the same line
+            ax.plot([altitude], 0, **marker_kwargs)
+
+
 def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
                                           non_stationary_context, uncertainty_method):
     # Compute ordered_return_level_uncertaines for a given massif_name, uncertainty methods, and non stationary context
@@ -160,8 +150,12 @@ def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitud
     valid_altitudes = list(np.array(altitudes)[not_nan_index])
     ordered_return_level_uncertainties = list(np.array(ordered_return_level_uncertainties)[not_nan_index])
     ci_method_name = str(uncertainty_method).split('.')[1].replace('_', ' ')
+    label_name = get_label_name(non_stationary_context, ci_method_name)
     ax.plot(valid_altitudes, mean, linestyle='--', marker='o', color=color,
-            label=get_label_name(non_stationary_context, ci_method_name))
+            label=label_name)
     lower_bound = [r.confidence_interval[0] for r in ordered_return_level_uncertainties]
     upper_bound = [r.confidence_interval[1] for r in ordered_return_level_uncertainties]
-    ax.fill_between(valid_altitudes, lower_bound, upper_bound, color=color, alpha=alpha)
+    confidence_interval_str = ' {}'.format(AbstractExtractEurocodeReturnLevel.percentage_confidence_interval)
+    confidence_interval_str += '\% confidence interval'
+    ax.fill_between(valid_altitudes, lower_bound, upper_bound, color=color, alpha=alpha, label=label_name + confidence_interval_str)
+    return valid_altitudes

experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -79,7 +79,10 @@ if __name__ == '__main__':
     # intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
     #                     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,
-                                             ConfidenceIntervalMethodFromExtremes.ci_bayes],
+    # intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
+    #                                          ConfidenceIntervalMethodFromExtremes.ci_bayes],
+    #                     non_stationary_uncertainty=[False, True])
+    intermediate_result(altitudes=[300, 600, 900], massif_names=None,
+                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
                         non_stationary_uncertainty=[False, True])

experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py

@@ -6,7 +6,7 @@ from typing import Dict, Tuple
 import numpy as np
 from cached_property import cached_property
 
-from experiment.eurocode_data.utils import EUROCODE_QUANTILE
+from experiment.eurocode_data.utils import EUROCODE_QUANTILE, EUROCODE_RETURN_LEVEL_STR
 from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map, get_colors
 from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
@@ -117,19 +117,19 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         self.study.visualize_study(massif_name_to_value=self.massif_name_to_tdrl_value,
                                    replace_blue_by_white=False, axis_off=False, show_label=False,
                                    add_colorbar=True,
-                                   massif_name_to_marker_style=self.massif_name_to_marker_style(),
+                                   massif_name_to_marker_style=self.massif_name_to_marker_style,
                                    massif_name_to_color=self.massif_name_to_tdrl_color,
                                    cmap=self.cmap,
                                    show=self.show,
                                    ticks_values_and_labels=self.ticks_values_and_labels,
-                                   label=self.label_trend_bar)
+                                   label=self.label_tdrl_bar + ' for {}'.format(EUROCODE_RETURN_LEVEL_STR))
         self.plot_name = 'tdlr_trends'
         self.show_or_save_to_file(add_classic_title=False, tight_layout=False, no_title=True)
         plt.close()
 
     @property
-    def label_trend_bar(self):
-        return 'Change in {} years for Eurocode quantile of SL (kN $m^-2$)'.format(AbstractGevTrendTest.nb_years_for_quantile_evolution)
+    def label_tdrl_bar(self):
+        return 'Change in {} years'.format(AbstractGevTrendTest.nb_years_for_quantile_evolution)
 
     @property
     def ticks_values_and_labels(self):
@@ -140,7 +140,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
             positive_ticks.append(round(tick, 1))
             tick += graduation
         all_ticks_labels = [-t for t in positive_ticks] + [0] + positive_ticks
-        ticks_values = [(t / 2 * self._max_abs_tdrl) + 0.5 for t in all_ticks_labels]
+        ticks_values = [((t / self._max_abs_tdrl) + 1) / 2 for t in all_ticks_labels]
         return ticks_values, all_ticks_labels
 
     @cached_property
@@ -157,12 +157,13 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         return {m: get_colors([v], self.cmap, -self._max_abs_tdrl, self._max_abs_tdrl)[0]
                 for m, v in self.massif_name_to_tdrl_value.items()}
 
-    def massif_name_to_marker_style(self, markersize=5):
+    @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():
             d[m] = {'marker': self.non_stationary_trend_test_to_marker[type(t)],
                     'color': 'k',
-                    'markersize': markersize,
+                    'markersize': 5,
                     'fillstyle': 'full' if t.is_significant else 'none'}
         return d
 
@@ -182,8 +183,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     def nb_uncertainty_method(self):
         return len(self.uncertainty_methods)
 
-    def massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_model) \
+    def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_model) \
             -> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]:
+        # Compute for the uncertainty massif names
         arguments = [
             [self.massif_name_to_non_null_years_and_maxima[m],
              self.massif_name_to_model_class(m, non_stationary_model),
@@ -195,17 +197,24 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                 res = p.starmap(compute_eurocode_confidence_interval, arguments)
         else:
             res = [compute_eurocode_confidence_interval(*argument) for argument in arguments]
-        massif_name_to_eurocode_return_level_uncertainty = OrderedDict(zip(self.uncertainty_massif_names, res))
+        massif_name_to_eurocode_return_level_uncertainty = dict(zip(self.uncertainty_massif_names, res))
+        # For the rest of the massif names. Create a Eurocode Return Level Uncertainty as nan
+        for massif_name in set(self.study.all_massif_names) - set(self.uncertainty_massif_names):
+            massif_name_to_eurocode_return_level_uncertainty[massif_name] = self.default_eurocode_uncertainty
         return massif_name_to_eurocode_return_level_uncertainty
 
+    @cached_property
+    def default_eurocode_uncertainty(self):
+        return EurocodeConfidenceIntervalFromExtremes(mean_estimate=np.nan, confidence_interval=(np.nan, np.nan))
+
     @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(
+                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, uncertainty_massif_name)] = eurocode_uncertainty
+                    d[(ci_method, non_stationary_uncertainty, massif_name)] = eurocode_uncertainty
         return d
 
     def model_name_to_uncertainty_method_to_ratio_above_eurocode(self):

extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py

@@ -11,6 +11,7 @@ from extreme_fit.model.margin_model.margin_function.linear_margin_function impor
 from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_bayesian_extremes import \
     ResultFromBayesianExtremes
 from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_mle_extremes import ResultFromMleExtremes
+from root_utils import classproperty
 
 
 class AbstractExtractEurocodeReturnLevel(object):
@@ -22,8 +23,11 @@ class AbstractExtractEurocodeReturnLevel(object):
         self.result_from_fit = self.estimator.result_from_model_fit
         self.temporal_covariate = temporal_covariate
         # Fixed Parameters
-        self.eurocode_quantile_level = quantile_level
-        self.alpha_for_confidence_interval = self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY
+        self.quantile_level = quantile_level
+
+    @classproperty
+    def percentage_confidence_interval(cls) -> int:
+        return int(100 * (1 - cls.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY))
 
     @property
     def mean_estimate(self) -> np.ndarray:
@@ -43,8 +47,8 @@ class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel)
 
     @cached_property
     def confidence_interval_method(self):
-        return self.result_from_fit.confidence_interval_method(self.eurocode_quantile_level,
-                                                               self.alpha_for_confidence_interval,
+        return self.result_from_fit.confidence_interval_method(self.quantile_level,
+                                                               self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY,
                                                                self.transformed_temporal_covariate,
                                                                self.ci_method)
 
@@ -81,7 +85,7 @@ class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeRe
     @cached_property
     def posterior_eurocode_return_level_samples_for_temporal_covariate(self) -> np.ndarray:
         return np.array(
-            [p.quantile(self.eurocode_quantile_level) for p in self.gev_params_from_fit_for_temporal_covariate])
+            [p.quantile(self.quantile_level) for p in self.gev_params_from_fit_for_temporal_covariate])
 
     @property
     def mean_estimate(self) -> np.ndarray:
@@ -91,7 +95,7 @@ class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeRe
     @property
     def confidence_interval(self):
         # Bottom and upper quantile correspond to the quantile
-        bottom_quantile = self.alpha_for_confidence_interval / 2
+        bottom_quantile = self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY / 2
         bottom_and_upper_quantile = (bottom_quantile, 1 - bottom_quantile)
         return [np.quantile(self.posterior_eurocode_return_level_samples_for_temporal_covariate, q=q)
                 for q in bottom_and_upper_quantile]