[exceeding] improve qqplot display

extreme_trend/abstract_gev_trend_test.py

@@ -14,7 +14,7 @@ from extreme_fit.distribution.gev.gev_params import GevParams
 from extreme_fit.model.margin_model.utils import \
     MarginFitMethod
 from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
-    StationaryTemporalModel
+    StationaryTemporalModel, GumbelTemporalModel
 from extreme_fit.model.utils import SafeRunException
 from root_utils import classproperty
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -50,6 +50,10 @@ class AbstractGevTrendTest(object):
         return fitted_linear_margin_estimator(self.unconstrained_model_class, self.coordinates, self.dataset,
                                               self.starting_year, self.fit_method)
 
+    @property
+    def unconstrained_model_is_stationary(self):
+        return self.unconstrained_model_class in [StationaryTemporalModel, GumbelTemporalModel]
+
     # Likelihood ratio test
 
     @property
@@ -168,6 +172,58 @@ class AbstractGevTrendTest(object):
 
     # Some display function
 
+    def intensity_plot_wrt_standard_gumbel_simple(self, massif_name, altitude, color='grey', psnow=1):
+        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', color='black')
+
+        ax_twiny = ax.twiny()
+
+        # Plot for the selected model with line break
+        # 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))
+
+        end_real_proba = 1 - (0.02 / psnow)
+        if self.unconstrained_model_is_stationary:
+            years = [None]
+        else:
+            color = None
+            years = [1959, 1989, 2019]
+        for j, year in enumerate(years):
+            year_suffix = '' if year is None else 'for the year t={}'.format(year)
+            if j == 0:
+                label = 'Selected model i.e. ${}$'.format(self.label)
+                if year is not None:
+                    label += '\n' + year_suffix
+            else:
+                label = year_suffix
+            print(label)
+            self.plot_model(ax, year, end_proba=end_real_proba, label=label, color=color)
+
+        ax_lim = [-1.5, 4]
+        ax.set_xlim(ax_lim)
+        ax_twiny.set_xlim(ax_lim)
+        ax.set_xticks([-1 + i for i in range(6)])
+        epsilon = 0.005
+        ax.set_ylim(bottom=-epsilon)
+        lalsize = 13
+        ax.tick_params(axis='both', which='major', labelsize=lalsize)
+        ax_twiny.tick_params(axis='both', which='major', labelsize=lalsize)
+
+        ax.yaxis.grid()
+
+        ax.set_xlabel("Standard Gumbel quantile", fontsize=size)
+        ax.set_ylabel("Non-zero annual maxima of GSL ({})".format(AbstractSnowLoadVariable.UNIT), fontsize=size)
+        ax.legend(prop={'size': 17})
+
+
     def intensity_plot_wrt_standard_gumbel(self, massif_name, altitude, psnow):
         ax = plt.gca()
         sorted_maxima = sorted(self.maxima)
@@ -238,13 +294,14 @@ class AbstractGevTrendTest(object):
         ax.legend(prop={'size': 17})
 
     def plot_model(self, ax, year, start_proba=0.02, end_proba=0.98, label='', color=None):
+        if year is None:
+            # for the stationary model
+            year = 2019
+
         standard_gumbel = GevParams(0, 1, 0)
         start_quantile = standard_gumbel.quantile(start_proba)
         end_quantile = standard_gumbel.quantile(end_proba)
         extended_quantiles = np.linspace(start_quantile, end_quantile, 500)
-        label = 'Y({})'.format(year) if year is not None else label
-        if year is None:
-            year = 2019
         gev_params_year = self.get_unconstrained_gev_params(year)
         extended_maxima = [gev_params_year.gumbel_inverse_standardization(q) for q in extended_quantiles]
 
@@ -301,21 +358,6 @@ class AbstractGevTrendTest(object):
         label_generic = '{} massif \nat {} m '.format(massif_name, altitude)
         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]
-            print(len(q), len(standard_gumbel_quantiles))
-            ax.plot(standard_gumbel_quantiles, q, linestyle='None',
-                    label=label_generic
-                          + '(discarded model is ${}$\n'.format('\mathcal{M}_{\zeta_0, \sigma_1}')
-                          + 'with $\zeta_0=0.84$)',
-                    marker='o')
 
         ax.set_xlabel("Standard Gumbel quantile", fontsize=size)
         ax.set_ylabel("Standard Empirical quantile", fontsize=size)
@@ -328,8 +370,6 @@ class AbstractGevTrendTest(object):
         ax.grid()
         ax.tick_params(labelsize=size)
 
-        plt.show()
-
     def get_standard_gumbel_quantiles(self):
         # Standard Gumbel quantiles
         standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
@@ -338,7 +378,6 @@ class AbstractGevTrendTest(object):
         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]

extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py

@@ -420,9 +420,10 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
     # Part 3 - QQPLOT
 
-    def intensity_plot(self, massif_name, psnow, color=None):
+    def intensity_plot(self, massif_name, color=None):
+        psnow = self.massif_name_to_psnow[massif_name]
         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)
+        trend_test.intensity_plot_wrt_standard_gumbel_simple(massif_name, self.altitude, color, psnow)
         self.plot_name = 'intensity_plot_{}_{}_{}_{}'.format(self.altitude, massif_name, psnow, trend_test.unconstrained_estimator_gev_params.shape)
         self.show_or_save_to_file(add_classic_title=False, no_title=True)
         plt.close()
@@ -430,6 +431,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     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)
+        self.plot_name = 'qpplot_plot_{}_{}_{}'.format(self.altitude, massif_name, trend_test.unconstrained_estimator_gev_params.shape)
+        self.show_or_save_to_file(add_classic_title=False, no_title=True)
+        plt.close()
 
     def return_level_plot(self, ax, ax2, massif_name, color=None):
         trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]

projects/exceeding_snow_loads/checks/qqplot/plot_qqplot.py

@@ -11,10 +11,10 @@ from extreme_fit.model.margin_model.utils import \
     MarginFitMethod
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
     AbstractExtractEurocodeReturnLevel
-from projects.exceeding_snow_loads.data.main_example_swe_total_plot import tuples_for_examples_paper1
 from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import \
     StudyVisualizerForNonStationaryTrends
 from extreme_fit.distribution.gev.gev_params import GevParams
+from projects.exceeding_snow_loads.section_data.main_example_swe_total_plot import tuples_for_examples_paper1
 
 
 def extract_time_serimes_with_worst_number_of_zeros(altitude_to_visualizer, nb_worst_examples):
@@ -49,6 +49,14 @@ def plot_intensity_against_gumbel_quantile_for_time_series_with_missing_zeros(
         v.intensity_plot(m, p)
 
 
+def plot_intensity_against_gumbel_quantile_for_3_examples(
+        altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+    for color, altitude, m in tuples_for_examples_paper1():
+        v = altitude_to_visualizer[altitude]
+        v.intensity_plot(m, color)
+        v.qqplot(m, color)
+
+
 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:

projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py

@@ -2,6 +2,8 @@ from multiprocessing.pool import Pool
 
 import matplotlib as mpl
 
+from projects.exceeding_snow_loads.checks.qqplot.plot_qqplot import \
+    plot_intensity_against_gumbel_quantile_for_3_examples
 from projects.exceeding_snow_loads.section_results.plot_selection_curves import plot_selection_curves
 from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map
 
@@ -69,8 +71,9 @@ def intermediate_result(altitudes, massif_names=None,
     # 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)
     # uncertainty_interval_size(altitude_to_visualizer)
+    plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)
 
 
 def major_result():
@@ -86,7 +89,7 @@ def major_result():
     model_subsets_for_uncertainty = None
     # study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
     altitudes = paper_altitudes
-    # altitudes = [900]
+    altitudes = [900, 1800, 2700]
     for study_class in study_classes:
         intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
                             uncertainty_methods, study_class,