diff --git a/experiment/meteo_france_SCM_study/abstract_study.py b/experiment/meteo_france_SCM_study/abstract_study.py
index 4420f433a9cbd12a0e363d88746fceeef15ae285..51b76c85a99079c0fa6b23b4a2a5ca15c24bad59 100644
--- a/experiment/meteo_france_SCM_study/abstract_study.py
+++ b/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]
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 df0115f1ba7e22452ff0aa64b0de14ad77bf6439..25ef29669f63a18955f4ced858c727f78ec492f3 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
@@ -1,4 +1,5 @@
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
+
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 2646e1ad01ebb246fad20a803fa71785de669bfd..07a1c6b9a22acfd59515e641783835c822fd2e6a 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
@@ -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()
diff --git a/experiment/trend_analysis/abstract_score.py b/experiment/trend_analysis/abstract_score.py
index ee653e54a6800b310ac85121cee532f6ac5128ff..d5506b51e605fa0cff15d889d4f032265f1575a0 100644
--- a/experiment/trend_analysis/abstract_score.py
+++ b/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
diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
index 3d1b29001ff4d4c772259b3a3e2ff443fc28e360..cc46b8a8713c05a6dee67d773ae4ae4bde91ecc7 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
+++ b/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:
diff --git a/experiment/trend_analysis/univariate_test/abstract_univariate_test.py b/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
index 2da67e6338c107bcb11aea0319d15985dcbef0ca..2bb3baf3b6d8c920f08a682662e813b728579086 100644
--- a/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
+++ b/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: