[contrasting] add general test for...

[contrasting] add general test for test_gev_spatio_temporal_polynomial_extremes_mle.py. refactor a bit

extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py

@@ -69,8 +69,8 @@ class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel)
 class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeReturnLevel):
     result_from_fit: ResultFromBayesianExtremes
 
-    def __init__(self, estimator: LinearMarginEstimator, ci_method, temporal_covariate, quantile_level=EUROCODE_QUANTILE):
-        super().__init__(estimator, ci_method, temporal_covariate, quantile_level)
+    def __init__(self, estimator: LinearMarginEstimator, ci_method, coordinate, quantile_level=EUROCODE_QUANTILE):
+        super().__init__(estimator, ci_method, coordinate, quantile_level)
         assert isinstance(self.result_from_fit, ResultFromBayesianExtremes)
 
     @property

extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py

@@ -61,15 +61,13 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
             'type': r.c("return.level")
         }
         if self.param_name_to_dim:
-            print("here")
             if isinstance(transformed_temporal_covariate, (int, float, np.int, np.float, np.int64)):
-                print("here2")
                 d = {GevParams.greek_letter_from_param_name_confidence_interval(param_name) + '1': r.c(transformed_temporal_covariate) for
                      param_name in self.param_name_to_dim.keys()}
             elif isinstance(transformed_temporal_covariate, np.ndarray):
                 d = OrderedDict()
                 linearity_in_shape = GevParams.SHAPE in self.param_name_to_dim
-                nb_calls = 2  # or 4 (1 and 3 did not work for the test)
+                nb_calls = 4  # or 4 (1 and 3 did not work for the test)
                 for param_name in GevParams.PARAM_NAMES:
                     suffix = '0' if param_name in self.param_name_to_dim else ''
                     covariate = np.array([1] * nb_calls)

extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py

+from typing import Union
+
+import numpy as np
+
 from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE
 from extreme_fit.estimator.margin_estimator.utils import fitted_linear_margin_estimator
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
@@ -29,12 +33,12 @@ class EurocodeConfidenceIntervalFromExtremes(object):
     @classmethod
     def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator,
                                 ci_method: ConfidenceIntervalMethodFromExtremes,
-                                temporal_covariate: int,
+                                coordinate: Union[int, np.ndarray],
                                 ):
         if ci_method == ConfidenceIntervalMethodFromExtremes.my_bayes:
-            extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, temporal_covariate, cls.quantile_level)
+            extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, coordinate, cls.quantile_level)
         else:
-            extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, temporal_covariate, cls.quantile_level)
+            extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, coordinate, cls.quantile_level)
         return cls(extractor.mean_estimate,  extractor.confidence_interval)
 
     @classmethod

extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py

@@ -23,7 +23,6 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
         d = self.get_python_dictionary(values)
         if 'par' in d:
             values = {i: param for i, param in enumerate(np.array(d['par']))}
-            print(values)
         else:
             values = {i: np.array(v)[0] for i, v in enumerate(d.values())}
         return get_margin_coef_ordered_dict(self.param_name_to_dim, values, self.type_for_mle,
@@ -38,7 +37,6 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
         res = r('ci.fevd.mle_fixed')(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
         if self.is_non_stationary:
             b = np.array(res)
-            print(b)
             a = b[0]
             lower, mean_estimate, upper, _ = a
         else:

extreme_fit/model/result_from_model_fit/utils.py

@@ -41,5 +41,4 @@ def get_margin_coef_ordered_dict(param_name_to_dims, mle_values, type_for_mle="G
                         coef_name = coef_template.format(k + offset)
                         coef_dict[coef_name] = mle_values[i]
                         i += 1
-    print(coef_dict)
     return coef_dict

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py

@@ -273,11 +273,10 @@ class OneFoldFit(object):
         has_not_opposite_sign = self.sign_of_change(estimator) * self.sign_of_change(new_estimator) >= 0
         return has_not_opposite_sign
 
-
     def sensitivity_of_fit_test_last_years(self, estimator: LinearMarginEstimator):
         # Build the dataset without the maxima for each altitude
         new_dataset = AbstractDataset.remove_last_maxima(self.dataset.observations,
-                                                        self.dataset.coordinates,
+                                                         self.dataset.coordinates,
                                                          nb_years=10)
         # Fit the new estimator
         model_class = type(estimator.margin_model)
@@ -288,7 +287,7 @@ class OneFoldFit(object):
         # Compare sign of change
         has_not_opposite_sign = self.sign_of_change(estimator) * self.sign_of_change(new_estimator) >= 0
         # if not has_not_opposite_sign:
-            # print('Last years', self.massif_name, model_class, self.sign_of_change(estimator), self.sign_of_change(new_estimator))
+        # print('Last years', self.massif_name, model_class, self.sign_of_change(estimator), self.sign_of_change(new_estimator))
         return has_not_opposite_sign
 
     def sign_of_change(self, estimator):
@@ -312,9 +311,9 @@ class OneFoldFit(object):
         standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)]
         return standard_gumbel_quantiles
 
-
-    # def best_confidence_interval(self):
-    #     EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(self.best_estimator,
-    #                                                                    ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
-    #                                                                    temporal_covariate=np.array([2019, self.altitude_plot]),)
-
+    @property
+    def best_confidence_interval(self) -> EurocodeConfidenceIntervalFromExtremes:
+        return EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(self.best_estimator,
+                                                                       ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
+                                                                       coordinate=np.array(
+                                                                           [2019, self.altitude_plot]), )

projects/others/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py

@@ -107,7 +107,7 @@ def get_return_level_bayesian_example(nb_iterations_for_bayesian_fit):
                                                       nb_iterations_for_bayesian_fit=nb_iterations_for_bayesian_fit)
     return_level_bayesian = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator=fitted_estimator,
                                                                              ci_method=ConfidenceIntervalMethodFromExtremes.my_bayes,
-                                                                             temporal_covariate=2019)
+                                                                             coordinate=2019)
     return return_level_bayesian
 
 

test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py

@@ -14,6 +14,8 @@ from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pari
 from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import \
     NonStationaryQuadraticLocationModel, \
     NonStationaryQuadraticScaleModel, NonStationaryQuadraticLocationGumbelModel, NonStationaryQuadraticScaleGumbelModel
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
+    ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
 from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
     ConfidenceIntervalMethodFromExtremes
 from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
@@ -76,45 +78,27 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
         self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
         self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
         # Assert we can compute the return level
-        covariate1_for_return_level = np.array([700, 0])
-        covariate2_for_return_level = np.array([700, 1000])
+        covariate1_for_return_level = np.array([500, 0])
+        covariate2_for_return_level = np.array([500, 50])
         covariate3_for_return_level = np.array([400, 0])
         coordinates = [covariate1_for_return_level, covariate2_for_return_level, covariate3_for_return_level]
         for coordinate in coordinates:
             EurocodeConfidenceIntervalFromExtremes.quantile_level = 0.98
             confidence_interval = EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator,
                                                                                                  ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
-                                                                                                 temporal_covariate=coordinate)
+                                                                                                 coordinate=coordinate)
             gev_params = estimator.function_from_fit.get_params(coordinate)
-            # print(gev_params, coordinate)
             return_level = gev_params.return_level(return_period=50)
-            # print("my return level", return_level)
-            # print("their return level", confidence_interval.mean_estimate)
-            # print("diff", confidence_interval.mean_estimate - return_level)
-            # print('\n\n')
-            self.assertAlmostEqual(return_level, confidence_interval.mean_estimate)
-            self.assertFalse(np.isnan(confidence_interval.confidence_interval[0]))
-            self.assertFalse(np.isnan(confidence_interval.confidence_interval[1]))
-
-    def test_gev_spatio_temporal_margin_fit_1(self):
-        self.function_test_gev_spatio_temporal_margin_fit_non_stationary(StationaryAltitudinal)
-
-    #
-    def test_gev_spatio_temporal_margin_fit_1_bis(self):
-        self.function_test_gev_spatio_temporal_margin_fit_non_stationary(AltitudinalShapeLinearTimeStationary)
-
-    def test_gev_spatio_temporal_margin_fit_2(self):
-        # first model with both a time and altitude non stationarity
-        self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
-            AltitudinalShapeConstantTimeLocationLinear)
-
-    # def test_gev_spatio_temporal_margin_fit_3(self):
-    #     self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
-    #         AltitudinalShapeLinearTimeLocationLinear)
-    #
-    # def test_gev_spatio_temporal_margin_fit_4(self):
-    #     self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
-    #         AltitudinalShapeLinearTimeLocScaleLinear)
+            if np.isnan(return_level) or np.isnan(confidence_interval.mean_estimate):
+                self.assertTrue(np.isnan(return_level) and np.isnan(confidence_interval.mean_estimate))
+            else:
+                self.assertAlmostEqual(return_level, confidence_interval.mean_estimate)
+                self.assertFalse(np.isnan(confidence_interval.confidence_interval[0]))
+                self.assertFalse(np.isnan(confidence_interval.confidence_interval[1]))
+
+    def test_gev_spatio_temporal_all(self):
+        for model_class in ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS:
+            self.function_test_gev_spatio_temporal_margin_fit_non_stationary(model_class)
 
 
 if __name__ == '__main__':