From d89b6e0329e51e5b85cb2188f9d2af739bee4fb6 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 18 Jun 2020 15:58:44 +0200
Subject: [PATCH] [contrasting] add plot selection for paper 2
---
extreme_trend/abstract_gev_trend_test.py | 12 ++-
.../main_snowfall_article.py | 12 ++-
.../plot_selection_curves_paper2.py | 73 +++++++++++++++++++
.../snowfall_plot.py | 10 +--
.../study_visualizer_for_mean_values.py | 30 +++++++-
5 files changed, 120 insertions(+), 17 deletions(-)
create mode 100644 projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/plot_selection_curves_paper2.py
diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index e205f35a..dc078514 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/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):
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 d46f3ec3..301ea5e6 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
@@ -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)
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/plot_selection_curves_paper2.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/plot_selection_curves_paper2.py
new file mode 100644
index 00000000..8e54d150
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/plot_selection_curves_paper2.py
@@ -0,0 +1,73 @@
+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()
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py
index fdc65301..e98023f8 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py
+++ b/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:
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
index 2c33eb43..d1a30c5c 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
+++ b/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):
--
GitLab