From 16f48fe8c55e8596e93419fe699dfcced7ee9123 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 5 Sep 2019 14:13:13 +0200
Subject: [PATCH] [POSTER EVAN] add hatch to represent the participation of the
mean & variance variation
---
.../scm_models_data/abstract_study.py | 14 ++++++
.../abstract_hypercube_visualizer.py | 8 ++++
.../altitude_hypercube_visualizer.py | 44 +++++++++++++++----
.../study_visualization/study_visualizer.py | 7 +--
.../poster_EVAN2019/main_poster_EVAN2019.py | 2 +-
.../abstract_gev_trend_test.py | 15 +++++++
.../gev_trend_test_one_parameter.py | 21 +++++++++
.../gev_trend_test_two_parameters.py | 24 ++++++++--
.../trend_analysis/univariate_test/utils.py | 9 +++-
9 files changed, 126 insertions(+), 18 deletions(-)
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 ea6d7a67..c4ee1b58 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -1,4 +1,5 @@
import datetime
+from matplotlib.patches import Polygon
import io
import os
import os.path as op
@@ -268,6 +269,7 @@ class AbstractStudy(object):
scaled=False,
fontsize=7,
axis_off=False,
+ massif_name_to_hatch_boolean_list=None
):
if ax is None:
ax = plt.gca()
@@ -292,8 +294,10 @@ class AbstractStudy(object):
# if j == 0:
# mask_outside_polygon(poly_verts=l, ax=ax)
# Plot the contour of the massif
+
coords_list = list(zip(*coords_list))
ax.plot(*coords_list, color='black')
+
# Potentially, fill the inside of the polygon with some color
if fill and coordinate_id in cls.coordinate_id_to_massif_name:
massif_name = cls.coordinate_id_to_massif_name[coordinate_id]
@@ -303,6 +307,16 @@ class AbstractStudy(object):
else:
ax.fill(*coords_list, **{'color': default_color_for_missing_massif})
+ # Add a hatch to visualize the mean & variance variation sign
+ hatch_list = ['//', '\\\\']
+ if massif_name_to_hatch_boolean_list is not None:
+ if massif_name in massif_name_to_hatch_boolean_list:
+ a = np.array(coords_list).transpose()
+ hatch_boolean_list = massif_name_to_hatch_boolean_list[massif_name]
+ for hatch, is_hatch in zip(hatch_list, hatch_boolean_list):
+ if is_hatch:
+ ax.add_patch(Polygon(xy=a, fill=False, hatch=hatch))
+
# Display the center of the massif
masssif_coordinate_for_display = cls.massifs_coordinates_for_display(massif_names)
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
index 728090c1..d5e2e7c2 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
@@ -102,6 +102,14 @@ class AbstractHypercubeVisualizer(object):
def df_hypercube_trend_constant_quantile(self) -> pd.DataFrame:
return self._df_hypercube_trend_meta(idx=3)
+ @cached_property
+ def df_hypercube_trend_mean_same_sign(self) -> pd.DataFrame:
+ return self._df_hypercube_trend_meta(idx=4)
+
+ @cached_property
+ def df_hypercube_trend_variance_same_sign(self) -> pd.DataFrame:
+ return self._df_hypercube_trend_meta(idx=5)
+
# Some properties
@property
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
index f288db92..76e38352 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
@@ -82,6 +82,10 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
df_constant = self.df_hypercube_trend_constant_quantile[df_bool]
s_trend_constant = reduction_function(df_constant)
series.extend([s_trend_strength, s_trend_constant])
+ # Add the mean and the variance anyway
+ s_trend_mean_sign = reduction_function(self.df_hypercube_trend_mean_same_sign[df_bool])
+ s_trend_variance_sign = reduction_function(self.df_hypercube_trend_variance_same_sign[df_bool])
+ series.extend([s_trend_mean_sign, s_trend_variance_sign])
return series
def subtitle_to_reduction_function(self, reduction_function, level=None, add_detailed_plot=False, subtitle=None):
@@ -294,6 +298,8 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
massif_to_year = {}
massif_to_strength = {}
massif_to_constant = {}
+ massif_to_mean_difference_same_sign = {}
+ massif_to_variance_difference_same_sign = {}
poster_trend_types = [AbstractUnivariateTest.SIGNIFICATIVE_POSITIVE_TREND,
AbstractUnivariateTest.SIGNIFICATIVE_NEGATIVE_TREND,
AbstractUnivariateTest.NEGATIVE_TREND,
@@ -314,32 +320,52 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
if k in massif_to_color_for_trend_type} for i in [1, 2]]
massif_to_strength.update(massif_to_value_for_trend_type[0])
massif_to_constant.update(massif_to_value_for_trend_type[1])
+ mean_idx, variance_idx = 2, 3
else:
+ mean_idx, variance_idx = 1, 2
+
massif_to_value_for_trend_type = {k: "$t_0=$" + str(int(v)) for k, v in
self.trend_type_to_series(reduction_function,
isin_parameters)[
- display_trend_type][1].items()
+ display_trend_type][3].items()
if k in massif_to_color_for_trend_type}
massif_to_year.update(massif_to_value_for_trend_type)
+
+ # Add the mean and variance sign anyway
+ massif_to_value_for_trend_type = [{k: v for k, v in
+ self.trend_type_to_series(reduction_function,
+ isin_parameters)[
+ display_trend_type][i].items()
+ if k in massif_to_color_for_trend_type} for i in
+ [mean_idx, variance_idx]]
+ massif_to_mean_difference_same_sign.update(massif_to_value_for_trend_type[0])
+ massif_to_variance_difference_same_sign.update(massif_to_value_for_trend_type[1])
+ # Compute massif to hatch boolean
+ massif_name_to_hatch_boolean_list = {
+ massif: [massif_to_mean_difference_same_sign[massif] == 1.0,
+ massif_to_variance_difference_same_sign[massif] == 1.0]
+ for massif in massif_to_color.keys()
+ }
+
# Compute massif_to_value
if self.reduce_strength_array:
massif_name_to_value = {m: "{} {}{}".format(
- int(massif_to_constant[m]),
- "+" if massif_to_strength[m] > 0 else "",
- round(massif_to_strength[m] * massif_to_constant[m], 1),
- AbstractGevTrendTest.nb_years_for_quantile_evolution)
- for m in massif_to_strength}
+ int(massif_to_constant[m]),
+ "+" if massif_to_strength[m] > 0 else "",
+ round(massif_to_strength[m] * massif_to_constant[m], 1),
+ AbstractGevTrendTest.nb_years_for_quantile_evolution)
+ for m in massif_to_strength}
else:
massif_name_to_value = massif_to_year
self.study.visualize_study(None, massif_name_to_color=massif_to_color, show=False,
show_label=False, scaled=True, add_text=write_text_on_massif,
massif_name_to_value=massif_name_to_value,
fontsize=4,
- axis_off=True)
+ axis_off=True,
+ massif_name_to_hatch_boolean_list=massif_name_to_hatch_boolean_list)
title = self.set_trend_test_reparition_title(subtitle, set=not poster_plot)
-
return title
def set_trend_test_reparition_title(self, subtitle, set=True):
@@ -464,7 +490,7 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
isin_parameters=None,
show_or_save_to_file=True,
poster_plot=False,
- write_text_on_massif=True):
+ write_text_on_massif=True):
last_title = ''
for subtitle, reduction_function in self.subtitle_to_reduction_function(self.index_reduction,
level=self.massif_index_level,
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
index d4fa910e..45a7526e 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
@@ -467,11 +467,12 @@ class StudyVisualizer(VisualizationParameters):
sample_one_massif_from_each_region)
for massif_name, gev_change_point_test_results in massif_name_to_gev_change_point_test_results.items():
trend_test_res, best_idxs = gev_change_point_test_results
- trend_test_res = [(a, b, c, d) if i in best_idxs else (np.nan, np.nan, c, np.nan)
- for i, (a, b, c, d, *_) in enumerate(trend_test_res)]
+ trend_test_res = [(a, b, c, d, e, f) if i in best_idxs else (np.nan, np.nan, c, np.nan, np.nan, np.nan)
+ for i, (a, b, c, d, e, f, *_) in enumerate(trend_test_res)]
massif_name_to_trend_res[massif_name] = list(zip(*trend_test_res))
nb_res = len(list(massif_name_to_trend_res.values())[0])
- assert nb_res == 4
+ assert nb_res == 6
+
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=self.starting_years_for_change_point_test).transpose()
diff --git a/experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py b/experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py
index 541c719a..db48f63a 100644
--- a/experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py
+++ b/experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py
@@ -9,7 +9,7 @@ POSTER_ALTITUDES = [900, 1800, 2700]
def main_poster_A_non_stationary_model_choice():
- nb = 1
+ nb = 3
for altitude in POSTER_ALTITUDES[:nb]:
for trend_test_class in [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest][-nb:]:
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
index 456390e3..80fd5214 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
@@ -131,6 +131,21 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
def _slope_strength(self):
raise NotImplementedError
+ @staticmethod
+ def same_sign(a, b):
+ return (a > 0 and b > 0) or (a < 0 and b < 0)
+
+ @property
+ def mean_difference_same_sign_as_slope_strenght(self) -> bool:
+ return False
+
+ @property
+ def variance_difference_same_sign_as_slope_strenght(self) -> bool:
+ return False
+
+ def mean_difference(self, zeta0, mu1=0.0, sigma1=0.0):
+ return GevParams(loc=mu1, scale=sigma1, shape=zeta0).mean
+
@property
def test_trend_constant_quantile(self):
if self.crashed:
diff --git a/experiment/trend_analysis/univariate_test/gev_trend_test_one_parameter.py b/experiment/trend_analysis/univariate_test/gev_trend_test_one_parameter.py
index b27320f7..f4582a5b 100644
--- a/experiment/trend_analysis/univariate_test/gev_trend_test_one_parameter.py
+++ b/experiment/trend_analysis/univariate_test/gev_trend_test_one_parameter.py
@@ -29,6 +29,16 @@ class GevLocationTrendTest(GevTrendTestOneParameter):
return self.non_stationary_constant_gev_params.quantile_strength_evolution(p=self.quantile_for_strength,
mu1=self.non_stationary_linear_coef)
+ @property
+ def mean_difference_same_sign_as_slope_strenght(self) -> bool:
+ zeta0 = self.non_stationary_constant_gev_params.shape
+ mean_difference = self.mean_difference(zeta0=zeta0, mu1=self.non_stationary_linear_coef)
+ return self.same_sign(mean_difference, self._slope_strength())
+
+ @property
+ def variance_difference_same_sign_as_slope_strenght(self) -> bool:
+ return False
+
class GevScaleTrendTest(GevTrendTestOneParameter):
@@ -41,6 +51,17 @@ class GevScaleTrendTest(GevTrendTestOneParameter):
p=self.quantile_for_strength,
sigma1=self.non_stationary_linear_coef)
+ @property
+ def mean_difference_same_sign_as_slope_strenght(self) -> bool:
+ zeta0 = self.non_stationary_constant_gev_params.shape
+ mean_difference = self.mean_difference(zeta0=zeta0, sigma1=self.non_stationary_linear_coef)
+ return self.same_sign(mean_difference, self._slope_strength())
+
+ @property
+ def variance_difference_same_sign_as_slope_strenght(self) -> bool:
+ sigma1 = self.non_stationary_linear_coef
+ return self.same_sign(sigma1, self._slope_strength())
+
class GevShapeTrendTest(GevTrendTestOneParameter):
diff --git a/experiment/trend_analysis/univariate_test/gev_trend_test_two_parameters.py b/experiment/trend_analysis/univariate_test/gev_trend_test_two_parameters.py
index bbd2c8ef..386ed310 100644
--- a/experiment/trend_analysis/univariate_test/gev_trend_test_two_parameters.py
+++ b/experiment/trend_analysis/univariate_test/gev_trend_test_two_parameters.py
@@ -17,9 +17,25 @@ class GevLocationAndScaleTrendTest(GevTrendTestTwoParameters):
super().__init__(years, maxima, starting_year,
NonStationaryLocationAndScaleModel)
+ @property
+ def mu1(self):
+ return self.get_non_stationary_linear_coef(gev_param_name=GevParams.LOC)
+
+ @property
+ def sigma1(self):
+ return self.get_non_stationary_linear_coef(gev_param_name=GevParams.SCALE)
+
def _slope_strength(self):
- mu1 = self.get_non_stationary_linear_coef(gev_param_name=GevParams.LOC)
- sigma1 = self.get_non_stationary_linear_coef(gev_param_name=GevParams.SCALE)
return self.non_stationary_constant_gev_params.quantile_strength_evolution(p=self.quantile_for_strength,
- mu1=mu1,
- sigma1=sigma1)
+ mu1=self.mu1,
+ sigma1=self.sigma1)
+
+ @property
+ def mean_difference_same_sign_as_slope_strenght(self) -> bool:
+ zeta0 = self.non_stationary_constant_gev_params.shape
+ mean_difference = self.mean_difference(zeta0=zeta0, mu1=self.mu1, sigma1=self.sigma1)
+ return self.same_sign(mean_difference, self._slope_strength())
+
+ @property
+ def variance_difference_same_sign_as_slope_strenght(self) -> bool:
+ return self.same_sign(self.sigma1, self._slope_strength())
diff --git a/experiment/trend_analysis/univariate_test/utils.py b/experiment/trend_analysis/univariate_test/utils.py
index 23da2440..18e3aafe 100644
--- a/experiment/trend_analysis/univariate_test/utils.py
+++ b/experiment/trend_analysis/univariate_test/utils.py
@@ -9,7 +9,14 @@ from utils import NB_CORES
def compute_gev_change_point_test_result(smooth_maxima, starting_year, trend_test_class, years):
trend_test = trend_test_class(years, smooth_maxima, starting_year) # type: AbstractGevTrendTest
assert isinstance(trend_test, AbstractGevTrendTest)
- return trend_test.test_trend_type, trend_test.test_trend_slope_strength, trend_test.non_stationary_nllh, trend_test.test_trend_constant_quantile, trend_test.non_stationary_deviance, trend_test.stationary_deviance
+ return trend_test.test_trend_type, \
+ trend_test.test_trend_slope_strength, \
+ trend_test.non_stationary_nllh, \
+ trend_test.test_trend_constant_quantile, \
+ trend_test.mean_difference_same_sign_as_slope_strenght, \
+ trend_test.variance_difference_same_sign_as_slope_strenght, \
+ trend_test.non_stationary_deviance, \
+ trend_test.stationary_deviance
def compute_gev_change_point_test_results(multiprocessing, maxima, starting_years, trend_test_class,
--
GitLab