Commit 9b0050b1 authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[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)
parent 937c09a0
No related merge requests found
Showing with 109 additions and 30 deletions
+109 -30
...@@ -17,6 +17,27 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel): ...@@ -17,6 +17,27 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
def param_name_to_list_dim_and_degree(self): def param_name_to_list_dim_and_degree(self):
raise NotImplementedError 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): class StationaryAltitudinal(AbstractAltitudinalModel):
......
...@@ -26,6 +26,8 @@ ALTITUDINAL_MODELS = [ ...@@ -26,6 +26,8 @@ ALTITUDINAL_MODELS = [
NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation,
][:] ][:]
MODELS_THAT_SHOULD_RAISE_AN_ASSERTION_ERROR = [NonStationaryLocationSpatioTemporalLinearityModelAssertError1, MODELS_THAT_SHOULD_RAISE_AN_ASSERTION_ERROR = [NonStationaryLocationSpatioTemporalLinearityModelAssertError1,
NonStationaryLocationSpatioTemporalLinearityModelAssertError2, NonStationaryLocationSpatioTemporalLinearityModelAssertError2,
NonStationaryLocationSpatioTemporalLinearityModelAssertError3] NonStationaryLocationSpatioTemporalLinearityModelAssertError3]
......
...@@ -45,24 +45,31 @@ class AltitudesStudies(object): ...@@ -45,24 +45,31 @@ class AltitudesStudies(object):
def spatio_temporal_dataset(self, massif_name, s_split_spatial: pd.Series = None, def spatio_temporal_dataset(self, massif_name, s_split_spatial: pd.Series = None,
s_split_temporal: pd.Series = None): s_split_temporal: pd.Series = None):
coordinate_values_to_maxima = {} coordinate_values_to_maxima = {}
massif_altitudes = []
for altitude in self.altitudes: for altitude in self.altitudes:
study = self.altitude_to_study[altitude] study = self.altitude_to_study[altitude]
if massif_name in study.study_massif_names: 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]): for year, maxima in zip(study.ordered_years, study.massif_name_to_annual_maxima[massif_name]):
coordinate_values_to_maxima[(altitude, year)] = [maxima] 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) observations = AnnualMaxima.from_coordinates(coordinates, coordinate_values_to_maxima)
return AbstractDataset(observations=observations, coordinates=coordinates) return AbstractDataset(observations=observations, coordinates=coordinates)
# Coordinates Loader # 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) slicer_class = get_slicer_class_from_s_splits(s_split_spatial, s_split_temporal)
return AbstractSpatioTemporalCoordinates(slicer_class=slicer_class, return AbstractSpatioTemporalCoordinates(slicer_class=slicer_class,
s_split_spatial=s_split_spatial, s_split_spatial=s_split_spatial,
s_split_temporal=s_split_temporal, s_split_temporal=s_split_temporal,
transformation_class=self.spatial_transformation_class, transformation_class=self.spatial_transformation_class,
spatial_coordinates=self.spatial_coordinates, spatial_coordinates=spatial_coordinates,
temporal_coordinates=self.temporal_coordinates) temporal_coordinates=self.temporal_coordinates)
@cached_property @cached_property
...@@ -73,7 +80,10 @@ class AltitudesStudies(object): ...@@ -73,7 +80,10 @@ class AltitudesStudies(object):
@cached_property @cached_property
def spatial_coordinates(self): 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) transformation_class=self.temporal_transformation_class)
@cached_property @cached_property
...@@ -124,12 +134,13 @@ class AltitudesStudies(object): ...@@ -124,12 +134,13 @@ class AltitudesStudies(object):
ax = plt.gca() ax = plt.gca()
self.run_for_each_massif(self._plot_mean_maxima_against_altitude, massif_names, ax=ax, std=std, change=change) 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) ax.legend(prop={'size': 7}, ncol=3)
moment = 'Mean' if not std else 'Std' moment = ''
if change is None: if change is None:
moment += ' relative change for' moment += ' Relative'
elif change: if change is True or change is None:
moment += ' change for' moment += ' change (between two block of 30 years) for'
plot_name = '{} annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class]) 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_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
ax.set_xlabel('altitudes', fontsize=15) ax.set_xlabel('altitudes', fontsize=15)
lim_down, lim_up = ax.get_ylim() lim_down, lim_up = ax.get_ylim()
......
...@@ -15,6 +15,7 @@ def plot_altitudinal_fit(studies, massif_names=None): ...@@ -15,6 +15,7 @@ def plot_altitudinal_fit(studies, massif_names=None):
show=False) show=False)
visualizer.plot_mean() visualizer.plot_mean()
visualizer.plot_relative_change() visualizer.plot_relative_change()
visualizer.plot_shape_map()
def plot_time_series(studies, massif_names=None): def plot_time_series(studies, massif_names=None):
...@@ -34,7 +35,9 @@ def main(): ...@@ -34,7 +35,9 @@ def main():
study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days, study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
SafranPrecipitation7Days][:] SafranPrecipitation7Days][:]
study_classes = [SafranPrecipitation1Day, SafranSnowfall1Day, SafranSnowfall3Days, SafranPrecipitation3Days][:1] 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: for study_class in study_classes:
studies = AltitudesStudies(study_class, altitudes, season=Season.winter_extended) studies = AltitudesStudies(study_class, altitudes, season=Season.winter_extended)
......
...@@ -10,6 +10,7 @@ from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils imp ...@@ -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_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 \ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
AbstractSpatioTemporalPolynomialModel 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.altitudes_studies import AltitudesStudies
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import \ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import \
OneFoldFit OneFoldFit
...@@ -20,16 +21,17 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): ...@@ -20,16 +21,17 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
def __init__(self, studies: AltitudesStudies, def __init__(self, studies: AltitudesStudies,
model_classes: List[AbstractSpatioTemporalPolynomialModel], model_classes: List[AbstractSpatioTemporalPolynomialModel],
show=False, show=False,
massif_names=None): massif_names=None,
study = studies.study 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.massif_names = massif_names if massif_names is not None else self.study.study_massif_names
self.studies = studies self.studies = studies
self.non_stationary_models = model_classes 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 = {} self.massif_name_to_one_fold_fit = {}
for massif_name in self.massif_names: for massif_name in self.massif_names:
dataset = studies.spatio_temporal_dataset(massif_name=massif_name) 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 self.massif_name_to_one_fold_fit[massif_name] = old_fold_fit
def plot_mean(self): def plot_mean(self):
...@@ -58,3 +60,24 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): ...@@ -58,3 +60,24 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
ax.tick_params(axis='both', which='major', labelsize=13) ax.tick_params(axis='both', which='major', labelsize=13)
self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show) self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
ax.clear() 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)
...@@ -7,7 +7,8 @@ from extreme_fit.model.margin_model.utils import MarginFitMethod ...@@ -7,7 +7,8 @@ from extreme_fit.model.margin_model.utils import MarginFitMethod
class OneFoldFit(object): 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.dataset = dataset
self.models_classes = models_classes self.models_classes = models_classes
self.fit_method = fit_method self.fit_method = fit_method
...@@ -18,28 +19,17 @@ class OneFoldFit(object): ...@@ -18,28 +19,17 @@ class OneFoldFit(object):
self.model_class_to_estimator[model_class] = fitted_linear_margin_estimator_short(model_class=model_class, self.model_class_to_estimator[model_class] = fitted_linear_margin_estimator_short(model_class=model_class,
dataset=self.dataset, dataset=self.dataset,
fit_method=self.fit_method) 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): def mean(self, altitudes, year=2019):
return [self.get_mean(altitude, year) for altitude in altitudes] return [self.get_mean(altitude, year) for altitude in altitudes]
def get_mean(self, altitude, year): 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]) coordinate = np.array([altitude, year])
gev_params = self.best_function_from_fit.get_params(coordinate, is_transformed=False) 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): def relative_changes_in_the_mean(self, altitudes, year=2019, nb_years=50):
relative_changes = [] relative_changes = []
...@@ -50,6 +40,35 @@ class OneFoldFit(object): ...@@ -50,6 +40,35 @@ class OneFoldFit(object):
relative_changes.append(relative_change) relative_changes.append(relative_change)
return relative_changes 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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment