From 60d5976ad4de223556181f6ae821547cca18a6a9 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 4 Mar 2020 12:53:22 +0100
Subject: [PATCH] [paper 1] fix main plots for the paper
---
.../shape/main_shape_repartition.py | 4 ++--
.../exceeding_snow_loads/paper_main_utils.py | 9 +++----
.../main_result_trends_and_return_levels.py | 14 +++++------
.../plot_selection_curves.py | 12 +++++-----
.../plot_trend_curves.py | 4 ++--
.../plot_uncertainty_curves.py | 2 +-
.../plot_uncertainty_histogram.py | 24 +++++++++++++++----
...dy_visualizer_for_non_stationary_trends.py | 10 ++++++--
8 files changed, 49 insertions(+), 30 deletions(-)
diff --git a/papers/exceeding_snow_loads/check_mle_convergence_for_trends/shape/main_shape_repartition.py b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/shape/main_shape_repartition.py
index 846d1aae..3ba5be62 100644
--- a/papers/exceeding_snow_loads/check_mle_convergence_for_trends/shape/main_shape_repartition.py
+++ b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/shape/main_shape_repartition.py
@@ -27,8 +27,8 @@ def main_shape_repartition(altitudes, massif_names=None,
if __name__ == '__main__':
# main_shape_repartition([900], save_to_file=False)
- # main_shape_repartition([900, 1800, 2700])
+ main_shape_repartition([900, 1800, 2700])
# main_shape_repartition([300, 600, 900, 1200, 1500, 1800, 2700])
- main_shape_repartition(paper_altitudes, study_visualizer_class=StudyVisualizerForShape, save_to_file=True)
+ # main_shape_repartition(paper_altitudes, study_visualizer_class=StudyVisualizerForShape, save_to_file=True)
# main_shape_repartition([300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200],
# study_visualizer_class=StudyVisualizerForShape, save_to_file=True)
diff --git a/papers/exceeding_snow_loads/paper_main_utils.py b/papers/exceeding_snow_loads/paper_main_utils.py
index 4b5fc27d..fabaea8a 100644
--- a/papers/exceeding_snow_loads/paper_main_utils.py
+++ b/papers/exceeding_snow_loads/paper_main_utils.py
@@ -11,16 +11,13 @@ def load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncer
study_visualizer_class=StudyVisualizerForNonStationaryTrends,
save_to_file=True):
fit_method = TemporalMarginFitMethod.extremes_fevd_mle
- select_only_acceptable_shape_parameter = True
- print('Fit method: {}, Select only acceptable shape parameter: {}'
- .format(fit_method, select_only_acceptable_shape_parameter))
altitude_to_visualizer = OrderedDict()
for altitude in altitudes:
altitude_to_visualizer[altitude] = study_visualizer_class(
study=study_class(altitude=altitude), multiprocessing=True, save_to_file=save_to_file,
uncertainty_massif_names=massif_names, uncertainty_methods=uncertainty_methods,
model_subsets_for_uncertainty=model_subsets_for_uncertainty, fit_method=fit_method,
- select_only_acceptable_shape_parameter=select_only_acceptable_shape_parameter)
+ 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)
return altitude_to_visualizer
-
-
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index 12e04ac9..8826ac1e 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -69,7 +69,7 @@ def intermediate_result(altitudes, massif_names=None,
plot_trend_map(altitude_to_visualizer)
# 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_uncertainty_histogram(altitude_to_visualizer)
# plot_selection_curves(altitude_to_visualizer)
@@ -77,7 +77,7 @@ def major_result():
uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
massif_names = None
- study_classes = paper_study_classes[:1]
+ study_classes = paper_study_classes[:2]
# model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
# ModelSubsetForUncertainty.stationary_gumbel_and_gev,
# ModelSubsetForUncertainty.non_stationary_gumbel,
@@ -90,11 +90,11 @@ def major_result():
if __name__ == '__main__':
- # major_result()
- intermediate_result(altitudes=[1800], massif_names=None,
- uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
- ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
- multiprocessing=True)
+ major_result()
+ # intermediate_result(altitudes=[1800], massif_names=None,
+ # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+ # ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+ # multiprocessing=True)
# intermediate_result(altitudes=[900, 1200], massif_names=['Maurienne'],
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
# ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py
index 454b0238..698222a9 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py
+++ b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py
@@ -31,12 +31,12 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
selected_and_signifcative_list = get_ordered_list_from_counter(selected_and_significative_counter, total_of_selected_models, visualizer, permutation)
labels = permute(['${}$'.format(t.label) for t in visualizer.non_stationary_trend_test], permutation)
- # print(l)
- # print(select_list)
- # print(selected_and_signifcative_list)
- # [(5, <class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters.GevLocationAgainstGumbel'> ), (6, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters.GevScaleAgainstGumbel' > ), (2, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelScaleTrendTest' > ), (1, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelLocationTrendTest' > ), (7, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_three_parameters.gev_trend_test_three_parameters.GevLocationAndScaleTrendTestAgainstGumbel' > ), (3, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gumbel_test_two_parameters.GumbelLocationAndScaleTrendTest' > ), (4, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GevStationaryVersusGumbel' > ), (0, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelVersusGumbel' > )]
- # [32.520325203252035, 25.609756097560975, 12.195121951219512, 9.34959349593496, 9.34959349593496, 4.471544715447155, 3.252032520325203, 3.252032520325203]
- # [0, 13.821138211382113, 7.317073170731708, 8.536585365853659, 5.691056910569106, 1.6260162601626016, 1.2195121951219512, 2.4390243902439024]
+ print(l)
+ print(select_list)
+ print(selected_and_signifcative_list)
+ # [(5, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters.GevLocationAgainstGumbel'> ), (6, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters.GevScaleAgainstGumbel' > ), (2, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelScaleTrendTest' > ), (1, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelLocationTrendTest' > ), (7, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_three_parameters.gev_trend_test_three_parameters.GevLocationAndScaleTrendTestAgainstGumbel' > ), (3, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gumbel_test_two_parameters.GumbelLocationAndScaleTrendTest' > ), (4, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GevStationaryVersusGumbel' > ), (0, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelVersusGumbel' > )]
+ # [32.64462809917355, 24.380165289256198, 12.396694214876034, 9.50413223140496, 9.090909090909092, 5.785123966942149, 3.71900826446281, 2.479338842975207]
+ # [0, 13.223140495867769, 7.851239669421488, 8.264462809917354, 4.958677685950414, 2.479338842975207, 0.8264462809917356, 2.0661157024793386]
# parameters
width = 5
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_trend_curves.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_trend_curves.py
index 890084bf..2d47750f 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_trend_curves.py
+++ b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_trend_curves.py
@@ -17,7 +17,7 @@ def plot_trend_map(altitude_to_visualizer):
for altitude, visualizer in altitude_to_visualizer.items():
if 900 <= altitude <= 4200:
- add_color = (visualizer.study.altitude - 1500) % 900 == 0
+ add_color = (visualizer.study.altitude - 1500) % 1200 == 0
visualizer.plot_trends(max_abs_tdrl_above_900, add_colorbar=add_color)
# Plot 2700 also with a colorbar
if altitude == 2700:
@@ -52,7 +52,7 @@ def plot_trend_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNonSta
ax.bar(altitudes, percent_decrease, width=width, color=color, edgecolor='blue', label='decreasing trend',
linewidth=linewidth)
ax.bar(altitudes, percent_decrease_signi, width=width, color=color, edgecolor='black',
- label='significative decreasing trend',
+ label='significant decreasing trend',
linewidth=linewidth)
ax.legend(loc='upper left', prop={'size': size})
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
index 83400351..c9e04aae 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
+++ b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
@@ -151,7 +151,7 @@ def add_title(ax, eurocode_region, massif_name, non_stationary_context):
def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, legend_size, fontsize):
visualizers = [v for a, v in altitude_to_visualizer.items()
- if a in valid_altitudes and massif_name in v.uncertainty_massif_names]
+ if a in valid_altitudes and massif_name in v.massif_names_fitted]
if len(visualizers) > 0:
tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in visualizers]
# Plot bars
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
index 76936a3d..3d9c660f 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
+++ b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
@@ -46,7 +46,6 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
width = 200
plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width=width)
-
ax.set_xticks(altitudes)
ax.tick_params(labelsize=fontsize_label)
if not (len(visualizer.uncertainty_methods) == 1
@@ -57,19 +56,36 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
ax.set_ylabel('Massifs exceeding French standards (\%)', fontsize=fontsize_label)
ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
ax.set_ylim([0, 100])
+ ax.yaxis.grid()
+
+ ax_twiny = ax.twiny()
+ ax_twiny.plot(altitudes, [0 for _ in altitudes], linewidth=0)
+ ax_twiny.tick_params(labelsize=fontsize_label)
+ ax_twiny.set_xlim(ax.get_xlim())
+ ax_twiny.set_xticks(altitudes)
+ nb_massif_names = [v.study.nb_study_massif_names for v in altitude_to_visualizer.values()]
+ print(nb_massif_names)
+ ax_twiny.set_xticklabels(nb_massif_names)
+ ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=fontsize_label)
+
ax.set_yticks([10 * i for i in range(11)])
- visualizer.plot_name = 'Percentages of exceedance with {}'.format(get_display_name_from_object_type(model_subset_for_uncertainty))
+ visualizer.plot_name = 'Percentages of exceedance with {}'.format(
+ get_display_name_from_object_type(model_subset_for_uncertainty))
# visualizer.show = True
visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
ax.clear()
+ ax_twiny.clear()
+ plt.close()
def plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width):
three_percentages_of_excess = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[:3] for v in
visualizers]
epsilon = 0.5
- three_percentages_of_excess = [(a, b, c) if a == b else (max(epsilon, a), b, c) for (a, b, c) in three_percentages_of_excess]
- three_percentages_of_excess = [(a, b, c) if b == c else (a, b, min(100 - epsilon, c)) for (a, b, c) in three_percentages_of_excess]
+ three_percentages_of_excess = [(a, b, c) if a == b else (max(epsilon, a), b, c) for (a, b, c) in
+ three_percentages_of_excess]
+ three_percentages_of_excess = [(a, b, c) if b == c else (a, b, min(100 - epsilon, c)) for (a, b, c) in
+ three_percentages_of_excess]
y = [d[1] for d in three_percentages_of_excess]
yerr = np.array([[d[1] - d[0], d[2] - d[1]] for d in three_percentages_of_excess]).transpose()
label = ci_method_to_label[ci_method]
diff --git a/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py b/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
index f33fe1ef..c2d1a3c2 100644
--- a/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -117,6 +117,10 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
else:
return {m: EUROCODE_QUANTILE for m in self.massif_name_to_psnow.keys()}
+ @property
+ def massif_names_fitted(self):
+ return list(self.massif_name_to_years_and_maxima_for_model_fitting.keys())
+
@cached_property
def massif_name_to_years_and_maxima_for_model_fitting(self):
if self.fit_gev_only_on_non_null_maxima:
@@ -162,6 +166,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
all_trend_test = [t for t in all_trend_test
if acceptable_shape_parameter(t.unconstrained_estimator_gev_params.shape)]
sorted_trend_test = sorted(all_trend_test, key=lambda t: t.aic)
+
# Extract the stationary or non-stationary model that minimized AIC
trend_test_that_minimized_aic = sorted_trend_test[0]
massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
@@ -175,6 +180,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
[GumbelVersusGumbel, GumbelLocationTrendTest, GumbelScaleTrendTest,
GumbelLocationAndScaleTrendTest]][0]
massif_name_to_gumbel_trend_test_that_minimized_aic[massif_name] = gumbel_trend_test_that_minimized_aic
+
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
# Part 1 - Trends
@@ -342,9 +348,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
res = p.starmap(compute_eurocode_confidence_interval, arguments)
else:
res = [compute_eurocode_confidence_interval(*argument) for argument in arguments]
- massif_name_to_eurocode_return_level_uncertainty = dict(zip(self.uncertainty_massif_names, res))
+ massif_name_to_eurocode_return_level_uncertainty = dict(zip(massifs_names, res))
# For the rest of the massif names. Create a Eurocode Return Level Uncertainty as nan
- for massif_name in set(self.study.all_massif_names) - set(self.uncertainty_massif_names):
+ for massif_name in set(self.study.all_massif_names) - set(massifs_names):
massif_name_to_eurocode_return_level_uncertainty[massif_name] = self.default_eurocode_uncertainty
return massif_name_to_eurocode_return_level_uncertainty
--
GitLab