[SCM][HYPERCUBE] add trend_strength to the hypercube_visualization.

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}
-
-
-
-

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):

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__':

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):

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):
         """

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):

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

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)]