[contrasting] add plot selection for paper 2

extreme_trend/abstract_gev_trend_test.py

@@ -70,14 +70,18 @@ class AbstractGevTrendTest(object):
 
     @property
     def goodness_of_fit_anderson_test(self):
-        assert self.SIGNIFICANCE_LEVEL == 0.05
-        # significance_level=array([25. , 10. ,  5. ,  2.5,  1. ]))
-        index_for_significance_level_5_percent = 2
+        # significance_level_to_index = dict(zip([0.25, 0.1, 0.05, 0.025, 0.01], list(range(5))))
+        # print(significance_level_to_index)
+        # assert self.SIGNIFICANCE_LEVEL in significance_level_to_index
+        # # significance_level=array([25. , 10. ,  5. ,  2.5,  1. ]))
+        # index_for_significance_level_5_percent = 2
         quantiles = self.compute_empirical_quantiles(estimator=self.unconstrained_estimator)
         # test_res = anderson(quantiles, dist='gumbel_r')  # type: AndersonResult
         # return test_res.statistic < test_res.critical_values[index_for_significance_level_5_percent]
+        # print(quantiles)
         test = cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution(quantiles)
-        return test[1] > self.SIGNIFICANCE_LEVEL
+        _, ander_darling_test_pvalue = test
+        return ander_darling_test_pvalue > self.SIGNIFICANCE_LEVEL
 
     @property
     def name(self):

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py

@@ -4,6 +4,9 @@ from multiprocessing.pool import Pool
 import matplotlib as mpl
 
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranPrecipitation1Day
+from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest
+from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.plot_selection_curves_paper2 import \
+    plot_selection_curves_paper2
 from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.shape_plot import shape_plot
 from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.snowfall_plot import \
     plot_snowfall_mean, plot_snowfall_change_mean
@@ -77,7 +80,7 @@ def intermediate_result(altitudes, massif_names=None,
     # validation_plot(altitude_to_visualizer, order_derivative=0)
     # validation_plot(altitude_to_visualizer, order_derivative=1)
     plot_snowfall_mean(altitude_to_visualizer)
-    plot_selection_curves(altitude_to_visualizer, paper1=False)
+    plot_selection_curves_paper2(altitude_to_visualizer)
     plot_snowfall_change_mean(altitude_to_visualizer)
     # shape_plot(altitude_to_visualizer)
 
@@ -86,15 +89,16 @@ def major_result():
     uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.ci_mle][:]
     # massif_names = ['Beaufortain', 'Vercors']
     massif_names = None
-    study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day][::-1]
+    study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day][:]
     model_subsets_for_uncertainty = None
     altitudes = paper_altitudes
-    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
+    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
     # altitudes = [900, 1200, 1500, 1800][:2]
     # altitudes = [1800, 2100, 2400, 2700][:2]
     # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
     # altitudes = draft_altitudes
-
+    # for significance_level in [0.1, 0.05][]:
+    AbstractGevTrendTest.SIGNIFICANCE_LEVEL = 0.1
     for study_class in study_classes:
         intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
                             uncertainty_methods, study_class, multiprocessing=False)

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/plot_selection_curves_paper2.py

+from typing import Dict
+import matplotlib.pyplot as plt
+
+from extreme_data.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest
+from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
+    StudyVisualizerForMeanValues
+from projects.exceeding_snow_loads.section_results.plot_selection_curves import merge_counter, \
+    get_ordered_list_from_counter, permute
+from projects.exceeding_snow_loads.utils import dpi_paper1_figure, get_trend_test_name
+from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
+
+
+def plot_selection_curves_paper2(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
+    visualizer = list(altitude_to_visualizer.values())[0]
+
+    ax = create_adjusted_axes(1, 1)
+
+    vs = [v for v in altitude_to_visualizer.values()]
+
+    selected_counter = merge_counter([v.selected_trend_test_class_counter for v in vs])
+    selected_and_anderson_counter = merge_counter([v.selected_and_anderson_trend_test_class_counter for v in vs])
+    selected_and_anderson_and_likelihood_counter = merge_counter(
+        [v.selected_and_anderson_and_likelihood_ratio_trend_test_class_counter() for v in vs])
+
+    total_of_selected_models = sum(selected_counter.values())
+    l = sorted(enumerate(visualizer.non_stationary_trend_test), key=lambda e: selected_counter[e[1]])
+    permutation = [i for i, v in l][::-1]
+
+    select_list = get_ordered_list_from_counter(selected_counter, total_of_selected_models, visualizer, permutation)
+    selected_and_anderson_list = get_ordered_list_from_counter(selected_and_anderson_counter, total_of_selected_models,
+                                                               visualizer, permutation)
+    selected_and_anderson_and_likelihood_list = get_ordered_list_from_counter(
+        selected_and_anderson_and_likelihood_counter, total_of_selected_models, visualizer, permutation)
+
+    labels = [get_trend_test_name(t) for t in visualizer.non_stationary_trend_test]
+    labels = permute(labels, permutation)
+    print(select_list, sum(select_list))
+
+    nb_selected_models = sum(select_list)
+    nb_selected_and_anderson_models = sum(selected_and_anderson_list)
+    nb_selected_and_anderson_and_likelihood_models = sum(selected_and_anderson_and_likelihood_list)
+    nb_selected_models_not_passing_any_test = nb_selected_models - nb_selected_and_anderson_models
+    nb_selected_models_just_passing_anderson = nb_selected_and_anderson_models - nb_selected_and_anderson_and_likelihood_models
+
+    # parameters
+    width = 5
+    size = 30
+    legend_fontsize = 30
+    labelsize = 15
+    linewidth = 3
+    x = [10 * (i + 1) for i in range(len(select_list))]
+    for l, color, label, nb in zip([select_list, selected_and_anderson_list, selected_and_anderson_and_likelihood_list],
+                               ['white', 'grey', 'black'],
+                               ['Not passing any test', 'Just passing anderson test ',
+                                'Passing both anderson and likelihood tests '],
+                               [nb_selected_models_not_passing_any_test, nb_selected_models_just_passing_anderson, nb_selected_and_anderson_and_likelihood_models]):
+        label += ' ({} \%)'.format(round(nb, 1))
+        ax.bar(x, l, width=width, color=color, edgecolor='black', label=label, linewidth=linewidth)
+
+    ax.legend(loc='upper right', prop={'size': size})
+    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)
+    ax.yaxis.grid()
+    ax.set_xticklabels(labels)
+
+    # Save plot
+    visualizer.plot_name = 'Selection curves with significance level = {} '.format(AbstractGevTrendTest.SIGNIFICANCE_LEVEL)
+    visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
+    plt.close()

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py

@@ -26,7 +26,8 @@ def plot_snowfall_change_mean(altitude_to_visualizer: Dict[int, StudyVisualizerF
 
     variables = ['mean annual maxima', '50 year return level']
     mean_indicators = [True, False]
-    for variable, mean_indicator in zip(variables, mean_indicators):
+    l = list(zip(variables, mean_indicators))
+    for variable, mean_indicator in l[:1]:
         for relative in [False, True]:
             if relative:
                 variable += str(' (relative)')
@@ -41,7 +42,8 @@ def plot_snowfall_change_mean(altitude_to_visualizer: Dict[int, StudyVisualizerF
             visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km,
                                           label='Slope for changes of {} of {}\n for every km of elevation ({})'.format(
                                               variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
-                                          add_x_label=False)
+                                          add_x_label=False,
+                                          add_text=True, massif_name_to_text=massif_to_r2_score_text)
             # Value at 2000 m
             massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
             visualizer.plot_abstract_fast(massif_name_to_mean_at_2000,
@@ -123,9 +125,7 @@ def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], d
                 indicator_str = 'mean'
             else:
                 indicator_str = '50 year return level'
-            res = [(a, t)
-                   for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))
-                   if not t.unconstrained_model_is_stationary]
+            res = [(a, t) for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))]
             if derivative:
                 if mean_indicator:
                     if relative:

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py

@@ -58,15 +58,37 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
     # Override the main dict massif_name_to_trend_test_that_minimized_aic
 
     @property
-    def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
+    def super_massif_name_to_trend_test_that_minimized_aic(self):
+        return super().massif_name_to_trend_test_that_minimized_aic
+
+    @property
+    def massif_name_to_trend_test_that_minimized_aic_after_anderson_test(self):
         return {m: t for m, t in super().massif_name_to_trend_test_that_minimized_aic.items()
                 if t.goodness_of_fit_anderson_test}
 
     @property
-    def super_massif_name_to_trend_test_that_minimized_aic(self):
-        return super().massif_name_to_trend_test_that_minimized_aic
+    def massif_name_to_trend_test_that_minimized_aic_after_anderson_test_and_likelihood_ratio_test(self):
+        return {m: t for m, t in self.massif_name_to_trend_test_that_minimized_aic_after_anderson_test.items()
+                if t.is_significant}
+
+    @property
+    def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
+        return self.massif_name_to_trend_test_that_minimized_aic_after_anderson_test
+
+    # Counter for the selection process
+
+    @cached_property
+    def selected_trend_test_class_counter(self):
+        return Counter([type(t) for t in self.super_massif_name_to_trend_test_that_minimized_aic.values()])
+
+    @cached_property
+    def selected_and_anderson_trend_test_class_counter(self):
+        return Counter([type(t) for t in self.massif_name_to_trend_test_that_minimized_aic_after_anderson_test.values()])
+
+    def selected_and_anderson_and_likelihood_ratio_trend_test_class_counter(self):
+        return Counter([type(t) for t in self.massif_name_to_trend_test_that_minimized_aic_after_anderson_test_and_likelihood_ratio_test.values()])
 
-    # Study of the mean
+# Study of the mean
 
     @cached_property
     def massif_name_to_empirical_mean(self):