diff --git a/experiment/eurocode_data/utils.py b/experiment/eurocode_data/utils.py
index 764ddbd093ce3f916bd9a0e59218433dc562ae42..7385d7ff77bc013a2d5707edc9669ecf6989a16e 100644
--- a/experiment/eurocode_data/utils.py
+++ b/experiment/eurocode_data/utils.py
@@ -1,4 +1,6 @@
+# Eurocode quantile str
+EUROCODE_RETURN_LEVEL_STR = '50-year return level of SL (kN $m^-2$)'
# Eurocode quantile correspond to a 50 year return period
EUROCODE_QUANTILE = 0.98
# Altitudes (between low and mid altitudes) < 2000m and should be > 200m
diff --git a/experiment/paper_past_snow_loads/data/main_eurocode_plot.py b/experiment/paper_past_snow_loads/data/main_eurocode_plot.py
index 693a400bb118d5678e22dda103c6c43c556eb877..5aee137de3c019b0ff247cf468de5d82f2b56355 100644
--- a/experiment/paper_past_snow_loads/data/main_eurocode_plot.py
+++ b/experiment/paper_past_snow_loads/data/main_eurocode_plot.py
@@ -3,10 +3,13 @@ import numpy as np
from experiment.eurocode_data.eurocode_region import C2, E, C1
from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
+from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from root_utils import get_display_name_from_object_type
+
+
def main_eurocode_norms(ax=None):
if ax is None:
ax = plt.gca()
@@ -20,7 +23,7 @@ def main_eurocode_norms(ax=None):
ax.legend()
ax.xaxis.set_ticks([250 * i for i in range(1, 9)])
ax.tick_params(axis='both', which='major', labelsize=13)
- ax.set_ylabel('50-year return level of SL (kN $m^-2$)')
+ ax.set_ylabel(EUROCODE_RETURN_LEVEL_STR)
ax.set_xlabel('Altitude (m)')
ax.set_ylim([0.0, 11.0])
ax.grid()
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/eurocode_visualizer.py
index 1c917c6ef131613c3dc39e336a56c7f2f65b93d0..1a01f1dbf4b1a89d2a51584e7b45202b4874a67a 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/eurocode_visualizer.py
@@ -3,8 +3,11 @@ from typing import Dict, List, Tuple
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.result_trends_and_return_levels.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+ AbstractExtractEurocodeReturnLevel
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
EurocodeConfidenceIntervalFromExtremes
from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
@@ -12,31 +15,11 @@ from experiment.meteo_france_data.scm_models_data.visualization.utils import cre
from root_utils import get_display_name_from_object_type
-def massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, massif_names,
- uncertainty_methods, temporal_covariate,
- non_stationary_model):
- massif_name_to_ordered_eurocode_level_uncertainty = {
- massif_name: {ci_method: [] for ci_method in uncertainty_methods} for massif_name in massif_names}
- for altitude, visualizer in altitude_to_visualizer.items():
- print('Processing altitude = {} '.format(altitude))
- for ci_method in uncertainty_methods:
- d = visualizer.massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(
- massif_names, ci_method,
- temporal_covariate, non_stationary_model)
- # Append the altitude one by one
- for massif_name, return_level_uncertainty in d.items():
- print(massif_name, return_level_uncertainty[0], return_level_uncertainty[1].confidence_interval,
- return_level_uncertainty[1].mean_estimate)
- massif_name_to_ordered_eurocode_level_uncertainty[massif_name][ci_method].append(
- return_level_uncertainty)
- return massif_name_to_ordered_eurocode_level_uncertainty
-
-
def plot_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
""" Plot several uncertainty plots
:return:
"""
- visualizer = list(altitude_to_visualizer.values())[0]
+ visualizer = list(altitude_to_visualizer.values())[-1]
# Subdivide massif names in group of 3
m = 1
uncertainty_massif_names = visualizer.uncertainty_massif_names
@@ -80,10 +63,9 @@ def plot_single_uncertainty_massif(altitude_to_visualizer: Dict[int, StudyVisual
altitude_to_visualizer)
-def get_label_name(model_name, ci_method_name: str):
- is_non_stationary = model_name == 'NonStationary'
- model_symbol = 'N' if is_non_stationary else '0'
- parameter = ', 2017' if is_non_stationary else ''
+def get_label_name(non_stationary_context, ci_method_name):
+ model_symbol = 'N' if non_stationary_context else '0'
+ parameter = ', 2017' if non_stationary_context else ''
model_name = ' $ \widehat{z_p}(\\boldsymbol{\\theta_{\mathcal{M}_'
model_name += model_symbol
model_name += '}}'
@@ -105,26 +87,18 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
# Display the return level from model class
for j, (color, uncertainty_method) in enumerate(zip(colors, visualizer.uncertainty_methods)):
- # Plot eurocode standards only for the first loop
if j == 0:
# Plot eurocode norm
eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax,
altitudes=altitudes)
- # Plot bars of TDRL only in the non stationary case
- if non_stationary_context:
- tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in altitude_to_visualizer.values()]
- # Plot bars
- colors = [v.massif_name_to_tdrl_color[massif_name] for v in altitude_to_visualizer.values()]
- ax.bar(altitudes, tdrl_values, width=150, color=colors)
- # Plot markers
- markers_kwargs = [v.massif_name_to_marker_style(markersize=4)[massif_name]
- for v, tdrl_value in zip(altitude_to_visualizer.values(), tdrl_values)]
- for altitude, marker_kwargs, value in zip(altitudes, markers_kwargs, tdrl_values):
- ax.plot([altitude], [value / 2], **marker_kwargs)
# Plot uncertainties
- plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
- non_stationary_context, uncertainty_method)
+ valid_altitudes = plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color,
+ massif_name, non_stationary_context, uncertainty_method)
+
+ # Plot bars of TDRL only in the non stationary case
+ if j == 0 and non_stationary_context:
+ plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes)
ax.legend(loc=2)
ax.set_ylim([-1, 16])
@@ -140,11 +114,27 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
non_stationary_context)
ax.set_title(title)
ax.set_xticks(altitudes)
- ax.set_ylabel('50-year return level of SL (kN $m^-2$)')
+ ax.set_ylabel(EUROCODE_RETURN_LEVEL_STR)
ax.set_xlabel('Altitude (m)')
ax.grid()
+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]
+ if len(visualizers) > 0:
+ tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in visualizers]
+ # Plot bars
+ colors = [v.massif_name_to_tdrl_color[massif_name] for v in visualizers]
+ ax.bar(valid_altitudes, tdrl_values, width=150, color=colors, label=visualizers[0].label_tdrl_bar,
+ edgecolor='black', hatch='//')
+ # Plot markers
+ markers_kwargs = [v.massif_name_to_marker_style[massif_name] for v in visualizers]
+ for altitude, marker_kwargs, value in zip(valid_altitudes, markers_kwargs, tdrl_values):
+ # ax.plot([altitude], [value / 2], **marker_kwargs)
+ # Better to plot all the markers on the same line
+ ax.plot([altitude], 0, **marker_kwargs)
+
+
def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
non_stationary_context, uncertainty_method):
# Compute ordered_return_level_uncertaines for a given massif_name, uncertainty methods, and non stationary context
@@ -160,8 +150,12 @@ def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitud
valid_altitudes = list(np.array(altitudes)[not_nan_index])
ordered_return_level_uncertainties = list(np.array(ordered_return_level_uncertainties)[not_nan_index])
ci_method_name = str(uncertainty_method).split('.')[1].replace('_', ' ')
+ label_name = get_label_name(non_stationary_context, ci_method_name)
ax.plot(valid_altitudes, mean, linestyle='--', marker='o', color=color,
- label=get_label_name(non_stationary_context, ci_method_name))
+ label=label_name)
lower_bound = [r.confidence_interval[0] for r in ordered_return_level_uncertainties]
upper_bound = [r.confidence_interval[1] for r in ordered_return_level_uncertainties]
- ax.fill_between(valid_altitudes, lower_bound, upper_bound, color=color, alpha=alpha)
+ 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)
+ return valid_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 70ad497e5548edf73876258cc328d883165ad8a3..41301d5233ce5fd5a8be0bf714b4c8282fde92fd 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
@@ -79,7 +79,10 @@ if __name__ == '__main__':
# intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
# non_stationary_uncertainty=[False])
- intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
- uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
- ConfidenceIntervalMethodFromExtremes.ci_bayes],
+ # intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
+ # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
+ # ConfidenceIntervalMethodFromExtremes.ci_bayes],
+ # non_stationary_uncertainty=[False, True])
+ intermediate_result(altitudes=[300, 600, 900], massif_names=None,
+ uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
non_stationary_uncertainty=[False, True])
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
index 0bcacfd00bad110ab15edbbf79fc383580df1b10..a418c2b50a1ce91f0e29773e4b9c677c138944a9 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
@@ -6,7 +6,7 @@ from typing import Dict, Tuple
import numpy as np
from cached_property import cached_property
-from experiment.eurocode_data.utils import EUROCODE_QUANTILE
+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
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
@@ -117,19 +117,19 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
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_marker_style=self.massif_name_to_marker_style,
massif_name_to_color=self.massif_name_to_tdrl_color,
cmap=self.cmap,
show=self.show,
ticks_values_and_labels=self.ticks_values_and_labels,
- label=self.label_trend_bar)
+ label=self.label_tdrl_bar + ' for {}'.format(EUROCODE_RETURN_LEVEL_STR))
self.plot_name = 'tdlr_trends'
self.show_or_save_to_file(add_classic_title=False, tight_layout=False, no_title=True)
plt.close()
@property
- def label_trend_bar(self):
- return 'Change in {} years for Eurocode quantile of SL (kN $m^-2$)'.format(AbstractGevTrendTest.nb_years_for_quantile_evolution)
+ def label_tdrl_bar(self):
+ return 'Change in {} years'.format(AbstractGevTrendTest.nb_years_for_quantile_evolution)
@property
def ticks_values_and_labels(self):
@@ -140,7 +140,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
positive_ticks.append(round(tick, 1))
tick += graduation
all_ticks_labels = [-t for t in positive_ticks] + [0] + positive_ticks
- ticks_values = [(t / 2 * self._max_abs_tdrl) + 0.5 for t in all_ticks_labels]
+ ticks_values = [((t / self._max_abs_tdrl) + 1) / 2 for t in all_ticks_labels]
return ticks_values, all_ticks_labels
@cached_property
@@ -157,12 +157,13 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
return {m: get_colors([v], self.cmap, -self._max_abs_tdrl, self._max_abs_tdrl)[0]
for m, v in self.massif_name_to_tdrl_value.items()}
- def massif_name_to_marker_style(self, markersize=5):
+ @cached_property
+ def massif_name_to_marker_style(self):
d = {}
for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items():
d[m] = {'marker': self.non_stationary_trend_test_to_marker[type(t)],
'color': 'k',
- 'markersize': markersize,
+ 'markersize': 5,
'fillstyle': 'full' if t.is_significant else 'none'}
return d
@@ -182,8 +183,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
def nb_uncertainty_method(self):
return len(self.uncertainty_methods)
- def massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_model) \
+ def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_model) \
-> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]:
+ # Compute for the uncertainty massif names
arguments = [
[self.massif_name_to_non_null_years_and_maxima[m],
self.massif_name_to_model_class(m, non_stationary_model),
@@ -195,17 +197,24 @@ 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 = OrderedDict(zip(self.uncertainty_massif_names, res))
+ massif_name_to_eurocode_return_level_uncertainty = dict(zip(self.uncertainty_massif_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):
+ massif_name_to_eurocode_return_level_uncertainty[massif_name] = self.default_eurocode_uncertainty
return massif_name_to_eurocode_return_level_uncertainty
+ @cached_property
+ def default_eurocode_uncertainty(self):
+ return EurocodeConfidenceIntervalFromExtremes(mean_estimate=np.nan, confidence_interval=(np.nan, np.nan))
+
@cached_property
def triplet_to_eurocode_uncertainty(self):
d = {}
for ci_method in self.uncertainty_methods:
for non_stationary_uncertainty in self.non_stationary_contexts:
- for uncertainty_massif_name, eurocode_uncertainty in self.massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(
+ for massif_name, eurocode_uncertainty in self.all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(
ci_method, non_stationary_uncertainty).items():
- d[(ci_method, non_stationary_uncertainty, uncertainty_massif_name)] = eurocode_uncertainty
+ d[(ci_method, non_stationary_uncertainty, massif_name)] = eurocode_uncertainty
return d
def model_name_to_uncertainty_method_to_ratio_above_eurocode(self):
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
index d02ef69da73761553f6f1bf2ebb9809a3815781a..559dae4f8528f162ee2c6ba38905fe911bd1f26b 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
@@ -11,6 +11,7 @@ from extreme_fit.model.margin_model.margin_function.linear_margin_function impor
from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_bayesian_extremes import \
ResultFromBayesianExtremes
from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_mle_extremes import ResultFromMleExtremes
+from root_utils import classproperty
class AbstractExtractEurocodeReturnLevel(object):
@@ -22,8 +23,11 @@ class AbstractExtractEurocodeReturnLevel(object):
self.result_from_fit = self.estimator.result_from_model_fit
self.temporal_covariate = temporal_covariate
# Fixed Parameters
- self.eurocode_quantile_level = quantile_level
- self.alpha_for_confidence_interval = self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY
+ self.quantile_level = quantile_level
+
+ @classproperty
+ def percentage_confidence_interval(cls) -> int:
+ return int(100 * (1 - cls.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY))
@property
def mean_estimate(self) -> np.ndarray:
@@ -43,8 +47,8 @@ class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel)
@cached_property
def confidence_interval_method(self):
- return self.result_from_fit.confidence_interval_method(self.eurocode_quantile_level,
- self.alpha_for_confidence_interval,
+ return self.result_from_fit.confidence_interval_method(self.quantile_level,
+ self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY,
self.transformed_temporal_covariate,
self.ci_method)
@@ -81,7 +85,7 @@ class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeRe
@cached_property
def posterior_eurocode_return_level_samples_for_temporal_covariate(self) -> np.ndarray:
return np.array(
- [p.quantile(self.eurocode_quantile_level) for p in self.gev_params_from_fit_for_temporal_covariate])
+ [p.quantile(self.quantile_level) for p in self.gev_params_from_fit_for_temporal_covariate])
@property
def mean_estimate(self) -> np.ndarray:
@@ -91,7 +95,7 @@ class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeRe
@property
def confidence_interval(self):
# Bottom and upper quantile correspond to the quantile
- bottom_quantile = self.alpha_for_confidence_interval / 2
+ bottom_quantile = self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY / 2
bottom_and_upper_quantile = (bottom_quantile, 1 - bottom_quantile)
return [np.quantile(self.posterior_eurocode_return_level_samples_for_temporal_covariate, q=q)
for q in bottom_and_upper_quantile]