diff --git a/experiment/meteo_france_SCM_study/visualization/studies_visualization/main_studies_visualizer.py b/experiment/meteo_france_SCM_study/visualization/studies_visualization/main_studies_visualizer.py index 19d45bc5db026bf820bfb630567722c7b5049be6..ea254cd671b23842e70d19a175d1e3b302ca34c7 100644 --- a/experiment/meteo_france_SCM_study/visualization/studies_visualization/main_studies_visualizer.py +++ b/experiment/meteo_france_SCM_study/visualization/studies_visualization/main_studies_visualizer.py @@ -1,5 +1,6 @@ from experiment.trend_analysis.abstract_score import MannKendall, WeigthedScore, MeanScore, MedianScore -from experiment.trend_analysis.univariate_trend_test.abstract_gev_trend_test import GevLocationTrendTest +from experiment.trend_analysis.univariate_trend_test.abstract_gev_trend_test import GevLocationTrendTest, \ + GevScaleTrendTest, GevShapeTrendTest from experiment.trend_analysis.univariate_trend_test.abstract_trend_test import MannKendallTrendTest from experiment.meteo_france_SCM_study.safran.safran import ExtendedSafranTotalPrecip from experiment.meteo_france_SCM_study.visualization.studies_visualization.studies import Studies @@ -45,7 +46,7 @@ def altitude_trends_significant(): # altitudes = ALL_ALTITUDES[3:5] # altitudes = ALL_ALTITUDES[2:4] for study_class in SCM_STUDIES[:1]: - trend_test_classes = [MannKendallTrendTest, GevLocationTrendTest][1:] + trend_test_classes = [MannKendallTrendTest, GevLocationTrendTest, GevScaleTrendTest, GevShapeTrendTest][2:] visualizers = [StudyVisualizer(study, temporal_non_stationarity=True, verbose=False) for study in study_iterator_global(study_classes=[study_class], only_first_one=only_first_one, altitudes=altitudes)] diff --git a/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py index 77248e32fba18652260fee0496faeb43398c271a..b5b3fd3e1c0c669d6d86ab192468fef257fc7881 100644 --- a/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py +++ b/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py @@ -5,7 +5,8 @@ from sklearn.preprocessing import normalize from experiment.trend_analysis.univariate_trend_test.abstract_trend_test import AbstractTrendTest from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator from extreme_estimator.extreme_models.margin_model.temporal_linear_margin_model import StationaryStationModel, \ - NonStationaryLocationStationModel + NonStationaryLocationStationModel, NonStationaryScaleStationModel, NonStationaryShapeStationModel +from extreme_estimator.margin_fits.gev.gev_params import GevParams from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \ AbstractTemporalCoordinates @@ -55,4 +56,27 @@ class GevLocationTrendTest(AbstractGevTrendTest): @property def test_sign(self) -> int: - return np.sign(self.non_stationary_estimator.margin_function_fitted.mu1_temporal_trend) + return np.sign(self.non_stationary_estimator.margin_function_fitted.get_coef(GevParams.LOC, + AbstractCoordinates.COORDINATE_T)) + + +class GevScaleTrendTest(AbstractGevTrendTest): + + def __init__(self, years_after_change_point, maxima_after_change_point): + super().__init__(years_after_change_point, maxima_after_change_point, NonStationaryScaleStationModel) + + @property + def test_sign(self) -> int: + return np.sign(self.non_stationary_estimator.margin_function_fitted.get_coef(GevParams.SCALE, + AbstractCoordinates.COORDINATE_T)) + + +class GevShapeTrendTest(AbstractGevTrendTest): + + def __init__(self, years_after_change_point, maxima_after_change_point): + super().__init__(years_after_change_point, maxima_after_change_point, NonStationaryShapeStationModel) + + @property + def test_sign(self) -> int: + return np.sign(self.non_stationary_estimator.margin_function_fitted.get_coef(GevParams.SHAPE, + AbstractCoordinates.COORDINATE_T)) diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py index ebbfd62fd4dd3205c58b41c3a59c16eaa1cddd4a..4c81dd5c38ad667a5e77071d1793b76b55ffa903 100644 --- a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py +++ b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py @@ -82,6 +82,11 @@ class LinearMarginFunction(ParametricMarginFunction): # Properties for the location parameter + def get_coef(self, gev_param_name, coef_name): + idx = 1 if coef_name in [AbstractCoordinates.COORDINATE_T] \ + else AbstractCoordinates.COORDINATE_SPATIAL_NAMES.index(coef_name) + 2 + return self.coef_dict[LinearCoef.coef_template_str(gev_param_name, coef_name).format(idx)] + @property def mu1_temporal_trend(self): return self.coef_dict[LinearCoef.coef_template_str(ExtremeParams.LOC, AbstractCoordinates.COORDINATE_T).format(1)] diff --git a/extreme_estimator/extreme_models/margin_model/temporal_linear_margin_model.py b/extreme_estimator/extreme_models/margin_model/temporal_linear_margin_model.py index aee5481c372ce97f46d8f58b38e2294b9d9c5370..48237de33f77c42c7709ac2004a08846d4ade4be 100644 --- a/extreme_estimator/extreme_models/margin_model/temporal_linear_margin_model.py +++ b/extreme_estimator/extreme_models/margin_model/temporal_linear_margin_model.py @@ -23,12 +23,21 @@ class TemporalLinearMarginModel(LinearMarginModel): # Gev Fit assert data.shape[1] == len(df_coordinates_temp.values) res = safe_run_r_estimator(function=r('gev.fit'), use_start=self.use_start_value, - xdat=ro.FloatVector(data[0]), y=df_coordinates_temp.values, mul=self.mul) + xdat=ro.FloatVector(data[0]), y=df_coordinates_temp.values, mul=self.mul, + sigl=self.sigl, shl=self.shl) return ResultFromIsmev(res, self.margin_function_start_fit.gev_param_name_to_dims) @property def mul(self): - return NotImplementedError + return get_null() + + @property + def sigl(self): + return get_null() + + @property + def shl(self): + return get_null() class StationaryStationModel(TemporalLinearMarginModel): @@ -36,10 +45,6 @@ class StationaryStationModel(TemporalLinearMarginModel): def load_margin_functions(self, gev_param_name_to_dims=None): super().load_margin_functions({}) - @property - def mul(self): - return get_null() - class NonStationaryLocationStationModel(TemporalLinearMarginModel): @@ -51,3 +56,21 @@ class NonStationaryLocationStationModel(TemporalLinearMarginModel): return 1 +class NonStationaryScaleStationModel(TemporalLinearMarginModel): + + def load_margin_functions(self, gev_param_name_to_dims=None): + super().load_margin_functions({GevParams.SCALE: [self.coordinates.idx_temporal_coordinates]}) + + @property + def sigl(self): + return 1 + + +class NonStationaryShapeStationModel(TemporalLinearMarginModel): + + def load_margin_functions(self, gev_param_name_to_dims=None): + super().load_margin_functions({GevParams.SHAPE: [self.coordinates.idx_temporal_coordinates]}) + + @property + def shl(self): + return 1 diff --git a/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_temporal.py b/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_temporal.py index 38636360b2a5d375a273a8211d66be3bc400bf4d..d8b1c71a98519ea9251ba3d64f79b12b482ffac1 100644 --- a/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_temporal.py +++ b/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_temporal.py @@ -13,6 +13,7 @@ from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_ from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \ AbstractSpatioTemporalObservations +from test.test_utils import load_non_stationary_temporal_margin_models class TestGevTemporal(unittest.TestCase): @@ -32,6 +33,8 @@ class TestGevTemporal(unittest.TestCase): df2 = pd.DataFrame(data=np.array(r['x_gev']), index=df.index) observations = AbstractSpatioTemporalObservations(df_maxima_gev=df2) self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates) + # + self.margin_models = load_non_stationary_temporal_margin_models(self.coordinates) def test_gev_temporal_margin_fit_stationary(self): # Create estimator @@ -46,14 +49,14 @@ class TestGevTemporal(unittest.TestCase): def test_gev_temporal_margin_fit_nonstationary(self): # Create estimator - margin_model = NonStationaryLocationStationModel(self.coordinates) - estimator = LinearMarginEstimator(self.dataset, margin_model) - estimator.fit() - self.assertNotEqual(estimator.margin_function_fitted.mu1_temporal_trend, 0.0) - # Checks that parameters returned are indeed different - mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(np.array([1])).to_dict() - mle_params_estimated_year3 = estimator.margin_function_fitted.get_gev_params(np.array([3])).to_dict() - self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3) + for margin_model in self.margin_models: + # margin_model = NonStationaryLocationStationModel(self.coordinates) + estimator = LinearMarginEstimator(self.dataset, margin_model) + estimator.fit() + # Checks that parameters returned are indeed different + mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(np.array([1])).to_dict() + mle_params_estimated_year3 = estimator.margin_function_fitted.get_gev_params(np.array([3])).to_dict() + self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3) def test_gev_temporal_margin_fit_nonstationary_with_start_point(self): # Create estimator @@ -80,7 +83,6 @@ class TestGevTemporal(unittest.TestCase): estimator2 = self.fit_non_stationary_estimator(starting_point=28) mu1_estimator1 = estimator1.margin_function_fitted.mu1_temporal_trend mu1_estimator2 = estimator2.margin_function_fitted.mu1_temporal_trend - print(mu1_estimator1, mu1_estimator2) self.assertNotEqual(mu1_estimator1, mu1_estimator2) diff --git a/test/test_utils.py b/test/test_utils.py index 072bd796382121dc76042e95bdb38bbb19e62a9f..e3eb691e16f7590aa0f00843aeedbca2613db220 100644 --- a/test/test_utils.py +++ b/test/test_utils.py @@ -9,6 +9,8 @@ from extreme_estimator.estimator.full_estimator.abstract_full_estimator import S from extreme_estimator.estimator.max_stable_estimator.abstract_max_stable_estimator import MaxStableEstimator from extreme_estimator.extreme_models.margin_model.linear_margin_model import LinearAllParametersAllDimsMarginModel, \ ConstantMarginModel +from extreme_estimator.extreme_models.margin_model.temporal_linear_margin_model import \ + NonStationaryLocationStationModel, NonStationaryScaleStationModel, NonStationaryShapeStationModel from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import \ AbstractMaxStableModelWithCovarianceFunction, CovarianceFunction from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith, BrownResnick, Schlather, \ @@ -37,10 +39,15 @@ TEST_3D_SPATIAL_COORDINATES = [AlpsStation3DCoordinatesWithAnisotropy] TEST_TEMPORAL_COORDINATES = [ConsecutiveTemporalCoordinates] TEST_SPATIO_TEMPORAL_COORDINATES = [UniformSpatioTemporalCoordinates, LinSpaceSpatial2DSpatioTemporalCoordinates] TEST_MARGIN_TYPES = [ConstantMarginModel, LinearAllParametersAllDimsMarginModel][:] +TEST_NON_STATIONARY_TEMPORAL_MARGIN_TYPES = [NonStationaryLocationStationModel, NonStationaryScaleStationModel, NonStationaryShapeStationModel] TEST_MAX_STABLE_ESTIMATOR = [MaxStableEstimator] TEST_FULL_ESTIMATORS = [SmoothMarginalsThenUnitaryMsp, FullEstimatorInASingleStepWithSmoothMargin][:] +def load_non_stationary_temporal_margin_models(coordinates): + return [margin_class(coordinates=coordinates) for margin_class in TEST_NON_STATIONARY_TEMPORAL_MARGIN_TYPES] + + def load_test_full_estimators(dataset, margin_model, max_stable_model): return [full_estimator(dataset=dataset, margin_model=margin_model, max_stable_model=max_stable_model) for full_estimator in TEST_FULL_ESTIMATORS] @@ -89,7 +96,8 @@ def load_test_temporal_coordinates(nb_steps, train_split_ratio=None): TEST_TEMPORAL_COORDINATES] -def load_test_spatiotemporal_coordinates(nb_points, nb_steps, train_split_ratio=None, transformation_class: type = None): +def load_test_spatiotemporal_coordinates(nb_points, nb_steps, train_split_ratio=None, + transformation_class: type = None): return [coordinate_class.from_nb_points_and_nb_steps(nb_points=nb_points, nb_steps=nb_steps, train_split_ratio=train_split_ratio, transformation_class=transformation_class)