[paper 1] intensity plots. fit ylabel pour selection curves.

experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py

@@ -6,6 +6,7 @@ from cached_property import cached_property
 from scipy.stats import chi2
 
 from experiment.eurocode_data.utils import EUROCODE_QUANTILE
+from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import AbstractSnowLoadVariable
 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
@@ -209,13 +210,87 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
 
     # Some display function
 
+    def intensity_plot_wrt_standard_gumbel(self, massif_name, altitude, psnow):
+        ax = plt.gca()
+        sorted_maxima = sorted(self.maxima)
+        label_generic = '{} massif \nat {} m '.format(massif_name, altitude)
+        size = 15
+
+        # Plot for the empirical
+        standard_gumbel_quantiles = self.get_standard_gumbel_quantiles()
+        ax.plot(standard_gumbel_quantiles, sorted_maxima, linestyle='None',
+                label='Empirical model', marker='o')
+
+        # Plot for the selected model
+        unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator)
+        ax.plot(unconstrained_empirical_quantiles, sorted_maxima,
+                label='Selected model, which is ${}$'.format(self.label))
+        print(unconstrained_empirical_quantiles)
+
+        ax_twiny = ax.twiny()
+        return_periods = [10, 25, 50, 100]
+        quantiles = self.get_standard_quantiles_for_return_periods(return_periods, psnow)
+        print(quantiles)
+        ax_twiny.plot(quantiles, [0 for _ in quantiles], linewidth=0)
+        ax_twiny.set_xticks(quantiles)
+        ax_twiny.set_xticklabels([return_period for return_period in return_periods])
+        ax_twiny.set_xlabel('Return period w.r.t all annual maxima of GSL (years)', fontsize=size)
+
+        # Plot for the discarded model
+        if 'Verdon' in massif_name and altitude == 300:
+            q = [-1.4541688117485054, -1.2811308174310914, -1.216589300814509, -0.7635793791201918, -0.6298883422064275,
+                 -0.5275954855697504, -0.4577268043676126, -0.4497570331795861, -0.1647955002136654,
+                 -0.14492222503785876, -0.139173823298689, -0.11945617994263039, -0.07303100174657867,
+                 -5.497308509286266e-05, 0.13906416388625908, 0.15274793441408543, 0.1717763342727519,
+                 0.17712605315013535, 0.17900143646245203, 0.371986176207554, 0.51640780422156, 0.7380550963951035,
+                 0.7783015252180445, 0.887836077295502, 0.917853338231094, 0.9832396811506262, 1.0359396416309927,
+                 1.1892663813729711, 1.2053261113817888, 1.5695111391491652, 2.3223652143938476, 2.674882764437432,
+                 2.6955728524900406, 2.8155882785356896, 3.282838470153471, 3.2885313947906765]
+            color = 'red'
+            ax.plot(q, sorted_maxima,
+                    label='Discarded model, which is ${}$\n'.format('\mathcal{M}_{\zeta_0, \sigma_1}')
+                          + 'with $\zeta_0=0.84$',
+                    color=color)
+
+            # Extend the curve linear a bit if the return period 50 is not in the quantiles
+
+            def compute_slope_intercept(x, y):
+                x1, x2 = x[-2:]
+                y1, y2 = y[-2:]
+                a = (y2 - y1) / (x2 - x1)
+                b = y1 - a * x1
+                return a, b
+
+            def compute_maxima_corresponding_to_return_period(return_period_quantiles, quantiles, model_maxima):
+                a, b = compute_slope_intercept(quantiles, model_maxima)
+                return a * return_period_quantiles + b
+
+            quantile_return_period_50 = quantiles[-1]
+            if max(q) < quantile_return_period_50:
+                maxima_extended = compute_maxima_corresponding_to_return_period(quantile_return_period_50,
+                                                                                q,
+                                                                                sorted_maxima)
+                ax.plot([q[-1], quantile_return_period_50],
+                        [sorted_maxima[-1], maxima_extended], linestyle='--',
+                        color=color)
+
+        all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles + quantiles
+        epsilon = 0.5
+        ax_lim = [min(all_quantiles) - epsilon, max(all_quantiles) + epsilon]
+        ax.set_xlim(ax_lim)
+        ax_twiny.set_xlim(ax_lim)
+
+        ax.legend()
+
+        ax.set_xlabel("Standard Gumbel quantile", fontsize=size)
+        ax.set_ylabel("Non-zero annual maxima of GSL ({})".format(AbstractSnowLoadVariable.UNIT), fontsize=size)
+        ax.legend(loc='lower right', prop={'size': 10})
+        plt.show()
+
     def qqplot_wrt_standard_gumbel(self, massif_name, altitude):
         ax = plt.gca()
         size = 15
-        # Standard Gumbel quantiles
-        standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
-        n = len(self.years)
-        standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)]
+        standard_gumbel_quantiles = self.get_standard_gumbel_quantiles()
         unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator)
         constrained_empirical_quantiles = self.compute_empirical_quantiles(self.constrained_estimator)
         all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles + constrained_empirical_quantiles
@@ -230,7 +305,14 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
         ax.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
                 label=label_generic + '(selected model is ${}$)'.format(self.label), marker='o')
         if 'Verdon' in massif_name and altitude == 300:
-            q = [-1.4541688117485054, -1.2811308174310914, -1.216589300814509, -0.7635793791201918, -0.6298883422064275, -0.5275954855697504, -0.4577268043676126, -0.4497570331795861, -0.1647955002136654, -0.14492222503785876, -0.139173823298689, -0.11945617994263039, -0.07303100174657867, -5.497308509286266e-05, 0.13906416388625908, 0.15274793441408543, 0.1717763342727519, 0.17712605315013535, 0.17900143646245203, 0.371986176207554, 0.51640780422156, 0.7380550963951035, 0.7783015252180445, 0.887836077295502, 0.917853338231094, 0.9832396811506262, 1.0359396416309927, 1.1892663813729711, 1.2053261113817888, 1.5695111391491652, 2.3223652143938476, 2.674882764437432, 2.6955728524900406, 2.8155882785356896, 3.282838470153471, 3.2885313947906765]
+            q = [-1.4541688117485054, -1.2811308174310914, -1.216589300814509, -0.7635793791201918, -0.6298883422064275,
+                 -0.5275954855697504, -0.4577268043676126, -0.4497570331795861, -0.1647955002136654,
+                 -0.14492222503785876, -0.139173823298689, -0.11945617994263039, -0.07303100174657867,
+                 -5.497308509286266e-05, 0.13906416388625908, 0.15274793441408543, 0.1717763342727519,
+                 0.17712605315013535, 0.17900143646245203, 0.371986176207554, 0.51640780422156, 0.7380550963951035,
+                 0.7783015252180445, 0.887836077295502, 0.917853338231094, 0.9832396811506262, 1.0359396416309927,
+                 1.1892663813729711, 1.2053261113817888, 1.5695111391491652, 2.3223652143938476, 2.674882764437432,
+                 2.6955728524900406, 2.8155882785356896, 3.282838470153471, 3.2885313947906765]
             print(len(q), len(standard_gumbel_quantiles))
             ax.plot(standard_gumbel_quantiles, q, linestyle='None',
                     label=label_generic
@@ -251,6 +333,20 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
 
         plt.show()
 
+    def get_standard_gumbel_quantiles(self):
+        # Standard Gumbel quantiles
+        standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
+        n = len(self.years)
+        standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)]
+        return standard_gumbel_quantiles
+
+    def get_standard_quantiles_for_return_periods(self, return_periods, psnow):
+        n = len(self.years)
+        p_list = [1 - ((1 / return_period) / psnow) for return_period in return_periods]
+        standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
+        corresponding_quantiles = [standard_gumbel_distribution.quantile(p) for p in p_list]
+        return corresponding_quantiles
+
     def compute_empirical_quantiles(self, estimator):
         empirical_quantiles = []
         for year, maximum in sorted(zip(self.years, self.maxima), key=lambda t: t[1]):

papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py

@@ -16,9 +16,7 @@ from papers.exceeding_snow_loads.study_visualizer_for_non_stationary_trends impo
 from extreme_fit.distribution.gev.gev_params import GevParams
 
 
-def plot_qqplot_for_time_series_with_missing_zeros(
-        altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
-        nb_worst_examples=3):
+def extract_time_serimes_with_worst_number_of_zeros(altitude_to_visualizer, nb_worst_examples):
     # Extract all the values
     l = []
     for a, v in altitude_to_visualizer.items():
@@ -29,9 +27,25 @@ def plot_qqplot_for_time_series_with_missing_zeros(
     for a, v, m, p in l:
         print(a, m, p)
         print('Last standard quantile (depends on the number of data):', last_quantile(p))
+    return l
+
+
+def plot_qqplot_for_time_series_with_missing_zeros(
+        altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
+        nb_worst_examples=3):
+    l = extract_time_serimes_with_worst_number_of_zeros(altitude_to_visualizer, nb_worst_examples)
+    for a, v, m, p in l:
         v.qqplot(m)
 
 
+def plot_intensity_against_gumbel_quantile_for_time_series_with_missing_zeros(
+        altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
+        nb_worst_examples=3):
+    l = extract_time_serimes_with_worst_number_of_zeros(altitude_to_visualizer, nb_worst_examples)
+    for a, v, m, p in l:
+        v.intensity_plot(m, p)
+
+
 def plot_qqplot_for_time_series_examples(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
     marker_altitude_massif_name_for_paper1 = tuples_for_examples_paper1()
     for color, a, m in marker_altitude_massif_name_for_paper1:
@@ -109,4 +123,4 @@ if __name__ == '__main__':
     # plot_qqplot_wrt_standard_gumbel(altitude_to_visualizer)
     # plot_hist_psnow(altitude_to_visualizer)
     # plot_qqplot_for_time_series_examples(altitude_to_visualizer)
-    plot_qqplot_for_time_series_with_missing_zeros(altitude_to_visualizer, nb_worst_examples=3)
+    plot_intensity_against_gumbel_quantile_for_time_series_with_missing_zeros(altitude_to_visualizer, nb_worst_examples=3)

papers/exceeding_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py

@@ -6,7 +6,7 @@ from cached_property import cached_property
 
 from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map
 from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
-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
 from experiment.trend_analysis.abstract_score import MeanScore
 

papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -68,10 +68,10 @@ def intermediate_result(altitudes, massif_names=None,
     # Plots
     # 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_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_selection_curves(altitude_to_visualizer)
+    plot_selection_curves(altitude_to_visualizer)
 
 
 def major_result():
@@ -91,11 +91,11 @@ def major_result():
 
 
 if __name__ == '__main__':
-    # major_result()
-    intermediate_result(altitudes=paper_altitudes, massif_names=['Beaufortain', 'Vercors'],
-                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
-                                             ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
-                        multiprocessing=True)
+    major_result()
+    # intermediate_result(altitudes=[900], massif_names=['Beaufortain', 'Vercors'],
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+    #                                          ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+    #                     multiprocessing=True)
     # intermediate_result(altitudes=[900, 1200], massif_names=['Maurienne'],
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
     #                                          ConfidenceIntervalMethodFromExtremes.ci_mle][1:],

papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py

@@ -51,8 +51,8 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
            label='Significative model',
            linewidth=linewidth)
     ax.legend(loc='upper right', prop={'size': size})
-    ax.set_ylabel('Time series where the model is selected\n'
-                  'i.e. where it minimizes the AIC (\%)', fontsize=legend_fontsize)
+    ax.set_ylabel(' Frequency of selected model w.r.t all time series \n '
+                  'i.e. for all massifs and altitudes (\%)', fontsize=legend_fontsize)
     ax.set_xlabel('Models', fontsize=legend_fontsize)
     ax.tick_params(axis='both', which='major', labelsize=labelsize)
     ax.set_xticks(x)

papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -392,6 +392,10 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
     # Part 3 - QQPLOT
 
+    def intensity_plot(self, massif_name, psnow, color=None):
+        trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
+        trend_test.intensity_plot_wrt_standard_gumbel(massif_name, self.altitude, psnow)
+
     def qqplot(self, massif_name, color=None):
         trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
         trend_test.qqplot_wrt_standard_gumbel(massif_name, self.altitude)