diff --git a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
index af0b69a423e2a17120c04d418000a120e673970a..36a4f22c044a3d449c88a682884d94664847395f 100644
--- a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
+++ b/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
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
index 11ad1ce45ccbf40455c1401371bcb2dbd8f517c0..4fe9ff3ae18cccca9cb44940730fa095d6e19bfe 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
@@ -1,10 +1,11 @@
-
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__':
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values_with_mean_aic.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values_with_mean_aic.py
new file mode 100644
index 0000000000000000000000000000000000000000..7317f7f6f3c65b3de5b0c89d76e1d9fc663c8421
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values_with_mean_aic.py
@@ -0,0 +1,57 @@
+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
+
diff --git a/projects/exceeding_snow_loads/utils.py b/projects/exceeding_snow_loads/utils.py
index 5cbcec03b375f33d5a03567ff9b706eef9adb78a..62104769a695a572439ff7eb53a43b210404d930 100644
--- a/projects/exceeding_snow_loads/utils.py
+++ b/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,
]