[contrasting] add 4 order model. add models with only scale linearity in the altitudes.

extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py

@@ -71,7 +71,7 @@ class LinearMarginEstimator(AbstractMarginEstimator):
 
     def aic(self, split=Split.all):
         aic = 2 * self.margin_model.nb_params + 2 * self.nllh(split=split)
-        npt.assert_almost_equal(self.result_from_model_fit.aic, aic, decimal=5)
+        npt.assert_almost_equal(self.result_from_model_fit.aic, aic, decimal=1)
         return aic
 
     def n(self, split=Split.all):

extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py

@@ -22,6 +22,8 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
         return [d for d, _ in self.param_name_to_list_dim_and_degree[param_name]]
 
     def dim_to_str_number(self, param_name, dim):
+        if param_name not in self.param_name_to_list_dim_and_degree:
+            return '0'
         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:
@@ -90,6 +92,15 @@ class NonStationaryAltitudinalLocationCubic(AbstractAltitudinalModel):
             GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
         }
 
+class NonStationaryAltitudinalLocationOrder4(AbstractAltitudinalModel):
+
+    @property
+    def param_name_to_list_dim_and_degree(self):
+        return {
+            GevParams.LOC: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 4)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
+        }
+
 
 class NonStationaryAltitudinalLocationLinearScaleLinear(AbstractAltitudinalModel):
 
@@ -118,9 +129,11 @@ class AbstractAddCrossTermForLocation(AbstractAltitudinalModel):
     @property
     def param_name_to_list_dim_and_degree_for_margin_function(self):
         d = self.param_name_to_list_dim_and_degree
-        assert 1 <= len(d[GevParams.LOC]) <= 2
-        assert self.coordinates.idx_x_coordinates == d[GevParams.LOC][0][0]
-        d[GevParams.LOC].insert(1, ((self.coordinates.idx_x_coordinates, self.coordinates.idx_temporal_coordinates), 1))
+        cross_term = ((self.coordinates.idx_x_coordinates, self.coordinates.idx_temporal_coordinates), 1)
+        if GevParams.LOC in d and self.coordinates.idx_x_coordinates == d[GevParams.LOC][0][0]:
+            d[GevParams.LOC].insert(1, cross_term)
+        else:
+            d[GevParams.LOC] = [cross_term]
         return d
 
 class AbstractAddCrossTermForScale(AbstractAltitudinalModel):
@@ -155,6 +168,11 @@ class NonStationaryAltitudinalLocationCubicCrossTermForLocation(AbstractAddCross
                                                                 ):
     pass
 
+class NonStationaryAltitudinalLocationOrder4CrossTermForLocation(AbstractAddCrossTermForLocation,
+                                                                NonStationaryAltitudinalLocationOrder4,
+                                                                ):
+    pass
+
 class NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation,
                                                                             NonStationaryAltitudinalLocationLinearScaleLinear):
     pass

extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_only_altitude_and_scale.py

+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import AbstractAltitudinalModel, \
+    AbstractAddCrossTermForLocation
+from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import PolynomialMarginModel
+from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
+    AbstractSpatioTemporalPolynomialModel
+
+
+class AltitudinalOnlyScale(AbstractAltitudinalModel):
+    pass
+
+class StationaryAltitudinalOnlyScale(AltitudinalOnlyScale):
+
+    @property
+    def param_name_to_list_dim_and_degree(self):
+        return {
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)],
+        }
+
+
+
+
+class NonStationaryAltitudinalOnlyScaleLocationLinear(AltitudinalOnlyScale):
+
+    @property
+    def param_name_to_list_dim_and_degree(self):
+        return {
+            GevParams.LOC: [(self.coordinates.idx_temporal_coordinates, 1)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
+        }
+
+
+class NonStationaryAltitudinalOnlyScaleLocationQuadratic(AltitudinalOnlyScale):
+
+    @property
+    def param_name_to_list_dim_and_degree(self):
+        return {
+            GevParams.LOC: [(self.coordinates.idx_temporal_coordinates, 2)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
+        }
+
+class NonStationaryAltitudinalOnlyScaleLocationCubic(AltitudinalOnlyScale):
+
+    @property
+    def param_name_to_list_dim_and_degree(self):
+        return {
+            GevParams.LOC: [(self.coordinates.idx_temporal_coordinates, 3)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
+        }
+
+class NonStationaryAltitudinalOnlyScaleLocationOrder4(AltitudinalOnlyScale):
+
+    @property
+    def param_name_to_list_dim_and_degree(self):
+        return {
+            GevParams.LOC: [(self.coordinates.idx_temporal_coordinates, 4)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
+        }
+
+# Add cross terms
+
+
+
+class NonStationaryOnlyScaleCrossTermForLocation(AbstractAddCrossTermForLocation, StationaryAltitudinalOnlyScale):
+    pass
+
+
+class NonStationaryAltitudinalOnlyScaleLocationLinearCrossTermForLocation(AbstractAddCrossTermForLocation,
+                                                                 NonStationaryAltitudinalOnlyScaleLocationLinear):
+    pass
+
+
+class NonStationaryAltitudinalOnlyScaleLocationQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation,
+                                                                    NonStationaryAltitudinalOnlyScaleLocationQuadratic):
+    pass
+
+class NonStationaryAltitudinalOnlyScaleLocationCubicCrossTermForLocation(AbstractAddCrossTermForLocation,
+                                                                NonStationaryAltitudinalOnlyScaleLocationCubic,
+                                                                ):
+    pass
+
+class NonStationaryAltitudinalOnlyScaleLocationOrder4CrossTermForLocation(AbstractAddCrossTermForLocation,
+                                                                NonStationaryAltitudinalOnlyScaleLocationOrder4,
+                                                                ):
+    pass
+

extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py

@@ -20,7 +20,7 @@ class AbstractGumbelAltitudinalModel(AbstractAltitudinalModel):
 
     def __init__(self, coordinates: AbstractCoordinates, params_user=None, starting_point=None,
                  fit_method=MarginFitMethod.extremes_fevd_mle, nb_iterations_for_bayesian_fit=5000,
-                 params_initial_fit_bayesian=None, params_class=GevParams, max_degree=3):
+                 params_initial_fit_bayesian=None, params_class=GevParams, max_degree=4):
         super().__init__(coordinates, params_user, starting_point, fit_method, nb_iterations_for_bayesian_fit,
                          params_initial_fit_bayesian, "Gumbel", params_class, max_degree)
 

extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py

@@ -8,7 +8,7 @@ class AbstractSpatioTemporalPolynomialModel(PolynomialMarginModel):
 
     def __init__(self, coordinates: AbstractCoordinates, params_user=None, starting_point=None,
                  fit_method=MarginFitMethod.extremes_fevd_mle, nb_iterations_for_bayesian_fit=5000,
-                 params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=3):
+                 params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=4):
         super().__init__(coordinates, params_user, starting_point, fit_method, nb_iterations_for_bayesian_fit,
                          params_initial_fit_bayesian, type_for_MLE, params_class, max_degree)
         self.drop_duplicates = False

extreme_fit/model/margin_model/polynomial_margin_model/utils.py

@@ -6,7 +6,8 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_mode
     NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation, \
     NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, NonStationaryAltitudinalScaleLinear, \
     NonStationaryAltitudinalScaleLinearCrossTermForLocation, NonStationaryAltitudinalLocationCubicCrossTermForLocation, \
-    NonStationaryAltitudinalLocationCubic
+    NonStationaryAltitudinalLocationCubic, NonStationaryAltitudinalLocationOrder4CrossTermForLocation, \
+    NonStationaryAltitudinalLocationOrder4
 from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_cross_term_in_scale import \
     NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForScale, \
     NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForScale, \
@@ -14,6 +15,14 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_mode
     NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForScale, \
     NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForScale, \
     NonStationaryAltitudinalScaleQuadraticCrossTermForScale
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_only_altitude_and_scale import \
+    StationaryAltitudinalOnlyScale, NonStationaryAltitudinalOnlyScaleLocationLinear, \
+    NonStationaryAltitudinalOnlyScaleLocationQuadratic, NonStationaryAltitudinalOnlyScaleLocationCubic, \
+    NonStationaryAltitudinalOnlyScaleLocationOrder4, NonStationaryOnlyScaleCrossTermForLocation, \
+    NonStationaryAltitudinalOnlyScaleLocationLinearCrossTermForLocation, \
+    NonStationaryAltitudinalOnlyScaleLocationQuadraticCrossTermForLocation, \
+    NonStationaryAltitudinalOnlyScaleLocationCubicCrossTermForLocation, \
+    NonStationaryAltitudinalOnlyScaleLocationOrder4CrossTermForLocation
 from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_with_scale import \
     NonStationaryAltitudinalScaleQuadraticCrossTermForLocation, NonStationaryAltitudinalScaleQuadratic, \
     NonStationaryAltitudinalLocationLinearScaleQuadratic, NonStationaryAltitudinalLocationQuadraticScaleQuadratic, \
@@ -41,14 +50,30 @@ ALTITUDINAL_GEV_MODELS_LOCATION = [
     NonStationaryAltitudinalLocationLinear,
     NonStationaryAltitudinalLocationQuadratic,
     NonStationaryAltitudinalLocationCubic,
+    NonStationaryAltitudinalLocationOrder4,
     # First order cross term
     NonStationaryCrossTermForLocation,
     NonStationaryAltitudinalLocationLinearCrossTermForLocation,
     NonStationaryAltitudinalLocationQuadraticCrossTermForLocation,
-    NonStationaryAltitudinalLocationCubicCrossTermForLocation
+    NonStationaryAltitudinalLocationCubicCrossTermForLocation,
+    NonStationaryAltitudinalLocationOrder4CrossTermForLocation,
 
 ]
 
+ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES = [
+StationaryAltitudinalOnlyScale,
+NonStationaryAltitudinalOnlyScaleLocationLinear,
+NonStationaryAltitudinalOnlyScaleLocationQuadratic,
+NonStationaryAltitudinalOnlyScaleLocationCubic,
+NonStationaryAltitudinalOnlyScaleLocationOrder4,
+# Cross terms
+NonStationaryOnlyScaleCrossTermForLocation,
+NonStationaryAltitudinalOnlyScaleLocationLinearCrossTermForLocation,
+NonStationaryAltitudinalOnlyScaleLocationQuadraticCrossTermForLocation,
+NonStationaryAltitudinalOnlyScaleLocationCubicCrossTermForLocation,
+NonStationaryAltitudinalOnlyScaleLocationOrder4CrossTermForLocation,
+]
+
 
 ALTITUDINAL_GEV_MODELS = [
                              StationaryAltitudinal,

projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py

@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \
-    ALTITUDINAL_GEV_MODELS_LOCATION
+    ALTITUDINAL_GEV_MODELS_LOCATION, ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic, \
     plot_individual_aic
 
@@ -34,6 +34,7 @@ def plot_moments(studies, massif_names=None):
 
 
 def plot_altitudinal_fit(studies, massif_names=None):
+    # model_classes = ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES
     model_classes = ALTITUDINAL_GEV_MODELS_LOCATION
     # model_classes = ALTITUDINAL_GEV_MODELS_QUADRATIC
     visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
@@ -41,7 +42,7 @@ def plot_altitudinal_fit(studies, massif_names=None):
                                                                   massif_names=massif_names,
                                                                   show=False)
     # Plot the results for the model that minimizes the individual aic
-    plot_individual_aic(visualizer)
+    # plot_individual_aic(visualizer)
     # Plot the results for the model that minimizes the total aic
     plot_total_aic(model_classes, visualizer)
 
@@ -51,15 +52,14 @@ def main():
     # todo: l ecart  pour les saisons de l automne ou de sprint
     #  vient probablement de certains zéros pas pris en compte pour le fit,
     # mais pris en compte pour le calcul de mon aic
-    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200][3:]
-    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
-    print(altitudes)
-    altitudes = [900, 1200, 1500][:]
+    # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
+    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200, 4500, 4800][:]
+    # altitudes = [1500, 1800][:]
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
     study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
                      SafranPrecipitation7Days][:]
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranPrecipitation1Day
-                     , SafranPrecipitation3Days][:2]
+                     , SafranPrecipitation3Days][:1]
     # seasons = [Season.automn, Season.winter, Season.spring][::-1]
     seasons = [Season.winter]
     # seasons = [Season.winter_extended]

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py

@@ -5,7 +5,7 @@ import numpy as np
 
 from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
-    SCM_STUDY_CLASS_TO_ABBREVIATION
+    SCM_STUDY_CLASS_TO_ABBREVIATION, ALL_ALTITUDES_WITHOUT_NAN
 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
@@ -143,25 +143,28 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         pass
 
     def plot_year_for_the_peak(self):
-        t_list = 1800 + np.arange(400)
+        t_list = 1959 + np.arange(141)
         for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
             ax = plt.gca()
             # One plot for each altitude
-            print(one_fold_fit.dataset.df_dataset.head())
-            altitudes = self.studies.altitudes
-            print(altitudes)
+            altitudes = np.arange(500, 3000, 500)
             for altitude in altitudes:
-
-                mean_list = []
-                for t in t_list:
-                    coordinate = np.array([altitude, t])
-                    gev_params = one_fold_fit.best_function_from_fit.get_params(coordinate, is_transformed=False)
-                    mean_list.append(gev_params.mean)
-                label = '{} m'.format(altitude)
-                ax.plot(t_list, mean_list, label=label)
+                i = 0
+                while self.studies.altitudes[i] < altitude:
+                    i += 1
+                nearest_altitude = self.studies.altitudes[i]
+                nearest_study = self.studies.altitude_to_study[nearest_altitude]
+                if massif_name in nearest_study.study_massif_names:
+                    mean_list = []
+                    for t in t_list:
+                        coordinate = np.array([altitude, t])
+                        gev_params = one_fold_fit.best_function_from_fit.get_params(coordinate, is_transformed=False)
+                        mean_list.append(gev_params.mean)
+                    label = '{} m'.format(altitude)
+                    ax.plot(t_list, mean_list, label=label)
             ax.legend()
             ax.set_xlabel('Year')
-            ax.set_xlabel('Mean annual maxima of ')
+            ax.set_ylabel('Mean annual maxima of {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)]))
             plot_name = '{}/{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, 'Peak year', massif_name.replace('_', ''))
             self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
             plt.close()

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py

@@ -8,6 +8,8 @@ from extreme_fit.distribution.gumbel.gumbel_gof import goodness_of_fit_anderson
 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.gev_altitudinal_models_only_altitude_and_scale import \
+    AltitudinalOnlyScale, StationaryAltitudinalOnlyScale
 from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_models import \
     StationaryGumbelAltitudinal, AbstractGumbelAltitudinalModel
 from extreme_fit.model.margin_model.utils import MarginFitMethod
@@ -141,6 +143,8 @@ class OneFoldFit(object):
     def stationary_estimator(self):
         if isinstance(self.best_estimator.margin_model, AbstractGumbelAltitudinalModel):
             return self.model_class_to_estimator_with_finite_aic[StationaryGumbelAltitudinal]
+        elif isinstance(self.best_estimator.margin_model, AltitudinalOnlyScale):
+            return self.model_class_to_estimator_with_finite_aic[StationaryAltitudinalOnlyScale]
         else:
             return self.model_class_to_estimator_with_finite_aic[StationaryAltitudinal]
 

projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py

@@ -40,6 +40,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
     legend_size = 15
     # Plot histogram
     ylim = (0, 110)
+    ticks = []
     for j, ci_method in enumerate(visualizer.uncertainty_methods):
         if len(visualizer.uncertainty_methods) == 2:
             offset = -50 if j == 0 else 50
@@ -52,7 +53,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
         plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width, legend_size)
 
         # Plot percentages of return level excess on the right axis
-        ylim = plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax.twinx(), fontsize_label, legend_size)
+        ylim, ticks = plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax.twinx(), fontsize_label, legend_size)
 
 
     ax.set_xticks(altitudes)
@@ -77,8 +78,8 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
     ax_twiny.set_xticklabels(nb_massif_names)
     ax_twiny.set_xlabel('Number of massifs at each altitude (for the percentage and the mean)', fontsize=fontsize_label)
 
-    nb_ticks =  1 + (ylim[1] - ylim[0]) // 20
-    ax.set_yticks([100 - 20 * i for i in range(nb_ticks)][::-1])
+
+    ax.set_yticks(ticks)
     visualizer.plot_name = 'Percentages of exceedance with {}'.format(
         get_display_name_from_object_type(model_subset_for_uncertainty))
     visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure, tight_layout=True)
@@ -131,12 +132,17 @@ def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_meth
                     label=label_confidence_interval)
 
     ax.tick_params(labelsize=fontsize_label)
-    ax.legend(loc='upper right', prop={'size': legend_size})
+    ax.legend(loc='lower right', prop={'size': legend_size})
     ax.set_ylabel(full_label_name, fontsize=fontsize_label)
 
     ylim = (-85, 110)
     ax.set_ylim(ylim)
-    return ylim
+
+    nb_ticks =  1 + (ylim[1] - ylim[0]) // 20
+    ticks = [100 - 20 * i for i in range(nb_ticks)][::-1]
+    ax.set_yticks(ticks)
+
+    return ylim, ticks
 
 
 def get_label_confidence_interval(label_name):