diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index dc078514207e2230959234a9f1312d0fa4bcb138..85df6b20bda42f18411aee377bcc33450fdc8c28 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/extreme_trend/abstract_gev_trend_test.py
@@ -372,11 +372,11 @@ class AbstractGevTrendTest(object):
size = 15
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
- epsilon = 0.5
+ # constrained_empirical_quantiles = self.compute_empirical_quantiles(self.constrained_estimator)
+ all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles
+ epsilon = 0.1
ax_lim = [min(all_quantiles) - epsilon, max(all_quantiles) + epsilon]
- ax.plot(standard_gumbel_quantiles, standard_gumbel_quantiles, color='k')
+
# ax.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x',
# label='Stationary Gumbel model $\mathcal{M}_0$')
@@ -390,10 +390,12 @@ class AbstractGevTrendTest(object):
ax.legend(loc='lower right', prop={'size': 10})
ax.set_xlim(ax_lim)
ax.set_ylim(ax_lim)
+
+ ax.plot(ax_lim, ax_lim, color='k')
ticks = [i for i in range(ceil(ax_lim[0]), floor(ax_lim[1]) + 1)]
ax.set_xticks(ticks)
ax.set_yticks(ticks)
- ax.grid()
+ # ax.grid()
ax.tick_params(labelsize=size)
def get_standard_gumbel_quantiles(self):
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 b07a5eca8e0167fd5898fdaefc5ebb87005c086f..7f99652def22204dd953989b6fd5a0800db13306 100644
--- a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
+++ b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
@@ -144,59 +144,21 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
d = {m: v for m, v in d.items() if self.massif_name_to_psnow[m] >= 0.9}
return d
- @property
- def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
- return self.massif_name_to_trend_test_tuple[0]
-
- @property
- def massif_name_to_stationary_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
- return self.massif_name_to_trend_test_tuple[1]
-
- @property
- def massif_name_to_gumbel_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
- return self.massif_name_to_trend_test_tuple[2]
-
@cached_property
- def massif_name_to_trend_test_tuple(self) -> Tuple[
- Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest]]:
- print('cached', self.altitude, id(self))
-
+ def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
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 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]
- if self.select_only_model_that_pass_anderson_test and \
- (not trend_test_that_minimized_aic.goodness_of_fit_anderson_test):
- print('not anderson')
- else:
- print('ok')
+ if (not self.select_only_model_that_pass_anderson_test) or \
+ trend_test_that_minimized_aic.goodness_of_fit_anderson_test:
massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
-
- # Extract the stationary model that minimized AIC
- stationary_trend_tests_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
- [GumbelVersusGumbel, GevStationaryVersusGumbel]]
- if len(stationary_trend_tests_that_minimized_aic) == 0:
- stationary_trend_test_that_minimized_aic = None
- else:
- stationary_trend_test_that_minimized_aic = stationary_trend_tests_that_minimized_aic[0]
- massif_name_to_stationary_trend_test_that_minimized_aic[
- massif_name] = stationary_trend_test_that_minimized_aic
- # Extract the Gumbel model that minimized AIC
- gumbel_trend_tests = [t for t in sorted_trend_test if
- type(t) in [GumbelVersusGumbel, GumbelLocationTrendTest, GumbelScaleTrendTest,
- GumbelLocationAndScaleTrendTest]]
- if len(gumbel_trend_tests) > 0:
- gumbel_trend_test_that_minimized_aic = gumbel_trend_tests[0]
else:
- gumbel_trend_test_that_minimized_aic = None
- massif_name_to_gumbel_trend_test_that_minimized_aic[massif_name] = gumbel_trend_test_that_minimized_aic
+ print('anderson test was not passed')
- 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
+ return massif_name_to_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]
@@ -360,12 +322,6 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
def massif_name_and_model_subset_to_model_class(self, massif_name, model_subset_for_uncertainty):
if model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gumbel:
return GumbelTemporalModel
- if model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gev:
- return StationaryTemporalModel
- elif model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gumbel_and_gev:
- return self.massif_name_to_stationary_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
- elif model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel:
- return self.massif_name_to_gumbel_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
elif model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel_and_gev:
return self.massif_name_to_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
else:
@@ -378,6 +334,10 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
# Compute for the uncertainty massif names
massifs_names = set(self.massif_name_to_years_and_maxima_for_model_fitting.keys()). \
intersection(self.uncertainty_massif_names)
+ # Update the name of the massif (because some massifs might have been removed by anderson test)
+ if model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel_and_gev\
+ and self.select_only_model_that_pass_anderson_test:
+ massifs_names = massifs_names.intersection(set(self.massif_name_to_trend_test_that_minimized_aic.keys()))
arguments = [
[self.massif_name_to_years_and_maxima_for_model_fitting[m],
self.massif_name_and_model_subset_to_model_class(m, model_subset_for_uncertainty),
diff --git a/projects/exceeding_snow_loads/section_data/main_example_swe_total_plot.py b/projects/exceeding_snow_loads/section_data/main_example_swe_total_plot.py
index bb78ef8dfe0f60caa629844f13daeaa043b8cba2..9b6a86ed4f3e4b5af803764a21caf99cc851b2a5 100644
--- a/projects/exceeding_snow_loads/section_data/main_example_swe_total_plot.py
+++ b/projects/exceeding_snow_loads/section_data/main_example_swe_total_plot.py
@@ -55,7 +55,7 @@ def max_graph_annual_maxima_poster():
last_plot = color == "magenta"
label = '{} massif at {}m'.format(massif_name, altitude)
tight_pad = {'h_pad': 0.2}
- snow_abbreviation = 'annual maxima of ' + snow_abbreviation
+ snow_abbreviation = 'annual maximum of ' + snow_abbreviation
study_visualizer.visualize_max_graphs_poster(massif_name, altitude, snow_abbreviation, color, label,
last_plot, ax, tight_pad=tight_pad,
dpi=dpi_paper1_figure,
diff --git a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
index 46f93e5223ea9a8262f0b88bbeee0f63e203478d..efc2f2d4ec54b4c0037ec39d78db2d771e0ca5f1 100644
--- a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
+++ b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
@@ -74,8 +74,10 @@ def intermediate_result(altitudes, massif_names=None,
plot_uncertainty_massifs(altitude_to_visualizer)
plot_uncertainty_histogram(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)
+
+ # Additional plots
+ # uncertainty_interval_size(altitude_to_visualizer)
# plot_full_diagnostic(altitude_to_visualizer)
@@ -85,13 +87,9 @@ def major_result():
# massif_names = ['Beaufortain', 'Vercors']
massif_names = None
study_classes = paper_study_classes[:1]
- # model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
- # ModelSubsetForUncertainty.stationary_gumbel_and_gev,
- # ModelSubsetForUncertainty.non_stationary_gumbel,
- # ModelSubsetForUncertainty.non_stationary_gumbel_and_gev]
model_subsets_for_uncertainty = None
# study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
- altitudes = [300, 600]
+ # altitudes = [300, 600, 900, 1800, 2700][:2]
altitudes = paper_altitudes
# altitudes = [900, 1800, 2700][:1]
for study_class in study_classes:
diff --git a/projects/exceeding_snow_loads/utils.py b/projects/exceeding_snow_loads/utils.py
index 62104769a695a572439ff7eb53a43b210404d930..c28ee205adc1f036a08d96263476ca947f309c5c 100644
--- a/projects/exceeding_snow_loads/utils.py
+++ b/projects/exceeding_snow_loads/utils.py
@@ -19,7 +19,7 @@ from extreme_trend.trend_test_two_parameters.gumbel_test_two_parameters import \
GumbelLocationAndScaleTrendTest, GumbelLocationQuadraticTrendTest, GumbelScaleQuadraticTrendTest
paper_altitudes = ALL_ALTITUDES_WITHOUT_NAN
-paper_study_classes = [CrocusSnowLoadTotal, CrocusSnowLoadEurocode, CrocusSnowLoad3Days][:2]
+paper_study_classes = [CrocusSnowLoadTotal, CrocusSnowLoadEurocode][:]
# dpi_paper1_figure = 700
dpi_paper1_figure = None
NON_STATIONARY_TREND_TEST_PAPER_1 = [GumbelVersusGumbel,