[SCM][HYPERCUBE] remove run time error and no trend from trend plot. add...

[SCM][HYPERCUBE] remove run time error and no trend from trend plot. add Altitude_Hypercube_Year_Visualizer as an extension of Altitude_Hypercube_Visualizer

experiment/meteo_france_SCM_study/abstract_study.py

@@ -237,7 +237,7 @@ class AbstractStudy(object):
     """ Visualization methods """
 
     def visualize_study(self, ax=None, massif_name_to_value=None, show=True, fill=True, replace_blue_by_white=True,
-                        label=None, add_text=False, cmap=None, vmax=100):
+                        label=None, add_text=False, cmap=None, vmax=100, vmin=0):
         if massif_name_to_value is None:
             massif_name_to_fill_kwargs = None
         else:
@@ -245,7 +245,7 @@ class AbstractStudy(object):
             if cmap is None:
                 colors = get_color_rbga_shifted(ax, replace_blue_by_white, values, label=label)
             else:
-                norm = Normalize(0, vmax)
+                norm = Normalize(vmin, vmax)
                 create_colorbase_axis(ax, label, cmap, norm)
                 m = cm.ScalarMappable(norm=norm, cmap=cmap)
                 colors = [m.to_rgba(value) for value in values]

experiment/meteo_france_SCM_study/visualization/hypercube_visualization/altitude_hypercube_visualizer.py

 import matplotlib.pyplot as plt
+import numpy as np
 import pandas as pd
 
 from experiment.meteo_france_SCM_study.visualization.hypercube_visualization.abstract_hypercube_visualizer import \
@@ -22,22 +23,28 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
     def trend_types(self):
         return self.trend_type_to_style.keys()
 
+    @property
+    def nb_axes(self):
+        return 2
+
     def trend_type_to_series(self, reduction_function):
         # Map each trend type to its serie with percentages
-        trend_type_to_series = {}
-        for trend_type in self.trend_types:
-            # Reduce df_bool df to a serie s_trend_type_percentage
-            df_bool = self.df_hypercube_trend_type.isin(AbstractUnivariateTest.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
+        return {trend_type: self.trend_type_reduction(reduction_function, trend_type)[0]
+                for trend_type in self.trend_types}
+
+    def trend_type_reduction(self, reduction_function, trend_type):
+        # Reduce df_bool df to a serie s_trend_type_percentage
+        df_bool = self.df_hypercube_trend_type.isin(AbstractUnivariateTest.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)
+        # Group result
+        series = [s_trend_type_percentage, s_trend_strength]
+        return series, df_bool
 
     def subtitle_to_reduction_function(self, reduction_function, level=None, add_detailed_plot=False, subtitle=None):
         def reduction_function_with_level(df_bool):
@@ -64,27 +71,21 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
     def visualize_trend_test_evolution(self, reduction_function, xlabel, xlabel_values, axes=None, marker='o',
                                        subtitle=''):
         if axes is None:
-            fig, axes = plt.subplots(2, 1, figsize=self.study_visualizer.figsize)
+            fig, axes = plt.subplots(self.nb_axes, 1, figsize=self.study_visualizer.figsize)
 
         trend_type_to_series = self.trend_type_to_series(reduction_function)
-        for i, ax in enumerate(axes):
+        for ax_idx, 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]
+                percentages_values = trend_type_to_series[trend_type][ax_idx]
                 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[AbstractUnivariateTest.SIGNIFICATIVE_NEGATIVE_TREND][i] \
-                #                        + trend_type_to_series[AbstractUnivariateTest.SIGNIFICATIVE_POSITIVE_TREND][i]
-                # ax.plot(xlabel_values, significative_values, 'y-' + marker,
-                #         label=AbstractUnivariateTest.SIGNIFICATIVE + ' trends')
-
+            if ax_idx == 0:
                 # Global information
                 ax.set_ylabel(self.get_title_plot(xlabel, ax_idx=0))
                 ax.set_yticks(list(range(0, 101, 10)))
             else:
-                ax.set_ylabel(self.get_title_plot(xlabel, ax_idx=1))
+                ax.set_ylabel(self.get_title_plot(xlabel, ax_idx=ax_idx))
 
             # Common function functions
             ax.set_xlabel(xlabel)
@@ -101,17 +102,18 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
     def visualize_trend_test_repartition(self, reduction_function, axes=None, subtitle=''):
         if axes is None:
             nb_trend_type = len(self.trend_test_class.trend_type_to_style())
-            fig, axes = plt.subplots(2, nb_trend_type, figsize=self.study_visualizer.figsize)
+            fig, axes = plt.subplots(self.nb_axes, nb_trend_type, figsize=self.study_visualizer.figsize)
 
         for i, axes_row in enumerate(axes):
             trend_type_to_serie = {k: v[i] for k, v in self.trend_type_to_series(reduction_function).items()}
             vmax = max([s.max() for s in trend_type_to_serie.values()])
+            vmin = min([s.min() for s in trend_type_to_serie.values()])
             vmax = max(vmax, 0.01)
             for ax, trend_type in zip(axes_row, self.trend_types):
                 s_percentages = trend_type_to_serie[trend_type]
                 massif_to_value = dict(s_percentages)
                 cmap = self.trend_test_class.get_cmap_from_trend_type(trend_type)
-                self.study.visualize_study(ax, massif_to_value, show=False, cmap=cmap, label=None, vmax=vmax)
+                self.study.visualize_study(ax, massif_to_value, show=False, cmap=cmap, label=None, vmax=vmax, vmin=vmin)
                 ax.set_title(trend_type)
             row_title = self.get_title_plot(xlabel='massifs', ax_idx=i)
             StudyVisualizer.clean_axes_write_title_on_the_left(axes_row, row_title, left_border=None)
@@ -131,9 +133,9 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
         return 1
 
     def visualize_year_trend_test(self, axes=None, marker='o', add_detailed_plots=False):
-        def year_reduction(df_bool):
+        def year_reduction(df):
             # Take the mean with respect to all the first axis indices
-            return df_bool.mean(axis=0)
+            return df.mean(axis=0)
 
         for subtitle, reduction_function in self.subtitle_to_reduction_function(year_reduction,
                                                                                 add_detailed_plot=add_detailed_plots).items():
@@ -141,14 +143,15 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
                                                 xlabel_values=self.starting_years, axes=axes, marker=marker,
                                                 subtitle=subtitle)
 
-    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.any(axis=1)
-            # Take the mean with respect the massifs
-            return df_bool.mean(level=level)
+    @staticmethod
+    def index_reduction(df, level):
+        # Take the sum with respect to the years, replace any missing data with np.nan
+        df = df.any(axis=1)
+        # Take the mean with respect to the level of interest
+        return df.mean(level=level)
 
-        for subtitle, reduction_function in self.subtitle_to_reduction_function(altitude_reduction,
+    def visualize_altitude_trend_test(self, axes=None, marker='o', add_detailed_plots=False):
+        for subtitle, reduction_function in self.subtitle_to_reduction_function(self.index_reduction,
                                                                                 level=self.altitude_index_level,
                                                                                 add_detailed_plot=add_detailed_plots).items():
             self.visualize_trend_test_evolution(reduction_function=reduction_function, xlabel='altitudes',
@@ -156,13 +159,37 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
                                                 subtitle=subtitle)
 
     def visualize_massif_trend_test(self, axes=None, add_detailed_plots=False):
-        def massif_reduction(df_bool, level):
-            # Take the mean with respect to the years
-            df_bool = df_bool.any(axis=1)
-            # Take the mean with respect the altitude
-            return df_bool.mean(level=level)
-
-        for subtitle, reduction_function in self.subtitle_to_reduction_function(massif_reduction,
+        for subtitle, reduction_function in self.subtitle_to_reduction_function(self.index_reduction,
                                                                                 level=self.massif_index_level,
                                                                                 add_detailed_plot=add_detailed_plots).items():
             self.visualize_trend_test_repartition(reduction_function, axes, subtitle=subtitle)
+
+
+class Altitude_Hypercube_Year_Visualizer(AltitudeHypercubeVisualizer):
+
+    def get_title_plot(self, xlabel, ax_idx=None):
+        if ax_idx == 2:
+            return 'mean starting year'
+        return super().get_title_plot(xlabel, ax_idx)
+
+    @property
+    def nb_axes(self):
+        return 3
+
+    @staticmethod
+    def index_reduction(df, level):
+        # Take the sum with respect to the years, replace any missing data with np.nan
+        df = df.sum(axis=1).replace(0.0, np.nan)
+        # Take the mean with respect to the level of interest
+        return df.mean(level=level)
+
+    def trend_type_reduction(self, reduction_function, trend_type):
+        series, df_bool = super().trend_type_reduction(reduction_function, trend_type)
+        # Create df argmax
+        df = df_bool.copy()
+        df = (df * df.columns)[df_bool]
+        # Reduce and append
+        serie = reduction_function(df)
+        series.append(serie)
+        return series, df_bool
+

experiment/meteo_france_SCM_study/visualization/hypercube_visualization/main_hypercube_visualization.py

@@ -3,7 +3,7 @@ from itertools import product
 from collections import OrderedDict
 
 from experiment.meteo_france_SCM_study.visualization.hypercube_visualization.altitude_hypercube_visualizer import \
-    AltitudeHypercubeVisualizer
+    AltitudeHypercubeVisualizer, Altitude_Hypercube_Year_Visualizer
 from experiment.meteo_france_SCM_study.visualization.hypercube_visualization.quantity_altitude_visualizer import \
     QuantityAltitudeHypercubeVisualizer
 from experiment.meteo_france_SCM_study.visualization.study_visualization.main_study_visualizer import ALL_ALTITUDES, \
@@ -66,6 +66,45 @@ def fast_altitude_hypercube():
             altitude_to_visualizer = OrderedDict(zip(altitudes, visualizers))
             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()
+
+
+def fast_altitude_year_hypercube():
+    save_to_file = False
+    only_first_one = False
+    fast = True
+    altitudes = ALL_ALTITUDES[2:4]
+    for study_class in SCM_STUDIES[:1]:
+        for trend_test_class in [GevLocationChangePointTest, GevScaleChangePointTest, GevShapeChangePointTest][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)]
+            altitude_to_visualizer = OrderedDict(zip(altitudes, visualizers))
+            visualizer = Altitude_Hypercube_Year_Visualizer(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()
+
+
+def full_altitude_year_hypercube():
+    save_to_file = False
+    only_first_one = False
+    fast = False
+    altitudes = ALL_ALTITUDES[3:-6]
+    for study_class in SCM_STUDIES[:1]:
+        for trend_test_class in [GevLocationChangePointTest, GevScaleChangePointTest,
+                                 GevShapeChangePointTest][:1]:
+            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)]
+            altitude_to_visualizer = OrderedDict(zip(altitudes, visualizers))
+            visualizer = Altitude_Hypercube_Year_Visualizer(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()
@@ -92,7 +131,9 @@ def fast_quantity_altitude_hypercube():
 
 
 def main_run():
-    fast_altitude_hypercube()
+    # fast_altitude_hypercube()
+    # fast_altitude_year_hypercube()
+    full_altitude_year_hypercube()
     # fast_quantity_altitude_hypercube()
     # full_quantity_altitude_hypercube()
 

experiment/trend_analysis/abstract_score.py

@@ -17,6 +17,7 @@ class AbstractTrendScore(object):
 
 
 class MannKendall(AbstractTrendScore):
+    # see here for the explanation: https://up-rs-esp.github.io/mkt/
 
     def get_detailed_score(self, years_after_change_point, maxima_after_change_point):
         score = 0.0

experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py

@@ -62,11 +62,11 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
 
     # Add a trend type that correspond to run that crashed
 
-    @classmethod
-    def trend_type_to_style(cls):
-        trend_type_to_style = super().trend_type_to_style()
-        trend_type_to_style[cls.RRunTimeError_TREND] = 'b:'
-        return trend_type_to_style
+    # @classmethod
+    # def trend_type_to_style(cls):
+    #     trend_type_to_style = super().trend_type_to_style()
+    #     trend_type_to_style[cls.RRunTimeError_TREND] = 'b:'
+    #     return trend_type_to_style
 
     @property
     def test_trend_type(self) -> str:

experiment/trend_analysis/univariate_test/abstract_univariate_test.py

@@ -15,8 +15,10 @@ class AbstractUnivariateTest(object):
     SIGNIFICATIVE = 'significative'
     # 5 possible types of trends
     NO_TREND = 'no trend'
+    ALL_TREND = 'all trend'
     POSITIVE_TREND = 'positive trend'
     NEGATIVE_TREND = 'negative trend'
+    SIGNIFICATIVE_ALL_TREND = SIGNIFICATIVE + ' ' + ALL_TREND
     SIGNIFICATIVE_POSITIVE_TREND = SIGNIFICATIVE + ' ' + POSITIVE_TREND
     SIGNIFICATIVE_NEGATIVE_TREND = SIGNIFICATIVE + ' ' + NEGATIVE_TREND
 
@@ -43,15 +45,21 @@ class AbstractUnivariateTest(object):
     @classmethod
     def trend_type_to_style(cls):
         d = OrderedDict()
+        d[cls.ALL_TREND] = 'k--'
         d[cls.POSITIVE_TREND] = 'g--'
         d[cls.NEGATIVE_TREND] = 'r--'
+        d[cls.SIGNIFICATIVE_ALL_TREND] = 'k-'
         d[cls.SIGNIFICATIVE_POSITIVE_TREND] = 'g-'
         d[cls.SIGNIFICATIVE_NEGATIVE_TREND] = 'r-'
-        d[cls.NO_TREND] = 'k--'
+        # d[cls.NO_TREND] = 'k--'
         return d
 
     @classmethod
     def get_trend_types(cls, trend_type):
+        if trend_type is cls.ALL_TREND:
+            return cls.get_trend_types(cls.POSITIVE_TREND) + cls.get_trend_types(cls.NEGATIVE_TREND)
+        elif trend_type is cls.SIGNIFICATIVE_ALL_TREND:
+            return [cls.SIGNIFICATIVE_POSITIVE_TREND, cls.SIGNIFICATIVE_NEGATIVE_TREND]
         if trend_type is cls.POSITIVE_TREND:
             return [cls.POSITIVE_TREND, cls.SIGNIFICATIVE_POSITIVE_TREND]
         elif trend_type is cls.NEGATIVE_TREND: