From 95efb390e95dce9cbba602f67d8c40faeea6e3ac Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 16 Jun 2020 17:38:47 +0200
Subject: [PATCH] [contrasting] add quadratic model for gumbel/gev and
 location/scale. modify code for the name of the models

---
 .../margin_model/polynomial_margin_model.py   | 25 +++++++++++++--
 extreme_trend/abstract_gev_trend_test.py      | 23 +++++--------
 .../gev_trend_test_three_parameters.py        | 32 +++++++++++++++++++
 .../gev_trend_test_two_parameters.py          | 22 +++++++++++++
 .../gumbel_test_two_parameters.py             | 31 ++++++++++++++++++
 projects/exceeding_snow_loads/utils.py        | 13 +++++---
 ...st_gev_temporal_polynomial_extremes_mle.py | 28 ++++++++++++----
 7 files changed, 146 insertions(+), 28 deletions(-)

diff --git a/extreme_fit/model/margin_model/polynomial_margin_model.py b/extreme_fit/model/margin_model/polynomial_margin_model.py
index 1f327961..c5066d72 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model.py
@@ -29,8 +29,9 @@ class PolynomialMarginModel(AbstractTemporalLinearMarginModel):
         return super().margin_function
 
     def load_margin_function(self, param_name_to_list_dim_and_degree=None):
-        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,
-                                                                                   param_name_and_dim_and_degree_to_default_coef=self.default_params)
+        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,
+            param_name_and_dim_and_degree_to_default_coef=self.default_params)
         return PolynomialMarginFunction(coordinates=self.coordinates,
                                         param_name_to_coef=param_name_to_polynomial_all_coef,
                                         param_name_to_dim_and_max_degree=param_name_to_list_dim_and_degree,
@@ -84,6 +85,24 @@ class NonStationaryQuadraticLocationModel(PolynomialMarginModel):
     def load_margin_function(self, param_name_to_dims=None):
         return super().load_margin_function({GevParams.LOC: [(self.coordinates.idx_temporal_coordinates, 2)]})
 
+    @property
+    def mul(self):
+        return 2
+
+
+class NonStationaryQuadraticScaleModel(PolynomialMarginModel):
+
+    def load_margin_function(self, param_name_to_dims=None):
+        return super().load_margin_function({GevParams.SCALE: [(self.coordinates.idx_temporal_coordinates, 2)]})
+
+    @property
+    def sigl(self):
+        return 2
+
+
+class NonStationaryQuadraticLocationGumbelModel(GumbelTemporalModel, NonStationaryQuadraticLocationModel):
+    pass
+
 
-class NonStationaryLocationGumbelModel(GumbelTemporalModel, NonStationaryQuadraticLocationModel):
+class NonStationaryQuadraticScaleGumbelModel(GumbelTemporalModel, NonStationaryQuadraticScaleModel):
     pass
diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index c50de998..13121a30 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/extreme_trend/abstract_gev_trend_test.py
@@ -17,6 +17,7 @@ from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_m
     StationaryTemporalModel, GumbelTemporalModel
 from extreme_fit.model.margin_model.utils import \
     MarginFitMethod
+from extreme_fit.model.utils import get_null
 from root_utils import classproperty
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.utils import load_temporal_coordinates_and_dataset
@@ -84,23 +85,15 @@ class AbstractGevTrendTest(object):
             family = 'Gum'
         else:
             family = 'Gev'
+            
+        def get_str_number(param_l):
+            nb = 0 if param_l == get_null() else param_l
+            return str(nb)
 
-        if self.unconstrained_estimator.margin_model.mul == 1:
-            family += '1'
-        else:
-            family += '0'
-
-        if self.unconstrained_estimator.margin_model.sigl == 1:
-            family += '1'
-        else:
-            family += '0'
-
+        family += get_str_number(self.unconstrained_estimator.margin_model.mul)
+        family += get_str_number(self.unconstrained_estimator.margin_model.sigl)
         if family.startswith('Gev'):
-            if self.unconstrained_estimator.margin_model.shl == 1:
-                family += '1'
-            else:
-                family += '0'
-
+            family += get_str_number(self.unconstrained_estimator.margin_model.shl)
         return family
 
     @property
diff --git a/extreme_trend/trend_test_three_parameters/gev_trend_test_three_parameters.py b/extreme_trend/trend_test_three_parameters/gev_trend_test_three_parameters.py
index 04d780b0..8af11133 100644
--- a/extreme_trend/trend_test_three_parameters/gev_trend_test_three_parameters.py
+++ b/extreme_trend/trend_test_three_parameters/gev_trend_test_three_parameters.py
@@ -1,3 +1,5 @@
+from extreme_fit.model.margin_model.polynomial_margin_model import NonStationaryQuadraticLocationModel, \
+    NonStationaryQuadraticScaleModel
 from extreme_fit.model.margin_model.utils import \
     MarginFitMethod
 from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE
@@ -65,3 +67,33 @@ class GevLocationAndScaleTrendTestAgainstGumbel(GevTrendTestThreeParameters):
     @classproperty
     def total_number_of_parameters_for_unconstrained_model(cls) -> int:
         return 5
+
+
+class GevLocationQuadraticTrendTestAgainstGumbel(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=NonStationaryQuadraticLocationModel,
+                         constrained_model_class=GumbelTemporalModel,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
+
+    @classproperty
+    def total_number_of_parameters_for_unconstrained_model(cls) -> int:
+        return 5
+
+
+class GevScaleQuadraticTrendTestAgainstGumbel(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=NonStationaryQuadraticScaleModel,
+                         constrained_model_class=GumbelTemporalModel,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
+
+    @classproperty
+    def total_number_of_parameters_for_unconstrained_model(cls) -> int:
+        return 5
\ No newline at end of file
diff --git a/extreme_trend/trend_test_two_parameters/gev_trend_test_two_parameters.py b/extreme_trend/trend_test_two_parameters/gev_trend_test_two_parameters.py
index a7f49d72..794e2b33 100644
--- a/extreme_trend/trend_test_two_parameters/gev_trend_test_two_parameters.py
+++ b/extreme_trend/trend_test_two_parameters/gev_trend_test_two_parameters.py
@@ -1,4 +1,6 @@
 from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE
+from extreme_fit.model.margin_model.polynomial_margin_model import NonStationaryQuadraticLocationModel, \
+    NonStationaryQuadraticScaleModel
 from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest
 from extreme_trend.trend_test_one_parameter.gev_trend_test_one_parameter import \
     GevLocationTrendTest, GevScaleTrendTest
@@ -45,6 +47,26 @@ class GevScaleAndShapeTrendTest(GevTrendTestTwoParametersAgainstGev):
                          constrained_model_class=constrained_model_class,
                          quantile_level=quantile_level,
                          fit_method=fit_method)
+        
+class GevQuadraticLocationTrendTest(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=NonStationaryQuadraticLocationModel,
+                         constrained_model_class=constrained_model_class,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
+
+class GevQuadraticScaleTrendTest(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=NonStationaryQuadraticScaleModel,
+                         constrained_model_class=constrained_model_class,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
 
 
 class GevLocationAndScaleTrendTest(GevTrendTestTwoParametersAgainstGev):
diff --git a/extreme_trend/trend_test_two_parameters/gumbel_test_two_parameters.py b/extreme_trend/trend_test_two_parameters/gumbel_test_two_parameters.py
index 426e593d..f4a9db11 100644
--- a/extreme_trend/trend_test_two_parameters/gumbel_test_two_parameters.py
+++ b/extreme_trend/trend_test_two_parameters/gumbel_test_two_parameters.py
@@ -1,4 +1,6 @@
 from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE
+from extreme_fit.model.margin_model.polynomial_margin_model import NonStationaryQuadraticLocationModel, \
+    NonStationaryQuadraticLocationGumbelModel, NonStationaryQuadraticScaleGumbelModel
 from extreme_trend.trend_test_two_parameters.gev_trend_test_two_parameters import \
     GevTrendTestTwoParameters
 from extreme_fit.distribution.gev.gev_params import GevParams
@@ -41,3 +43,32 @@ class GumbelLocationAndScaleTrendTest(GevTrendTestTwoParameters):
     @classproperty
     def marker(self):
         return 'd'
+
+class GumbelLocationQuadraticTrendTest(GevTrendTestTwoParameters):
+
+    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=NonStationaryQuadraticLocationGumbelModel,
+                         constrained_model_class=GumbelTemporalModel,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
+
+    @classproperty
+    def total_number_of_parameters_for_unconstrained_model(cls) -> int:
+        return 4
+
+class GumbelScaleQuadraticTrendTest(GevTrendTestTwoParameters):
+
+    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=NonStationaryQuadraticScaleGumbelModel,
+                         constrained_model_class=GumbelTemporalModel,
+                         quantile_level=quantile_level,
+                         fit_method=fit_method)
+
+    @classproperty
+    def total_number_of_parameters_for_unconstrained_model(cls) -> int:
+        return 4
+
diff --git a/projects/exceeding_snow_loads/utils.py b/projects/exceeding_snow_loads/utils.py
index 2594b132..477660cd 100644
--- a/projects/exceeding_snow_loads/utils.py
+++ b/projects/exceeding_snow_loads/utils.py
@@ -9,12 +9,14 @@ from extreme_trend.trend_test_one_parameter.gev_trend_test_one_parameter import
 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, GevLocationAndScaleAndShapeTrendTest
+    GevLocationAndScaleTrendTestAgainstGumbel, GevLocationAndScaleAndShapeTrendTest, \
+    GevLocationQuadraticTrendTestAgainstGumbel, \
+    GevScaleQuadraticTrendTestAgainstGumbel
 from extreme_trend.trend_test_two_parameters.gev_trend_test_two_parameters import \
     GevLocationAgainstGumbel, GevScaleAgainstGumbel, GevLocationAndScaleTrendTest, GevScaleAndShapeTrendTest, \
-    GevLocationAndShapeTrendTest
+    GevLocationAndShapeTrendTest, GevQuadraticLocationTrendTest, GevQuadraticScaleTrendTest
 from extreme_trend.trend_test_two_parameters.gumbel_test_two_parameters import \
-    GumbelLocationAndScaleTrendTest
+    GumbelLocationAndScaleTrendTest, GumbelLocationQuadraticTrendTest, GumbelScaleQuadraticTrendTest
 
 paper_altitudes = ALL_ALTITUDES_WITHOUT_NAN
 paper_study_classes = [CrocusSnowLoadTotal, CrocusSnowLoadEurocode, CrocusSnowLoad3Days][:2]
@@ -33,7 +35,10 @@ NON_STATIONARY_TREND_TEST_PAPER_2 = [
     # GEV models with constant shape
     GevVersusGev, GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest,
     # GEV models with linear shape
-    GevShapeTrendTest, GevLocationAndShapeTrendTest, GevScaleAndShapeTrendTest, GevLocationAndScaleAndShapeTrendTest
+    GevShapeTrendTest, GevLocationAndShapeTrendTest, GevScaleAndShapeTrendTest, GevLocationAndScaleAndShapeTrendTest,
+    # Quadratic model for the Gev/Gumbel and for the location/scale
+    GevQuadraticLocationTrendTest, GevQuadraticScaleTrendTest, GumbelLocationQuadraticTrendTest,
+    GumbelScaleQuadraticTrendTest,
 ]
 
 
diff --git a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_polynomial_extremes_mle.py b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_polynomial_extremes_mle.py
index 54817ce0..a5b09acf 100644
--- a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_polynomial_extremes_mle.py
+++ b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_polynomial_extremes_mle.py
@@ -4,7 +4,8 @@ import numpy as np
 import pandas as pd
 
 from extreme_fit.distribution.gev.gev_params import GevParams
-from extreme_fit.model.margin_model.polynomial_margin_model import NonStationaryQuadraticLocationModel
+from extreme_fit.model.margin_model.polynomial_margin_model import NonStationaryQuadraticLocationModel, \
+    NonStationaryQuadraticScaleModel, NonStationaryQuadraticLocationGumbelModel, NonStationaryQuadraticScaleGumbelModel
 from extreme_trend.abstract_gev_trend_test import fitted_linear_margin_estimator
 from extreme_fit.model.margin_model.utils import \
     MarginFitMethod
@@ -38,9 +39,9 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
         self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
         self.fit_method = MarginFitMethod.extremes_fevd_mle
 
-    def test_gev_temporal_margin_fit_non_stationary_quadratic_location(self):
+    def function_test_gev_temporal_margin_fit_non_stationary_quadratic(self, model_class, quadratic_param):
         # Create estimator
-        estimator = fitted_linear_margin_estimator(NonStationaryQuadraticLocationModel,
+        estimator = fitted_linear_margin_estimator(model_class,
                                                    self.coordinates, self.dataset,
                                                    starting_year=0,
                                                    fit_method=self.fit_method)
@@ -51,15 +52,30 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
         self.assertNotEqual(mle_params_estimated_year3, mle_params_estimated_year5)
         # Assert the relationship for the location is indeed quadratic
-        location_gev_params = GevParams.LOC
-        diff1 = mle_params_estimated_year1[location_gev_params] - mle_params_estimated_year3[location_gev_params]
-        diff2 = mle_params_estimated_year3[location_gev_params] - mle_params_estimated_year5[location_gev_params]
+        diff1 = mle_params_estimated_year1[quadratic_param] - mle_params_estimated_year3[quadratic_param]
+        diff2 = mle_params_estimated_year3[quadratic_param] - mle_params_estimated_year5[quadratic_param]
         self.assertNotAlmostEqual(diff1, diff2)
         # Assert that indicators are correctly computed
         self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh())
         self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
         self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
 
+    def test_gev_temporal_margin_fit_non_stationary_quadratic_location(self):
+        self.function_test_gev_temporal_margin_fit_non_stationary_quadratic(NonStationaryQuadraticLocationModel,
+                                                                            quadratic_param=GevParams.LOC)
+
+    def test_gev_temporal_margin_fit_non_stationary_quadratic_scale(self):
+        self.function_test_gev_temporal_margin_fit_non_stationary_quadratic(NonStationaryQuadraticScaleModel,
+                                                                            quadratic_param=GevParams.SCALE)
+
+    def test_gumbel_temporal_margin_fit_non_stationary_quadratic_location(self):
+        self.function_test_gev_temporal_margin_fit_non_stationary_quadratic(NonStationaryQuadraticLocationGumbelModel,
+                                                                            quadratic_param=GevParams.LOC)
+
+    def test_gumbel_temporal_margin_fit_non_stationary_quadratic_scale(self):
+        self.function_test_gev_temporal_margin_fit_non_stationary_quadratic(NonStationaryQuadraticScaleGumbelModel,
+                                                                            quadratic_param=GevParams.SCALE)
+
 
 if __name__ == '__main__':
     unittest.main()
-- 
GitLab