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 6008f2074117cf4020a0bf8b474ac7d02bb69ae5..fceadfceb8331cc2fcd611b1a24246c4ff67fadc 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 78346c36ca15ec4e3635998b3fd9d1d8f80bcef3..a788af48c84bc346c6da28303c18a42ac94fbd6f 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 9bbadcf2bf62d4a62c664d8ee936a98665f9299b..acd7f9d5c010670fe8a4b752ca330bb9b89fc8ad 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 5b3e651af99b9dd2491da77cc6b922d03950cbb8..8c800b5dd67949bd31c9b85a8a363eb0ffc6d764 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 0ea2b5e9bb366e6f74ffc26841e918b170f95d86..8b7aeb81ac60daff2100fa569ffc7707412940d2 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 312caafb98d14da434a8beeeb3c2cbe5436276b3..1ed8d6795523d4ec762fed46b0f3b142d16c02c7 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