[contrasting] add NON_STATIONARY_TREND_TEST_PAPER_2 and associated model and trend test.

extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py

@@ -57,7 +57,7 @@ class NonStationaryLocationAndScaleTemporalModel(AbstractTemporalLinearMarginMod
 
     def load_margin_function(self, param_name_to_dims=None):
         return super().load_margin_function({GevParams.LOC: [self.coordinates.idx_temporal_coordinates],
-                                      GevParams.SCALE: [self.coordinates.idx_temporal_coordinates]})
+                                             GevParams.SCALE: [self.coordinates.idx_temporal_coordinates]})
 
     @property
     def mul(self):
@@ -68,6 +68,56 @@ class NonStationaryLocationAndScaleTemporalModel(AbstractTemporalLinearMarginMod
         return 1
 
 
+class NonStationaryLocationAndShapeTemporalModel(AbstractTemporalLinearMarginModel):
+
+    def load_margin_function(self, param_name_to_dims=None):
+        return super().load_margin_function({GevParams.LOC: [self.coordinates.idx_temporal_coordinates],
+                                             GevParams.SHAPE: [self.coordinates.idx_temporal_coordinates]})
+
+    @property
+    def mul(self):
+        return 1
+
+    @property
+    def shl(self):
+        return 1
+
+
+class NonStationaryScaleAndShapeTemporalModel(AbstractTemporalLinearMarginModel):
+
+    def load_margin_function(self, param_name_to_dims=None):
+        return super().load_margin_function({GevParams.SHAPE: [self.coordinates.idx_temporal_coordinates],
+                                             GevParams.SCALE: [self.coordinates.idx_temporal_coordinates]})
+
+    @property
+    def shl(self):
+        return 1
+
+    @property
+    def sigl(self):
+        return 1
+
+
+class NonStationaryLocationAndScaleAndShapeTemporalModel(AbstractTemporalLinearMarginModel):
+
+    def load_margin_function(self, param_name_to_dims=None):
+        return super().load_margin_function({GevParams.LOC: [self.coordinates.idx_temporal_coordinates],
+                                             GevParams.SCALE: [self.coordinates.idx_temporal_coordinates],
+                                             GevParams.SHAPE: [self.coordinates.idx_temporal_coordinates]})
+
+    @property
+    def mul(self):
+        return 1
+
+    @property
+    def sigl(self):
+        return 1
+
+    @property
+    def shl(self):
+        return 1
+
+
 class GumbelTemporalModel(StationaryTemporalModel):
 
     def __init__(self, coordinates: AbstractCoordinates,

extreme_trend/trend_test_one_parameter/gev_trend_test_one_parameter.py

@@ -6,6 +6,7 @@ from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_m
     NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel, NonStationaryShapeTemporalModel, \
     StationaryTemporalModel
 from extreme_fit.distribution.gev.gev_params import GevParams
+from root_utils import classproperty
 
 
 class GevTrendTestOneParameter(AbstractGevTrendTest):
@@ -15,6 +16,28 @@ class GevTrendTestOneParameter(AbstractGevTrendTest):
         return 1
 
 
+class GevVersusGev(GevTrendTestOneParameter):
+
+    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE,
+                 fit_method=MarginFitMethod.extremes_fevd_mle):
+        super().__init__(years, maxima, starting_year,
+                         unconstrained_model_class=StationaryTemporalModel,
+                         constrained_model_class=StationaryTemporalModel,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
+
+    @property
+    def is_significant(self) -> bool:
+        return False
+
+    @classproperty
+    def total_number_of_parameters_for_unconstrained_model(cls) -> int:
+        return 3
+
+    def _slope_strength(self):
+        return 0.0
+
+
 class GevTrendTestOneParameterAgainstStationary(GevTrendTestOneParameter):
 
     def __init__(self, years, maxima, starting_year, unconstrained_model_class, param_name,
@@ -32,10 +55,15 @@ class GevTrendTestOneParameterAgainstStationary(GevTrendTestOneParameter):
     def non_stationary_linear_coef(self):
         return self.get_non_stationary_linear_coef(param_name=self.param_name)
 
+    @classproperty
+    def total_number_of_parameters_for_unconstrained_model(cls) -> int:
+        return 4
+
 
 class GevLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
 
-    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE, constrained_model_class=StationaryTemporalModel,
+    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE,
+                 constrained_model_class=StationaryTemporalModel,
                  fit_method=MarginFitMethod.extremes_fevd_mle):
         super().__init__(years, maxima, starting_year,
                          unconstrained_model_class=NonStationaryLocationTemporalModel,
@@ -58,9 +86,11 @@ class GevLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
     def variance_difference_same_sign_as_slope_strenght(self) -> bool:
         return False
 
+
 class GevScaleTrendTest(GevTrendTestOneParameterAgainstStationary):
 
-    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE, constrained_model_class=StationaryTemporalModel,
+    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE,
+                 constrained_model_class=StationaryTemporalModel,
                  fit_method=MarginFitMethod.extremes_fevd_mle):
         super().__init__(years, maxima, starting_year,
                          unconstrained_model_class=NonStationaryScaleTemporalModel,
@@ -88,11 +118,10 @@ class GevScaleTrendTest(GevTrendTestOneParameterAgainstStationary):
 
 class GevShapeTrendTest(GevTrendTestOneParameterAgainstStationary):
 
-    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE, fit_method=MarginFitMethod.extremes_fevd_mle):
+    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE,
+                 fit_method=MarginFitMethod.extremes_fevd_mle):
         super().__init__(years, maxima, starting_year,
                          unconstrained_model_class=NonStationaryShapeTemporalModel,
                          param_name=GevParams.SHAPE,
                          quantile_level=quantile_level,
                          fit_method=fit_method)
-
-

extreme_trend/trend_test_three_parameters/gev_trend_test_three_parameters.py

@@ -3,7 +3,8 @@ from extreme_fit.model.margin_model.utils import \
 from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE
 from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest
 from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
-    NonStationaryLocationAndScaleTemporalModel, GumbelTemporalModel
+    NonStationaryLocationAndScaleTemporalModel, GumbelTemporalModel, StationaryTemporalModel, \
+    NonStationaryLocationAndScaleAndShapeTemporalModel
 from extreme_fit.distribution.gev.gev_params import GevParams
 from root_utils import classproperty
 
@@ -15,9 +16,25 @@ class GevTrendTestThreeParameters(AbstractGevTrendTest):
         return 3
 
 
+class GevLocationAndScaleAndShapeTrendTest(GevTrendTestThreeParameters):
+
+    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE,
+                 fit_method=MarginFitMethod.extremes_fevd_mle):
+        super().__init__(years, maxima, starting_year,
+                         unconstrained_model_class=StationaryTemporalModel,
+                         constrained_model_class=NonStationaryLocationAndScaleAndShapeTemporalModel,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
+
+    @classproperty
+    def total_number_of_parameters_for_unconstrained_model(cls) -> int:
+        return 6
+
+
 class GevLocationAndScaleTrendTestAgainstGumbel(GevTrendTestThreeParameters):
 
-    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE, fit_method=MarginFitMethod.extremes_fevd_mle):
+    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE,
+                 fit_method=MarginFitMethod.extremes_fevd_mle):
         super().__init__(years, maxima, starting_year,
                          unconstrained_model_class=NonStationaryLocationAndScaleTemporalModel,
                          constrained_model_class=GumbelTemporalModel,
@@ -47,4 +64,4 @@ class GevLocationAndScaleTrendTestAgainstGumbel(GevTrendTestThreeParameters):
 
     @classproperty
     def total_number_of_parameters_for_unconstrained_model(cls) -> int:
-        return 5
\ No newline at end of file
+        return 5

extreme_trend/trend_test_two_parameters/gev_trend_test_two_parameters.py

@@ -5,7 +5,8 @@ from extreme_trend.trend_test_one_parameter.gev_trend_test_one_parameter import
 from extreme_fit.model.margin_model.utils import \
     MarginFitMethod
 from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
-    NonStationaryLocationAndScaleTemporalModel, StationaryTemporalModel, GumbelTemporalModel
+    NonStationaryLocationAndScaleTemporalModel, StationaryTemporalModel, GumbelTemporalModel, \
+    NonStationaryScaleAndShapeTemporalModel, NonStationaryLocationAndShapeTemporalModel
 from extreme_fit.distribution.gev.gev_params import GevParams
 from root_utils import classproperty
 
@@ -16,8 +17,37 @@ class GevTrendTestTwoParameters(AbstractGevTrendTest):
     def degree_freedom_chi2(self) -> int:
         return 2
 
+class GevTrendTestTwoParametersAgainstGev(AbstractGevTrendTest):
 
-class GevLocationAndScaleTrendTest(GevTrendTestTwoParameters):
+    @classproperty
+    def total_number_of_parameters_for_unconstrained_model(cls) -> int:
+        return 5
+
+
+class GevLocationAndShapeTrendTest(GevTrendTestTwoParametersAgainstGev):
+    
+    def __init__(self, years, maxima, starting_year, constrained_model_class=StationaryTemporalModel,
+                 quantile_level=EUROCODE_QUANTILE, fit_method=MarginFitMethod.extremes_fevd_mle):
+        super().__init__(years, maxima, starting_year,
+                         unconstrained_model_class=NonStationaryLocationAndShapeTemporalModel,
+                         constrained_model_class=constrained_model_class,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
+
+
+class GevScaleAndShapeTrendTest(GevTrendTestTwoParametersAgainstGev):
+
+
+    def __init__(self, years, maxima, starting_year, constrained_model_class=StationaryTemporalModel,
+                 quantile_level=EUROCODE_QUANTILE, fit_method=MarginFitMethod.extremes_fevd_mle):
+        super().__init__(years, maxima, starting_year,
+                         unconstrained_model_class=NonStationaryScaleAndShapeTemporalModel,
+                         constrained_model_class=constrained_model_class,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
+
+
+class GevLocationAndScaleTrendTest(GevTrendTestTwoParametersAgainstGev):
 
     def __init__(self, years, maxima, starting_year, constrained_model_class=StationaryTemporalModel,
                  quantile_level=EUROCODE_QUANTILE, fit_method=MarginFitMethod.extremes_fevd_mle):

extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py

@@ -17,7 +17,7 @@ from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study impo
 from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \
     StudyVisualizer
-from projects.exceeding_snow_loads.utils import NON_STATIONARY_TREND_TEST_PAPER
+from projects.exceeding_snow_loads.utils import NON_STATIONARY_TREND_TEST_PAPER_1
 from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest
 from extreme_trend.trend_test_one_parameter.gumbel_trend_test_one_parameter import \
     GumbelLocationTrendTest, GevStationaryVersusGumbel, GumbelScaleTrendTest, GumbelVersusGumbel
@@ -87,7 +87,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
             assert set(self.uncertainty_massif_names).issubset(set(self.study.study_massif_names))
         if self.non_stationary_trend_test_to_marker is None:
             # Assign default argument for the non stationary trends
-            self.non_stationary_trend_test = NON_STATIONARY_TREND_TEST_PAPER
+            self.non_stationary_trend_test = NON_STATIONARY_TREND_TEST_PAPER_1
             self.non_stationary_trend_test_to_marker = {t: t.marker for t in self.non_stationary_trend_test}
         else:
             self.non_stationary_trend_test = list(self.non_stationary_trend_test_to_marker.keys())

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/data.py

@@ -63,7 +63,7 @@ def max_graph_annual_maxima_poster():
             tight_pad = {'h_pad': 0.2}
             snow_abbreviation = 'max ' + snow_abbreviation
             for color, _, massif_name in marker_altitude_massif_name_for_paper1[::-1]:
-                last_plot = color == "magenta" and altitude == second_altitude
+                last_plot = massif_name == 'Ubaye' and altitude == second_altitude
                 label = '{} massif at {}m'.format(massif_name, altitude)
                 study_visualizer.visualize_max_graphs_poster(massif_name, altitude, snow_abbreviation, color, label,
                                                              last_plot, ax, tight_pad=tight_pad,

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py

@@ -72,10 +72,10 @@ def intermediate_result(altitudes, massif_names=None,
             _ = compute_minimized_aic(visualizer)
 
     # Plots
-    validation_plot(altitude_to_visualizer, order_derivative=0)
-    validation_plot(altitude_to_visualizer, order_derivative=1)
+    # validation_plot(altitude_to_visualizer, order_derivative=0)
+    # validation_plot(altitude_to_visualizer, order_derivative=1)
     plot_snowfall_mean(altitude_to_visualizer)
-    plot_snowfall_time_derivative_mean(altitude_to_visualizer)
+    # plot_snowfall_time_derivative_mean(altitude_to_visualizer)
 
 
 def major_result():
@@ -86,7 +86,7 @@ def major_result():
     model_subsets_for_uncertainty = None
     altitudes = paper_altitudes
     altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
-    # altitudes = [900, 1200]
+    altitudes = [900, 1200]
 
     for study_class in study_classes:
         intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,

projects/exceeding_snow_loads/utils.py

@@ -4,12 +4,15 @@ from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusS
     CrocusSnowLoad3Days
 from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
     ALL_ALTITUDES_WITHOUT_NAN
+from extreme_trend.trend_test_one_parameter.gev_trend_test_one_parameter import GevVersusGev, GevScaleTrendTest, \
+    GevLocationTrendTest, GevShapeTrendTest
 from extreme_trend.trend_test_one_parameter.gumbel_trend_test_one_parameter import \
     GumbelVersusGumbel, GumbelLocationTrendTest, GumbelScaleTrendTest, GevStationaryVersusGumbel
 from extreme_trend.trend_test_three_parameters.gev_trend_test_three_parameters import \
-    GevLocationAndScaleTrendTestAgainstGumbel
+    GevLocationAndScaleTrendTestAgainstGumbel, GevLocationAndScaleAndShapeTrendTest
 from extreme_trend.trend_test_two_parameters.gev_trend_test_two_parameters import \
-    GevLocationAgainstGumbel, GevScaleAgainstGumbel
+    GevLocationAgainstGumbel, GevScaleAgainstGumbel, GevLocationAndScaleTrendTest, GevScaleAndShapeTrendTest, \
+    GevLocationAndShapeTrendTest
 from extreme_trend.trend_test_two_parameters.gumbel_test_two_parameters import \
     GumbelLocationAndScaleTrendTest
 
@@ -17,12 +20,20 @@ paper_altitudes = ALL_ALTITUDES_WITHOUT_NAN
 paper_study_classes = [CrocusSnowLoadTotal, CrocusSnowLoadEurocode, CrocusSnowLoad3Days][:2]
 # dpi_paper1_figure = 700
 dpi_paper1_figure = None
-NON_STATIONARY_TREND_TEST_PAPER = [GumbelVersusGumbel,
-                                   GumbelLocationTrendTest, GumbelScaleTrendTest,
-                                   GumbelLocationAndScaleTrendTest,
-                                   GevStationaryVersusGumbel,
-                                   GevLocationAgainstGumbel, GevScaleAgainstGumbel,
-                                   GevLocationAndScaleTrendTestAgainstGumbel]
+NON_STATIONARY_TREND_TEST_PAPER_1 = [GumbelVersusGumbel,
+                                     GumbelLocationTrendTest, GumbelScaleTrendTest,
+                                     GumbelLocationAndScaleTrendTest,
+                                     GevStationaryVersusGumbel,
+                                     GevLocationAgainstGumbel, GevScaleAgainstGumbel,
+                                     GevLocationAndScaleTrendTestAgainstGumbel]
 
+NON_STATIONARY_TREND_TEST_PAPER_2 = [
+    # Gumbel models
+    GumbelVersusGumbel, GumbelLocationTrendTest, GumbelScaleTrendTest, GumbelLocationAndScaleTrendTest,
+    # GEV models with constant shape
+    GevVersusGev, GevLocationTrendTest, GevScaleTrendTest,  GevLocationAndScaleTrendTest,
+    # GEV models with linear shape
+    GevShapeTrendTest, GevLocationAndShapeTrendTest, GevScaleAndShapeTrendTest, GevLocationAndScaleAndShapeTrendTest
+]
 
 

test/test_extreme_trend/test_extreme_trend.py

 import unittest
 
-from projects.exceeding_snow_loads.utils import NON_STATIONARY_TREND_TEST_PAPER
+from projects.exceeding_snow_loads.utils import NON_STATIONARY_TREND_TEST_PAPER_1, NON_STATIONARY_TREND_TEST_PAPER_2
 
 
 class TestTrendAnalysis(unittest.TestCase):
 
-    def test_nb_parameters(self):
-        trend_test_classes = NON_STATIONARY_TREND_TEST_PAPER
+    def test_nb_parameters_paper1(self):
+        trend_test_classes = NON_STATIONARY_TREND_TEST_PAPER_1
         nb_expected = [2, 3, 3, 4, 3, 4, 4, 5]
         for trend_test_class, nb in zip(trend_test_classes, nb_expected):
             self.assertEqual(trend_test_class.total_number_of_parameters_for_unconstrained_model, nb)
 
+    def test_nb_parameters_paper2(self):
+        trend_test_classes = NON_STATIONARY_TREND_TEST_PAPER_2
+        nb_expected = [2, 3, 3, 4,
+                       3, 4, 4, 5,
+                       4, 5, 5, 6]
+        for trend_test_class, nb in zip(trend_test_classes, nb_expected):
+            self.assertEqual(trend_test_class.total_number_of_parameters_for_unconstrained_model, nb)
+
 
 if __name__ == '__main__':
     unittest.main()