From 9473d0f3bdc645a0488e548efade5abfc4ebd0d1 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 29 Mar 2021 12:15:11 +0200
Subject: [PATCH] [projection snowfall] refactor first_year and last_year in
OneFoldFit. and apply the code to check the difference in spatial patterns
between safran 2020 reanalysis and ADAMont v2
---
.../ensemble_fit/abstract_ensemble_fit.py | 4 --
.../one_fold_fit_merge.py | 4 +-
.../visualizer_merge.py | 4 +-
...es_visualizer_for_non_stationary_models.py | 57 +++----------------
extreme_trend/one_fold_fit/one_fold_fit.py | 12 ++--
.../main_altitudes_studies.py | 20 +++----
.../results/main_projections_ensemble.py | 18 +++---
test/test_extreme_trend/test_one_fold_fit.py | 24 ++++++--
8 files changed, 60 insertions(+), 83 deletions(-)
diff --git a/extreme_trend/ensemble_fit/abstract_ensemble_fit.py b/extreme_trend/ensemble_fit/abstract_ensemble_fit.py
index b9348041..ecaa0db4 100644
--- a/extreme_trend/ensemble_fit/abstract_ensemble_fit.py
+++ b/extreme_trend/ensemble_fit/abstract_ensemble_fit.py
@@ -27,10 +27,6 @@ class AbstractEnsembleFit(object):
self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test
self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
- # Set appropriate setting
- OneFoldFit.last_year = 2100
- OneFoldFit.nb_years = OneFoldFit.last_year - 2005
-
@property
def altitudes(self):
raise self.visualizer_list.studies.altitudes
diff --git a/extreme_trend/ensemble_fit/independent_ensemble_fit/one_fold_fit_merge.py b/extreme_trend/ensemble_fit/independent_ensemble_fit/one_fold_fit_merge.py
index 4db04016..2345fca5 100644
--- a/extreme_trend/ensemble_fit/independent_ensemble_fit/one_fold_fit_merge.py
+++ b/extreme_trend/ensemble_fit/independent_ensemble_fit/one_fold_fit_merge.py
@@ -8,13 +8,15 @@ from extreme_trend.one_fold_fit.one_fold_fit import OneFoldFit
class OneFoldFitMerge(OneFoldFit):
def __init__(self, one_fold_fit_list: List[OneFoldFit], massif_name, altitude_class, temporal_covariate_for_fit,
- merge_function=np.median):
+ first_year, last_year, merge_function=np.median):
assert len(one_fold_fit_list) > 0
self.one_fold_fit_list = one_fold_fit_list
self.altitude_group = altitude_class()
self.massif_name = massif_name
self.temporal_covariate_for_fit = temporal_covariate_for_fit
self.merge_function = merge_function
+ self.first_year = first_year
+ self.last_year = last_year
def get_moment(self, altitude, temporal_covariate, order=1):
return self.merge_function([o.get_moment(altitude, temporal_covariate, order) for o in self.one_fold_fit_list])
diff --git a/extreme_trend/ensemble_fit/independent_ensemble_fit/visualizer_merge.py b/extreme_trend/ensemble_fit/independent_ensemble_fit/visualizer_merge.py
index 1f2361fb..5ed3fcda 100644
--- a/extreme_trend/ensemble_fit/independent_ensemble_fit/visualizer_merge.py
+++ b/extreme_trend/ensemble_fit/independent_ensemble_fit/visualizer_merge.py
@@ -37,7 +37,9 @@ class VisualizerMerge(AltitudesStudiesVisualizerForNonStationaryModels):
if massif_name in v.massif_name_to_one_fold_fit]
if len(one_fold_fit_list) > 0:
one_fold_fit_merge = OneFoldFitMerge(one_fold_fit_list, massif_name,
- type(self.altitude_group), self.temporal_covariate_for_fit)
+ type(self.altitude_group), self.temporal_covariate_for_fit,
+ first_year=self.study.year_min, last_year=self.study.year_max,
+ merge_function=self.merge_function)
self._massif_name_to_one_fold_fit[massif_name] = one_fold_fit_merge
@property
diff --git a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
index f856a926..7706a557 100644
--- a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -80,7 +80,10 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self.massif_name_to_massif_altitudes[massif_name] = massif_altitudes
# Load dataset
dataset = self.get_dataset(massif_altitudes, massif_name)
- old_fold_fit = OneFoldFit(massif_name, dataset, self.model_classes, self.fit_method,
+ old_fold_fit = OneFoldFit(massif_name, dataset, self.model_classes,
+ self.study.year_min,
+ self.study.year_max,
+ self.fit_method,
self.temporal_covariate_for_fit,
type(self.altitude_group),
self.display_only_model_that_pass_test,
@@ -163,7 +166,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
return self._method_name_and_order_to_massif_name_to_value[c]
def ratio_groups(self):
- return [self.ratio_uncertainty_interval_size(altitude, OneFoldFit.last_year) for altitude in
+ return [self.ratio_uncertainty_interval_size(altitude, self.first_one_fold_fit.last_year) for altitude in
self.studies.altitudes]
def ratio_uncertainty_interval_size(self, altitude, year):
@@ -264,7 +267,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
ax.legend(prop={'size': 7}, ncol=3)
moment = ' '.join(method_name.split('_'))
moment = moment.replace('moment',
- '{} in {}'.format(OneFoldFit.get_moment_str(order=order), OneFoldFit.last_year))
+ '{} in {}'.format(OneFoldFit.get_moment_str(order=order), self.first_one_fold_fit.last_year))
plot_name = 'Model {} annual maxima of {}'.format(moment,
SCM_STUDY_CLASS_TO_ABBREVIATION[self.studies.study_class])
ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
@@ -273,13 +276,6 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True)
ax.clear()
- # def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm, add_x_label=True,
- # negative_and_positive_values=True, massif_name_to_text=None):
- # plot_name = '{}{}'.format(OneFoldFit.folder_for_plots, label)
- # self.plot_map(cmap, self.fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label,
- # negative_and_positive_values,
- # massif_name_to_text)
-
@property
def massif_name_to_shape(self):
return {massif_name: one_fold_fit.best_shape
@@ -290,48 +286,9 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
return {massif_name: one_fold_fit.best_name
for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items()}
- def plot_best_coef_maps(self):
- for param_name in GevParams.PARAM_NAMES:
- coordinate_names = [AbstractCoordinates.COORDINATE_X, AbstractCoordinates.COORDINATE_T]
- dim_to_coordinate_name = dict(zip([0, 1], coordinate_names))
- for dim in [0, 1, (0, 1)]:
- coordinate_name = LinearCoef.coefficient_name(dim, dim_to_coordinate_name)
- for degree in range(4):
- coef_name = ' '.join([param_name + coordinate_name + str(degree)])
- massif_name_to_best_coef = {}
- for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
- best_coef = one_fold_fit.best_coef(param_name, dim, degree)
- if best_coef is not None:
- massif_name_to_best_coef[massif_name] = best_coef
-
- if len(massif_name_to_best_coef) > 0:
- for evaluate_coordinate in [False, True][:1]:
- if evaluate_coordinate:
- coef_name += 'evaluated at coordinates'
- for massif_name in massif_name_to_best_coef.keys():
- if AbstractCoordinates.COORDINATE_X in coordinate_name:
- massif_name_to_best_coef[massif_name] *= np.power(1000, degree)
- if AbstractCoordinates.COORDINATE_T in coordinate_name:
- massif_name_to_best_coef[massif_name] *= np.power(OneFoldFit.last_year, degree)
- self.plot_best_coef_map(coef_name.replace('_', ''), massif_name_to_best_coef)
-
- def plot_best_coef_map(self, coef_name, massif_name_to_best_coef):
- values = list(massif_name_to_best_coef.values())
- graduation = (max(values) - min(values)) / 6
- print(coef_name, graduation, max(values), min(values))
- negative_and_positive_values = (max(values) > 0) and (min(values) < 0)
- raise NotImplementedError
- self.plot_map(massif_name_to_value=massif_name_to_best_coef,
- label='Coef/{} plot for {} {}'.format(coef_name,
- SCM_STUDY_CLASS_TO_ABBREVIATION[
- type(self.study)],
- self.study.variable_unit),
- add_x_label=False, graduation=graduation, massif_name_to_text=self.massif_name_to_best_name,
- negative_and_positive_values=negative_and_positive_values)
-
def plot_shape_map(self):
- label = 'Shape parameter in {} (no unit)'.format(OneFoldFit.last_year)
+ label = 'Shape parameter in {} (no unit)'.format(self.first_one_fold_fit.last_year)
max_abs_change = self._max_abs_for_shape + 0.05
self.plot_map(massif_name_to_value=self.massif_name_to_shape,
label=label,
diff --git a/extreme_trend/one_fold_fit/one_fold_fit.py b/extreme_trend/one_fold_fit/one_fold_fit.py
index e8a59383..7cdfc50f 100644
--- a/extreme_trend/one_fold_fit/one_fold_fit.py
+++ b/extreme_trend/one_fold_fit/one_fold_fit.py
@@ -42,11 +42,9 @@ class OneFoldFit(object):
SIGNIFICANCE_LEVEL = 0.05
return_period = 100
quantile_level = 1 - (1 / return_period)
- nb_years = 60
- last_year = 2019
- last_anomaly = 2
def __init__(self, massif_name: str, dataset: AbstractDataset, models_classes,
+ first_year, last_year,
fit_method=MarginFitMethod.extremes_fevd_mle,
temporal_covariate_for_fit=None,
altitude_class=DefaultAltitudeGroup,
@@ -54,6 +52,9 @@ class OneFoldFit(object):
confidence_interval_based_on_delta_method=False,
remove_physically_implausible_models=False,
):
+ self.first_year = first_year
+ self.last_year = last_year
+ self.years_of_difference = last_year - first_year
self.remove_physically_implausible_models = remove_physically_implausible_models
self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test
@@ -64,6 +65,7 @@ class OneFoldFit(object):
self.fit_method = fit_method
self.temporal_covariate_for_fit = temporal_covariate_for_fit
+
# Fit Estimators
self.model_class_to_estimator = {}
for model_class in models_classes:
@@ -133,10 +135,10 @@ class OneFoldFit(object):
@property
def _covariate_before_and_after(self):
if self.temporal_covariate_for_fit in [None, TimeTemporalCovariate]:
- return self.last_year - self.nb_years, self.last_year
+ return self.first_year, self.last_year
elif self.temporal_covariate_for_fit is AnomalyTemperatureWithSplineTemporalCovariate:
# In 2020, we are roughly at 1 degree. Thus it natural to see the augmentation from 1 to 2 degree.
- return 1, self.last_anomaly
+ return 1, 2
else:
raise NotImplementedError
diff --git a/projects/past_extreme_snowfall/section_data_and_results/main_altitudes_studies.py b/projects/past_extreme_snowfall/section_data_and_results/main_altitudes_studies.py
index 7c389bd6..d04f3860 100644
--- a/projects/past_extreme_snowfall/section_data_and_results/main_altitudes_studies.py
+++ b/projects/past_extreme_snowfall/section_data_and_results/main_altitudes_studies.py
@@ -43,13 +43,14 @@ def main():
SafranSnowfall5Days, SafranSnowfall7Days][:1]
study_classes = [SafranSnowfall2020, SafranSnowfall2019, SafranSnowfallCenterOnDay1day,
SafranSnowfallNotCenterOnDay1day,
- SafranSnowfallCenterOnDay1dayMeanRate, SafranSnowfall1Day][1:3]
- study_classes = [SafranSnowfallNotCenterOnDay1day, SafranSnowfall2019]
+ SafranSnowfallCenterOnDay1dayMeanRate, SafranSnowfall1Day][:1]
+ # study_classes = [SafranSnowfallNotCenterOnDay1day, SafranSnowfall2019]
seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
set_seed_for_test()
- model_must_pass_the_test = False
+ model_must_pass_the_test = True
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
+ year_max = 2005
fast = False
if fast is None:
@@ -58,28 +59,27 @@ def main():
altitudes_list = altitudes_for_groups[2:3]
elif fast:
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
- massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
+ massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:]
altitudes_list = altitudes_for_groups[2:3]
else:
massif_names = None
altitudes_list = altitudes_for_groups[:]
start = time.time()
- main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_pass_the_test)
+ main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_pass_the_test, year_max=year_max)
end = time.time()
duration = str(datetime.timedelta(seconds=end - start))
print('Total duration', duration)
-def main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_pass_the_test):
+def main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_pass_the_test, year_max=None):
assert isinstance(altitudes_list, List)
assert isinstance(altitudes_list[0], List)
for season in seasons:
for study_class in study_classes:
print('Inner loop', season, study_class)
visualizer_list = load_visualizer_list(season, study_class, altitudes_list, massif_names,
- model_must_pass_the_test,
- year_max=2019)
+ model_must_pass_the_test, year_max=year_max)
plot_visualizers(massif_names, visualizer_list)
for visualizer in visualizer_list:
plot_visualizer(massif_names, visualizer)
@@ -89,10 +89,10 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_p
def plot_visualizers(massif_names, visualizer_list):
# plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
- plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list, with_significance=True)
+ plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list, with_significance=False)
# plot_shoe_plot_ratio_interval_size_against_altitude(massif_names, visualizer_list)
for relative in [True, False]:
- plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative, with_significance=True)
+ plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative, with_significance=False)
# plot_coherence_curves(['Vanoise'], visualizer_list)
pass
diff --git a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py b/projects/projected_extreme_snowfall/results/main_projections_ensemble.py
index c44d3cbf..49ea7796 100644
--- a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py
+++ b/projects/projected_extreme_snowfall/results/main_projections_ensemble.py
@@ -24,7 +24,7 @@ from extreme_fit.model.utils import set_seed_for_test
from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
- rcp_scenarios, rcm_scenarios_extended
+ rcp_scenarios, rcm_scenarios_extended, get_gcm_list
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
TimeTemporalCovariate
@@ -41,14 +41,19 @@ def main():
study_class = AdamontSnowfall
ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][1:]
temporal_covariate_for_fit = [TimeTemporalCovariate,
- AnomalyTemperatureWithSplineTemporalCovariate][1]
+ AnomalyTemperatureWithSplineTemporalCovariate][0]
set_seed_for_test()
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
- fast = True
+ fast = False
scenarios = rcp_scenarios[::-1] if fast is False else [AdamontScenario.rcp85]
scenarios = rcm_scenarios_extended[::-1]
+ scenarios = [AdamontScenario.histo]
+ gcm_to_year_min_and_year_max = {
+ gcm: (1959, 2005) for gcm in get_gcm_list(adamont_version=2)
+ }
+
for scenario in scenarios:
gcm_rcm_couples = get_gcm_rcm_couples(scenario)
if fast is None:
@@ -62,8 +67,8 @@ def main():
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
altitudes_list = altitudes_for_groups[:1]
else:
- massif_names = ['Vanoise']
- altitudes_list = altitudes_for_groups[3:]
+ massif_names = None
+ altitudes_list = altitudes_for_groups[:]
assert isinstance(gcm_rcm_couples, list)
@@ -73,7 +78,6 @@ def main():
print('Covariate is {}'.format(temporal_covariate_for_fit))
model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
- assert scenario is not AdamontScenario.histo
visualizer = VisualizerForProjectionEnsemble(
altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
@@ -82,7 +86,7 @@ def main():
massif_names=massif_names,
temporal_covariate_for_fit=temporal_covariate_for_fit,
remove_physically_implausible_models=True,
- gcm_to_year_min_and_year_max=None,
+ gcm_to_year_min_and_year_max=gcm_to_year_min_and_year_max,
)
visualizer.plot()
diff --git a/test/test_extreme_trend/test_one_fold_fit.py b/test/test_extreme_trend/test_one_fold_fit.py
index 30829f58..80e3e706 100644
--- a/test/test_extreme_trend/test_one_fold_fit.py
+++ b/test/test_extreme_trend/test_one_fold_fit.py
@@ -41,7 +41,10 @@ class TestOneFoldFit(unittest.TestCase):
dataset = self.load_dataset(study_class)
one_fold_fit = OneFoldFit(self.massif_name, dataset,
models_classes=self.model_classes, temporal_covariate_for_fit=None,
- only_models_that_pass_goodness_of_fit_test=False)
+ only_models_that_pass_goodness_of_fit_test=False,
+ first_year=1959,
+ last_year=2019
+ )
_ = one_fold_fit.best_estimator.margin_model
self.assertTrue(True)
@@ -51,7 +54,10 @@ class TestOneFoldFit(unittest.TestCase):
one_fold_fit = OneFoldFit(self.massif_name, dataset,
models_classes=self.model_classes,
temporal_covariate_for_fit=TimeTemporalCovariate,
- only_models_that_pass_goodness_of_fit_test=False)
+ only_models_that_pass_goodness_of_fit_test=False,
+ first_year=1959,
+ last_year=2019
+ )
_ = one_fold_fit.best_estimator.margin_model
self.assertTrue(True)
@@ -61,7 +67,10 @@ class TestOneFoldFit(unittest.TestCase):
one_fold_fit = OneFoldFit(self.massif_name, dataset,
models_classes=self.model_classes,
temporal_covariate_for_fit=AnomalyTemperatureWithSplineTemporalCovariate,
- only_models_that_pass_goodness_of_fit_test=False)
+ only_models_that_pass_goodness_of_fit_test=False,
+ first_year=1959,
+ last_year=2019
+ )
_ = one_fold_fit.best_estimator.margin_model
self.assertTrue(True)
@@ -75,7 +84,10 @@ class TestOneFoldFit(unittest.TestCase):
models_classes=self.model_classes,
temporal_covariate_for_fit=None,
only_models_that_pass_goodness_of_fit_test=False,
- remove_physically_implausible_models=True)
+ remove_physically_implausible_models=True,
+ first_year=1959,
+ last_year=2019
+ )
self.assertFalse(one_fold_fit.has_at_least_one_valid_model)
def test_assertion_error_for_a_specific_case(self):
@@ -94,7 +106,9 @@ class TestOneFoldFit(unittest.TestCase):
temporal_covariate_for_fit=AnomalyTemperatureWithSplineTemporalCovariate,
altitude_class=VeyHighAltitudeGroup,
only_models_that_pass_goodness_of_fit_test=False,
- remove_physically_implausible_models=True)
+ remove_physically_implausible_models=True,
+ first_year=1959,
+ last_year=2019)
self.assertTrue(one_fold_fit.has_at_least_one_valid_model)
if __name__ == '__main__':
--
GitLab