From e939952dc7df04663c8e5d3bb6f4cd2c68416a09 Mon Sep 17 00:00:00 2001 From: Le Roux Erwan <erwan.le-roux@irstea.fr> Date: Mon, 29 Jul 2019 11:12:15 +0200 Subject: [PATCH] [NON STATIONARY] fix constant gev params. also some refactoring for the regions, the colors, and functions. --- .../departementalpesfrancaises.py | 9 ++------- .../scm_models_data/abstract_study.py | 2 +- .../main_study_visualizer.py | 2 +- .../comparisons_visualization.py | 19 ------------------- .../main1_good_stationary_gev_fit.py | 5 +++-- ...main3_non_stationary_strength_evolution.py | 6 +++--- .../abstract_gev_trend_test.py | 6 +++++- .../abstract_univariate_test.py | 4 ++-- .../temporal_linear_margin_model.py | 14 ++++++++++++++ .../extreme_models/result_from_fit.py | 10 ---------- .../test_gev/test_gev_temporal.py | 5 ++--- 11 files changed, 33 insertions(+), 49 deletions(-) diff --git a/experiment/eurocode_data/departementalpesfrancaises.py b/experiment/eurocode_data/departementalpesfrancaises.py index d36d9f41..4b30db05 100644 --- a/experiment/eurocode_data/departementalpesfrancaises.py +++ b/experiment/eurocode_data/departementalpesfrancaises.py @@ -51,11 +51,6 @@ class Drome(AbstractDepartementAlpesFrancaises): super().__init__(C2) -""" -Quand c'est à cheval, je mets les deux massifs -Quand juste un bout du massif est dans un autre departement -(Chartreuse, Belledonne sont un peu en Savoie -""" massif_name_to_departements = { 'Chablais': [HauteSavoie], 'Aravis': [HauteSavoie, Savoie], @@ -63,8 +58,8 @@ massif_name_to_departements = { 'Bauges': [HauteSavoie, Savoie], 'Beaufortain': [HauteSavoie, Savoie], 'Haute-Tarentaise': [Savoie], - 'Chartreuse': [Isere], - 'Belledonne': [Isere], + 'Chartreuse': [Isere, Savoie], + 'Belledonne': [Isere, Savoie], 'Maurienne': [Savoie], 'Vanoise': [Savoie], 'Haute-Maurienne': [Savoie], diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py index 121eec14..424398e0 100644 --- a/experiment/meteo_france_data/scm_models_data/abstract_study.py +++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py @@ -229,7 +229,7 @@ class AbstractStudy(object): def visualize_study(cls, ax=None, massif_name_to_value: Union[None, Dict[str, float]] = None, show=True, fill=True, replace_blue_by_white=True, label=None, add_text=False, cmap=None, vmax=100, vmin=0, - default_color_for_missing_massif='grey', + default_color_for_missing_massif='gainsboro', default_color_for_nan_values='w', massif_name_to_color=None, show_label=True, diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py index 4bc1d621..11fb476c 100644 --- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py +++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py @@ -24,7 +24,7 @@ SCM_STUDY_NAME_TO_SCM_STUDY = dict(zip(SCM_STUDIES_NAMES, SCM_STUDIES)) SCM_STUDY_CLASS_TO_ABBREVIATION = { SafranSnowfall: 'SF3', CrocusTotalSwe: 'TSWE', - CrocusRecentSwe: 'RSWE', + CrocusRecentSwe: 'SWE3', CrocusDepth: 'SD', } diff --git a/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py b/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py index 4e1b2427..5ee81c4f 100644 --- a/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py +++ b/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py @@ -14,7 +14,6 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat VisualizationParameters from experiment.meteo_france_data.stations_data.comparison_analysis import ComparisonAnalysis, MASSIF_COLUMN_NAME, \ REANALYSE_STR, ALTITUDE_COLUMN_NAME, STATION_COLUMN_NAME -from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import GevLocationChangePointTest from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest from experiment.trend_analysis.univariate_test.utils import compute_gev_change_point_test_results from extreme_estimator.extreme_models.result_from_fit import ResultFromIsmev @@ -259,21 +258,3 @@ class ComparisonsVisualization(VisualizationParameters): ax.grid() ax.plot(years, maxima, label=label, color=plot_color) return ordered_dict - - def visualize_gev(self): - return self._visualize_main(self.plot_gev) - - def plot_gev(self, ax, ax2, years, maxima, label, plot_color): - # todo should I normalize here ? - # fit gev - data = maxima - res = safe_run_r_estimator(function=r('gev.fit'), xdat=ro.FloatVector(data), - use_start=True) - res = ResultFromIsmev(res, {}) - gev_params = res.constant_gev_params - - lim = 1.5 * max(data) - x = np.linspace(0, lim, 1000) - y = gev_params.density(x) - # display the gev distribution that was obtained - ax.plot(x, y, label=label, color=plot_color) diff --git a/experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py b/experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py index 0cfb9bcb..935d017c 100644 --- a/experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py +++ b/experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py @@ -8,14 +8,15 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat def maxima_analysis(): save_to_file = False only_first_one = False - durand_altitude = [900, 1500, 2100, 2700] + durand_altitude = [900, 1500, 1800, 2100, 2700][2:-2] altitudes = durand_altitude study_classes = [CrocusRecentSwe][:] for study in study_iterator_global(study_classes, only_first_one=only_first_one, altitudes=altitudes): study_visualizer = StudyVisualizer(study, save_to_file=save_to_file, verbose=True, multiprocessing=True) - study_visualizer.visualize_summary_of_annual_values_and_stationary_gev_fit() + # study_visualizer.visualize_summary_of_annual_values_and_stationary_gev_fit() + study_visualizer.visualize_all_mean_and_max_graphs() if __name__ == '__main__': diff --git a/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py b/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py index 5ec783b1..29713ba1 100644 --- a/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py +++ b/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py @@ -13,7 +13,7 @@ from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALT def main_fast_spatial_risk_evolution(): - for altitude in FULL_ALTITUDES[-1:]: + for altitude in [1800]: vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude, exact_starting_year=1958, reduce_strength_array=True, trend_test_class=GevLocationAndScaleTrendTest) @@ -31,8 +31,8 @@ def main_full_spatial_risk_evolution(): def main_run(): - main_full_spatial_risk_evolution() - # main_fast_spatial_risk_evolution() + # main_full_spatial_risk_evolution() + main_fast_spatial_risk_evolution() if __name__ == '__main__': diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py index e3584e95..456390e3 100644 --- a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py +++ b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py @@ -113,14 +113,18 @@ class AbstractGevTrendTest(AbstractUnivariateTest): @property def non_stationary_constant_gev_params(self) -> GevParams: - return self.non_stationary_estimator.result_from_fit.constant_gev_params + # Constant parameters correspond to the gev params in 1958 + return self.non_stationary_estimator.margin_function_fitted.get_gev_params(coordinate=np.array([1958]), + is_transformed=False) @property def test_trend_slope_strength(self): if self.crashed: return 0.0 else: + # Compute the slope strength slope = self._slope_strength() + # Delta T must in the same unit as were the parameter of slope mu1 and sigma1 slope *= self.nb_years_for_quantile_evolution * self.coordinates.transformed_distance_between_two_successive_years[0] return slope diff --git a/experiment/trend_analysis/univariate_test/abstract_univariate_test.py b/experiment/trend_analysis/univariate_test/abstract_univariate_test.py index 7381e71c..7cbbf78e 100644 --- a/experiment/trend_analysis/univariate_test/abstract_univariate_test.py +++ b/experiment/trend_analysis/univariate_test/abstract_univariate_test.py @@ -61,8 +61,8 @@ class AbstractUnivariateTest(object): d[cls.ALL_TREND] = 'k-' d[cls.NON_SIGNIFICATIVE_TREND] = 'b-' # d[cls.SIGNIFICATIVE_ALL_TREND] = 'k-' - d[cls.SIGNIFICATIVE_POSITIVE_TREND] = 'g-' - d[cls.SIGNIFICATIVE_NEGATIVE_TREND] = 'r-' + d[cls.SIGNIFICATIVE_POSITIVE_TREND] = 'darkgreen-' + d[cls.SIGNIFICATIVE_NEGATIVE_TREND] = 'darkred-' # d[cls.NO_TREND] = 'k--' return d 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 90d95ba1..156558e1 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 @@ -37,6 +37,10 @@ class TemporalLinearMarginModel(LinearMarginModel): def shl(self): return get_null() + @property + def siglink(self): + return r('identity') + class StationaryStationModel(TemporalLinearMarginModel): @@ -64,6 +68,16 @@ class NonStationaryScaleStationModel(TemporalLinearMarginModel): return 1 +class NonStationaryLogScaleStationModel(NonStationaryScaleStationModel): + + def load_margin_functions(self, gev_param_name_to_dims=None): + super().load_margin_functions({GevParams.SCALE: [self.coordinates.idx_temporal_coordinates]}) + + @property + def siglink(self): + return r('exp') + + class NonStationaryShapeStationModel(TemporalLinearMarginModel): def load_margin_functions(self, gev_param_name_to_dims=None): diff --git a/extreme_estimator/extreme_models/result_from_fit.py b/extreme_estimator/extreme_models/result_from_fit.py index 78d9d83a..c3170e18 100644 --- a/extreme_estimator/extreme_models/result_from_fit.py +++ b/extreme_estimator/extreme_models/result_from_fit.py @@ -48,10 +48,6 @@ class ResultFromFit(object): def convergence(self) -> str: raise NotImplementedError - @property - def constant_gev_params(self) -> GevParams: - raise NotImplementedError - class ResultFromIsmev(ResultFromFit): @@ -79,12 +75,6 @@ class ResultFromIsmev(ResultFromFit): i += 1 return coef_dict - @property - def constant_gev_params(self) -> GevParams: - params = {k.split('Coeff1')[0]: v for k, v in self.margin_coef_dict.items() - if 'Coeff1' in k and 'temp' not in k} - return GevParams.from_dict(params) - @property def all_parameters(self): return self.margin_coef_dict 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 8ffe3513..f0a17fb3 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 @@ -33,8 +33,6 @@ 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 @@ -49,7 +47,8 @@ class TestGevTemporal(unittest.TestCase): def test_gev_temporal_margin_fit_nonstationary(self): # Create estimator - for margin_model in self.margin_models: + margin_models = load_non_stationary_temporal_margin_models(self.coordinates) + for margin_model in margin_models: # margin_model = NonStationaryLocationStationModel(self.coordinates) estimator = LinearMarginEstimator(self.dataset, margin_model) estimator.fit() -- GitLab