[PAPER 1] improve all uncertainty plots

experiment/eurocode_data/eurocode_region.py

@@ -44,7 +44,9 @@ class AbstractEurocodeRegion(object):
     def lois_de_variation_1000_and_2000(self):
         return 3.5, -2.45
 
-    def plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(self, ax, altitudes, label='Eurocode standards', linestyle=None):
+    def plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(self, ax, altitudes,
+                                                                               label='French standards',
+                                                                               linestyle=None):
         ax.plot(altitudes, [self.valeur_caracteristique(altitude) for altitude in altitudes],
                 label=label, color=self.eurocode_color, linewidth=5, linestyle=linestyle)
 

experiment/meteo_france_data/scm_models_data/abstract_study.py

@@ -401,22 +401,26 @@ class AbstractStudy(object):
 
         # Add legend for the marker
         if massif_name_to_marker_style is not None:
-            labels = ['\mathcal{M}_{\mu_1}', '\mathcal{M}_{\sigma_1}', '\mathcal{M}_{\mu_1, \sigma_1}']
-            markers = ["s", "^", "D"]
-            legend_elements = [
-                Line2D([0], [0], marker=marker, color='w', label='${}$'.format(label),
-                           markerfacecolor='w', markeredgecolor='k', markersize=8)
-                for label, marker in zip(labels, markers)
-            ]
+            legend_elements = cls.get_legend_for_model_symbol(markersize=8)
             ax.legend(handles=legend_elements, bbox_to_anchor=(0.01, 0.03), loc='lower left')
-            ax.annotate("Filled symbol = significant trend ", xy=(0.05, 0.015), xycoords='axes fraction',
-                         fontsize=7)
+            ax.annotate("Filled symbol = significant trend ", xy=(0.05, 0.015), xycoords='axes fraction', fontsize=7)
 
         if show:
             plt.show()
 
         return ax
 
+    @classmethod
+    def get_legend_for_model_symbol(cls, markersize):
+        labels = ['\mathcal{M}_{\mu_1}', '\mathcal{M}_{\sigma_1}', '\mathcal{M}_{\mu_1, \sigma_1}']
+        markers = ["s", "^", "D"]
+        legend_elements = [
+            Line2D([0], [0], marker=marker, color='w', label='${}$'.format(label),
+                   markerfacecolor='w', markeredgecolor='k', markersize=markersize)
+            for label, marker in zip(labels, markers)
+        ]
+        return legend_elements
+
     """ 
     CLASS ATTRIBUTES COMMON TO ALL OBJECTS 
     (written as object attributes/methods for simplicity)

experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py

@@ -34,7 +34,7 @@ SCM_STUDY_CLASS_TO_ABBREVIATION = {
     SafranSnowfall: 'SF3',
     CrocusSweTotal: 'SWE',
     CrocusSwe3Days: 'SWE3',
-    CrocusSnowLoadEurocode: 'GSL_Eurocode',
+    CrocusSnowLoadEurocode: 'GSL from annual maximum of HS and {}'.format(eurocode_snow_density),
     CrocusDepth: 'SD',
     CrocusSnowLoadTotal: 'max GSL',
     CrocusSnowLoad3Days: 'GSL3',

experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -7,7 +7,7 @@ from experiment.paper_past_snow_loads.paper_utils import paper_study_classes, pa
     load_altitude_to_visualizer
 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.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoadEurocode
 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 \
@@ -58,15 +58,17 @@ def major_result():
                            ConfidenceIntervalMethodFromExtremes.ci_mle][:]
     massif_names = None
     non_stationary_uncertainty = [False, True][:]
-    for study_class in paper_study_classes[2:]:
+    for study_class in paper_study_classes[:2]:
         intermediate_result(paper_altitudes, massif_names, non_stationary_uncertainty, uncertainty_methods, study_class)
 
 
 if __name__ == '__main__':
-    # major_result()
-    intermediate_result(altitudes=[900, 1200], massif_names=['Chartreuse', 'Vanoise', 'Mercantour'],
-                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
-                           ConfidenceIntervalMethodFromExtremes.ci_mle])
+    major_result()
+    # intermediate_result(altitudes=paper_altitudes, massif_names=['Chartreuse', 'Vanoise', 'Mercantour'],
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+    #                        ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+    #                     non_stationary_uncertainty=[False, True][1:],
+    #                     study_class=CrocusSnowLoadEurocode)
     # intermediate_result(altitudes=[900, 1200], massif_names=None)
     # intermediate_result(ALL_ALTITUDES_WITHOUT_NAN)
     # intermediate_result(paper_altitudes)

experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py

@@ -4,6 +4,9 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_ALTITUDES
+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.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 \
@@ -102,20 +105,20 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
             plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes)
 
     ax.legend(loc=2)
-    ax.set_ylim([-1, 16])
+    # ax.set_ylim([-1, 16])
     massif_name_str = massif_name.replace('_', ' ')
     eurocode_region_str = get_display_name_from_object_type(type(eurocode_region))
     is_non_stationary_model = non_stationary_context if isinstance(non_stationary_context,
                                                                    bool) else 'Non' in non_stationary_context
     if is_non_stationary_model:
-        non_stationary_context = 'non-stationary'
+        non_stationary_context = 'selected non-stationary models'
     else:
-        non_stationary_context = 'stationary'
-    title = '{} ({} Eurocodes area) with a {} model'.format(massif_name_str, eurocode_region_str,
-                                                            non_stationary_context)
+        non_stationary_context = 'the stationary model'
+    title = '{} massif with {}'.format(massif_name_str,  non_stationary_context)
     ax.set_title(title)
     ax.set_xticks(altitudes)
-    ax.set_ylabel(EUROCODE_RETURN_LEVEL_STR)
+    ylabel = EUROCODE_RETURN_LEVEL_STR.replace('GSL', SCM_STUDY_CLASS_TO_ABBREVIATION[type(visualizer.study)])
+    ax.set_ylabel(ylabel)
     ax.set_xlabel('Altitude (m)')
     ax.grid()
 
@@ -131,10 +134,20 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes):
                edgecolor='black', hatch='//')
         # Plot markers
         markers_kwargs = [v.massif_name_to_marker_style[massif_name] for v in visualizers]
+        for k in markers_kwargs:
+            k['markersize'] = 7
         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)
+    # Add a legend plot
+    legend_elements = AbstractStudy.get_legend_for_model_symbol(markersize=9)
+    ax2 = ax.twinx()
+    # ax2.legend(handles=legend_elements, bbox_to_anchor=(0.93, 0.7), loc='upper right')
+    # ax2.annotate("Filled symbol = significant trend ", xy=(0.85, 0.5), xycoords='axes fraction', fontsize=7)
+    ax2.legend(handles=legend_elements, loc='upper right')
+    ax2.annotate("Filled symbol = significant trend ", xy=(0.75, 0.97), xycoords='axes fraction', fontsize=10)
+    ax2.set_yticks([])
 
 
 def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,

experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py

@@ -2,16 +2,18 @@ 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.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_ALTITUDES
 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
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color, \
+    ci_method_to_label
 
 
 def plot_uncertainty_histogram(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
     """ Plot one graph for each non-stationary context
     :return:
     """
+    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)
@@ -33,24 +35,27 @@ def plot_histogram(altitude_to_visualizer, non_stationary_context):
         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.legend()
     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.plot_name = 'Percentages of exceedance with non_stationary={}'.format(non_stationary_context)
     # visualizer.show = True
-    # visualizer.show_or_save_to_file(no_title=True)
-    plt.show()
+    visualizer.show_or_save_to_file(no_title=True)
 
 
 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]
+    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]
+    three_percentages_of_excess = [(a, b, c) if b == c else (a, b, min(100 - epsilon, c)) for (a, b, c) in three_percentages_of_excess]
     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)
+    label = ci_method_to_label[ci_method]
+    color = ci_method_to_color[ci_method]
+    ax.bar(bincenters, y, width=width, color=color, yerr=yerr, label=label, ecolor='black', capsize=5)

extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py

@@ -17,7 +17,18 @@ ci_method_to_method_name = {
     ConfidenceIntervalMethodFromExtremes.ci_proflik: 'proflik',
 }
 
+# ci_method_to_color = {
+#     ConfidenceIntervalMethodFromExtremes.ci_mle: 'tab:brown',
+#     ConfidenceIntervalMethodFromExtremes.my_bayes: 'tab:green'
+# }
+
 ci_method_to_color = {
-    ConfidenceIntervalMethodFromExtremes.ci_mle: 'tab:brown',
-    ConfidenceIntervalMethodFromExtremes.my_bayes: 'tab:green'
+    ConfidenceIntervalMethodFromExtremes.ci_mle: 'yellowgreen',
+    ConfidenceIntervalMethodFromExtremes.my_bayes: 'darkgreen'
+}
+
+common_part_and_uncertainty = 'Return level and its uncertainty with the'
+ci_method_to_label = {
+    ConfidenceIntervalMethodFromExtremes.ci_mle: 'Bayesian procedure',
+    ConfidenceIntervalMethodFromExtremes.my_bayes: 'Delta method '
 }
\ No newline at end of file