diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
index 43e1b0e5bdef88f21253d0acdb88e2f16cae8813..57128a0c42464e25f940f00960a0f8630afda259 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
@@ -116,17 +116,17 @@ def complete_analysis(only_first_one=False):
def trend_analysis():
- save_to_file = False
- only_first_one = True
+ save_to_file = True
+ only_first_one = False
short_altitudes = [300, 1200, 2100, 3000][:1]
- full_altitude_with_at_least_2_stations = [0, 300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200][-1:]
+ full_altitude_with_at_least_2_stations = [0, 300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200][-5:]
altitudes = full_altitude_with_at_least_2_stations
normalization_class = [None, BetweenMinusOneAndOneNormalization, BetweenZeroAndOneNormalization][-1]
- study_classes = [CrocusSwe, CrocusDepth, SafranSnowfall, SafranRainfall, SafranTemperature][:]
+ study_classes = [CrocusSwe, CrocusDepth, SafranSnowfall, SafranRainfall, SafranTemperature][:1]
for study in study_iterator_global(study_classes, only_first_one=only_first_one, altitudes=altitudes):
study_visualizer = StudyVisualizer(study, save_to_file=save_to_file,
transformation_class=normalization_class)
- study_visualizer.visualize_temporal_trend_relevance(complete_analysis=True, multiprocessing=True)
+ study_visualizer.visualize_temporal_trend_relevance(complete_analysis=False, multiprocessing=True)
if __name__ == '__main__':
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py
index 89be0f2c2fdbbc6214fa7fd41613fd752b86276f..b805e09a72412775158e50e3467ec660d553bb4d 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py
@@ -2,6 +2,8 @@ import time
from multiprocessing import Pool
from typing import Union
+import pandas as pd
+
from experiment.meteo_france_SCM_study.visualization.utils import align_yaxis_on_zero
from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
from scipy.stats import chi2
@@ -32,89 +34,124 @@ class AbstractNonStationaryTrendTest(object):
self.non_stationary_margin_model_class = non_stationary_margin_model_class
# Compute a dictionary that maps couple (margin model class, starting point)
# to the corresponding fitted estimator
- self._margin_model_class_and_starting_point_to_estimator = {}
+ self._starting_point_to_estimator = {}
# parallelization arguments
self.multiprocessing = multiprocessing
self.nb_cores = 7
- def get_estimator(self, margin_model_class, starting_point) -> Union[
+ def get_estimator(self, starting_point):
+ if starting_point not in self._starting_point_to_estimator:
+ estimator = self.load_estimator(starting_point)
+ self._starting_point_to_estimator[starting_point] = estimator
+ return self._starting_point_to_estimator[starting_point]
+
+
+ def load_estimator(self, starting_point) -> Union[
AbstractFullEstimator, AbstractMarginEstimator]:
- if (margin_model_class, starting_point) not in self._margin_model_class_and_starting_point_to_estimator:
- margin_model = margin_model_class(coordinates=self.dataset.coordinates, starting_point=starting_point)
- estimator = self._load_estimator(margin_model)
- start = time.time()
- estimator.fit()
- duration = time.time() - start
- if self.verbose:
- if self.verbose:
- estimator_name = get_display_name_from_object_type(estimator)
- margin_model_name = get_display_name_from_object_type(margin_model)
- text = 'Fittig {} with margin: {} for starting_point={}\n'.format(estimator_name,
- margin_model_name,
- starting_point)
- text += 'Fit took {}s and was {}'.format(round(duration, 1), estimator.result_from_fit.convergence)
- print(text)
- self._margin_model_class_and_starting_point_to_estimator[(margin_model_class, starting_point)] = estimator
- return self._margin_model_class_and_starting_point_to_estimator[(margin_model_class, starting_point)]
+ margin_model_class = self.stationary_margin_model_class if starting_point is None else self.non_stationary_margin_model_class
+ assert starting_point not in self._starting_point_to_estimator
+ margin_model = margin_model_class(coordinates=self.dataset.coordinates, starting_point=starting_point)
+ estimator = self._load_estimator(margin_model)
+ start = time.time()
+ estimator.fit()
+ duration = time.time() - start
+ if self.verbose:
+ estimator_name = get_display_name_from_object_type(estimator)
+ margin_model_name = get_display_name_from_object_type(margin_model)
+ text = 'Fittig {} with margin: {} for starting_point={}\n'.format(estimator_name,
+ margin_model_name,
+ starting_point)
+ text += 'Fit took {}s and was {}'.format(round(duration, 1), estimator.result_from_fit.convergence)
+ print(text)
+ return estimator
def _load_estimator(self, margin_model) -> Union[AbstractFullEstimator, AbstractMarginEstimator]:
return self.estimator_class(self.dataset, margin_model)
def get_metric(self, starting_point):
- margin_model_class = self.stationary_margin_model_class if starting_point is None else self.non_stationary_margin_model_class
- estimator = self.get_estimator(margin_model_class, starting_point)
+ estimator = self.get_estimator(starting_point)
metric = estimator.result_from_fit.__getattribute__(self.RESULT_ATTRIBUTE_METRIC)
assert isinstance(metric, float)
return metric
- def get_mu1(self, starting_point):
+ def get_mu_coefs(self, starting_point):
# for the non stationary model gives the mu1 parameters that was fitted
- estimator = self.get_estimator(self.non_stationary_margin_model_class, starting_point)
+ estimator = self.get_estimator(starting_point)
margin_function = estimator.margin_function_fitted # type: LinearMarginFunction
assert isinstance(margin_function, LinearMarginFunction)
- return margin_function.mu1_temporal_trend
+ mu_coefs = [margin_function.mu_intercept, margin_function.mu_longitude_trend, margin_function.mu_latitude_trend,
+ margin_function.mu1_temporal_trend]
+ return dict(zip(self.mu_coef_names, mu_coefs))
+
+ @property
+ def mu_coef_names(self):
+ return ['mu_intercept', 'mu_longitude', 'mu_latitude', 'mu_temporal']
+
+ @property
+ def mu_coef_colors(self):
+ return ['b', 'g', 'y', 'c']
def visualize(self, ax, complete_analysis=True):
years = self.years(complete_analysis)
- # Compute metric with parallelization or not
+ # Load the estimator only once
if self.multiprocessing:
with Pool(self.nb_cores) as p:
- stationary_metric, *non_stationary_metrics = p.map(self.get_metric, [None] + years)
+ stationary_estimator, *non_stationary_estimators = p.map(self.load_estimator, [None] + years)
else:
- stationary_metric = self.get_metric(starting_point=None)
- non_stationary_metrics = [self.get_metric(year) for year in years]
+ stationary_estimator = self.load_estimator(None)
+ non_stationary_estimators = [self.load_estimator(year) for year in years]
+ self._starting_point_to_estimator[None] = stationary_estimator
+ for year, non_stationary_estimator in zip(years, non_stationary_estimators):
+ self._starting_point_to_estimator[year] = non_stationary_estimator
# Plot differences
+ stationary_metric, *non_stationary_metrics = [self.get_metric(starting_point) for starting_point in
+ [None] + years]
difference = [m - stationary_metric for m in non_stationary_metrics]
- color_difference = 'b'
- ax.plot(years, difference, color_difference + 'o-')
- ax.set_ylabel(self.RESULT_ATTRIBUTE_METRIC + ' difference', color=color_difference)
+ color_difference = 'r'
+ label_difference = self.RESULT_ATTRIBUTE_METRIC + ' difference'
+ ax.plot(years, difference, color_difference + 'x-', label=label_difference)
+ ax.set_ylabel(label_difference, color=color_difference, )
# Plot zero line
# years_line = [years[0] -10, years[-1] + 10]
- ax.plot(years, [0 for _ in years], 'k-', label='zero line')
+ # ax.plot(years, [0 for _ in years], 'kx-', label='zero line')
# Plot significative line corresponding to 0.05 relevance
alpha = 0.05
significative_deviance = chi2.ppf(q=1 - alpha, df=1)
- ax.plot(years, [significative_deviance for _ in years], 'g-', label='significative line')
+ ax.plot(years, [significative_deviance for _ in years], 'mx-', label='significative line')
+
+ all_deviance_data = [significative_deviance] + difference
+ min_deviance_data, max_deviance_data = min(all_deviance_data), max(all_deviance_data)
# Plot the mu1 parameter
- mu1_trends = [self.get_mu1(starting_point=year) for year in years]
+ mu_trends = [self.get_mu_coefs(starting_point=year) for year in years]
+ mus = {mu_coef_name: [mu_trend[mu_coef_name] for mu_trend in mu_trends] for mu_coef_name in self.mu_coef_names}
+
ax2 = ax.twinx()
- color_mu1 = 'c'
- if self.verbose:
- print('mu1 trends:', mu1_trends, '\n')
- ax2.plot(years, mu1_trends, color_mu1 + 'o-')
- ax2.set_ylabel('mu1 parameter', color=color_mu1)
+ for j, (mu_coef_name, mu_coef_values) in enumerate(mus.items()):
+ color_mu_coef = self.mu_coef_colors[j]
+ if self.verbose:
+ print(mu_coef_name, mu_coef_values)
+ ax2.plot(years, mu_coef_values, color_mu_coef + 'o-', label=mu_coef_name)
+ # ax2.set_ylabel(mu_coef_name, color=color_mu_coef)
+
+ df_mus = pd.DataFrame(mus)
+ min_mus, max_mus = df_mus.min().min(), df_mus.max().max()
+ min_global, max_global = min(min_deviance_data, min_mus), max(max_deviance_data, max_mus)
+ ax2.set_ylim(min_global, max_global)
+
+ ax.set_xlabel('starting year for the linear trend of {}'.format(self.mu_coef_names[-1]))
+ if min_mus < 0.0 < max_mus:
+ align_yaxis_on_zero(ax2, ax)
- ax.set_xlabel('starting year for the linear trend of mu1')
- if min(mu1_trends) < 0.0 < max(mu1_trends):
- align_yaxis_on_zero(ax, ax2)
title = self.display_name
ax.set_title(title)
- ax.legend()
+ ax.grid()
+ ax.legend(loc=6)
+ ax2.legend(loc=7)
def years(self, complete_analysis=True):
# Define the year_min and year_max for the starting point
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 35a0428f9e5020c9c52063e740a5e08a9578b004..8e5062e518ef3b03d3c1581035b00c83a70c14b0 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
@@ -177,11 +177,11 @@ class StudyVisualizer(object):
trend_tests = [ConditionalIndedendenceLocationTrendTest(dataset=self.dataset, verbose=verbose,
multiprocessing=multiprocessing)]
- max_stable_models = load_test_max_stable_models(default_covariance_function=self.default_covariance_function)
- for max_stable_model in [max_stable_models[1], max_stable_models[-2]]:
- trend_tests.append(MaxStableLocationTrendTest(dataset=self.dataset,
- max_stable_model=max_stable_model, verbose=verbose,
- multiprocessing=multiprocessing))
+ # max_stable_models = load_test_max_stable_models(default_covariance_function=self.default_covariance_function)
+ # for max_stable_model in [max_stable_models[1], max_stable_models[-2]]:
+ # trend_tests.append(MaxStableLocationTrendTest(dataset=self.dataset,
+ # max_stable_model=max_stable_model, verbose=verbose,
+ # multiprocessing=multiprocessing))
nb_trend_tests = len(trend_tests)
fig, axes = plt.subplots(1, nb_trend_tests, figsize=self.figsize)
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 4c72f5faaf076bdcce3c541ab052d1312bb5c196..ebbfd62fd4dd3205c58b41c3a59c16eaa1cddd4a 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
@@ -59,14 +59,6 @@ class LinearMarginFunction(ParametricMarginFunction):
coef_dict.update(coef.coef_dict(dims, self.idx_to_coefficient_name(self.coordinates)))
return coef_dict
- @property
- def mu1_temporal_trend(self):
- return self.coef_dict[LinearCoef.coef_template_str(ExtremeParams.LOC, AbstractCoordinates.COORDINATE_T).format(1)]
-
- @property
- def mu0(self):
- return self.coef_dict[LinearCoef.coef_template_str(ExtremeParams.LOC, LinearCoef.INTERCEPT_NAME).format(1)]
-
@property
def form_dict(self) -> Dict[str, str]:
form_dict = {}
@@ -87,3 +79,21 @@ class LinearMarginFunction(ParametricMarginFunction):
assert not any(['1' in formula for formula in temporal_form.values()])
form_dict.update(temporal_form)
return form_dict
+
+ # Properties for the location parameter
+
+ @property
+ def mu1_temporal_trend(self):
+ return self.coef_dict[LinearCoef.coef_template_str(ExtremeParams.LOC, AbstractCoordinates.COORDINATE_T).format(1)]
+
+ @property
+ def mu_intercept(self):
+ return self.coef_dict[LinearCoef.coef_template_str(ExtremeParams.LOC, LinearCoef.INTERCEPT_NAME).format(1)]
+
+ @property
+ def mu_longitude_trend(self):
+ return self.coef_dict[LinearCoef.coef_template_str(ExtremeParams.LOC, AbstractCoordinates.COORDINATE_X).format(2)]
+
+ @property
+ def mu_latitude_trend(self):
+ return self.coef_dict[LinearCoef.coef_template_str(ExtremeParams.LOC, AbstractCoordinates.COORDINATE_Y).format(3)]
\ No newline at end of file
diff --git a/spatio_temporal_dataset/spatio_temporal_observations/abstract_spatio_temporal_observations.py b/spatio_temporal_dataset/spatio_temporal_observations/abstract_spatio_temporal_observations.py
index e65c488e289f21786eb1a3b1ab6db475abfbb35e..7ff8ecb67d5889d166ec4ad134c00876f0ba6ed8 100644
--- a/spatio_temporal_dataset/spatio_temporal_observations/abstract_spatio_temporal_observations.py
+++ b/spatio_temporal_dataset/spatio_temporal_observations/abstract_spatio_temporal_observations.py
@@ -105,6 +105,15 @@ class AbstractSpatioTemporalObservations(object):
def __str__(self) -> str:
return self._df_maxima.__str__()
+ def print_summary(self):
+ # Write a summary of observations
+ df = self.df_maxima_gev
+ print('Observations summary:', ' ', end='')
+ print('Mean value:', df.mean().mean(), ' ', end='')
+ print('Min value:', df.min().min(), ' ', end='')
+ print('Max value:', df.max().max(), ' ', end='')
+ print('# of zero values:', df.size - np.count_nonzero(df.values), '\n')
+
@_df_maxima.setter
def _df_maxima(self, value):
self.__df_maxima = value
diff --git a/test/test_experiment/test_coordinate_sensitivity.py b/test/test_experiment/test_coordinate_sensitivity.py
index 2672d7a6768e68ec9b52f1afca76e4bf5f780cf1..24e24d44d6021a50eda419a97892575eca6711ec 100644
--- a/test/test_experiment/test_coordinate_sensitivity.py
+++ b/test/test_experiment/test_coordinate_sensitivity.py
@@ -17,7 +17,8 @@ class TestCoordinateSensitivity(unittest.TestCase):
def test_coordinate_normalization_sensitivity(self):
altitudes = [300, 600, 900, 1200, 2100, 3000][-1:]
- transformation_classes = [None, BetweenZeroAndOneNormalization, BetweenZeroAndOneNormalizationMinEpsilon, BetweenZeroAndOneNormalizationMaxEpsilon][1:2]
+ transformation_classes = [None, BetweenZeroAndOneNormalization, BetweenZeroAndOneNormalizationMinEpsilon,
+ BetweenZeroAndOneNormalizationMaxEpsilon][1:2]
study_classes = [CrocusSwe]
for study in study_iterator_global(study_classes, altitudes=altitudes, verbose=False):
@@ -28,19 +29,19 @@ class TestCoordinateSensitivity(unittest.TestCase):
study_visualizer.temporal_non_stationarity = True
trend_test = ConditionalIndedendenceLocationTrendTest(study_visualizer.dataset)
years = [1960, 1990]
- mu1s = [trend_test.get_mu1(year) for year in years]
+ mu1s = [trend_test.get_mu_coefs(year)['mu_temporal'] for year in years]
if self.DISPLAY:
print('Stationary')
- stationary_est = trend_test.get_estimator(trend_test.stationary_margin_model_class,
- starting_point=None)
+ stationary_est = trend_test.load_estimator(trend_test.stationary_margin_model_class,
+ starting_point=None)
print(stationary_est.result_from_fit.convergence)
print(stationary_est.margin_function_fitted.coef_dict)
print('Non Stationary')
- non_stationary_est = trend_test.get_estimator(trend_test.non_stationary_margin_model_class,
- starting_point=1960)
+ non_stationary_est = trend_test.load_estimator(trend_test.non_stationary_margin_model_class,
+ starting_point=1960)
print(non_stationary_est.result_from_fit.convergence)
- non_stationary_est = trend_test.get_estimator(trend_test.non_stationary_margin_model_class,
- starting_point=1990)
+ non_stationary_est = trend_test.load_estimator(trend_test.non_stationary_margin_model_class,
+ starting_point=1990)
print(non_stationary_est.result_from_fit.convergence)
print(non_stationary_est.margin_function_fitted.coef_dict)
print(get_display_name_from_object_type(transformation_class), 'mu1s: ', mu1s)