[PAPER 1] add histogram plot

experiment/eurocode_data/main_eurocode_drawing.py

@@ -7,7 +7,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
     StudyVisualizer
 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 \
+from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import \
     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

experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -5,9 +5,11 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
     ALL_ALTITUDES_WITHOUT_NAN
 from experiment.paper_past_snow_loads.paper_utils import paper_study_classes, paper_altitudes, \
     load_altitude_to_visualizer
-from experiment.paper_past_snow_loads.result_trends_and_return_levels.eurocode_visualizer import \
+from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import \
     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.plot_uncertainty_histogram import \
+    plot_uncertainty_histogram
 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.confidence_interval_method import \
@@ -47,8 +49,8 @@ def intermediate_result(altitudes, massif_names=None,
         visualizer.plot_trends(max_abs_tdrl)
     # Plot graph
     plot_uncertainty_massifs(altitude_to_visualizer)
-    return altitude_to_visualizer
-
+    # Plot histogram
+    plot_uncertainty_histogram(altitude_to_visualizer)
 
 
 def major_result():
@@ -62,8 +64,12 @@ def major_result():
 
 if __name__ == '__main__':
     # major_result()
+    intermediate_result(altitudes=[900, 1200], massif_names=['Chartreuse', 'Vanoise', 'Mercantour'],
+                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+                           ConfidenceIntervalMethodFromExtremes.ci_mle])
+    # intermediate_result(altitudes=[900, 1200], massif_names=None)
     # intermediate_result(ALL_ALTITUDES_WITHOUT_NAN)
-    intermediate_result(paper_altitudes)
+    # intermediate_result(paper_altitudes)
     # minor_result(altitude=600)
     # intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,

experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py → experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py

 from typing import Dict
+import matplotlib.pyplot as plt
 
 import numpy as np
 
@@ -9,6 +10,7 @@ from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extra
     AbstractExtractEurocodeReturnLevel
 from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
 from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color
 from root_utils import get_display_name_from_object_type
 
 
@@ -16,7 +18,7 @@ def plot_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerFo
     """ Plot several uncertainty plots
     :return:
     """
-    altitude_to_visualizer = {a:v for a,v in altitude_to_visualizer.items() if a in EUROCODE_ALTITUDES}
+    altitude_to_visualizer = {a: v for a, v in altitude_to_visualizer.items() if a in EUROCODE_ALTITUDES}
     visualizer = list(altitude_to_visualizer.values())[-1]
     # Subdivide massif names in group of 3
     m = 1
@@ -49,6 +51,7 @@ def plot_subgroup_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVis
     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.close()
 
 
 def plot_single_uncertainty_massif(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
@@ -78,19 +81,19 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
     """ 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:brown'][::-1]
     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, uncertainty_method) in enumerate(zip(colors, visualizer.uncertainty_methods)):
+    for j, uncertainty_method in enumerate(visualizer.uncertainty_methods):
         if j == 0:
             # Plot eurocode norm
             eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax,
                                                                                                    altitudes=altitudes)
 
         # 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)
 
@@ -118,7 +121,8 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
 
 
 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]
+    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
@@ -155,5 +159,6 @@ def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitud
     upper_bound = [r.confidence_interval[1] for r in ordered_return_level_uncertainties]
     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)
+    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/plot_uncertainty_histogram.py

+from typing import Dict
+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.study_visualizer_for_non_stationary_trends import \
+    StudyVisualizerForNonStationaryTrends
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color
+
+
+def plot_uncertainty_histogram(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+    """ Plot one graph for each non-stationary context
+    :return:
+    """
+    visualizer = list(altitude_to_visualizer.values())[0]
+    for non_stationary_context in visualizer.non_stationary_contexts:
+        plot_histogram(altitude_to_visualizer, non_stationary_context)
+
+
+def plot_histogram(altitude_to_visualizer, non_stationary_context):
+    """
+    Plot a single graph for potentially several confidence interval method
+    :param altitude_to_visualizer:
+    :param non_stationary_context:
+    :return:
+    """
+    visualizers = list(altitude_to_visualizer.values())
+    visualizer = visualizers[0]
+    ax = plt.gca()
+    altitudes = np.array(list(altitude_to_visualizer.keys()))
+    bincenters = altitudes
+    for j, ci_method in enumerate(visualizer.uncertainty_methods):
+        if len(visualizer.uncertainty_methods) == 2:
+            offset = -50 if j == 0 else 50
+            bincenters = altitudes + offset
+
+        width = 100
+        plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width=width)
+    ax.set_xticks(altitudes)
+    ax.set_ylabel('Massifs exceeding French standards (\%)')
+    ax.set_xlabel('Altitude (m)')
+    ax.set_ylim([0, 100])
+    ax.set_yticks([10 * i for i in range(11)])
+    # visualizer.plot_name = 'Percentage '
+    # visualizer.show = True
+    # visualizer.show_or_save_to_file(no_title=True)
+    plt.show()
+
+
+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
+                                   visualizers]
+    y = [d[1] for d in three_percentages_of_excess]
+    yerr = np.array([[d[1] - d[0], d[2] - d[1]] for d in three_percentages_of_excess]).transpose()
+    print(y, three_percentages_of_excess)
+    ax.bar(bincenters, y, width=width, color=ci_method_to_color[ci_method], yerr=yerr)

experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -6,6 +6,7 @@ from typing import Dict, Tuple
 import numpy as np
 from cached_property import cached_property
 
+from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
 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
@@ -117,15 +118,15 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         if max_abs_tdrl is not None:
             self.global_max_abs_tdrl = max_abs_tdrl
         ax = 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_color=self.massif_name_to_tdrl_color,
-                                   cmap=self.cmap,
-                                   show=False,
-                                   ticks_values_and_labels=self.ticks_values_and_labels,
-                                   label=self.label_tdrl_bar + ' for {}'.format(EUROCODE_RETURN_LEVEL_STR))
+                                        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_color=self.massif_name_to_tdrl_color,
+                                        cmap=self.cmap,
+                                        show=False,
+                                        ticks_values_and_labels=self.ticks_values_and_labels,
+                                        label=self.label_tdrl_bar + ' for {}'.format(EUROCODE_RETURN_LEVEL_STR))
         ax.get_xaxis().set_visible(True)
         ax.set_xticks([])
         ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=12)
@@ -187,7 +188,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         return len(self.non_stationary_contexts)
 
     def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_context) \
-            -> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]:
+            -> Dict[str, EurocodeConfidenceIntervalFromExtremes]:
         # Compute for the uncertainty massif names
         arguments = [
             [self.massif_name_to_non_null_years_and_maxima[m],
@@ -212,6 +213,7 @@ 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:
@@ -234,3 +236,22 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         else:
             res = [compute_gelman_convergence_value(*argument) for argument in arguments]
         return dict(zip(self.uncertainty_massif_names, res))
+
+    # Some values for the histogram
+
+    @cached_property
+    def massif_name_to_eurocode_values(self):
+        """Eurocode values for the altitude"""
+        return {m: r().lois_de_variation_de_la_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):
+        eurocode_and_uncertainties = [(self.massif_name_to_eurocode_values[massif_name],
+                                       self.triplet_to_eurocode_uncertainty[
+                                           (ci_method, non_stationary_context, 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)

extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py

@@ -16,3 +16,8 @@ ci_method_to_method_name = {
     ConfidenceIntervalMethodFromExtremes.ci_boot: 'boot',
     ConfidenceIntervalMethodFromExtremes.ci_proflik: 'proflik',
 }
+
+ci_method_to_color = {
+    ConfidenceIntervalMethodFromExtremes.ci_mle: 'tab:brown',
+    ConfidenceIntervalMethodFromExtremes.my_bayes: 'tab:green'
+}
\ No newline at end of file