From 9b0050b19e8f5271e1f9795b8c7c8509692c43fe Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 25 Jun 2020 16:13:59 +0200
Subject: [PATCH] [contrasting]add shape plot with the name of the selected
models. the selection for the model with minimal aic only takes into account
models with with finite aic (there are often several models with infinite
aic)
---
.../altitudinal_models.py | 21 ++++++++
.../polynomial_margin_model/utils.py | 2 +
.../altitudes_fit/altitudes_studies.py | 29 +++++++----
.../altitudes_fit/main_altitudes_studies.py | 5 +-
...es_visualizer_for_non_stationary_models.py | 31 +++++++++--
.../one_fold_analysis/one_fold_fit.py | 51 +++++++++++++------
6 files changed, 109 insertions(+), 30 deletions(-)
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/altitudinal_models.py b/extreme_fit/model/margin_model/polynomial_margin_model/altitudinal_models.py
index 6008f207..fceadfce 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/altitudinal_models.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/altitudinal_models.py
@@ -17,6 +17,27 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
def param_name_to_list_dim_and_degree(self):
raise NotImplementedError
+ def dims(self, param_name):
+ return [d for d, _ in self.param_name_to_list_dim_and_degree[param_name]]
+
+ def dim_to_str_number(self, param_name, dim):
+ list_dim_and_degree = self.param_name_to_list_dim_and_degree[param_name]
+ dims = [d for d, _ in list_dim_and_degree]
+ if dim not in dims:
+ return '0'
+ else:
+ idx = dims.index(dim)
+ return str(list_dim_and_degree[idx][1])
+
+ @property
+ def name_str(self):
+ name = 'Gev'
+ name += self.dim_to_str_number(GevParams.LOC, self.coordinates.idx_temporal_coordinates)
+ name += self.dim_to_str_number(GevParams.SCALE, self.coordinates.idx_temporal_coordinates)
+ if isinstance(self, AbstractAddCrossTermForLocation):
+ name += 'x'
+ return name
+
class StationaryAltitudinal(AbstractAltitudinalModel):
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
index 78346c36..a788af48 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
@@ -26,6 +26,8 @@ ALTITUDINAL_MODELS = [
NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation,
][:]
+
+
MODELS_THAT_SHOULD_RAISE_AN_ASSERTION_ERROR = [NonStationaryLocationSpatioTemporalLinearityModelAssertError1,
NonStationaryLocationSpatioTemporalLinearityModelAssertError2,
NonStationaryLocationSpatioTemporalLinearityModelAssertError3]
diff --git a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
index 9bbadcf2..acd7f9d5 100644
--- a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
@@ -45,24 +45,31 @@ class AltitudesStudies(object):
def spatio_temporal_dataset(self, massif_name, s_split_spatial: pd.Series = None,
s_split_temporal: pd.Series = None):
coordinate_values_to_maxima = {}
+ massif_altitudes = []
for altitude in self.altitudes:
study = self.altitude_to_study[altitude]
if massif_name in study.study_massif_names:
+ massif_altitudes.append(altitude)
for year, maxima in zip(study.ordered_years, study.massif_name_to_annual_maxima[massif_name]):
coordinate_values_to_maxima[(altitude, year)] = [maxima]
- coordinates = self.spatio_temporal_coordinates(s_split_spatial, s_split_temporal)
+ coordinates = self.spatio_temporal_coordinates(s_split_spatial, s_split_temporal, massif_altitudes)
observations = AnnualMaxima.from_coordinates(coordinates, coordinate_values_to_maxima)
return AbstractDataset(observations=observations, coordinates=coordinates)
# Coordinates Loader
- def spatio_temporal_coordinates(self, s_split_spatial: pd.Series = None, s_split_temporal: pd.Series = None):
+ def spatio_temporal_coordinates(self, s_split_spatial: pd.Series = None, s_split_temporal: pd.Series = None,
+ massif_altitudes=None):
+ if massif_altitudes is None or set(massif_altitudes) == set(self.altitudes):
+ spatial_coordinates = self.spatial_coordinates
+ else:
+ spatial_coordinates = self.spatial_coordinates_for_altitudes(massif_altitudes)
slicer_class = get_slicer_class_from_s_splits(s_split_spatial, s_split_temporal)
return AbstractSpatioTemporalCoordinates(slicer_class=slicer_class,
s_split_spatial=s_split_spatial,
s_split_temporal=s_split_temporal,
transformation_class=self.spatial_transformation_class,
- spatial_coordinates=self.spatial_coordinates,
+ spatial_coordinates=spatial_coordinates,
temporal_coordinates=self.temporal_coordinates)
@cached_property
@@ -73,7 +80,10 @@ class AltitudesStudies(object):
@cached_property
def spatial_coordinates(self):
- return AbstractSpatialCoordinates.from_list_x_coordinates(x_coordinates=self.altitudes,
+ return self.spatial_coordinates_for_altitudes(self.altitudes)
+
+ def spatial_coordinates_for_altitudes(self, altitudes):
+ return AbstractSpatialCoordinates.from_list_x_coordinates(x_coordinates=altitudes,
transformation_class=self.temporal_transformation_class)
@cached_property
@@ -124,12 +134,13 @@ class AltitudesStudies(object):
ax = plt.gca()
self.run_for_each_massif(self._plot_mean_maxima_against_altitude, massif_names, ax=ax, std=std, change=change)
ax.legend(prop={'size': 7}, ncol=3)
- moment = 'Mean' if not std else 'Std'
+ moment = ''
if change is None:
- moment += ' relative change for'
- elif change:
- moment += ' change for'
- plot_name = '{} annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class])
+ moment += ' Relative'
+ if change is True or change is None:
+ moment += ' change (between two block of 30 years) for'
+ moment = 'mean' if not std else 'std'
+ plot_name = '{} of annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class])
ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
ax.set_xlabel('altitudes', fontsize=15)
lim_down, lim_up = ax.get_ylim()
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index 5b3e651a..8c800b5d 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -15,6 +15,7 @@ def plot_altitudinal_fit(studies, massif_names=None):
show=False)
visualizer.plot_mean()
visualizer.plot_relative_change()
+ visualizer.plot_shape_map()
def plot_time_series(studies, massif_names=None):
@@ -34,7 +35,9 @@ def main():
study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
SafranPrecipitation7Days][:]
study_classes = [SafranPrecipitation1Day, SafranSnowfall1Day, SafranSnowfall3Days, SafranPrecipitation3Days][:1]
- massif_names = ['Belledonne']
+ massif_names = None
+ massif_names = ['Aravis']
+ # massif_names = ['Chartreuse', 'Belledonne']
for study_class in study_classes:
studies = AltitudesStudies(study_class, altitudes, season=Season.winter_extended)
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 0ea2b5e9..8b7aeb81 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -10,6 +10,7 @@ from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils imp
from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer
from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
AbstractSpatioTemporalPolynomialModel
+from extreme_fit.model.margin_model.utils import MarginFitMethod
from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import \
OneFoldFit
@@ -20,16 +21,17 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
def __init__(self, studies: AltitudesStudies,
model_classes: List[AbstractSpatioTemporalPolynomialModel],
show=False,
- massif_names=None):
- study = studies.study
+ massif_names=None,
+ fit_method=MarginFitMethod.extremes_fevd_mle):
+ super().__init__(studies.study, show=show, save_to_file=not show)
self.massif_names = massif_names if massif_names is not None else self.study.study_massif_names
self.studies = studies
self.non_stationary_models = model_classes
- super().__init__(study, show=show, save_to_file=not show)
+ self.fit_method = fit_method
self.massif_name_to_one_fold_fit = {}
for massif_name in self.massif_names:
dataset = studies.spatio_temporal_dataset(massif_name=massif_name)
- old_fold_fit = OneFoldFit(dataset, model_classes)
+ old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method)
self.massif_name_to_one_fold_fit[massif_name] = old_fold_fit
def plot_mean(self):
@@ -58,3 +60,24 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
ax.tick_params(axis='both', which='major', labelsize=13)
self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
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):
+ super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap, 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
+ for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items()}
+
+ @property
+ def massif_name_to_name(self):
+ 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_shape_map(self):
+ self.plot_abstract_fast(self.massif_name_to_shape,
+ label='Shape plot for {} {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)],
+ self.study.variable_unit),
+ add_x_label=False, graduation=0.1, massif_name_to_text=self.massif_name_to_name)
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index 312caafb..1ed8d679 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -7,7 +7,8 @@ from extreme_fit.model.margin_model.utils import MarginFitMethod
class OneFoldFit(object):
- def __init__(self, dataset, models_classes, fit_method=MarginFitMethod.extremes_fevd_mle):
+ def __init__(self, massif_name, dataset, models_classes, fit_method=MarginFitMethod.extremes_fevd_mle):
+ self.massif_name = massif_name
self.dataset = dataset
self.models_classes = models_classes
self.fit_method = fit_method
@@ -18,28 +19,17 @@ class OneFoldFit(object):
self.model_class_to_estimator[model_class] = fitted_linear_margin_estimator_short(model_class=model_class,
dataset=self.dataset,
fit_method=self.fit_method)
- # Some display
- for estimator in self.model_class_to_estimator.values():
- print(estimator.result_from_model_fit.aic)
-
- @cached_property
- def best_estimator(self):
- sorted_estimators = sorted([estimator for estimator in self.model_class_to_estimator.values()],
- key=lambda e: e.result_from_model_fit.aic)
- estimator_that_minimizes_aic = sorted_estimators[0]
- return estimator_that_minimizes_aic
-
- @property
- def best_function_from_fit(self):
- return self.best_estimator.function_from_fit
def mean(self, altitudes, year=2019):
return [self.get_mean(altitude, year) for altitude in altitudes]
def get_mean(self, altitude, year):
+ return self.get_gev_params(altitude, year).mean
+
+ def get_gev_params(self, altitude, year):
coordinate = np.array([altitude, year])
gev_params = self.best_function_from_fit.get_params(coordinate, is_transformed=False)
- return gev_params.mean
+ return gev_params
def relative_changes_in_the_mean(self, altitudes, year=2019, nb_years=50):
relative_changes = []
@@ -50,6 +40,35 @@ class OneFoldFit(object):
relative_changes.append(relative_change)
return relative_changes
+ # Minimizing the AIC and some properties
+
+ @cached_property
+ def sorted_estimators_with_finite_aic(self):
+ estimators = list(self.model_class_to_estimator.values())
+ estimators_with_finite_aic = []
+ for estimator in estimators:
+ try:
+ aic = estimator.aic()
+ print(self.massif_name, estimator.margin_model.name_str, aic)
+ estimators_with_finite_aic.append(estimator)
+ except AssertionError:
+ print(self.massif_name, estimator.margin_model.name_str, 'infinite aic')
+ sorted_estimators_with_finite_aic = sorted([estimator for estimator in estimators_with_finite_aic],
+ key=lambda e: e.aic())
+ return sorted_estimators_with_finite_aic
+ @cached_property
+ def best_estimator(self):
+ return self.sorted_estimators_with_finite_aic[0]
+ @property
+ def best_shape(self):
+ return self.get_gev_params(1000, year=2019).shape
+ @property
+ def best_function_from_fit(self):
+ return self.best_estimator.function_from_fit
+
+ @property
+ def best_name(self):
+ return self.best_estimator.margin_model.name_str
--
GitLab