[paper 1] add plot with the difference in return level (and relative...

[paper 1] add plot with the difference in return level (and relative difference) for the big shapes. we conclude that there is big difference indeed

experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py

-import numpy as np
 from math import ceil, floor
+
 import matplotlib.pyplot as plt
-import pandas as pd
+import numpy as np
 from cached_property import cached_property
 from scipy.stats import chi2
 
@@ -9,13 +9,12 @@ from experiment.eurocode_data.utils import EUROCODE_QUANTILE
 from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
 from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
     fitted_linear_margin_estimator
-from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
+from extreme_fit.distribution.gev.gev_params import GevParams
 from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
-    AbstractTemporalLinearMarginModel, TemporalMarginFitMethod
+    TemporalMarginFitMethod
 from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
     StationaryTemporalModel
 from extreme_fit.model.utils import SafeRunException
-from extreme_fit.distribution.gev.gev_params import GevParams
 from root_utils import classproperty
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
@@ -240,25 +239,58 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
 
         plt.show()
 
+    def compute_empirical_quantiles(self, estimator):
+        empirical_quantiles = []
+        for year, maximum in sorted(zip(self.years, self.maxima), key=lambda t: t[1]):
+            gev_param = estimator.margin_function_from_fit.get_gev_params_with_big_shape_and_correct_shape(
+                coordinate=np.array([year]),
+                is_transformed=False)
+            maximum_standardized = gev_param.gumbel_standardization(maximum)
+            empirical_quantiles.append(maximum_standardized)
+        return empirical_quantiles
+
+    # For some visualizations
+
     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, gev_params_with_corrected_shape = self.get_gev_params_with_big_shape_and_correct_shape()
+        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 return_level_plot_difference(self, ax, ax2, label, color=None):
+        gev_params, gev_params_with_corrected_shape = self.get_gev_params_with_big_shape_and_correct_shape()
+        return_periods = list(range(2, 61))
+        quantile_with_big_shape = gev_params.get_return_level(return_periods)
+        quantile_with_small_shape = gev_params_with_corrected_shape.get_return_level(return_periods)
+        difference_quantile = quantile_with_big_shape - quantile_with_small_shape
+        relative_difference = 100 * difference_quantile / quantile_with_small_shape
+        # Plot the difference on the two axis
+
+        ax.vlines(50, 0, np.max(difference_quantile))
+        ax.plot(return_periods, difference_quantile, color=color, linestyle='-', label=label)
+        ax.legend(loc='upper left')
+        difference_ylabel = 'difference return level in 2017'
+        ax.set_ylabel(difference_ylabel + ' (kPa)')
+
+        ax2.vlines(50, 0, np.max(relative_difference))
+        ax2.plot(return_periods, relative_difference, color=color, linestyle='--', label=label)
+        ax2.legend(loc='upper right')
+        ax2.yaxis.grid()
+        ax2.set_ylabel('relative ' + difference_ylabel + ' (%)')
+
+        ax.set_xlabel('Return period')
+        ax.set_xticks([10 * i for i in range(1, 7)])
+        plt.gca().set_ylim(bottom=0)
+
+    def get_gev_params_with_big_shape_and_correct_shape(self):
         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]):
-            gev_param = estimator.margin_function_from_fit.get_gev_params(
-                coordinate=np.array([year]),
-                is_transformed=False)
-            maximum_standardized = gev_param.gumbel_standardization(maximum)
-            empirical_quantiles.append(maximum_standardized)
-        return empirical_quantiles
+        return gev_params, gev_params_with_corrected_shape

extreme_fit/distribution/gev/gev_params.py

@@ -157,7 +157,7 @@ class GevParams(AbstractParams):
             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]
+        quantiles = self.get_return_level(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)
@@ -169,3 +169,6 @@ class GevParams(AbstractParams):
         plt.gca().set_ylim(bottom=0)
         if show:
             plt.show()
+
+    def get_return_level(self, return_periods):
+        return np.array([self.quantile(1 - 1 / return_period) for return_period in return_periods])

papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/main_qqplot_for_big_shapes.py

@@ -45,10 +45,11 @@ def plot_return_level_for_time_series_with_big_shape_parameters(tuple_ordered_by
     l = tuple_ordered_by_shape
     print('Highest examples:')
     ax = plt.gca()
+    ax2 = ax.twinx()
     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)
+        v.return_level_plot(ax, ax2, m, color)
     plt.show()
 
 

papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -397,11 +397,11 @@ 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):
+    def return_level_plot(self, ax, ax2, 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)
+        trend_test.return_level_plot_difference(ax, ax2, label, color)
 
     # Part 4 - Trend plot