[paper 2] add plot diagnosis risk. refactor excess metric in...

[paper 2] add plot diagnosis risk. refactor excess metric in study_visualizer_for_non_stationary_trends.py.

experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -8,6 +8,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
     ALL_ALTITUDES_WITHOUT_NAN
 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, ModelSubsetForUncertainty
+from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_diagnosis_risk import plot_diagnosis_risk
 from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_trend_curves import plot_trend_curves, \
     plot_trend_map
 from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import \
@@ -68,10 +69,11 @@ def intermediate_result(altitudes, massif_names=None,
             _ = compute_minimized_aic(visualizer)
 
     # Plots
-    plot_trend_map(altitude_to_visualizer)
+    # plot_trend_map(altitude_to_visualizer)
+    plot_diagnosis_risk(altitude_to_visualizer)
     # plot_trend_curves(altitude_to_visualizer={a: v for a, v in altitude_to_visualizer.items() if a >= 900})
     # plot_uncertainty_massifs(altitude_to_visualizer)
-    # plot_uncertainty_histogram(altitude_to_visualizer)
+    plot_uncertainty_histogram(altitude_to_visualizer)
 
 
 def major_result():
@@ -92,7 +94,7 @@ def major_result():
 
 if __name__ == '__main__':
     # major_result()
-    intermediate_result(altitudes=ALL_ALTITUDES_WITHOUT_NAN[:2], massif_names=None,
+    intermediate_result(altitudes=[1500, 1800], massif_names=None,
                         uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
                                              ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
                         multiprocessing=True)

experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_diagnosis_risk.py

+import matplotlib.pyplot as plt
+
+from experiment.paper_past_snow_loads.paper_utils import ModelSubsetForUncertainty, dpi_paper1_figure
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+    ConfidenceIntervalMethodFromExtremes
+
+
+def plot_diagnosis_risk(altitude_to_visualizer):
+    ax = plt.gca()
+    altitudes = list(altitude_to_visualizer.keys())
+    visualizers = list(altitude_to_visualizer.values())
+    visualizer = visualizers[0]
+    model_subset_for_uncertainty = ModelSubsetForUncertainty.non_stationary_gumbel_and_gev
+    ci_method = ConfidenceIntervalMethodFromExtremes.ci_mle
+
+    plot_mean_exceedance(visualizers, altitudes, ax, model_subset_for_uncertainty, ci_method)
+
+    visualizer.plot_name = 'Diagnosis plot'
+    visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
+    plt.close()
+
+    plt.show()
+
+
+def plot_mean_exceedance(visualizers, altitudes, ax, model_subset_for_uncertainty, ci_method):
+    l = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[3:] for v in visualizers]
+    diff, diff_c, diff_e = zip(*l)
+    ax.plot(altitudes, diff, marker='o', color='red', label='diff')
+    ax.plot(altitudes, diff_c, marker='o', color='yellow', label='diff c')
+    ax.plot(altitudes, diff_e, marker='o', color='orange', label='diff e')
+
+    ax.set_ylabel("Mean exceedance (kN)", fontsize=15)
+    ax.set_xlabel("Altitude", fontsize=15)
+    ax.yaxis.grid()
+
+    ax.legend()

experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py

@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_ALTITUDES
-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.confidence_interval_method import ci_method_to_color, \
@@ -33,6 +33,10 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
     ax = plt.gca()
     altitudes = np.array(list(altitude_to_visualizer.keys()))
     bincenters = altitudes
+
+    fontsize_label = 15
+    legend_size = 15
+    # Plot histogram
     for j, ci_method in enumerate(visualizer.uncertainty_methods):
         if len(visualizer.uncertainty_methods) == 2:
             offset = -50 if j == 0 else 50
@@ -41,13 +45,15 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
         else:
             width = 200
         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)
     ax.tick_params(labelsize=fontsize_label)
     if not (len(visualizer.uncertainty_methods) == 1
             and visualizer.uncertainty_methods[0] == ConfidenceIntervalMethodFromExtremes.ci_mle):
         ax.legend(loc='upper left', prop={'size': legend_size})
+    # ax.set_ylabel('Massifs whose 50-year return level\n'
+    #               'exceeds French standards (\%)', fontsize=fontsize_label)
     ax.set_ylabel('Massifs exceeding French standards (\%)', fontsize=fontsize_label)
     ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
     ax.set_ylim([0, 100])
@@ -59,7 +65,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
 
 
 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
+    three_percentages_of_excess = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[:3] 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

@@ -6,6 +6,7 @@ import matplotlib.pyplot as plt
 import numpy as np
 from cached_property import cached_property
 
+from experiment.eurocode_data.eurocode_region import C2, C1, E
 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
@@ -359,16 +360,31 @@ 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, model_subset_for_uncertainty):
-        eurocode_and_uncertainties = [(self.massif_name_to_eurocode_values[massif_name],
-                                       self.triplet_to_eurocode_uncertainty[
-                                           (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)
+    def excess_metrics(self, ci_method, model_subset_for_uncertainty):
+        triplet = [(massif_name_to_eurocode_region[massif_name],
+                    self.massif_name_to_eurocode_values[massif_name],
+                    self.triplet_to_eurocode_uncertainty[(ci_method, model_subset_for_uncertainty, massif_name)])
+                   for massif_name in self.uncertainty_massif_names]
+        # First array for histogram
+        a = 100 * np.array([(uncertainty.confidence_interval[0] > eurocode,
+                             uncertainty.mean_estimate > eurocode,
+                             uncertainty.confidence_interval[1] > eurocode)
+                            for _, eurocode, uncertainty in triplet])
+        a = np.mean(a, axis=0)
+        # Second array for curve
+        b = [[] for _ in range(3)]
+        for eurocode_region, eurocode, uncertainty in triplet:
+            diff = uncertainty.mean_estimate - eurocode
+            b[0].append(diff)
+            if eurocode_region in [C1, C2]:
+                b[1].append(diff)
+            if eurocode_region in [E]:
+                b[2].append(diff)
+
+        b = [np.mean(np.array(e)) for e in b]
+        # Return the concatenated results
+        concatenated_result = list(a) + list(b)
+        return concatenated_result
 
     # Part 3 - QQPLOT