[contrasting] harmonize gumbel altitudinal models. create plot_total_aic.py....

[contrasting] harmonize gumbel altitudinal models. create plot_total_aic.py. remove exception for model with infinite aic.

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
 
 

projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py

 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']

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)],

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:

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py

+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()