[contrasting] add analysis with mean aic

extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py

@@ -155,24 +155,13 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     @cached_property
     def massif_name_to_trend_test_tuple(self) -> Tuple[
         Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest]]:
-        starting_year = None
+
         massif_name_to_trend_test_that_minimized_aic = {}
         massif_name_to_stationary_trend_test_that_minimized_aic = {}
         massif_name_to_gumbel_trend_test_that_minimized_aic = {}
-        for massif_name, (x, y) in self.massif_name_to_years_and_maxima_for_model_fitting.items():
-            quantile_level = self.massif_name_to_eurocode_quantile_level_in_practice[massif_name]
-            all_trend_test = [
-                t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level,
-                  fit_method=self.fit_method)
-                for t in self.non_stationary_trend_test]  # type: List[AbstractGevTrendTest]
-            # Exclude GEV models whose shape parameter is not in the support of the prior distribution for GMLE
-            if self.select_only_acceptable_shape_parameter:
-                acceptable_shape_parameter = lambda s: -0.5 <= s <= 0.5  # physically acceptable prior
-                all_trend_test = [t for t in all_trend_test
-                                  if (acceptable_shape_parameter(t.unconstrained_estimator_gev_params_last_year.shape)
-                                      and acceptable_shape_parameter(
-                                t.unconstrained_estimator_gev_params_first_year.shape))]
-            sorted_trend_test = sorted(all_trend_test, key=lambda t: t.aic)
+        for massif_name in self.massif_name_to_years_and_maxima_for_model_fitting.keys():
+            # Compute sorted trend test
+            sorted_trend_test = self.get_sorted_trend_test(massif_name)
 
             # Extract the stationary or non-stationary model that minimized AIC
             trend_test_that_minimized_aic = sorted_trend_test[0]
@@ -198,6 +187,24 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
         return massif_name_to_trend_test_that_minimized_aic, massif_name_to_stationary_trend_test_that_minimized_aic, massif_name_to_gumbel_trend_test_that_minimized_aic
 
+    def get_sorted_trend_test(self, massif_name):
+        x, y = self.massif_name_to_years_and_maxima_for_model_fitting[massif_name]
+        starting_year = None
+        quantile_level = self.massif_name_to_eurocode_quantile_level_in_practice[massif_name]
+        all_trend_test = [
+            t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level,
+              fit_method=self.fit_method)
+            for t in self.non_stationary_trend_test]  # type: List[AbstractGevTrendTest]
+        # Exclude GEV models whose shape parameter is not in the support of the prior distribution for GMLE
+        if self.select_only_acceptable_shape_parameter:
+            acceptable_shape_parameter = lambda s: -0.5 <= s <= 0.5  # physically acceptable prior
+            all_trend_test = [t for t in all_trend_test
+                              if (acceptable_shape_parameter(t.unconstrained_estimator_gev_params_last_year.shape)
+                                  and acceptable_shape_parameter(
+                            t.unconstrained_estimator_gev_params_first_year.shape))]
+        sorted_trend_test = sorted(all_trend_test, key=lambda t: t.aic)
+        return sorted_trend_test
+
     # Part 1 - Trends
 
     @property

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py

-
 from multiprocessing.pool import Pool
 
 import matplotlib as mpl
+import numpy as np
 
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranPrecipitation1Day, \
     SafranPrecipitation3Days, SafranSnowfall3Days
+from extreme_data.meteo_france_data.scm_models_data.utils import Season
 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
@@ -13,7 +14,10 @@ from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and
     plot_snowfall_mean, plot_snowfall_change_mean
 from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
     StudyVisualizerForMeanValues
-from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.validation_plot import validation_plot
+from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values_with_mean_aic import \
+    StudyVisualizerForMeanValuesWithMeanAic
+from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.validation_plot import \
+    validation_plot
 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
 
@@ -30,7 +34,6 @@ from projects.exceeding_snow_loads.section_results.plot_uncertainty_curves impor
 from projects.exceeding_snow_loads.utils import paper_study_classes, paper_altitudes
 from root_utils import NB_CORES
 
-
 import matplotlib.pyplot as plt
 
 
@@ -43,14 +46,19 @@ def minor_result(altitude):
 
 
 def compute_minimized_aic(visualizer):
-    _ = visualizer.massif_name_to_trend_test_that_minimized_aic
+    if isinstance(visualizer, StudyVisualizerForMeanValuesWithMeanAic):
+        _ = visualizer.massif_name_and_trend_test_class_to_trend_test
+    else:
+        _ = visualizer.massif_name_to_trend_test_that_minimized_aic
     return True
 
 
 def intermediate_result(altitudes, massif_names=None,
-                        model_subsets_for_uncertainty=None, uncertainty_methods=None,
+                        model_subsets_for_uncertainty=None,
+                        uncertainty_methods=None,
                         study_class=SafranSnowfall1Day,
-                        multiprocessing=False):
+                        multiprocessing=False, study_visualizer_class=StudyVisualizerForMeanValues,
+                        season=Season.annual):
     """
     Plot all the trends for all altitudes
     And enable to plot uncertainty plot for some specific massif_names, uncertainty methods to be fast
@@ -64,7 +72,8 @@ def intermediate_result(altitudes, massif_names=None,
     # Load altitude to visualizer
     altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncertainty,
                                                          study_class, uncertainty_methods,
-                                                         study_visualizer_class=StudyVisualizerForMeanValues)
+                                                         study_visualizer_class=study_visualizer_class,
+                                                         season=season)
     # Load variable object efficiently
     for v in altitude_to_visualizer.values():
         _ = v.study.year_to_variable_object
@@ -76,34 +85,69 @@ def intermediate_result(altitudes, massif_names=None,
     else:
         for visualizer in visualizers:
             _ = compute_minimized_aic(visualizer)
+    # Aggregate the choice for the minimizer
+    aggregate(visualizers)
 
     # Plots
-    # validation_plot(altitude_to_visualizer, order_derivative=0)
+    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_paper2(altitude_to_visualizer)
+    # plot_selection_curves_paper2(altitude_to_visualizer)
     plot_snowfall_change_mean(altitude_to_visualizer)
-    # shape_plot(altitude_to_visualizer)
+    shape_plot(altitude_to_visualizer)
 
 
 def major_result():
     uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.ci_mle][:]
     # massif_names = ['Beaufortain', 'Vercors']
     massif_names = None
-    study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day][::-1][1:]
+    study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day][::-1][:]
     # study_classes = [SafranSnowfall3Days, SafranPrecipitation3Days][::-1]
     model_subsets_for_uncertainty = None
     altitudes = paper_altitudes
     altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
     # altitudes = [900, 1200, 1500, 1800][:2]
-    # altitudes = [1800, 2100, 2400, 2700][:2]
+    # altitudes = [1800, 2100, 2400, 2700][:3]
     # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
     # altitudes = draft_altitudes
     # for significance_level in [0.1, 0.05][]:
     AbstractGevTrendTest.SIGNIFICANCE_LEVEL = 0.05
+    study_visualizer_class = [StudyVisualizerForMeanValues, StudyVisualizerForMeanValuesWithMeanAic][0]
+    season = [Season.annual, Season.winter_extended][1]
     for study_class in study_classes:
         intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
-                            uncertainty_methods, study_class, multiprocessing=False)
+                            uncertainty_methods, study_class, multiprocessing=False,
+                            study_visualizer_class=study_visualizer_class,
+                            season=season)
+
+
+def aggregate(visualizers):
+    visualizer = visualizers[0]
+    if not isinstance(visualizer, StudyVisualizerForMeanValuesWithMeanAic):
+        return
+    massif_names = set.union(*[set(v.massif_names) for v in visualizers])
+    massif_names = list(massif_names)
+    trend_tests = visualizer.non_stationary_trend_test
+
+    massif_name_to_trend_test_to_aic_list = {m: {t: [] for t in trend_tests}
+                                             for m in massif_names}
+    for v in visualizers:
+        for (m, t), t2 in v.massif_name_and_trend_test_class_to_trend_test.items():
+            massif_name_to_trend_test_to_aic_list[m][t].append(t2.aic)
+
+    massif_name_to_trend_test_with_minimial_mean_aic = {}
+    for m in massif_names:
+        trend_test_and_mean_aic = [(t, np.array(massif_name_to_trend_test_to_aic_list[m][t]).mean())
+                                   for t in trend_tests]
+        sorted_l = sorted(trend_test_and_mean_aic, key=lambda e: e[1])
+        trend_test_with_minimial_mean_aic = sorted_l[0][0]
+        massif_name_to_trend_test_with_minimial_mean_aic[m] = trend_test_with_minimial_mean_aic
+    print(massif_name_to_trend_test_with_minimial_mean_aic)
+    for v in visualizers:
+        v.massif_name_to_trend_test_with_minimial_mean_aic = {}
+        for m, t in massif_name_to_trend_test_with_minimial_mean_aic.items():
+            if (m, t) in v.massif_name_and_trend_test_class_to_trend_test:
+                v.massif_name_to_trend_test_with_minimial_mean_aic[m] = v.massif_name_and_trend_test_class_to_trend_test[(m, t)]
 
 
 if __name__ == '__main__':

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values_with_mean_aic.py

+from collections import Counter
+from typing import Dict
+
+import matplotlib.pyplot as plt
+import numpy as np
+from cached_property import cached_property
+
+from extreme_data.eurocode_data.utils import YEAR_OF_INTEREST_FOR_RETURN_LEVEL
+from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest
+from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
+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.utils import NON_STATIONARY_TREND_TEST_PAPER_2
+
+
+class StudyVisualizerForMeanValuesWithMeanAic(StudyVisualizerForMeanValues):
+
+    def __init__(self, study: AbstractStudy, show=True, save_to_file=False, only_one_graph=False, only_first_row=False,
+                 vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False,
+                 temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False,
+                 complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True,
+                 uncertainty_methods=None, model_subsets_for_uncertainty=None, uncertainty_massif_names=None,
+                 effective_temporal_covariate=YEAR_OF_INTEREST_FOR_RETURN_LEVEL, relative_change_trend_plot=True,
+                 non_stationary_trend_test_to_marker=None, fit_method=MarginFitMethod.extremes_fevd_mle,
+                 select_only_acceptable_shape_parameter=True, fit_gev_only_on_non_null_maxima=False,
+                 fit_only_time_series_with_ninety_percent_of_non_null_values=True):
+        super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
+                         year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
+                         transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis,
+                         normalization_under_one_observations, uncertainty_methods, model_subsets_for_uncertainty,
+                         uncertainty_massif_names, effective_temporal_covariate, relative_change_trend_plot,
+                         non_stationary_trend_test_to_marker, fit_method, select_only_acceptable_shape_parameter,
+                         fit_gev_only_on_non_null_maxima, fit_only_time_series_with_ninety_percent_of_non_null_values)
+        self.massif_name_to_trend_test_with_minimial_mean_aic = None
+
+    @property
+    def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
+        if self.massif_name_to_trend_test_with_minimial_mean_aic is None:
+            raise NotImplementedError('Aggregation must be run first')
+        else:
+            return self.massif_name_to_trend_test_with_minimial_mean_aic
+
+    @property
+    def massif_names(self):
+        return [m for m, _ in self.massif_name_and_trend_test_class_to_trend_test.keys()]
+
+    @cached_property
+    def massif_name_and_trend_test_class_to_trend_test(self):
+        d = {}
+        for massif_name in self.massif_name_to_years_and_maxima_for_model_fitting.keys():
+            trend_tests = self.get_sorted_trend_test(massif_name)
+            for trend_test in trend_tests:
+                d[(massif_name, type(trend_test))] = trend_test
+        return d
+

projects/exceeding_snow_loads/utils.py

@@ -31,20 +31,20 @@ NON_STATIONARY_TREND_TEST_PAPER_1 = [GumbelVersusGumbel,
 
 NON_STATIONARY_TREND_TEST_PAPER_2 = [
     # Gumbel models
-    #GumbelVersusGumbel,
-    #GumbelLocationTrendTest,
-    #GumbelScaleTrendTest,
-    #GumbelLocationAndScaleTrendTest,
+    GumbelVersusGumbel,
+    GumbelLocationTrendTest,
+    GumbelScaleTrendTest,
+    GumbelLocationAndScaleTrendTest,
     # GEV models with constant shape
-    #GevVersusGev,
+    GevVersusGev,
     GevLocationTrendTest,
-    # GevScaleTrendTest,
-    #GevLocationAndScaleTrendTest,
+    GevScaleTrendTest,
+    GevLocationAndScaleTrendTest,
     # GEV models with linear shape
     #GevShapeTrendTest,
     #GevLocationAndShapeTrendTest, GevScaleAndShapeTrendTest, GevLocationAndScaleAndShapeTrendTest,
     # Quadratic model for the Gev/Gumbel and for the location/scale
-    #GevQuadraticLocationTrendTest, GevQuadraticScaleTrendTest, GumbelLocationQuadraticTrendTest, GumbelScaleQuadraticTrendTest,
+    GevQuadraticLocationTrendTest, GevQuadraticScaleTrendTest, GumbelLocationQuadraticTrendTest, GumbelScaleQuadraticTrendTest,
 ]