From 7574a7481b69ad259da6cde34508fb93ed1b334a Mon Sep 17 00:00:00 2001 From: Le Roux Erwan <erwan.le-roux@irstea.fr> Date: Sat, 23 Nov 2019 19:02:36 +0100 Subject: [PATCH] [PAPER 1] major refactor for drawing uncertainty --- .../eurocode_data/main_eurocode_drawing.py | 4 +- .../scm_models_data/abstract_study.py | 2 + .../eurocode_visualizer.py | 141 +++++++++--------- .../main_result_trends_and_return_levels.py | 54 +++++-- ...dy_visualizer_for_non_stationary_trends.py | 56 +++++-- 5 files changed, 158 insertions(+), 99 deletions(-) diff --git a/experiment/eurocode_data/main_eurocode_drawing.py b/experiment/eurocode_data/main_eurocode_drawing.py index c3db7810..4b19537e 100644 --- a/experiment/eurocode_data/main_eurocode_drawing.py +++ b/experiment/eurocode_data/main_eurocode_drawing.py @@ -8,7 +8,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat 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 \ - plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict, get_model_name + 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 from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days, CrocusSweTotal @@ -73,7 +73,7 @@ def plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, tem model_name_to_massif_name_to_ordered_return_level.keys()} massif_name_to_model_name_to_ordered_return_level_uncertainties[massif_name] = d2 # Plot graph - plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict( + plot_uncertainty_massifs( massif_name_to_model_name_to_ordered_return_level_uncertainties, nb_massif_names=len(massif_names), nb_model_names=len(model_class_and_last_year)) if show: diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py index 792dd816..e194878b 100644 --- a/experiment/meteo_france_data/scm_models_data/abstract_study.py +++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py @@ -225,6 +225,7 @@ class AbstractStudy(object): @property def study_massif_names(self) -> List[str]: + # Massif names that are present in the current study (i.e. for the current altitude) return self.altitude_to_massif_names[self.altitude] @property @@ -490,6 +491,7 @@ class AbstractStudy(object): for massif_name in self.massif_name_to_altitudes.keys(): for altitude in self.massif_name_to_altitudes[massif_name]: altitude_to_massif_names[altitude].append(massif_name) + # massif_names are ordered in the same way as all_massif_names return altitude_to_massif_names """ Visualization methods """ 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 df17fadc..6b887c9c 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,6 +3,8 @@ from typing import Dict, List, Tuple import matplotlib.pyplot as plt import numpy as np +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.eurocode_return_level_uncertainties import \ EurocodeConfidenceIntervalFromExtremes from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region @@ -10,21 +12,6 @@ from experiment.meteo_france_data.scm_models_data.visualization.utils import cre from root_utils import get_display_name_from_object_type -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 '' - model_name = ' $ \widehat{z_p}(\\boldsymbol{\\theta_{\mathcal{M}_' - model_name += model_symbol - model_name += '}}' - model_name += parameter - model_name += ')_{ \\textrm{' + ci_method_name.upper().split(' ')[1] + '}} $ ' - return model_name - - -def get_model_name(model_class): - return get_display_name_from_object_type(model_class).split('Stationary')[0] + 'Stationary' - def massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, massif_names, uncertainty_methods, temporal_covariate, non_stationary_model): @@ -33,7 +20,7 @@ def massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, ma for altitude, visualizer in altitude_to_visualizer.items(): print('Processing altitude = {} '.format(altitude)) for ci_method in uncertainty_methods: - d = visualizer.massif_name_to_altitude_and_eurocode_level_uncertainty_for_minimized_aic_model_class( + 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 @@ -45,78 +32,88 @@ def massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, ma return massif_name_to_ordered_eurocode_level_uncertainty -def plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(altitude_to_visualizer, - massif_names, - non_stationary_models_for_uncertainty, - uncertainty_methods): - """ - Rows correspond to massif names - Columns correspond to stationary/non stationary model name for a given date - Uncertainty result_trends_and_return_levels correpsond to the different plot on the graph +def plot_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]): + """ Plot several uncertainty plots :return: """ - # Compute the dictionary of interest - # Plot uncertainties - model_name_to_massif_name_to_ordered_return_level = {} - for non_stationary_model in non_stationary_models_for_uncertainty: - d = massif_name_to_ordered_return_level_uncertainties(altitude_to_visualizer, massif_names, - uncertainty_methods, - temporal_covariate=2017, - non_stationary_model=non_stationary_model) - model_name_to_massif_name_to_ordered_return_level[non_stationary_model] = d - - # Transform the dictionary into the desired format - d = {} - for massif_name in massif_names: - d2 = {model_name: model_name_to_massif_name_to_ordered_return_level[model_name][massif_name] for model_name - in - model_name_to_massif_name_to_ordered_return_level.keys()} - d[massif_name] = d2 - + visualizer = list(altitude_to_visualizer.values())[0] + # Subdivide massif names in group of 5 + m = 4 + uncertainty_massif_names = visualizer.uncertainty_massif_names + n = (len(uncertainty_massif_names) // m) + 1 + for i in list(range(n))[:]: + massif_names = uncertainty_massif_names[m * i: m * (i + 1)] + print(massif_names) + plot_subgroup_uncertainty_massifs(altitude_to_visualizer, massif_names) + + +def plot_subgroup_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends], + massif_names): + """Create a plot with a maximum of 4 massif names + We will save the plot at this level + """ + visualizer = list(altitude_to_visualizer.values())[0] nb_massif_names = len(massif_names) - nb_model_names = len(non_stationary_models_for_uncertainty) - axes = create_adjusted_axes(nb_massif_names, nb_model_names) + assert nb_massif_names <= 5 + axes = create_adjusted_axes(nb_massif_names, visualizer.nb_contexts) if nb_massif_names == 1: axes = [axes] - for ax, (massif_name, model_name_to_uncertainty_level) in zip(axes, d.items()): - plot_model_name_to_uncertainty_method_to_ordered_dict(model_name_to_uncertainty_level, - massif_name, ax) + for ax, massif_name in zip(axes, massif_names): + plot_single_uncertainty_massif(altitude_to_visualizer, + massif_name, ax) # Save plot - visualizer = list(altitude_to_visualizer.values())[0] massif_names_str = '_'.join(massif_names) - model_names_str = 'NonStationarity=' + '_'.join([str(e) for e in non_stationary_models_for_uncertainty]) + 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.suptitle('50-year return levels of extreme snow loads in France for several confiance interval methods.') -def plot_model_name_to_uncertainty_method_to_ordered_dict(d, massif_name, axes): - if len(d) == 1: +def plot_single_uncertainty_massif(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends], + massif_name, axes): + visualizer = list(altitude_to_visualizer.values())[0] + if visualizer.nb_contexts == 1: axes = [axes] - for ax, (model_name, uncertainty_method_to_ordered_dict) in zip(axes, d.items()): - plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name, - uncertainty_method_to_ordered_dict) + for ax, non_stationary_context in zip(axes, visualizer.non_stationary_contexts): + plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, non_stationary_context, + altitude_to_visualizer) -def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name, - label_to_ordered_return_level_uncertainties: - Dict[str, List[ - EurocodeConfidenceIntervalFromExtremes]]): +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 '' + model_name = ' $ \widehat{z_p}(\\boldsymbol{\\theta_{\mathcal{M}_' + model_name += model_symbol + model_name += '}}' + model_name += parameter + model_name += ')_{ \\textrm{' + ci_method_name.upper().split(' ')[1] + '}} $ ' + return model_name + + +def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, non_stationary_context, + altitude_to_visualizer: Dict[ + int, StudyVisualizerForNonStationaryTrends]): """ 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:olive'] 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, (label, l)) in enumerate(zip(colors, label_to_ordered_return_level_uncertainties.items())): - l = list(zip(*l)) - altitudes = l[0] - ordered_return_level_uncertaines = l[1] # type: List[EurocodeConfidenceIntervalFromExtremes] + for j, (color, uncertainty_method) in enumerate(zip(colors, visualizer.uncertainty_methods)): # Plot eurocode standards only for the first loop if j == 0: - eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax, altitudes=altitudes) + eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax, + altitudes=altitudes) + # Compute ordered_return_level_uncertaines for a given massif_name, uncertainty methods, and non stationary context + ordered_return_level_uncertaines = [] + for visualizer in altitude_to_visualizer.values(): + u = visualizer.triplet_to_eurocode_uncertainty[(uncertainty_method, non_stationary_context, massif_name)] + ordered_return_level_uncertaines.append(u) + # Display mean = [r.mean_estimate for r in ordered_return_level_uncertaines] # Filter and keep only non nan values not_nan_index = [not np.isnan(m) for m in mean] @@ -124,8 +121,9 @@ def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name altitudes = list(np.array(altitudes)[not_nan_index]) ordered_return_level_uncertaines = list(np.array(ordered_return_level_uncertaines)[not_nan_index]) - ci_method_name = str(label).split('.')[1].replace('_', ' ') - ax.plot(altitudes, mean, linestyle='--', marker='o', color=color, label=get_label_name(model_name, ci_method_name)) + ci_method_name = str(uncertainty_method).split('.')[1].replace('_', ' ') + ax.plot(altitudes, mean, linestyle='--', marker='o', color=color, + label=get_label_name(non_stationary_context, ci_method_name)) lower_bound = [r.confidence_interval[0] for r in ordered_return_level_uncertaines] upper_bound = [r.confidence_interval[1] for r in ordered_return_level_uncertaines] ax.fill_between(altitudes, lower_bound, upper_bound, color=color, alpha=alpha) @@ -133,14 +131,15 @@ def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name ax.set_ylim([0.0, 16]) massif_name_str = massif_name.replace('_', ' ') eurocode_region_str = get_display_name_from_object_type(type(eurocode_region)) - is_non_stationary_model = model_name if isinstance(model_name, bool) else 'Non' in model_name + is_non_stationary_model = non_stationary_context if isinstance(non_stationary_context, + bool) else 'Non' in non_stationary_context if is_non_stationary_model: - model_name = 'non-stationary' + non_stationary_context = 'non-stationary' else: - model_name = 'stationary' - title = '{} ({} Eurocodes area) with a {} model'.format(massif_name_str, eurocode_region_str, model_name) + non_stationary_context = 'stationary' + title = '{} ({} Eurocodes area) with a {} model'.format(massif_name_str, eurocode_region_str, + non_stationary_context) ax.set_title(title) ax.set_ylabel('50-year return level of SL (kN $m^-2$)') ax.set_xlabel('Altitude (m)') ax.grid() - 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 98e7beb7..172b0396 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 @@ -7,7 +7,7 @@ import matplotlib.pyplot as plt from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \ StudyVisualizer from experiment.paper_past_snow_loads.result_trends_and_return_levels.eurocode_visualizer import \ - plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict + 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.study_visualizer_for_non_stationary_trends import \ StudyVisualizerForNonStationaryTrends @@ -19,33 +19,57 @@ mpl.rcParams['text.usetex'] = True mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}'] -def draw_snow_load_map(altitude): +def minor_result(altitude): + """Plot trends for a single altitude to be fast""" visualizer = StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude), multiprocessing=True) visualizer.plot_trends() -def main_results(): - altitudes = [[1500, 1800]][0] - uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes, - ConfidenceIntervalMethodFromExtremes.ci_mle][1:] - massif_names = ['Chartreuse'] - non_stationary_models_for_uncertainty = [False, True][:1] +def intermediate_result(altitudes, massif_names=None, + non_stationary_uncertainty=None, uncertainty_methods=None, + study_class=CrocusSnowLoadTotal): + """ + Plot all the trends for all altitudes + And enable to plot uncertainty plot for some specific massif_names, uncertainty methods to be fast + :param altitudes: + :param massif_names: + :param non_stationary_uncertainty: + :param uncertainty_methods: + :param study_class: + :return: + """ # Load altitude to visualizer altitude_to_visualizer = OrderedDict() for altitude in altitudes: altitude_to_visualizer[altitude] = StudyVisualizerForNonStationaryTrends( - study=CrocusSnowLoadTotal(altitude=altitude), multiprocessing=True, save_to_file=True) + study=study_class(altitude=altitude), multiprocessing=True, save_to_file=True, + uncertainty_massif_names=massif_names, uncertainty_methods=uncertainty_methods, + non_stationary_contexts=non_stationary_uncertainty) # Plot trends max_abs_tdrl = max([visualizer.max_abs_tdrl for visualizer in altitude_to_visualizer.values()]) for visualizer in altitude_to_visualizer.values(): visualizer.plot_trends(max_abs_tdrl) # Plot graph - plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(altitude_to_visualizer, - massif_names, - non_stationary_models_for_uncertainty, - uncertainty_methods) + plot_uncertainty_massifs(altitude_to_visualizer) + return altitude_to_visualizer + + +def major_result(): + altitudes = [[1500, 1800]][0] + uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes, + ConfidenceIntervalMethodFromExtremes.ci_mle][1:] + massif_names = ['Chartreuse'] + non_stationary_models_for_uncertainty = [False, True][:1] + # + # altitudes + # study_class if __name__ == '__main__': - # draw_snow_load_map(altitude=1800) - main_results() + # minor_result(altitude=1800) + intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'], + 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], + non_stationary_uncertainty=[False]) 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 429ef994..e4f21054 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 @@ -15,6 +15,8 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter impo GevLocationTrendTest from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel +from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \ + ConfidenceIntervalMethodFromExtremes from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \ compute_eurocode_confidence_interval, EurocodeConfidenceIntervalFromExtremes from root_utils import NB_CORES @@ -26,11 +28,29 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer): vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False, temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False, complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True, - score_class=MeanScore): + score_class=MeanScore, + uncertainty_methods=None, + non_stationary_contexts=None, + uncertainty_massif_names=None, + effective_temporal_covariate=2017): super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot, year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity, transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis, normalization_under_one_observations, score_class) + # Add some attributes + self.effective_temporal_covariate = effective_temporal_covariate + self.non_stationary_contexts = non_stationary_contexts + self.uncertainty_methods = uncertainty_methods + self.uncertainty_massif_names = uncertainty_massif_names + # Assign some default arguments + if self.non_stationary_contexts is None: + self.non_stationary_contexts = [False, True][:1] + if self.uncertainty_methods is None: + self.uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes, + ConfidenceIntervalMethodFromExtremes.ci_mle][1:] + if self.uncertainty_massif_names is None: + self.uncertainty_massif_names = self.study.study_massif_names + # Assign default argument for the non stationary trends self.non_stationary_trend_test = [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest] self.non_stationary_trend_test_to_marker = dict(zip(self.non_stationary_trend_test, ["s", "^", "D"])) @@ -103,7 +123,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer): 'fillstyle': 'full' if t.is_significant else 'none'} return d - # Part 1 - Uncertainty return level plot + # Part 2 - Uncertainty return level plot def massif_name_to_model_class(self, massif_name, non_stationary_model): if not non_stationary_model: @@ -111,23 +131,37 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer): else: return self.massif_name_to_minimized_aic_non_stationary_trend_test[massif_name].unconstrained_model_class - def massif_name_to_uncertainty(self): - pass + @property + def nb_contexts(self): + return len(self.non_stationary_contexts) + + @property + def nb_uncertainty_method(self): + return len(self.uncertainty_methods) - def massif_name_to_altitude_and_eurocode_level_uncertainty_for_minimized_aic_model_class(self, massif_names, - ci_method, - effective_temporal_covariate, - non_stationary_model) \ + def massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_model) \ -> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]: arguments = [ [self.massif_name_to_non_null_years_and_maxima[m], self.massif_name_to_model_class(m, non_stationary_model), - ci_method, effective_temporal_covariate] for m in massif_names] + ci_method, self.effective_temporal_covariate] for m in self.uncertainty_massif_names] if self.multiprocessing: with Pool(NB_CORES) as p: res = p.starmap(compute_eurocode_confidence_interval, arguments) else: res = [compute_eurocode_confidence_interval(*argument) for argument in arguments] - altitudes_and_res = [(self.study.altitude, r) for r in res] - massif_name_to_eurocode_return_level_uncertainty = OrderedDict(zip(massif_names, altitudes_and_res)) + massif_name_to_eurocode_return_level_uncertainty = OrderedDict(zip(self.uncertainty_massif_names, res)) return massif_name_to_eurocode_return_level_uncertainty + + @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( + ci_method, non_stationary_uncertainty).items(): + d[(ci_method, non_stationary_uncertainty, uncertainty_massif_name)] = eurocode_uncertainty + return d + + def model_name_to_uncertainty_method_to_ratio_above_eurocode(self): + assert self.uncertainty_massif_names == self.study.study_massif_names -- GitLab