diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py b/extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py
index 35da1544b76b397403af77b766230dd7f9ad8836..6ff8e3b6cc71a9bc65e0a5a879c092be0abdf197 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py
@@ -60,12 +60,13 @@ class NonStationaryGumbelAltitudinalLocationQuadraticScaleLinear(AbstractGumbelA
# Add cross terms
-class NonStationaryGumbelCrossTermForLocation(AbstractAddCrossTermForLocation, AbstractGumbelAltitudinalModel, StationaryAltitudinal):
+class NonStationaryGumbelCrossTermForLocation(AbstractGumbelAltitudinalModel,
+ NonStationaryCrossTermForLocation):
pass
-class NonStationaryGumbelAltitudinalLocationLinearCrossTermForLocation(AbstractAddCrossTermForLocation, AbstractGumbelAltitudinalModel,
- NonStationaryAltitudinalLocationLinear):
+class NonStationaryGumbelAltitudinalLocationLinearCrossTermForLocation(AbstractGumbelAltitudinalModel,
+ NonStationaryAltitudinalLocationLinearCrossTermForLocation):
pass
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 baa704c6c71a173de606dec5b7fcf092c32ea27c..46242d6c64a0e900efcadb23603c603668ec2e1b 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -1,11 +1,16 @@
import matplotlib as mpl
+import matplotlib.pyplot as plt
import numpy as np
-from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+from extreme_data.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
+from projects.exceeding_snow_loads.utils import dpi_paper1_figure
+
from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
SafranSnowfall5Days, SafranSnowfall7Days, SafranPrecipitation1Day, SafranPrecipitation3Days, \
SafranPrecipitation5Days, SafranPrecipitation7Days
@@ -19,35 +24,17 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_s
def plot_altitudinal_fit(studies, massif_names=None):
model_classes = ALTITUDINAL_GEV_MODELS + ALTITUDINAL_GUMBEL_MODELS
- # model_classes = ALTITUDINAL_GUMBEL_MODELS
visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
model_classes=model_classes,
massif_names=massif_names,
show=False)
# Plot the results for the model that minimizes the individual aic
- OneFoldFit.best_estimator_minimizes_mean_aic = False
+ OneFoldFit.best_estimator_minimizes_total_aic = False
visualizer.plot_moments()
+ # visualizer.plot_best_coef_maps()
visualizer.plot_shape_map()
- # Compute the mean AIC for each model_class
- OneFoldFit.best_estimator_minimizes_mean_aic = True
- model_class_to_aic_scores = {model_class: [] for model_class in model_classes}
- for one_fold_fit in visualizer.massif_name_to_one_fold_fit.values():
- for model_class, estimator in one_fold_fit.model_class_to_estimator_with_finite_aic.items():
- aic_score_normalized = estimator.aic() / estimator.n()
- model_class_to_aic_scores[model_class].append(aic_score_normalized)
- model_class_to_mean_aic_score = {model_class: np.array(aic_scores).mean()
- for model_class, aic_scores in model_class_to_aic_scores.items()}
- print(model_class_to_mean_aic_score)
- sorted_model_class = sorted(model_classes, key=lambda m: model_class_to_mean_aic_score[m])
- best_model_class_for_mean_aic = sorted_model_class[0]
- print(best_model_class_for_mean_aic)
- for one_fold_fit in visualizer.massif_name_to_one_fold_fit.values():
- one_fold_fit.best_estimator_class_for_mean_aic = best_model_class_for_mean_aic
-
- # Plot the results for the model that minimizes the mean aic
- visualizer.plot_moments()
- visualizer.plot_shape_map()
+ plot_total_aic(model_classes, visualizer)
def plot_time_series(studies, massif_names=None):
@@ -61,13 +48,14 @@ def plot_moments(studies, massif_names=None):
def main():
- # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][4:7]
+ altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][4:7]
# altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
- altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
+ # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
SafranPrecipitation7Days][:]
- study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day, SafranSnowfall3Days, SafranPrecipitation3Days][:]
+ study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day,
+ SafranSnowfall3Days, SafranPrecipitation3Days][:1]
massif_names = None
# massif_names = ['Aravis']
# massif_names = ['Chartreuse', 'Belledonne']
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 e378f3841fcbd47e5fb73cf88aece49fee5a8b8a..63b186b625fbc5d7195a082d251e6c831f49906d 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
@@ -8,12 +8,15 @@ from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_vis
SCM_STUDY_CLASS_TO_ABBREVIATION
from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import plot_against_altitude
from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.function.param_function.linear_coef import LinearCoef
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
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
@@ -48,9 +51,8 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
ind = (min(massif_altitudes) <= altitudes_plot) & (altitudes_plot <= max(massif_altitudes))
massif_altitudes_plot = altitudes_plot[ind]
one_fold_fit = self.massif_name_to_one_fold_fit[massif_name]
- if one_fold_fit.best_estimator_has_finite_aic:
- values = one_fold_fit.__getattribute__(method_name)(massif_altitudes_plot, order=order)
- plot_against_altitude(massif_altitudes_plot, ax, massif_id, massif_name, values)
+ values = one_fold_fit.__getattribute__(method_name)(massif_altitudes_plot, order=order)
+ plot_against_altitude(massif_altitudes_plot, ax, massif_id, massif_name, values)
# Plot settings
ax.legend(prop={'size': 7}, ncol=3)
moment = ' '.join(method_name.split('_'))
@@ -76,14 +78,51 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
@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()
- if one_fold_fit.best_estimator_has_finite_aic}
+ 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_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(3):
+ 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]:
+ if evaluate_coordinate:
+ coef_name += 'evaluated at coordinates'
+ for m in massif_name_to_best_coef.values():
+ if AbstractCoordinates.COORDINATE_X in coordinate_name:
+ massif_name_to_best_coef[m] *= np.power(1000, degree)
+ if AbstractCoordinates.COORDINATE_T in coordinate_name:
+ massif_name_to_best_coef[m] *= np.power(1000, 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)
+ self.plot_abstract_fast(massif_name_to_best_coef,
+ label='{}/Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots,
+ 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_name,
+ negative_and_positive_values=negative_and_positive_values)
+
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)],
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 cc348b81518cb9d5c09c57275e4277ef5c4d0de6..d288900c7603cd60962b0b4cf6e67f2c6f938b18 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
@@ -5,6 +5,7 @@ from cached_property import cached_property
from scipy.stats import chi2
from extreme_fit.estimator.margin_estimator.utils import fitted_linear_margin_estimator_short
+from extreme_fit.function.param_function.polynomial_coef import PolynomialAllCoef, PolynomialCoef
from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import StationaryAltitudinal
from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_models import \
StationaryGumbelAltitudinal, AbstractGumbelAltitudinalModel
@@ -14,7 +15,7 @@ from root_utils import classproperty
class OneFoldFit(object):
SIGNIFICANCE_LEVEL = 0.05
- best_estimator_minimizes_mean_aic = False
+ best_estimator_minimizes_total_aic = False
def __init__(self, massif_name, dataset, models_classes, fit_method=MarginFitMethod.extremes_fevd_mle):
self.massif_name = massif_name
@@ -30,11 +31,11 @@ class OneFoldFit(object):
fit_method=self.fit_method)
# Best estimator definition
- self.best_estimator_class_for_mean_aic = None
+ self.best_estimator_class_for_total_aic = None
@classproperty
def folder_for_plots(cls):
- return 'Mean aic' if cls.best_estimator_minimizes_mean_aic else 'Individual aic'
+ return 'Total aic' if cls.best_estimator_minimizes_total_aic else 'Individual aic'
@classmethod
def get_moment_str(cls, order):
@@ -85,49 +86,56 @@ class OneFoldFit(object):
# Minimizing the AIC and some properties
@cached_property
- def sorted_estimators_with_finite_aic(self):
+ def sorted_estimators(self):
estimators = list(self.model_class_to_estimator.values())
- estimators_with_finite_aic = []
- for estimator in estimators:
- try:
- aic = estimator.aic()
- npt.assert_almost_equal(estimator.result_from_model_fit.aic, aic, decimal=5)
- print(self.massif_name, estimator.margin_model.name_str, aic)
- estimators_with_finite_aic.append(estimator)
- except (AssertionError, rpy2.rinterface.RRuntimeError):
- print(self.massif_name, estimator.margin_model.name_str, 'infinite aic')
- print('Summary {}: {}/{} fitted'.format(self.massif_name, len(estimators_with_finite_aic), len(estimators)))
- 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
+ # estimators_with_finite_aic = []
+ # for estimator in estimators:
+ # # try:
+ # aic = estimator.aic()
+ # npt.assert_almost_equal(estimator.result_from_model_fit.aic, aic, decimal=5)
+ # print(self.massif_name, estimator.margin_model.name_str, aic)
+ # estimators_with_finite_aic.append(estimator)
+ # # except (AssertionError, rpy2.rinterface.RRuntimeError) as e:
+ # # raise e
+ # # print(self.massif_name, estimator.margin_model.name_str, 'infinite aic')
+ # # print(e)
+ # print('Summary {}: {}/{} fitted'.format(self.massif_name, len(estimators_with_finite_aic), len(estimators)))
+ sorted_estimators = sorted([estimator for estimator in estimators],
+ key=lambda e: e.aic())
+ return sorted_estimators
@property
def model_class_to_estimator_with_finite_aic(self):
- return {type(estimator.margin_model): estimator for estimator in self.sorted_estimators_with_finite_aic}
-
- @cached_property
- def set_estimators_with_finite_aic(self):
- return set(self.sorted_estimators_with_finite_aic)
+ return {type(estimator.margin_model): estimator for estimator in self.sorted_estimators}
@property
def best_estimator(self):
- if self.best_estimator_minimizes_mean_aic and self.best_estimator_class_for_mean_aic is not None:
- return self.model_class_to_estimator[self.best_estimator_class_for_mean_aic]
+ if self.best_estimator_minimizes_total_aic and self.best_estimator_class_for_total_aic is not None:
+ return self.model_class_to_estimator[self.best_estimator_class_for_total_aic]
else:
- return self.sorted_estimators_with_finite_aic[0]
+ return self.sorted_estimators[0]
@property
- def best_estimator_has_finite_aic(self):
- return self.best_estimator in self.set_estimators_with_finite_aic
+ def best_margin_model(self):
+ return self.best_estimator.margin_model
+
+ @property
+ def best_function_from_fit(self):
+ return self.best_estimator.function_from_fit
@property
def best_shape(self):
# We take any altitude (altitude=1000 for instance) as the shape is constant w.r.t the altitude
return self.get_gev_params(altitude=1000, year=2019).shape
- @property
- def best_function_from_fit(self):
- return self.best_estimator.function_from_fit
+ def best_coef(self, param_name, dim, degree):
+ try:
+ coef = self.best_function_from_fit.param_name_to_coef[param_name] # type: PolynomialAllCoef
+ coef = coef.dim_to_polynomial_coef[dim] # type: PolynomialCoef
+ coef = coef.idx_to_coef[degree]
+ return coef
+ except (TypeError, KeyError):
+ return None
@property
def best_name(self):
@@ -155,7 +163,7 @@ class OneFoldFit(object):
@property
def is_significant(self) -> bool:
stationary_model_classes = [StationaryAltitudinal, StationaryGumbelAltitudinal]
- if any([isinstance(self.best_estimator.margin_model, c)
+ if any([isinstance(self.best_estimator.margin_model, c)
for c in stationary_model_classes]):
return False
else:
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
new file mode 100644
index 0000000000000000000000000000000000000000..bd28d16ec76cb96ed838f5a622e793959d1bd35a
--- /dev/null
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
@@ -0,0 +1,77 @@
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
+from extreme_data.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
+from projects.exceeding_snow_loads.utils import dpi_paper1_figure
+
+
+def plot_total_aic(model_classes, visualizer):
+ # Compute the mean AIC for each model_class
+ OneFoldFit.best_estimator_minimizes_total_aic = True
+ model_class_to_total_aic = {model_class: 0 for model_class in model_classes}
+ model_class_to_name_str = {}
+ # model_class_to_aic_scores = {model_class: [] for model_class in model_classes}
+ # model_class_to_n = {model_class: [] for model_class in model_classes}
+ for one_fold_fit in visualizer.massif_name_to_one_fold_fit.values():
+ for model_class, estimator in one_fold_fit.model_class_to_estimator_with_finite_aic.items():
+ model_class_to_total_aic[model_class] += estimator.aic()
+ model_class_to_name_str[model_class] = estimator.margin_model.name_str
+ # model_class_to_aic_scores[model_class].append(estimator.aic())
+ # model_class_to_n[model_class].append(estimator.n())
+ # model_class_to_mean_aic_score = {model_class: np.array(aic_scores).sum() / np.array(model_class_to_n[model_class]).sum()
+ # for model_class, aic_scores in model_class_to_aic_scores.items()}
+ # print(model_class_to_mean_aic_score)
+ sorted_model_class = sorted(model_classes, key=lambda m: model_class_to_total_aic[m])
+ sorted_scores = [model_class_to_total_aic[model_class] for model_class in sorted_model_class]
+ sorted_labels = [model_class_to_name_str[model_class] for model_class in sorted_model_class]
+ print(sorted_model_class)
+ print(sorted_scores)
+ print(sorted_labels)
+ best_model_class_for_total_aic = sorted_model_class[0]
+ for one_fold_fit in visualizer.massif_name_to_one_fold_fit.values():
+ one_fold_fit.best_estimator_class_for_total_aic = best_model_class_for_total_aic
+ # Plot the ranking of the
+ plot_total_aic_repartition(visualizer, sorted_labels, sorted_scores)
+ # Plot the results for the model that minimizes the mean aic
+ visualizer.plot_moments()
+ visualizer.plot_shape_map()
+ visualizer.plot_best_coef_maps()
+
+
+
+def plot_total_aic_repartition(visualizer, labels, scores):
+ """
+ Plot a single trend curves
+ :return:
+ """
+ scores = np.array(scores)
+ ax = create_adjusted_axes(1, 1)
+
+ # parameters
+ width = 3
+ size = 30
+ legend_fontsize = 30
+ labelsize = 10
+ linewidth = 3
+ x = [2 * width * (i + 1) for i in range(len(labels))]
+ ax.bar(x, scores, width=width, color='grey', edgecolor='grey', label='Total Aic',
+ linewidth=linewidth)
+ ax.legend(loc='upper right', prop={'size': size})
+ ax.set_ylabel(' Total AIC score \n '
+ 'i.e. sum of AIC score for all massifs ', fontsize=legend_fontsize)
+ ax.set_xlabel('Models', fontsize=legend_fontsize)
+ ax.tick_params(axis='both', which='major', labelsize=labelsize)
+ ax.set_xticks(x)
+ ax.set_ylim([scores.min(), scores.max()])
+ ax.yaxis.grid()
+ ax.set_xticklabels(labels)
+
+ # Save plot
+ visualizer.plot_name = 'Total AIC ranking'
+ visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
+ plt.close()