diff --git a/extreme_data/meteo_france_data/scm_models_data/altitudes_studies.py b/extreme_data/meteo_france_data/scm_models_data/altitudes_studies.py index d094282551aa9acf8fdfbf676c5796626fbadd56..730fa5eab8e83cea3dd99c7273ccf8063fb0f5ac 100644 --- a/extreme_data/meteo_france_data/scm_models_data/altitudes_studies.py +++ b/extreme_data/meteo_france_data/scm_models_data/altitudes_studies.py @@ -18,6 +18,8 @@ from spatio_temporal_dataset.coordinates.spatio_temporal_coordinates.abstract_sp AbstractSpatioTemporalCoordinates from spatio_temporal_dataset.coordinates.spatio_temporal_coordinates.spatio_temporal_coordinates_for_climate_models import \ SpatioTemporalCoordinatesForClimateModels +from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \ + AbstractTemporalCoordinates from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal_coordinates import \ ConsecutiveTemporalCoordinates from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset @@ -57,8 +59,17 @@ class AltitudesStudies(object): for altitude in massif_altitudes: study = self.altitude_to_study[altitude] for year, maxima in zip(study.ordered_years, study.massif_name_to_annual_maxima[massif_name]): - coordinate_values_to_maxima[(altitude, year)] = [maxima] + if len(massif_altitudes) == 1: + coordinate_values_to_maxima[year] = [maxima] + else: + coordinate_values_to_maxima[(altitude, year)] = [maxima] + coordinates = self.spatio_temporal_coordinates(s_split_spatial, s_split_temporal, massif_altitudes) + # Remove the spatial coordinate if we only have one altitude + if len(massif_altitudes) == 1: + df = pd.concat([coordinates.df_temporal_coordinates(), coordinates.df_coordinate_climate_model], axis=1) + coordinates = AbstractTemporalCoordinates.from_df(df) + observations = AnnualMaxima.from_coordinates(coordinates, coordinate_values_to_maxima) return AbstractDataset(observations=observations, coordinates=coordinates) diff --git a/extreme_fit/function/param_function/spline_coef.py b/extreme_fit/function/param_function/spline_coef.py index 204a8f68bc411505b23c99d81b052a64b7f3cb16..0c62b8befdd8e69f671c95c7061993060e3143c2 100644 --- a/extreme_fit/function/param_function/spline_coef.py +++ b/extreme_fit/function/param_function/spline_coef.py @@ -27,6 +27,9 @@ class SplineCoef(AbstractCoef): def nb_coefficients(self): return len(self.coefficients) + @property + def nb_params(self): + return self.nb_knots + self.nb_coefficients class SplineAllCoef(LinearCoef): @@ -45,3 +48,7 @@ class SplineAllCoef(LinearCoef): spline_coef.nb_coefficients)) formula_str = ' '.join(formula_list) return {self.param_name + '.form': self.param_name + ' ~ ' + formula_str} + + @property + def nb_params(self): + return sum([spline_coef.nb_params for spline_coef in self.dim_to_spline_coef.values()]) \ No newline at end of file diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py index f852fdd6fcef500d60e46bc45dfa51ff48b96a65..70238d0f39d198f440a5d707091d17a21ba27144 100644 --- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py +++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py @@ -53,7 +53,7 @@ class ResultFromEvgam(AbstractResultFromExtremes): @property def nb_parameters(self): - return len(np.array(self.name_to_value['coefficients'])) + return len(np.array(self.name_to_value['coefficients'])) + self.nb_knots @property def aic(self): @@ -107,14 +107,20 @@ class ResultFromEvgam(AbstractResultFromExtremes): raise NotImplementedError else: dim, max_degree = dims[0] - d = self.get_python_dictionary(self.name_to_value[r_param_name]) - smooth = list(d['smooth'])[0] - knots = np.array(self.get_python_dictionary(smooth)['knots']) + knots = self.load_knots(r_param_name) x = np.array(self.name_to_value["data"])[1] y = np.array(self.get_python_dictionary(self.name_to_value[r_param_name])['fitted']) assert len(knots) == 5 x_for_interpolation = [int(knots[1]+1), int((knots[1] + knots[3])/2), int(knots[3]-1)] - index = [i for i, e in enumerate(x) if e in x_for_interpolation][:len(x_for_interpolation)] + + # For the time covariate, the distance will be zero for the closer year + # For the temperature covariate, the distance will be minimal for the closer covariate + index = [] + for x_to_find in x_for_interpolation: + distances = np.power(x - x_to_find, 2) + closer_index = np.argmin(distances) + index.append(closer_index) + x = [x[i] for i in index] y = [y[i] for i in index] spline = make_interp_spline(x, y, k=1, t=knots) @@ -123,6 +129,19 @@ class ResultFromEvgam(AbstractResultFromExtremes): dim_knots_and_coefficient[dim] = (knots, coefficients) return dim_knots_and_coefficient + def load_knots(self, r_param_name): + try: + d = self.get_python_dictionary(self.name_to_value[r_param_name]) + smooth = list(d['smooth'])[0] + knots = np.array(self.get_python_dictionary(smooth)['knots']) + except (IndexError, KeyError): + knots = [] + return knots + + @property + def nb_knots(self): + return sum([len(self.load_knots(r_param_name)) for r_param_name in self.r_param_name_to_param_name.keys()]) + @property def r_param_name_to_param_name(self): return { diff --git a/extreme_trend/ensemble_fit/independent_ensemble_fit/one_fold_fit_merge.py b/extreme_trend/ensemble_fit/independent_ensemble_fit/one_fold_fit_merge.py index d749c0811cfbc795adb9ebe0064ec5533bb91014..48e012b6813f356150e4cf25e78f6f9559b03e90 100644 --- a/extreme_trend/ensemble_fit/independent_ensemble_fit/one_fold_fit_merge.py +++ b/extreme_trend/ensemble_fit/independent_ensemble_fit/one_fold_fit_merge.py @@ -7,11 +7,12 @@ from extreme_trend.one_fold_fit.one_fold_fit import OneFoldFit class OneFoldFitMerge(OneFoldFit): - def __init__(self, one_fold_fit_list: List[OneFoldFit], massif_name, altitude_class, temporal_covariate_for_fit, + def __init__(self, one_fold_fit_list: List[OneFoldFit], massif_name, + altitude_group, temporal_covariate_for_fit, first_year, last_year, merge_function=np.median): assert len(one_fold_fit_list) > 0 self.one_fold_fit_list = one_fold_fit_list - self.altitude_group = altitude_class() + self.altitude_group = altitude_group self.massif_name = massif_name self.temporal_covariate_for_fit = temporal_covariate_for_fit self.merge_function = merge_function diff --git a/extreme_trend/ensemble_fit/independent_ensemble_fit/visualizer_merge.py b/extreme_trend/ensemble_fit/independent_ensemble_fit/visualizer_merge.py index 5ed3fcda9452a796ef17399bd88e4ed846e70627..0c3c12fff32cad0a606fe5c336459d85d7b9614c 100644 --- a/extreme_trend/ensemble_fit/independent_ensemble_fit/visualizer_merge.py +++ b/extreme_trend/ensemble_fit/independent_ensemble_fit/visualizer_merge.py @@ -36,8 +36,8 @@ class VisualizerMerge(AltitudesStudiesVisualizerForNonStationaryModels): one_fold_fit_list = [v.massif_name_to_one_fold_fit[massif_name] for v in self.visualizers if massif_name in v.massif_name_to_one_fold_fit] if len(one_fold_fit_list) > 0: - one_fold_fit_merge = OneFoldFitMerge(one_fold_fit_list, massif_name, - type(self.altitude_group), self.temporal_covariate_for_fit, + one_fold_fit_merge = OneFoldFitMerge(one_fold_fit_list, massif_name, self.altitude_group, + self.temporal_covariate_for_fit, first_year=self.study.year_min, last_year=self.study.year_max, merge_function=self.merge_function) self._massif_name_to_one_fold_fit[massif_name] = one_fold_fit_merge diff --git a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py index 096e5e9ec9f7b65f780fdeeb2ec90b9ef6d4e5b5..dfd0943656d83525bd78d0433af37f75e121be5f 100644 --- a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py +++ b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py @@ -100,7 +100,7 @@ class VisualizerForProjectionEnsemble(object): for ensemble_fit_class in self.ensemble_fit_classes: for ensemble_fit in self.ensemble_fits(ensemble_fit_class): visualizer_list.extend(ensemble_fit.visualizer_list) - compute_and_assign_max_abs(visualizer_list) + # compute_and_assign_max_abs(visualizer_list) # Plot if IndependentEnsembleFit in self.ensemble_fit_classes: self.plot_independent() diff --git a/extreme_trend/one_fold_fit/altitude_group.py b/extreme_trend/one_fold_fit/altitude_group.py index 60924af08cbf243e402c5adf75d287b054d94ab4..c1d01f14611a3517cde3017fa9fdf87b4c9b818e 100644 --- a/extreme_trend/one_fold_fit/altitude_group.py +++ b/extreme_trend/one_fold_fit/altitude_group.py @@ -136,13 +136,17 @@ class VeyHighAltitudeGroup(AbstractAltitudeGroup): class DefaultAltitudeGroup(AbstractAltitudeGroup): + def __init__(self, altitudes): + assert len(altitudes) == 1 + self.altitude = list(altitudes)[0] + @property def name(self): - return 'default' + return str(self.altitude) @property def reference_altitude(self): - return 500 + return self.altitude def get_altitude_class_from_altitudes(altitudes): @@ -160,7 +164,7 @@ def get_altitude_group_from_altitudes(altitudes): elif s == set(altitudes_for_groups[3]): return VeyHighAltitudeGroup() else: - return DefaultAltitudeGroup() + return DefaultAltitudeGroup(altitudes) def get_linestyle_for_altitude_class(altitude_class): diff --git a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py index 7706a557c0ca297064cb2f5f3451e52643e98fe5..6b1a375d3f1c907b1e831a6893ac166ce5587941 100644 --- a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py +++ b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py @@ -30,6 +30,7 @@ from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covari class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): + consider_at_least_two_altitudes = True def __init__(self, studies: AltitudesStudies, model_classes: List[AbstractSpatioTemporalPolynomialModel], @@ -85,7 +86,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): self.study.year_max, self.fit_method, self.temporal_covariate_for_fit, - type(self.altitude_group), + self.altitude_group, self.display_only_model_that_pass_test, self.confidence_interval_based_on_delta_method, self.remove_physically_implausible_models) @@ -120,9 +121,10 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): def load_condition(self, massif_altitudes): # At least two altitudes for the estimated # reference_altitude_is_in_altitudes = (self.altitude_group.reference_altitude in massif_altitudes) - at_least_two_altitudes = (len(massif_altitudes) >= 2) - # return reference_altitude_is_in_altitudes and at_least_two_altitudes - return at_least_two_altitudes + if self.consider_at_least_two_altitudes: + return len(massif_altitudes) >= 2 + else: + return True @property def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]: diff --git a/extreme_trend/one_fold_fit/one_fold_fit.py b/extreme_trend/one_fold_fit/one_fold_fit.py index 7cdfc50f48677405d7ba64cf024f4b0866da6e81..7fe04a0a7d0e061b44c798e5402e1628389d23d5 100644 --- a/extreme_trend/one_fold_fit/one_fold_fit.py +++ b/extreme_trend/one_fold_fit/one_fold_fit.py @@ -47,7 +47,7 @@ class OneFoldFit(object): first_year, last_year, fit_method=MarginFitMethod.extremes_fevd_mle, temporal_covariate_for_fit=None, - altitude_class=DefaultAltitudeGroup, + altitude_group=None, only_models_that_pass_goodness_of_fit_test=True, confidence_interval_based_on_delta_method=False, remove_physically_implausible_models=False, @@ -58,7 +58,7 @@ class OneFoldFit(object): self.remove_physically_implausible_models = remove_physically_implausible_models self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test - self.altitude_group = altitude_class() + self.altitude_group = altitude_group self.massif_name = massif_name self.dataset = dataset self.models_classes = models_classes @@ -197,7 +197,11 @@ class OneFoldFit(object): return sorted_estimators def _compute_shape_for_reference_altitude(self, estimator): - coordinate = np.array([self.altitude_plot, self.last_year]) + if isinstance(self.altitude_group, DefaultAltitudeGroup): + coordinate = np.array([self.last_year]) + else: + coordinate = np.array([self.altitude_plot, self.last_year]) + print(coordinate) gev_params = estimator.function_from_fit.get_params(coordinate, is_transformed=False) shape = gev_params.shape return shape diff --git a/extreme_trend/one_fold_fit/utils_altitude_studies_visualizer.py b/extreme_trend/one_fold_fit/utils_altitude_studies_visualizer.py index a6d6ae62bcada070bb355c2328546aca733703ac..9f36764eae275c53c260b80b4f912d70ac043635 100644 --- a/extreme_trend/one_fold_fit/utils_altitude_studies_visualizer.py +++ b/extreme_trend/one_fold_fit/utils_altitude_studies_visualizer.py @@ -5,7 +5,8 @@ from extreme_trend.one_fold_fit.altitudes_studies_visualizer_for_non_stationary_ AltitudesStudiesVisualizerForNonStationaryModels -def load_visualizer_list(season, study_class, altitudes_list, massif_names, model_must_pass_the_test=True, **kwargs_study): +def load_visualizer_list(season, study_class, altitudes_list, massif_names, model_must_pass_the_test=True, + do_compute_and_assign_max_abs=True, **kwargs_study): model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS visualizer_list = [] # Load all studies @@ -20,7 +21,8 @@ def load_visualizer_list(season, study_class, altitudes_list, massif_names, mode display_only_model_that_pass_anderson_test=model_must_pass_the_test ) visualizer_list.append(visualizer) - compute_and_assign_max_abs(visualizer_list) + if do_compute_and_assign_max_abs: + compute_and_assign_max_abs(visualizer_list) return visualizer_list diff --git a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py b/projects/projected_extreme_snowfall/results/main_projections_ensemble.py index 864f6f95f88a9a06a432d5378ed7d5903879eaf1..2069ed808f519b32c04d49d8d7d963ef9b3958f0 100644 --- a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py +++ b/projects/projected_extreme_snowfall/results/main_projections_ensemble.py @@ -3,7 +3,11 @@ import time from typing import List import matplotlib +from extreme_fit.model.margin_model.utils import MarginFitMethod from extreme_trend.ensemble_fit.together_ensemble_fit.together_ensemble_fit import TogetherEnsembleFit +from extreme_trend.one_fold_fit.altitudes_studies_visualizer_for_non_stationary_models import \ + AltitudesStudiesVisualizerForNonStationaryModels +from projects.projected_extreme_snowfall.results.utils import SPLINE_MODELS_FOR_PROJECTION_ONE_ALTITUDE matplotlib.use('Agg') import matplotlib as mpl @@ -41,49 +45,49 @@ def main(): study_class = AdamontSnowfall ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][:1] temporal_covariate_for_fit = [TimeTemporalCovariate, - AnomalyTemperatureWithSplineTemporalCovariate][0] + AnomalyTemperatureWithSplineTemporalCovariate][1] set_seed_for_test() AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2 + scenarios = [AdamontScenario.rcp85_extended] - fast = False - scenarios = rcp_scenarios[::-1] if fast is False else [AdamontScenario.rcp85] - scenarios = rcm_scenarios_extended[::-1] - - scenarios = [AdamontScenario.histo] - gcm_to_year_min_and_year_max = { - gcm: (1959, 2005) for gcm in get_gcm_list(adamont_version=2) - } - + fast = True for scenario in scenarios: gcm_rcm_couples = get_gcm_rcm_couples(scenario) if fast is None: - massif_names = None gcm_rcm_couples = gcm_rcm_couples[:2] AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10 - altitudes_list = altitudes_for_groups[3:] + altitudes_list = [1800, 2100] elif fast: - massif_names = ['Vanoise', 'Haute-Maurienne'] gcm_rcm_couples = gcm_rcm_couples[:2] AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10 - altitudes_list = altitudes_for_groups[:1] + altitudes_list = [2400] else: - massif_names = None - altitudes_list = altitudes_for_groups[:] + altitudes_list = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600] assert isinstance(gcm_rcm_couples, list) + altitudes_list = [[a] for a in altitudes_list] assert isinstance(altitudes_list, List) assert isinstance(altitudes_list[0], List) + for altitudes in altitudes_list: + assert len(altitudes) == 1 + AltitudesStudiesVisualizerForNonStationaryModels.consider_at_least_two_altitudes = False + print('Scenario is', scenario) print('Covariate is {}'.format(temporal_covariate_for_fit)) - model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS + # Default parameters + gcm_to_year_min_and_year_max = None + massif_names = ['Vanoise'] + model_classes = SPLINE_MODELS_FOR_PROJECTION_ONE_ALTITUDE + visualizer = VisualizerForProjectionEnsemble( altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario, model_classes=model_classes, ensemble_fit_classes=ensemble_fit_classes, massif_names=massif_names, + fit_method=MarginFitMethod.evgam, temporal_covariate_for_fit=temporal_covariate_for_fit, remove_physically_implausible_models=True, gcm_to_year_min_and_year_max=gcm_to_year_min_and_year_max, diff --git a/projects/projected_extreme_snowfall/results/utils.py b/projects/projected_extreme_snowfall/results/utils.py new file mode 100644 index 0000000000000000000000000000000000000000..7c711d3314fe292783858c0122396cb88b6e2e03 --- /dev/null +++ b/projects/projected_extreme_snowfall/results/utils.py @@ -0,0 +1,16 @@ +from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \ + NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel +from extreme_fit.model.margin_model.spline_margin_model.temporal_spline_model_degree_1 import \ + NonStationaryTwoLinearLocationModel, NonStationaryTwoLinearScaleModel + +SPLINE_MODELS_FOR_PROJECTION_ONE_ALTITUDE = [ + StationaryTemporalModel, + + # Location only non-stationarity + NonStationaryLocationTemporalModel, + NonStationaryTwoLinearLocationModel, + + # Scale only non-stationarity + NonStationaryScaleTemporalModel, + NonStationaryTwoLinearScaleModel, +] diff --git a/spatio_temporal_dataset/spatio_temporal_observations/annual_maxima_observations.py b/spatio_temporal_dataset/spatio_temporal_observations/annual_maxima_observations.py index c8a889b7ad3f74f64710db04c1498a6ca8b41ad5..abd109d1f96f74df3afff1e706a4b795af3fa5a3 100644 --- a/spatio_temporal_dataset/spatio_temporal_observations/annual_maxima_observations.py +++ b/spatio_temporal_dataset/spatio_temporal_observations/annual_maxima_observations.py @@ -20,7 +20,10 @@ class AnnualMaxima(AbstractSpatioTemporalObservations): def from_coordinates(cls, coordinates: AbstractCoordinates, coordinate_values_to_maxima): index_to_maxima = {} for i, coordinate_values in coordinates.df_all_coordinates.iterrows(): - coordinate_values = tuple([int(v) for v in coordinate_values]) + if len(coordinate_values) == 1: + coordinate_values = coordinate_values[0] + else: + coordinate_values = tuple([int(v) for v in coordinate_values]) index_to_maxima[i] = coordinate_values_to_maxima[coordinate_values] df = pd.DataFrame(index_to_maxima).transpose() df.index = coordinates.index diff --git a/test/test_extreme_trend/test_one_fold_fit.py b/test/test_extreme_trend/test_one_fold_fit.py index 80e3e7064192bf5d514e21acf2813a83b7e3e86f..49432249cfd67c68bc1cf01449f37e3eda2497c0 100644 --- a/test/test_extreme_trend/test_one_fold_fit.py +++ b/test/test_extreme_trend/test_one_fold_fit.py @@ -11,7 +11,7 @@ from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pari from extreme_fit.model.margin_model.polynomial_margin_model.utils import \ ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies -from extreme_trend.one_fold_fit.altitude_group import VeyHighAltitudeGroup +from extreme_trend.one_fold_fit.altitude_group import VeyHighAltitudeGroup, LowAltitudeGroup from extreme_trend.one_fold_fit.one_fold_fit import OneFoldFit from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \ TimeTemporalCovariate @@ -85,6 +85,7 @@ class TestOneFoldFit(unittest.TestCase): temporal_covariate_for_fit=None, only_models_that_pass_goodness_of_fit_test=False, remove_physically_implausible_models=True, + altitude_group=LowAltitudeGroup(), first_year=1959, last_year=2019 ) @@ -104,7 +105,7 @@ class TestOneFoldFit(unittest.TestCase): one_fold_fit = OneFoldFit(self.massif_name, dataset, models_classes=self.model_classes, temporal_covariate_for_fit=AnomalyTemperatureWithSplineTemporalCovariate, - altitude_class=VeyHighAltitudeGroup, + altitude_group=VeyHighAltitudeGroup(), only_models_that_pass_goodness_of_fit_test=False, remove_physically_implausible_models=True, first_year=1959,