From b88a67cbaa011eefbd5c460fad20798603efd419 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 2 Dec 2019 12:02:25 +0100
Subject: [PATCH] [PAPER 1] add histogram plot
---
.../eurocode_data/main_eurocode_drawing.py | 2 +-
.../main_result_trends_and_return_levels.py | 14 +++--
...sualizer.py => plot_uncertainty_curves.py} | 15 +++--
.../plot_uncertainty_histogram.py | 56 +++++++++++++++++++
...dy_visualizer_for_non_stationary_trends.py | 41 ++++++++++----
.../confidence_interval_method.py | 5 ++
6 files changed, 113 insertions(+), 20 deletions(-)
rename experiment/paper_past_snow_loads/result_trends_and_return_levels/{eurocode_visualizer.py => plot_uncertainty_curves.py} (93%)
create mode 100644 experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
diff --git a/experiment/eurocode_data/main_eurocode_drawing.py b/experiment/eurocode_data/main_eurocode_drawing.py
index 4b19537e..ca616dcd 100644
--- a/experiment/eurocode_data/main_eurocode_drawing.py
+++ b/experiment/eurocode_data/main_eurocode_drawing.py
@@ -7,7 +7,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
StudyVisualizer
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
ConfidenceIntervalMethodFromExtremes
-from experiment.paper_past_snow_loads.result_trends_and_return_levels.eurocode_visualizer import \
+from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import \
plot_uncertainty_massifs, get_model_name
from experiment.eurocode_data.massif_name_to_departement import MASSIF_NAMES_ALPS
from experiment.eurocode_data.utils import EUROCODE_ALTITUDES
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index 93bd053f..3c6df869 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -5,9 +5,11 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
ALL_ALTITUDES_WITHOUT_NAN
from experiment.paper_past_snow_loads.paper_utils import paper_study_classes, paper_altitudes, \
load_altitude_to_visualizer
-from experiment.paper_past_snow_loads.result_trends_and_return_levels.eurocode_visualizer import \
+from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import \
plot_uncertainty_massifs
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_histogram import \
+ plot_uncertainty_histogram
from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
@@ -47,8 +49,8 @@ def intermediate_result(altitudes, massif_names=None,
visualizer.plot_trends(max_abs_tdrl)
# Plot graph
plot_uncertainty_massifs(altitude_to_visualizer)
- return altitude_to_visualizer
-
+ # Plot histogram
+ plot_uncertainty_histogram(altitude_to_visualizer)
def major_result():
@@ -62,8 +64,12 @@ def major_result():
if __name__ == '__main__':
# major_result()
+ intermediate_result(altitudes=[900, 1200], massif_names=['Chartreuse', 'Vanoise', 'Mercantour'],
+ uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+ ConfidenceIntervalMethodFromExtremes.ci_mle])
+ # intermediate_result(altitudes=[900, 1200], massif_names=None)
# intermediate_result(ALL_ALTITUDES_WITHOUT_NAN)
- intermediate_result(paper_altitudes)
+ # intermediate_result(paper_altitudes)
# minor_result(altitude=600)
# intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
similarity index 93%
rename from experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
rename to experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
index 5d6bbedf..bd9723a1 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
@@ -1,4 +1,5 @@
from typing import Dict
+import matplotlib.pyplot as plt
import numpy as np
@@ -9,6 +10,7 @@ from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extra
AbstractExtractEurocodeReturnLevel
from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color
from root_utils import get_display_name_from_object_type
@@ -16,7 +18,7 @@ def plot_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerFo
""" Plot several uncertainty plots
:return:
"""
- altitude_to_visualizer = {a:v for a,v in altitude_to_visualizer.items() if a in EUROCODE_ALTITUDES}
+ altitude_to_visualizer = {a: v for a, v in altitude_to_visualizer.items() if a in EUROCODE_ALTITUDES}
visualizer = list(altitude_to_visualizer.values())[-1]
# Subdivide massif names in group of 3
m = 1
@@ -49,6 +51,7 @@ def plot_subgroup_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVis
model_names_str = 'NonStationarity=' + '_'.join([str(e) for e in visualizer.non_stationary_contexts])
visualizer.plot_name = model_names_str + '_' + massif_names_str
visualizer.show_or_save_to_file(no_title=True)
+ plt.close()
def plot_single_uncertainty_massif(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
@@ -78,19 +81,19 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
""" Generic function that might be used by many other more global functions"""
altitudes = list(altitude_to_visualizer.keys())
visualizer = list(altitude_to_visualizer.values())[0]
- colors = ['tab:green', 'tab:brown'][::-1]
alpha = 0.2
# Display the EUROCODE return level
eurocode_region = massif_name_to_eurocode_region[massif_name]()
# Display the return level from model class
- for j, (color, uncertainty_method) in enumerate(zip(colors, visualizer.uncertainty_methods)):
+ for j, uncertainty_method in enumerate(visualizer.uncertainty_methods):
if j == 0:
# Plot eurocode norm
eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax,
altitudes=altitudes)
# Plot uncertainties
+ color = ci_method_to_color[uncertainty_method]
valid_altitudes = plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color,
massif_name, non_stationary_context, uncertainty_method)
@@ -118,7 +121,8 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes):
- visualizers = [v for a, v in altitude_to_visualizer.items() if a in valid_altitudes and massif_name in v.uncertainty_massif_names]
+ visualizers = [v for a, v in altitude_to_visualizer.items() if
+ a in valid_altitudes and massif_name in v.uncertainty_massif_names]
if len(visualizers) > 0:
tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in visualizers]
# Plot bars
@@ -155,5 +159,6 @@ def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitud
upper_bound = [r.confidence_interval[1] for r in ordered_return_level_uncertainties]
confidence_interval_str = ' {}'.format(AbstractExtractEurocodeReturnLevel.percentage_confidence_interval)
confidence_interval_str += '\% confidence interval'
- ax.fill_between(valid_altitudes, lower_bound, upper_bound, color=color, alpha=alpha, label=label_name + confidence_interval_str)
+ ax.fill_between(valid_altitudes, lower_bound, upper_bound, color=color, alpha=alpha,
+ label=label_name + confidence_interval_str)
return valid_altitudes
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
new file mode 100644
index 00000000..dbbfc221
--- /dev/null
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
@@ -0,0 +1,56 @@
+from typing import Dict
+import matplotlib.pyplot as plt
+import numpy as np
+
+from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR
+from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
+ StudyVisualizerForNonStationaryTrends
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color
+
+
+def plot_uncertainty_histogram(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+ """ Plot one graph for each non-stationary context
+ :return:
+ """
+ visualizer = list(altitude_to_visualizer.values())[0]
+ for non_stationary_context in visualizer.non_stationary_contexts:
+ plot_histogram(altitude_to_visualizer, non_stationary_context)
+
+
+def plot_histogram(altitude_to_visualizer, non_stationary_context):
+ """
+ Plot a single graph for potentially several confidence interval method
+ :param altitude_to_visualizer:
+ :param non_stationary_context:
+ :return:
+ """
+ visualizers = list(altitude_to_visualizer.values())
+ visualizer = visualizers[0]
+ ax = plt.gca()
+ altitudes = np.array(list(altitude_to_visualizer.keys()))
+ bincenters = altitudes
+ for j, ci_method in enumerate(visualizer.uncertainty_methods):
+ if len(visualizer.uncertainty_methods) == 2:
+ offset = -50 if j == 0 else 50
+ bincenters = altitudes + offset
+
+ width = 100
+ plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width=width)
+ ax.set_xticks(altitudes)
+ ax.set_ylabel('Massifs exceeding French standards (\%)')
+ ax.set_xlabel('Altitude (m)')
+ ax.set_ylim([0, 100])
+ ax.set_yticks([10 * i for i in range(11)])
+ # visualizer.plot_name = 'Percentage '
+ # visualizer.show = True
+ # visualizer.show_or_save_to_file(no_title=True)
+ plt.show()
+
+
+def plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width):
+ three_percentages_of_excess = [v.three_percentages_of_excess(ci_method, non_stationary_context) for v in
+ visualizers]
+ 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()
+ print(y, three_percentages_of_excess)
+ ax.bar(bincenters, y, width=width, color=ci_method_to_color[ci_method], yerr=yerr)
diff --git a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
index 39bc7ca2..dbc85bf4 100644
--- a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -6,6 +6,7 @@ from typing import Dict, Tuple
import numpy as np
from cached_property import cached_property
+from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
from experiment.eurocode_data.utils import EUROCODE_QUANTILE, EUROCODE_RETURN_LEVEL_STR
from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map, get_colors
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
@@ -117,15 +118,15 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
if max_abs_tdrl is not None:
self.global_max_abs_tdrl = max_abs_tdrl
ax = self.study.visualize_study(massif_name_to_value=self.massif_name_to_tdrl_value,
- replace_blue_by_white=False,
- axis_off=False, show_label=False,
- add_colorbar=True,
- massif_name_to_marker_style=self.massif_name_to_marker_style,
- massif_name_to_color=self.massif_name_to_tdrl_color,
- cmap=self.cmap,
- show=False,
- ticks_values_and_labels=self.ticks_values_and_labels,
- label=self.label_tdrl_bar + ' for {}'.format(EUROCODE_RETURN_LEVEL_STR))
+ replace_blue_by_white=False,
+ axis_off=False, show_label=False,
+ add_colorbar=True,
+ massif_name_to_marker_style=self.massif_name_to_marker_style,
+ massif_name_to_color=self.massif_name_to_tdrl_color,
+ cmap=self.cmap,
+ show=False,
+ ticks_values_and_labels=self.ticks_values_and_labels,
+ label=self.label_tdrl_bar + ' for {}'.format(EUROCODE_RETURN_LEVEL_STR))
ax.get_xaxis().set_visible(True)
ax.set_xticks([])
ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=12)
@@ -187,7 +188,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
return len(self.non_stationary_contexts)
def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_context) \
- -> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]:
+ -> Dict[str, EurocodeConfidenceIntervalFromExtremes]:
# Compute for the uncertainty massif names
arguments = [
[self.massif_name_to_non_null_years_and_maxima[m],
@@ -212,6 +213,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
@cached_property
def triplet_to_eurocode_uncertainty(self):
+ # -> Dict[(str, bool, str), EurocodeConfidenceIntervalFromExtremes]
d = {}
for ci_method in self.uncertainty_methods:
for non_stationary_uncertainty in self.non_stationary_contexts:
@@ -234,3 +236,22 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
else:
res = [compute_gelman_convergence_value(*argument) for argument in arguments]
return dict(zip(self.uncertainty_massif_names, res))
+
+ # Some values for the histogram
+
+ @cached_property
+ def massif_name_to_eurocode_values(self):
+ """Eurocode values for the altitude"""
+ return {m: r().lois_de_variation_de_la_valeur_caracteristique(altitude=self.study.altitude)
+ for m, r in massif_name_to_eurocode_region.items() if m in self.uncertainty_massif_names}
+
+ def three_percentages_of_excess(self, ci_method, non_stationary_context):
+ eurocode_and_uncertainties = [(self.massif_name_to_eurocode_values[massif_name],
+ self.triplet_to_eurocode_uncertainty[
+ (ci_method, non_stationary_context, massif_name)])
+ for massif_name in self.uncertainty_massif_names]
+ a = np.array([(uncertainty.confidence_interval[0] > eurocode,
+ uncertainty.mean_estimate > eurocode,
+ uncertainty.confidence_interval[1] > eurocode)
+ for eurocode, uncertainty in eurocode_and_uncertainties])
+ return 100 * np.mean(a, axis=0)
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
index 958dae9f..e28d30fb 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
@@ -16,3 +16,8 @@ ci_method_to_method_name = {
ConfidenceIntervalMethodFromExtremes.ci_boot: 'boot',
ConfidenceIntervalMethodFromExtremes.ci_proflik: 'proflik',
}
+
+ci_method_to_color = {
+ ConfidenceIntervalMethodFromExtremes.ci_mle: 'tab:brown',
+ ConfidenceIntervalMethodFromExtremes.my_bayes: 'tab:green'
+}
\ No newline at end of file
--
GitLab