From 07091db861e41f48cb8b50b3ec5da137493dcfa9 Mon Sep 17 00:00:00 2001 From: Le Roux Erwan <erwan.le-roux@irstea.fr> Date: Wed, 22 May 2019 16:22:51 +0200 Subject: [PATCH] [SCM][HYPERCUBE] add trend_strength to the hypercube_visualization. --- .../abstract_hypercube_visualizer.py | 22 ++-- .../altitude_hypercube_visualizer.py | 111 +++++++++--------- .../main_hypercube_visualization.py | 19 +-- .../quantity_altitude_visualizer.py | 5 +- .../study_visualization/study_visualizer.py | 29 +++-- .../abstract_gev_trend_test.py | 19 +-- .../abstract_trend_test.py | 13 ++ .../margin_function/linear_margin_function.py | 2 +- 8 files changed, 121 insertions(+), 99 deletions(-) diff --git a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/abstract_hypercube_visualizer.py b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/abstract_hypercube_visualizer.py index 13524c33..32ff017e 100644 --- a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/abstract_hypercube_visualizer.py +++ b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/abstract_hypercube_visualizer.py @@ -41,21 +41,21 @@ class AbstractHypercubeVisualizer(object): def tuple_values(self, idx): return sorted(set([t[idx] if isinstance(t, tuple) else t for t in self.tuple_to_study_visualizer.keys()])) + @cached_property + def df_trends_spatio_temporal(self): + return [study_visualizer.df_trend_spatio_temporal(self.trend_test_class, self.starting_years, + self.nb_data_for_fast_mode) + for study_visualizer in self.tuple_to_study_visualizer.values()] + @cached_property def df_hypercube_trend_type(self) -> pd.DataFrame: - df_spatio_temporal_trend_types = [ - study_visualizer.df_trend_spatio_temporal(self.trend_test_class, self.starting_years, - self.nb_data_for_fast_mode) - for study_visualizer in self.tuple_to_study_visualizer.values()] + df_spatio_temporal_trend_types = [e[0] for e in self.df_trends_spatio_temporal] return pd.concat(df_spatio_temporal_trend_types, keys=list(self.tuple_to_study_visualizer.keys()), axis=0) @cached_property def df_hypercube_trend_strength(self) -> pd.DataFrame: - df_spatio_temporal_trend_types = [ - study_visualizer.df_trend_spatio_temporal(self.trend_test_class, self.starting_years, - self.nb_data_for_fast_mode) - for study_visualizer in self.tuple_to_study_visualizer.values()] - return pd.concat(df_spatio_temporal_trend_types, keys=list(self.tuple_to_study_visualizer.keys()), axis=0) + df_spatio_temporal_trend_strength = [e[1] for e in self.df_trends_spatio_temporal] + return pd.concat(df_spatio_temporal_trend_strength, keys=list(self.tuple_to_study_visualizer.keys()), axis=0) # Some properties @@ -89,7 +89,3 @@ class AbstractHypercubeVisualizer(object): # Load uniform weights by default uniform_weight = 1 / len(self.starting_years) return {year: uniform_weight for year in self.starting_years} - - - - diff --git a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/altitude_hypercube_visualizer.py b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/altitude_hypercube_visualizer.py index 5c692a91..89b9943a 100644 --- a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/altitude_hypercube_visualizer.py +++ b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/altitude_hypercube_visualizer.py @@ -22,64 +22,66 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer): def trend_types(self): return self.trend_type_to_style.keys() - def trend_type_to_s_percentages(self, reduction_function): + def trend_type_to_series(self, reduction_function): # Map each trend type to its serie with percentages - trend_type_to_s_percentages = {} + trend_type_to_series = {} for trend_type in self.trend_types: - df_bool = (self.df_hypercube_trend_type == trend_type) - # Reduce the entire dataframe to a serie - s_percentages = reduction_function(df_bool) - assert isinstance(s_percentages, pd.Series) - assert not isinstance(s_percentages.index, pd.MultiIndex) - s_percentages *= 100 - trend_type_to_s_percentages[trend_type] = s_percentages - # Post processing - Add the significant trend into the count of normal trend - trend_type_to_s_percentages[AbstractTrendTest.POSITIVE_TREND] += trend_type_to_s_percentages[ - AbstractTrendTest.SIGNIFICATIVE_POSITIVE_TREND] - trend_type_to_s_percentages[AbstractTrendTest.NEGATIVE_TREND] += trend_type_to_s_percentages[ - AbstractTrendTest.SIGNIFICATIVE_NEGATIVE_TREND] - return trend_type_to_s_percentages - - def subtitle_to_reduction_function(self, reduction_function, level=None, add_detailed_plot=False): + # Reduce df_bool df to a serie s_trend_type_percentage + df_bool = self.df_hypercube_trend_type.isin(AbstractTrendTest.get_trend_types(trend_type)) + s_trend_type_percentage = reduction_function(df_bool) + assert isinstance(s_trend_type_percentage, pd.Series) + assert not isinstance(s_trend_type_percentage.index, pd.MultiIndex) + s_trend_type_percentage *= 100 + # Reduce df_strength to a serie s_trend_strength + df_strength = self.df_hypercube_trend_strength[df_bool] + s_trend_strength = reduction_function(df_strength) + # Store results + trend_type_to_series[trend_type] = (s_trend_type_percentage, s_trend_strength) + return trend_type_to_series + + def subtitle_to_reduction_function(self, reduction_function, level=None, add_detailed_plot=False, subtitle=None): def reduction_function_with_level(df_bool): return reduction_function(df_bool) if level is None else reduction_function(df_bool, level) + if subtitle is None: + subtitle = self.study.variable_name + return {subtitle: reduction_function_with_level} - return {'global': reduction_function_with_level} - - def visualize_trend_test_evolution(self, reduction_function, xlabel, xlabel_values, ax=None, marker='o', + def visualize_trend_test_evolution(self, reduction_function, xlabel, xlabel_values, axes=None, marker='o', subtitle=''): - if ax is None: - fig, ax = plt.subplots(1, 1, figsize=self.study_visualizer.figsize) - - trend_type_to_percentages_values = {k: s.values for k, s in - self.trend_type_to_s_percentages(reduction_function).items()} - for trend_type in self.trend_types: - style = self.trend_type_to_style[trend_type] - percentages_values = trend_type_to_percentages_values[trend_type] - ax.plot(xlabel_values, percentages_values, style + marker, label=trend_type) - - # Plot the total value of significative values - significative_values = trend_type_to_percentages_values[AbstractTrendTest.SIGNIFICATIVE_NEGATIVE_TREND] \ - + trend_type_to_percentages_values[AbstractTrendTest.SIGNIFICATIVE_POSITIVE_TREND] - ax.plot(xlabel_values, significative_values, 'y-' + marker, label=AbstractTrendTest.SIGNIFICATIVE + ' trends') + if axes is None: + fig, axes = plt.subplots(2, 1, figsize=self.study_visualizer.figsize) + + trend_type_to_series = self.trend_type_to_series(reduction_function) + for i, ax in enumerate(axes): + for trend_type in self.trend_types: + style = self.trend_type_to_style[trend_type] + percentages_values = trend_type_to_series[trend_type][i] + ax.plot(xlabel_values, percentages_values, style + marker, label=trend_type) + + if i == 0: + # Plot the total value of significative values + significative_values = trend_type_to_series[AbstractTrendTest.SIGNIFICATIVE_NEGATIVE_TREND][i] \ + + trend_type_to_series[AbstractTrendTest.SIGNIFICATIVE_POSITIVE_TREND][i] + ax.plot(xlabel_values, significative_values, 'y-' + marker, label=AbstractTrendTest.SIGNIFICATIVE + ' trends') + + # Global information + added_str = 'weighted ' + ylabel = '% averaged on massifs & {}averaged on {}'.format(added_str, xlabel) + ylabel += ' (with uniform weights)' + ax.set_ylabel(ylabel) + + ax.set_yticks(list(range(0, 101, 10))) + + # Common function functions + ax.set_xlabel(xlabel) + ax.grid() + ax.set_xticks(xlabel_values) + ax.legend() - # Global information - added_str = 'weighted ' - ylabel = '% averaged on massifs & {}averaged on {}'.format(added_str, xlabel) - ylabel += ' (with uniform weights)' - ax.set_ylabel(ylabel) - ax.set_xlabel(xlabel) - ax.set_xticks(xlabel_values) - ax.set_yticks(list(range(0, 101, 10))) - ax.grid() - ax.legend() - - variable_name = self.study.variable_class.NAME name = get_display_name_from_object_type(self.trend_test_class) - title = 'Evolution of {} trends (significative or not) wrt to the {} with {}'.format(variable_name, xlabel, + title = 'Evolution of {} trends (significative or not) wrt to the {} with {}'.format(subtitle, xlabel, name) - title += 'with {} data'.format(subtitle) - ax.set_title(title) + plt.suptitle(title) self.show_or_save_to_file(specific_title=title) def visualize_trend_test_repartition(self, reduction_function, axes=None, subtitle=''): @@ -88,7 +90,8 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer): fig, axes = plt.subplots(1, nb_trend_type, figsize=self.study_visualizer.figsize) # Plot weighted percentages over the years - trend_type_to_s_percentages = self.trend_type_to_s_percentages(reduction_function) + # todo: implement strength plot spatially + trend_type_to_s_percentages = {k: v[0] for k, v in self.trend_type_to_series(reduction_function).items()} for ax, trend_type in zip(axes, self.trend_types): s_percentages = trend_type_to_s_percentages[trend_type] massif_to_value = dict(s_percentages) @@ -113,7 +116,7 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer): def massif_index_level(self): return 1 - def visualize_year_trend_test(self, ax=None, marker='o', add_detailed_plots=False): + def visualize_year_trend_test(self, axes=None, marker='o', add_detailed_plots=False): def year_reduction(df_bool): # Take the mean with respect to all the first axis indices return df_bool.mean(axis=0) @@ -121,10 +124,10 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer): for subtitle, reduction_function in self.subtitle_to_reduction_function(year_reduction, add_detailed_plot=add_detailed_plots).items(): self.visualize_trend_test_evolution(reduction_function=reduction_function, xlabel='starting years', - xlabel_values=self.starting_years, ax=ax, marker=marker, + xlabel_values=self.starting_years, axes=axes, marker=marker, subtitle=subtitle) - def visualize_altitude_trend_test(self, ax=None, marker='o', add_detailed_plots=False): + def visualize_altitude_trend_test(self, axes=None, marker='o', add_detailed_plots=False): def altitude_reduction(df_bool, level): # Take the mean with respect to the years df_bool = df_bool.mean(axis=1) @@ -135,7 +138,7 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer): level=self.altitude_index_level, add_detailed_plot=add_detailed_plots).items(): self.visualize_trend_test_evolution(reduction_function=reduction_function, xlabel='altitude', - xlabel_values=self.altitudes, ax=ax, marker=marker, + xlabel_values=self.altitudes, axes=axes, marker=marker, subtitle=subtitle) def visualize_massif_trend_test(self, axes=None, add_detailed_plots=False): diff --git a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/main_hypercube_visualization.py b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/main_hypercube_visualization.py index 28e8e3bc..534541c8 100644 --- a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/main_hypercube_visualization.py +++ b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/main_hypercube_visualization.py @@ -37,9 +37,10 @@ def full_trends_with_quantity_altitude_hypercube(): save_to_file = True only_first_one = False fast = False + add_detailed_plots = False altitudes = ALL_ALTITUDES[3:-6] study_classes = SCM_STUDIES - for trend_test_class in [MannKendallTrendTest, GevLocationTrendTest, GevScaleTrendTest, GevShapeTrendTest][:]: + for trend_test_class in [MannKendallTrendTest, GevLocationTrendTest, GevScaleTrendTest, GevShapeTrendTest][1:2]: visualizers = [StudyVisualizer(study, temporal_non_stationarity=True, verbose=False, multiprocessing=True) for study in study_iterator_global(study_classes=study_classes, only_first_one=only_first_one, altitudes=altitudes)] @@ -48,9 +49,9 @@ def full_trends_with_quantity_altitude_hypercube(): quantity_altitude_to_visualizer = OrderedDict(zip(quantity_altitude_tuples, visualizers)) visualizer = QuantityAltitudeHypercubeVisualizer(quantity_altitude_to_visualizer, save_to_file=save_to_file, trend_test_class=trend_test_class, fast=fast) - visualizer.visualize_year_trend_test(add_detailed_plots=True) - visualizer.visualize_massif_trend_test(add_detailed_plots=True) - visualizer.visualize_altitude_trend_test(add_detailed_plots=True) + visualizer.visualize_year_trend_test(add_detailed_plots=add_detailed_plots) + visualizer.visualize_massif_trend_test(add_detailed_plots=add_detailed_plots) + visualizer.visualize_altitude_trend_test(add_detailed_plots=add_detailed_plots) def fast_trends_with_altitude_hypercube(): @@ -59,7 +60,7 @@ def fast_trends_with_altitude_hypercube(): fast = True altitudes = ALL_ALTITUDES[2:4] for study_class in SCM_STUDIES[:1]: - for trend_test_class in [MannKendallTrendTest, GevLocationTrendTest, GevScaleTrendTest, GevShapeTrendTest][:1]: + for trend_test_class in [MannKendallTrendTest, GevLocationTrendTest, GevScaleTrendTest, GevShapeTrendTest][1:2]: visualizers = [StudyVisualizer(study, temporal_non_stationarity=True, verbose=False, multiprocessing=True) for study in study_iterator(study_class=study_class, only_first_one=only_first_one, altitudes=altitudes)] @@ -67,8 +68,8 @@ def fast_trends_with_altitude_hypercube(): visualizer = AltitudeHypercubeVisualizer(altitude_to_visualizer, save_to_file=save_to_file, trend_test_class=trend_test_class, fast=fast) visualizer.visualize_year_trend_test() - # visualizer.visualize_massif_trend_test() - # visualizer.visualize_altitude_trend_test() + visualizer.visualize_massif_trend_test() + visualizer.visualize_altitude_trend_test() def fast_trends_with_quantity_altitude_hypercube(): @@ -92,9 +93,9 @@ def fast_trends_with_quantity_altitude_hypercube(): def main_run(): - fast_trends_with_altitude_hypercube() + # fast_trends_with_altitude_hypercube() # fast_trends_with_quantity_altitude_hypercube() - # full_trends_with_quantity_altitude_hypercube() + full_trends_with_quantity_altitude_hypercube() if __name__ == '__main__': diff --git a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/quantity_altitude_visualizer.py b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/quantity_altitude_visualizer.py index 92d20988..f4e29b8b 100644 --- a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/quantity_altitude_visualizer.py +++ b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/quantity_altitude_visualizer.py @@ -10,9 +10,10 @@ class QuantityAltitudeHypercubeVisualizer(AltitudeHypercubeVisualizer): def study_title(self): return 'Quantity Altitude Study' - def subtitle_to_reduction_function(self, reduction_function, level=None, add_detailed_plot=False): + def subtitle_to_reduction_function(self, reduction_function, level=None, add_detailed_plot=False, subtitle=None): subtitle_to_reduction_function = super().subtitle_to_reduction_function(reduction_function, - level, add_detailed_plot) + level, add_detailed_plot, + 'global') def get_function_from_tuple(tuple_for_axis_0): def f(df_bool: pd.DataFrame): diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py index 98f0aa6f..fc46ae26 100644 --- a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py +++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py @@ -376,7 +376,8 @@ class StudyVisualizer(object): trend_type_and_weight = [] years, smooth_maxima = self.smooth_maxima_x_y(massif_id) for starting_year, weight in starting_year_to_weight.items(): - test_trend_type = self.compute_trend_test_type(smooth_maxima, starting_year, trend_test_class, years) + test_trend_res = self.compute_trend_test_result(smooth_maxima, starting_year, trend_test_class, years) + test_trend_type = test_trend_res[0] trend_type_and_weight.append((test_trend_type, weight)) df = pd.DataFrame(trend_type_and_weight, columns=['trend type', 'weight']) massif_name_to_df_trend_type[massif_name] = df @@ -391,29 +392,35 @@ class StudyVisualizer(object): :param starting_year_to_weight: :return: """ - massif_name_to_trend_types = {} + massif_name_to_trend_res = {} massif_names = self.study.study_massif_names if nb_massif_for_fast_mode is not None: massif_names = massif_names[:nb_massif_for_fast_mode] for massif_id, massif_name in enumerate(massif_names): years, smooth_maxima = self.smooth_maxima_x_y(massif_id) if self.multiprocessing: - list_args = [(smooth_maxima, starting_year, trend_test_class, years) for starting_year in starting_years] + list_args = [(smooth_maxima, starting_year, trend_test_class, years) for starting_year in + starting_years] with Pool(NB_CORES) as p: - trend_types = p.starmap(self.compute_trend_test_type, list_args) + trend_test_res = p.starmap(self.compute_trend_test_result, list_args) else: - trend_types = [self.compute_trend_test_type(smooth_maxima, starting_year, trend_test_class, years) - for starting_year in starting_years] - massif_name_to_trend_types[massif_name] = trend_types - df = pd.DataFrame(massif_name_to_trend_types, index=starting_years) - return df.transpose() + trend_test_res = [self.compute_trend_test_result(smooth_maxima, starting_year, trend_test_class, years) + for starting_year in starting_years] + massif_name_to_trend_res[massif_name] = list(zip(*trend_test_res)) + # Build one DataFrame per trend test res + nb_res = len(list(massif_name_to_trend_res.values())[0]) + assert nb_res == 2 + all_massif_name_to_res = [{k: v[idx_res] for k, v in massif_name_to_trend_res.items()} + for idx_res in range(nb_res)] + return [pd.DataFrame(massif_name_to_res, index=starting_years).transpose() + for massif_name_to_res in all_massif_name_to_res] @staticmethod - def compute_trend_test_type(smooth_maxima, starting_year, trend_test_class, years): + def compute_trend_test_result(smooth_maxima, starting_year, trend_test_class, years): idx = years.index(starting_year) # assert years[0] == starting_year, "{} {}".format(years[0], starting_year) trend_test = trend_test_class(years[:][idx:], smooth_maxima[:][idx:]) # type: AbstractTrendTest - return trend_test.test_trend_type + return trend_test.test_trend_type, trend_test.test_trend_strength def df_trend_test_count(self, trend_test_class, starting_year_to_weight): """ diff --git a/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py index 42ae3c49..c9811d9e 100644 --- a/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py +++ b/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py @@ -4,6 +4,7 @@ from scipy.stats import chi2 from experiment.trend_analysis.univariate_trend_test.abstract_trend_test import AbstractTrendTest from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator +from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef from extreme_estimator.extreme_models.margin_model.temporal_linear_margin_model import StationaryStationModel, \ NonStationaryLocationStationModel, NonStationaryScaleStationModel, NonStationaryShapeStationModel from extreme_estimator.extreme_models.utils import SafeRunException @@ -67,24 +68,24 @@ class AbstractGevTrendTest(AbstractTrendTest): else: return super().test_trend_type - def get_coef(self, estimator): - return estimator.margin_function_fitted.get_coef(self.gev_param_name, AbstractCoordinates.COORDINATE_T) + def get_coef(self, estimator, coef_name): + return estimator.margin_function_fitted.get_coef(self.gev_param_name, coef_name) @property - def stationary_coef(self): - return self.get_coef(self.stationary_estimator) + def non_stationary_intercept_coef(self): + return self.get_coef(self.non_stationary_estimator, LinearCoef.INTERCEPT_NAME) @property - def non_stationary_coef(self): - return self.get_coef(self.non_stationary_estimator) + def non_stationary_linear_coef(self): + return self.get_coef(self.non_stationary_estimator, AbstractCoordinates.COORDINATE_T) @property - def ratio_non_stationary_coef(self): - pass + def test_trend_strength(self): + return 100 * np.abs(self.non_stationary_linear_coef) / np.abs(self.non_stationary_intercept_coef) @property def test_sign(self) -> int: - return np.sign(self.non_stationary_coef) + return np.sign(self.non_stationary_linear_coef) class GevLocationTrendTest(AbstractGevTrendTest): diff --git a/experiment/trend_analysis/univariate_trend_test/abstract_trend_test.py b/experiment/trend_analysis/univariate_trend_test/abstract_trend_test.py index 07466f6c..1462db94 100644 --- a/experiment/trend_analysis/univariate_trend_test/abstract_trend_test.py +++ b/experiment/trend_analysis/univariate_trend_test/abstract_trend_test.py @@ -31,6 +31,15 @@ class AbstractTrendTest(object): d[cls.NO_TREND] = 'k--' return d + @classmethod + def get_trend_types(cls, trend_type): + if trend_type is cls.POSITIVE_TREND: + return [cls.POSITIVE_TREND, cls.SIGNIFICATIVE_POSITIVE_TREND] + elif trend_type is cls.NEGATIVE_TREND: + return [cls.NEGATIVE_TREND, cls.SIGNIFICATIVE_NEGATIVE_TREND] + else: + return [trend_type] + @classmethod def get_cmap_from_trend_type(cls, trend_type): if 'positive' in trend_type: @@ -49,6 +58,10 @@ class AbstractTrendTest(object): def n(self): return len(self.years_after_change_point) + @property + def test_trend_strength(self): + return 0.0 + @property def test_trend_type(self) -> str: test_sign = self.test_sign diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py index 4c81dd5c..4c1b0f2a 100644 --- a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py +++ b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py @@ -83,7 +83,7 @@ class LinearMarginFunction(ParametricMarginFunction): # Properties for the location parameter def get_coef(self, gev_param_name, coef_name): - idx = 1 if coef_name in [AbstractCoordinates.COORDINATE_T] \ + idx = 1 if coef_name in [AbstractCoordinates.COORDINATE_T, LinearCoef.INTERCEPT_NAME] \ else AbstractCoordinates.COORDINATE_SPATIAL_NAMES.index(coef_name) + 2 return self.coef_dict[LinearCoef.coef_template_str(gev_param_name, coef_name).format(idx)] -- GitLab