[exceeding] add plot year for the peak. add many altitudinal models.

extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py

@@ -36,7 +36,9 @@ 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 += 'x'
+            name += 'l'
+        if isinstance(self, AbstractAddCrossTermForScale):
+            name += 's'
         return name
 
 
@@ -112,6 +114,19 @@ class AbstractAddCrossTermForLocation(AbstractAltitudinalModel):
         d[GevParams.LOC].insert(1, ((self.coordinates.idx_x_coordinates, self.coordinates.idx_temporal_coordinates), 1))
         return d
 
+class AbstractAddCrossTermForScale(AbstractAltitudinalModel):
+
+    @property
+    def param_name_to_list_dim_and_degree_for_margin_function(self):
+        d = self.param_name_to_list_dim_and_degree
+        # The two insert below enable to check that the insert_index should indeed be 1
+        assert 1 <= len(d[GevParams.SCALE]) <= 2
+        assert self.coordinates.idx_x_coordinates == d[GevParams.SCALE][0][0]
+        insert_index = 1
+        d[GevParams.SCALE].insert(insert_index, ((self.coordinates.idx_x_coordinates, self.coordinates.idx_temporal_coordinates), 1))
+        return d
+
+
 
 class NonStationaryCrossTermForLocation(AbstractAddCrossTermForLocation, StationaryAltitudinal):
     pass

extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_cross_term_in_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, \
+    NonStationaryAltitudinalScaleLinear, NonStationaryAltitudinalLocationLinearScaleLinear, \
+    NonStationaryAltitudinalLocationQuadraticScaleLinear, AbstractAddCrossTermForScale
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_with_scale import \
+    NonStationaryAltitudinalLocationQuadraticScaleQuadratic, NonStationaryAltitudinalScaleQuadratic, \
+    NonStationaryAltitudinalLocationLinearScaleQuadratic
+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 NonStationaryAltitudinalScaleLinearCrossTermForScale(AbstractAddCrossTermForScale,
+                                                              NonStationaryAltitudinalScaleLinear):
+    pass
+
+
+class NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForScale(AbstractAddCrossTermForScale,
+                                                                            NonStationaryAltitudinalLocationLinearScaleLinear):
+    pass
+
+
+class NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForScale(AbstractAddCrossTermForScale,
+                                                                               NonStationaryAltitudinalLocationQuadraticScaleLinear):
+    pass
+
+class NonStationaryAltitudinalScaleQuadraticCrossTermForScale(AbstractAddCrossTermForScale, NonStationaryAltitudinalScaleQuadratic):
+    pass
+
+
+class NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForScale(AbstractAddCrossTermForScale,
+                                                                 NonStationaryAltitudinalLocationLinearScaleQuadratic):
+    pass
+
+
+class NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForScale(AbstractAddCrossTermForScale,
+                                                                    NonStationaryAltitudinalLocationQuadraticScaleQuadratic):
+    pass
+
+
+
+
+

extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_with_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 NonStationaryAltitudinalScaleQuadratic(AbstractAltitudinalModel):
+
+    @property
+    def param_name_to_list_dim_and_degree(self):
+        return {
+            GevParams.LOC: [(self.coordinates.idx_x_coordinates, 1)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 2)],
+        }
+
+class NonStationaryAltitudinalLocationLinearScaleQuadratic(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, 1)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 2)],
+        }
+
+class NonStationaryAltitudinalLocationQuadraticScaleQuadratic(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, 2)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 2)],
+        }
+
+
+# Add cross terms
+
+
+
+class NonStationaryAltitudinalScaleQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation, NonStationaryAltitudinalScaleQuadratic):
+    pass
+
+
+class NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation,
+                                                                 NonStationaryAltitudinalLocationLinearScaleQuadratic):
+    pass
+
+
+class NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation,
+                                                                    NonStationaryAltitudinalLocationQuadraticScaleQuadratic):
+    pass
+

extreme_fit/model/margin_model/polynomial_margin_model/utils.py

@@ -6,6 +6,18 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_mode
     NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation, \
     NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, NonStationaryAltitudinalScaleLinear, \
     NonStationaryAltitudinalScaleLinearCrossTermForLocation
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_cross_term_in_scale import \
+    NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForScale, \
+    NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForScale, \
+    NonStationaryAltitudinalScaleLinearCrossTermForScale, \
+    NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForScale, \
+    NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForScale, \
+    NonStationaryAltitudinalScaleQuadraticCrossTermForScale
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_with_scale import \
+    NonStationaryAltitudinalScaleQuadraticCrossTermForLocation, NonStationaryAltitudinalScaleQuadratic, \
+    NonStationaryAltitudinalLocationLinearScaleQuadratic, NonStationaryAltitudinalLocationQuadraticScaleQuadratic, \
+    NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForLocation, \
+    NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForLocation
 from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_models import \
     StationaryGumbelAltitudinal, NonStationaryGumbelAltitudinalScaleLinear, \
     NonStationaryGumbelAltitudinalLocationLinear, NonStationaryGumbelAltitudinalLocationQuadratic, \
@@ -22,6 +34,7 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
     NonStationaryLocationSpatioTemporalLinearityModelAssertError2, \
     NonStationaryLocationSpatioTemporalLinearityModelAssertError3, NonStationaryLocationSpatioTemporalLinearityModel6
 
+
 ALTITUDINAL_GEV_MODELS = [
                              StationaryAltitudinal,
                              NonStationaryAltitudinalScaleLinear,
@@ -36,8 +49,26 @@ ALTITUDINAL_GEV_MODELS = [
                              NonStationaryAltitudinalLocationQuadraticCrossTermForLocation,
                              NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation,
                              NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation,
+
+                             # Quadratic in the scale
+                             NonStationaryAltitudinalScaleQuadratic,
+                             NonStationaryAltitudinalLocationLinearScaleQuadratic,
+                             NonStationaryAltitudinalLocationQuadraticScaleQuadratic,
+                             NonStationaryAltitudinalScaleQuadraticCrossTermForLocation,
+                             NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForLocation,
+                             NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForLocation,
+
+                             # Cross term for the scale
+                             NonStationaryAltitudinalScaleLinearCrossTermForScale,
+                             NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForScale,
+                             NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForScale,
+                             NonStationaryAltitudinalScaleQuadraticCrossTermForScale,
+                             NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForScale,
+                             NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForScale,
+
                          ][:]
 
+
 ALTITUDINAL_GUMBEL_MODELS = [
                                 StationaryGumbelAltitudinal,
                                 NonStationaryGumbelAltitudinalScaleLinear,

projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py

@@ -2,6 +2,7 @@ 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 projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic, \
     plot_individual_aic
 
@@ -16,8 +17,6 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranS
     SafranSnowfall5Days, SafranSnowfall7Days, SafranPrecipitation1Day, SafranPrecipitation3Days, \
     SafranPrecipitation5Days, SafranPrecipitation7Days
 from extreme_data.meteo_france_data.scm_models_data.utils import Season
-from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \
-    ALTITUDINAL_GUMBEL_MODELS
 from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
     AltitudesStudiesVisualizerForNonStationaryModels
@@ -34,24 +33,23 @@ def plot_moments(studies, massif_names=None):
 
 
 def plot_altitudinal_fit(studies, massif_names=None):
-    model_classes = ALTITUDINAL_GEV_MODELS + ALTITUDINAL_GUMBEL_MODELS
+    model_classes = ALTITUDINAL_GEV_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
-    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)
 
 
 def main():
-    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][4:6]
     # 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 = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900][:]
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
     study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
                      SafranPrecipitation7Days][:]
@@ -59,6 +57,7 @@ def main():
                      SafranSnowfall3Days, SafranPrecipitation3Days][:1]
     # seasons = [Season.automn, Season.winter, Season.spring][::-1]
     seasons = [Season.winter]
+    # seasons = [Season.winter_extended]
 
     massif_names = None
     # massif_names = ['Aravis']

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py

@@ -26,7 +26,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                  show=False,
                  massif_names=None,
                  fit_method=MarginFitMethod.extremes_fevd_mle,
-                 display_only_model_that_pass_anderson_test=True):
+                 display_only_model_that_pass_anderson_test=False):
         super().__init__(studies.study, show=show, save_to_file=not show)
         self.studies = studies
         self.non_stationary_models = model_classes
@@ -127,8 +127,9 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         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),
+                                                                         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)
 
@@ -138,5 +139,30 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                                                                     self.study.variable_unit),
                                 add_x_label=False, graduation=0.1, massif_name_to_text=self.massif_name_to_name)
 
-    def plot_altitude_change_of_sign(self):
+    def plot_altitude_for_the_peak(self):
         pass
+
+    def plot_year_for_the_peak(self):
+        t_list = 1900 + np.arange(200)
+        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)
+            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)
+            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)
+            plt.close()
+

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py

@@ -163,7 +163,9 @@ class OneFoldFit(object):
         quantiles = self.compute_empirical_quantiles(estimator=self.best_estimator)
         goodness_of_fit_anderson_test = goodness_of_fit_anderson(quantiles, self.SIGNIFICANCE_LEVEL)
         if not goodness_of_fit_anderson_test:
-            print('{} with {} does not pass the anderson test'.format(self.massif_name, self.folder_for_plots))
+            print('{} with {} does not pass the anderson test for model {}'.format(self.massif_name,
+                                                                                   self.folder_for_plots,
+                                                                                   type(self.best_margin_model)))
         return goodness_of_fit_anderson_test
 
     def compute_empirical_quantiles(self, estimator):

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py

@@ -10,11 +10,16 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fi
 from projects.exceeding_snow_loads.utils import dpi_paper1_figure
 
 
-def plot_individual_aic(visualizer):
-    OneFoldFit.best_estimator_minimizes_total_aic = False
+def plots(visualizer):
     visualizer.plot_moments()
-    # visualizer.plot_best_coef_maps()
+    visualizer.plot_best_coef_maps()
     visualizer.plot_shape_map()
+    visualizer.plot_year_for_the_peak()
+
+
+def plot_individual_aic(visualizer):
+    OneFoldFit.best_estimator_minimizes_total_aic = False
+    plots(visualizer)
 
 
 def plot_total_aic(model_classes, visualizer):
@@ -45,9 +50,7 @@ def plot_total_aic(model_classes, visualizer):
     # Plot the ranking of the model based on their total aic
     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()
+    plots(visualizer)
 
 
 def plot_total_aic_repartition(visualizer, labels, scores):