[paper 1] add return level plot

experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py

@@ -223,9 +223,10 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
         epsilon = 0.5
         ax_lim = [min(all_quantiles) - epsilon, max(all_quantiles) + epsilon]
         ax.plot(standard_gumbel_quantiles, standard_gumbel_quantiles, color='k')
-        ax.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x', label='Stationary Gumbel model $\mathcal{M}_0$')
+        ax.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x',
+                label='Stationary Gumbel model $\mathcal{M}_0$')
         ax.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
-                 label='Selected model $\mathcal{M}_N$', **marker)
+                label='Selected model $\mathcal{M}_N$', **marker)
         ax.set_xlabel("Standard Gumbel quantile", fontsize=size)
         ax.set_ylabel("Standard Empirical quantile", fontsize=size)
         ax.legend(loc='upper left', prop={'size': 10})
@@ -239,6 +240,19 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
 
         plt.show()
 
+    def return_level_plot_comparison(self, ax, label, color=None):
+        # ax = plt.gca()
+        size = 15
+        # Load Gev parameter in 2017 for the unconstrained estimator
+        gev_params = self.unconstrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([2017]),
+                                                                                          is_transformed=False)  # type: GevParams
+        gev_params_with_corrected_shape = GevParams(loc=gev_params.location,
+                                                    scale=gev_params.scale,
+                                                    shape=0.5)
+        suffix = 'in 2017'
+        gev_params.return_level_plot_against_return_period(ax, color, linestyle='-', label=label, suffix_return_level_label=suffix)
+        gev_params_with_corrected_shape.return_level_plot_against_return_period(ax, color=color, linestyle='--', suffix_return_level_label=suffix)
+
     def compute_empirical_quantiles(self, estimator):
         empirical_quantiles = []
         for year, maximum in sorted(zip(self.years, self.maxima), key=lambda t: t[1]):

extreme_fit/distribution/gev/gev_params.py

 from collections import OrderedDict
+import matplotlib.pyplot as plt
 from typing import List
 
 from cached_property import cached_property
@@ -148,4 +149,23 @@ class GevParams(AbstractParams):
         if self.shape <= 0:
             return np.inf
         else:
-            return self.bound
\ No newline at end of file
+            return self.bound
+
+    def return_level_plot_against_return_period(self, ax=None, color=None, linestyle=None, label=None, show=False,
+                                                suffix_return_level_label=''):
+        if ax is None:
+            ax = plt.gca()
+        # Plot return level against return period
+        return_periods = list(range(2, 61))
+        quantiles = [self.quantile(1 - 1 / return_period) for return_period in return_periods]
+        return_period_to_quantile = dict(zip(return_periods, quantiles))
+        ax.vlines(50, 0, return_period_to_quantile[50])
+        ax.plot(return_periods, quantiles, color=color, linestyle=linestyle, label=label)
+        ax.set_xlabel('Return period')
+        ax.legend()
+
+        ax.set_xticks([10 * i for i in range(1, 7)])
+        ax.set_ylabel('Return level {}'.format(suffix_return_level_label))
+        plt.gca().set_ylim(bottom=0)
+        if show:
+            plt.show()

papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/main_qqplot_for_big_shapes.py

 from typing import Dict
+import matplotlib.pyplot as plt
 
 from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
     ALL_ALTITUDES_WITHOUT_NAN
-from experiment.exceeding_snow_loads.study_visualizer_for_non_stationary_trends import \
+from papers.exceeding_snow_loads.study_visualizer_for_non_stationary_trends import \
     StudyVisualizerForNonStationaryTrends
 
 
-def qqplots_for_biggest_shape_parameters_before_selection():
-    altitudes = ALL_ALTITUDES_WITHOUT_NAN
+
+def get_tuple_ordered_by_shape(fast=False):
+    if fast:
+        altitudes = [300]
+    else:
+        altitudes = ALL_ALTITUDES_WITHOUT_NAN
     altitude_to_visualizer = {altitude: StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude),
                                                                               select_only_acceptable_shape_parameter=False,
                                                                               multiprocessing=True)
                               for altitude in altitudes}
-    plot_qqplot_for_time_series_with_worst_shape_parameters(altitude_to_visualizer, nb_worst_examples=5)
-
-
-def plot_qqplot_for_time_series_with_worst_shape_parameters(
-        altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
-        nb_worst_examples=3):
     # Extract all the values
     l = []
     for a, v in altitude_to_visualizer.items():
-        l.extend([(a, v, m, t.unconstrained_estimator_gev_params.shape) for m, t in v.massif_name_to_trend_test_that_minimized_aic.items()])
+        l.extend([(a, v, m, t.unconstrained_estimator_gev_params.shape) for m, t in
+                  v.massif_name_to_trend_test_that_minimized_aic.items()])
     # Sort them and keep the highest examples
     l = sorted(l, key=lambda t: t[-1])
+    return l
+
+
+def plot_qqplot_for_time_series_with_worst_shape_parameters(tuple_ordered_by_shape, nb_worst_examples=5):
+    l = tuple_ordered_by_shape
     print('Highest examples:')
     for a, v, m, shape in l[-nb_worst_examples:][::-1]:
         print(a, m, shape)
@@ -34,6 +39,19 @@ def plot_qqplot_for_time_series_with_worst_shape_parameters(
         print(a, m, shape)
         v.qqplot(m)
 
+
+def plot_return_level_for_time_series_with_big_shape_parameters(tuple_ordered_by_shape, nb_worst_examples=5):
+    # Extract all the values
+    l = tuple_ordered_by_shape
+    print('Highest examples:')
+    ax = plt.gca()
+    colors = ['orange', 'red', 'blue', 'green', 'yellow']
+    for (a, v, m, shape), color in zip(l[-nb_worst_examples:][::-1], colors):
+        print(a, m, shape, color)
+        v.return_level_plot(ax, m, color)
+    plt.show()
+
+
 """
 10 Worst examples:
 300 Oisans 1.1070477650581747
@@ -50,7 +68,6 @@ def plot_qqplot_for_time_series_with_worst_shape_parameters(
 300 Bauges 0.39291367817905504
 """
 
-
 """
 for the worst example for -shape
 
@@ -66,6 +83,8 @@ for the worst example for -shape
 2700 Mont-Blanc 0.2712596135797868
 """
 
-
 if __name__ == '__main__':
-    qqplots_for_biggest_shape_parameters_before_selection()
+    fast = False
+    nb = 1 if fast else 5
+    tuple_ordered_by_shape = get_tuple_ordered_by_shape(fast=fast)
+    plot_return_level_for_time_series_with_big_shape_parameters(tuple_ordered_by_shape, nb_worst_examples=nb)

papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -2,7 +2,7 @@ from multiprocessing.pool import Pool
 
 import matplotlib as mpl
 
-from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days
 from papers.exceeding_snow_loads.paper_main_utils import load_altitude_to_visualizer
 from papers.exceeding_snow_loads.paper_utils import paper_study_classes, paper_altitudes
 from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_diagnosis_risk import plot_diagnosis_risk

papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -401,6 +401,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         marker = self.massif_name_to_marker_style[massif_name]
         trend_test.qqplot_wrt_standard_gumbel(marker, color)
 
+    def return_level_plot(self, ax, massif_name, color=None):
+        trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
+        model_name = '${}$'.format(trend_test.label)
+        label = '{} at {}m with {}'.format(massif_name, self.altitude, model_name)
+        trend_test.return_level_plot_comparison(ax, label, color)
+
     # Part 4 - Trend plot
 
     @property

test/test_extreme_fit/test_distribution/test_gev/test_gev_params.py

@@ -54,6 +54,10 @@ class TestGevParams(unittest.TestCase):
         gev_params = GevParams(loc=mu, shape=chi, scale=sigma)
         self.assertEqual(gev_params.variance, ((sigma / chi) ** 2) * (gamma(1 - 2 * chi) - (gamma(1 - chi) ** 2)))
 
+    def test_return_level_plot(self):
+        gev_params = GevParams(loc=0.0, shape=0.0, scale=1.0)
+        gev_params.return_level_plot_against_return_period(show=False)
+
 
 if __name__ == '__main__':
     unittest.main()