[contrasting] add cubic models

extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py

@@ -36,7 +36,7 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
         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 += 'l'
+            name += 'x'
         if isinstance(self, AbstractAddCrossTermForScale):
             name += 's'
         return name
@@ -81,6 +81,15 @@ class NonStationaryAltitudinalLocationQuadratic(AbstractAltitudinalModel):
             GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
         }
 
+class NonStationaryAltitudinalLocationCubic(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, 3)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
+        }
+
 
 class NonStationaryAltitudinalLocationLinearScaleLinear(AbstractAltitudinalModel):
 
@@ -141,6 +150,10 @@ class NonStationaryAltitudinalLocationQuadraticCrossTermForLocation(AbstractAddC
                                                                     NonStationaryAltitudinalLocationQuadratic):
     pass
 
+class NonStationaryAltitudinalLocationCubicCrossTermForLocation(AbstractAddCrossTermForLocation,
+                                                                NonStationaryAltitudinalLocationCubic,
+                                                                ):
+    pass
 
 class NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation,
                                                                             NonStationaryAltitudinalLocationLinearScaleLinear):
@@ -155,3 +168,4 @@ class NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation(A
 class NonStationaryAltitudinalScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation,
                                                               NonStationaryAltitudinalScaleLinear):
     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=2):
+                 params_initial_fit_bayesian=None, params_class=GevParams, max_degree=3):
         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/polynomial_margin_model.py

@@ -40,6 +40,10 @@ class PolynomialMarginModel(AbstractTemporalLinearMarginModel):
                 assert dims.index(self.coordinates.idx_x_coordinates) == 0
             if self.coordinates.has_temporal_coordinates and self.coordinates.idx_temporal_coordinates in dims:
                 assert dims.index(self.coordinates.idx_temporal_coordinates) == len(dims) - 1
+        # Assert that the degree are inferior to the max degree
+        for list_dim_and_degree in param_name_to_list_dim_and_degree.values():
+            for _, max_degree in list_dim_and_degree:
+                assert max_degree <= self.max_degree, 'Max degree (={}) specified is too high'.format(max_degree)
         # Load param_name_to_polynomial_all_coef
         param_name_to_polynomial_all_coef = self.param_name_to_polynomial_all_coef(
             param_name_to_list_dim_and_degree=param_name_to_list_dim_and_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=2):
+                 params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=3):
         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

@@ -5,7 +5,8 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_mode
     NonStationaryAltitudinalLocationQuadraticCrossTermForLocation, \
     NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation, \
     NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, NonStationaryAltitudinalScaleLinear, \
-    NonStationaryAltitudinalScaleLinearCrossTermForLocation
+    NonStationaryAltitudinalScaleLinearCrossTermForLocation, NonStationaryAltitudinalLocationCubicCrossTermForLocation, \
+    NonStationaryAltitudinalLocationCubic
 from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_cross_term_in_scale import \
     NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForScale, \
     NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForScale, \
@@ -35,6 +36,20 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
     NonStationaryLocationSpatioTemporalLinearityModelAssertError3, NonStationaryLocationSpatioTemporalLinearityModel6
 
 
+ALTITUDINAL_GEV_MODELS_LOCATION = [
+    StationaryAltitudinal,
+    NonStationaryAltitudinalLocationLinear,
+    NonStationaryAltitudinalLocationQuadratic,
+    NonStationaryAltitudinalLocationCubic,
+    # First order cross term
+    NonStationaryCrossTermForLocation,
+    NonStationaryAltitudinalLocationLinearCrossTermForLocation,
+    NonStationaryAltitudinalLocationQuadraticCrossTermForLocation,
+    NonStationaryAltitudinalLocationCubicCrossTermForLocation
+
+]
+
+
 ALTITUDINAL_GEV_MODELS = [
                              StationaryAltitudinal,
                              NonStationaryAltitudinalScaleLinear,

extreme_fit/model/result_from_model_fit/utils.py

@@ -33,9 +33,6 @@ def get_margin_coef_ordered_dict(param_name_to_dims, mle_values, type_for_mle="G
                 coef_dict[coef_name] = mle_values[i]
                 i += 1
             else:
-                # We found (thanks to the test) that time was the first parameter when len(param_name_to_dims) == 1
-                # otherwise time is the second parameter in the order of the mle parameters
-                # inverted_dims = dims[::-1] if len(param_name_to_dims) == 1 else dims
                 for dim, max_degree in dims[:]:
                     coefficient_name = LinearCoef.coefficient_name(dim, dim_to_coordinate_name)
                     coef_template = LinearCoef.coef_template_str(param_name, coefficient_name)

projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py

@@ -2,7 +2,8 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
-from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \
+    ALTITUDINAL_GEV_MODELS_LOCATION
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic, \
     plot_individual_aic
 
@@ -33,7 +34,7 @@ def plot_moments(studies, massif_names=None):
 
 
 def plot_altitudinal_fit(studies, massif_names=None):
-    model_classes = ALTITUDINAL_GEV_MODELS
+    model_classes = ALTITUDINAL_GEV_MODELS_LOCATION
     # model_classes = ALTITUDINAL_GEV_MODELS_QUADRATIC
     visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
                                                                   model_classes=model_classes,
@@ -42,7 +43,7 @@ def plot_altitudinal_fit(studies, massif_names=None):
     # Plot the results for the model that minimizes the individual aic
     plot_individual_aic(visualizer)
     # Plot the results for the model that minimizes the total aic
-    # plot_total_aic(model_classes, visualizer)
+    plot_total_aic(model_classes, visualizer)
 
 
 def main():
@@ -50,12 +51,15 @@ 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 = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900][:]
+    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][:]
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
     study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
                      SafranPrecipitation7Days][:]
-    study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day,
-                     SafranSnowfall3Days, SafranPrecipitation3Days][:1]
+    study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranPrecipitation1Day
+                     , SafranPrecipitation3Days][:2]
     # 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

@@ -32,7 +32,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         self.non_stationary_models = model_classes
         self.fit_method = fit_method
         self.display_only_model_that_pass_anderson_test = display_only_model_that_pass_anderson_test
-        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.all_massif_names()
         self.massif_name_to_massif_id = {m: i for i, m in enumerate(self.massif_names)}
         # Load one fold fit
         self._massif_name_to_one_fold_fit = {}
@@ -143,7 +143,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         pass
 
     def plot_year_for_the_peak(self):
-        t_list = 1900 + np.arange(200)
+        t_list = 1800 + np.arange(400)
         for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
             ax = plt.gca()
             # One plot for each altitude
@@ -162,7 +162,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
             ax.legend()
             ax.set_xlabel('Year')
             ax.set_xlabel('Mean annual maxima of ')
-            self.plot_name = 'Peak year for {}'.format(massif_name)
-            self.show_or_save_to_file(no_title=True)
+            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

@@ -116,7 +116,8 @@ class OneFoldFit(object):
 
     @property
     def best_shape(self):
-        # We take any altitude (altitude=1000 for instance) as the shape is constant w.r.t the altitude
+        # We take any altitude (altitude=1000 for instance) and any year
+        # as the shape is constant w.r.t the altitude and the year
         return self.get_gev_params(altitude=1000, year=2019).shape
 
     def best_coef(self, param_name, dim, degree):

projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py

@@ -74,7 +74,7 @@ def intermediate_result(altitudes, massif_names=None,
         # Plots
         # plot_trend_map(altitude_to_visualizer)
         # plot_trend_curves(altitude_to_visualizer={a: v for a, v in altitude_to_visualizer.items() if a >= 900})
-        plot_uncertainty_massifs(altitude_to_visualizer)
+        # plot_uncertainty_massifs(altitude_to_visualizer)
         plot_uncertainty_histogram(altitude_to_visualizer)
         # plot_selection_curves(altitude_to_visualizer)
         # plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)

projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py

-from typing import Dict
+from typing import Dict, Tuple
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -39,6 +39,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
     fontsize_label = 15
     legend_size = 15
     # Plot histogram
+    ylim = (0, 110)
     for j, ci_method in enumerate(visualizer.uncertainty_methods):
         if len(visualizer.uncertainty_methods) == 2:
             offset = -50 if j == 0 else 50
@@ -51,7 +52,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
-        plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax.twinx(), fontsize_label, legend_size)
+        ylim = plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax.twinx(), fontsize_label, legend_size)
 
 
     ax.set_xticks(altitudes)
@@ -63,7 +64,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
     #               'exceeds French standards (\%)', fontsize=fontsize_label)
     ax.set_ylabel('Massifs exceeding French standards (\%)', fontsize=fontsize_label)
     ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
-    ax.set_ylim([0, 110])
+    ax.set_ylim(ylim)
     ax.yaxis.grid()
 
     ax_twiny = ax.twiny()
@@ -76,7 +77,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)
 
-    ax.set_yticks([20 * i for i in range(6)])
+    nb_ticks =  1 + (ylim[1] - ylim[0]) // 20
+    ax.set_yticks([100 - 20 * i for i in range(nb_ticks)][::-1])
     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)
@@ -109,7 +111,7 @@ def plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_metho
     ax.legend(loc='upper left', prop={'size': legend_size})
 
 
-def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax, fontsize_label, legend_size=10):
+def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax, fontsize_label, legend_size=10) -> Tuple[int, int]:
     l = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[2] for v in visualizers]
     lower_bound, mean, upper_bound = list(zip(*l))
     other_mean = [e[1] for e in l]
@@ -121,7 +123,6 @@ def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_meth
                       '50-year return levels and French standards (\%)'
 
     label_name = 'Mean relative difference'
-    print('mean relative difference', mean)
     ax.plot(altitudes, mean, linestyle='--', marker='o', color=color,
             label=label_name)
 
@@ -132,7 +133,10 @@ def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_meth
     ax.tick_params(labelsize=fontsize_label)
     ax.legend(loc='upper right', prop={'size': legend_size})
     ax.set_ylabel(full_label_name, fontsize=fontsize_label)
-    ax.set_ylim([0, 110])
+
+    ylim = (-85, 110)
+    ax.set_ylim(ylim)
+    return ylim
 
 
 def get_label_confidence_interval(label_name):