diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py b/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py index 604a3a1227a0b420e08f6de623f6f2dcf6cd46f6..e9003913572326e256b33e7ed1de648124ba861c 100644 --- a/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py +++ b/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py @@ -2,6 +2,8 @@ from enum import Enum from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import get_gcm_rcm_couple_adamont_to_full_name +SEPARATOR_STR = ' / ' + class AdamontScenario(Enum): histo = 0 @@ -136,4 +138,4 @@ def scenario_to_real_scenarios(adamont_scenario): def gcm_rcm_couple_to_str(gcm_rcm_couple): - return ' / '.join(gcm_rcm_couple) + return SEPARATOR_STR.join(gcm_rcm_couple) diff --git a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py index 6c71b8c53a0bef12871d86bf9c9da2cb31e8de24..47aa7b4f5e4978ea6fe65a1b1445f24ecde46e19 100644 --- a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py +++ b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py @@ -65,18 +65,9 @@ class LinearMarginEstimator(AbstractMarginEstimator): return pd.concat([self.df_coordinates_spat(split=split), self.df_coordinates_temp(split=split)], axis=1).values def nllh(self, split=Split.all): - nllh = 0 maxima_values = self.dataset.maxima_gev(split=split) coordinate_values = self.coordinates_for_nllh(split=split) - for maximum, coordinate in zip(maxima_values, coordinate_values): - assert len(maximum) == 1, \ - 'So far, only one observation for each coordinate, but code would be easy to change' - maximum = maximum[0] - gev_params = self.function_from_fit.get_params(coordinate, is_transformed=True) - p = gev_params.density(maximum) - nllh -= np.log(p) - assert not np.isinf(nllh) - return nllh + return compute_nllh(coordinate_values, maxima_values, self.function_from_fit) def sorted_empirical_standard_gumbel_quantiles(self, split=Split.all, coordinate_for_filter=None): sorted_empirical_quantiles = [] @@ -121,3 +112,18 @@ class LinearMarginEstimator(AbstractMarginEstimator): def bic(self, split=Split.all): return np.log(self.n(split=split)) * self.margin_model.nb_params + 2 * self.nllh(split=split) + + +def compute_nllh(coordinate_values, maxima_values, function_from_fit, maximum_from_obs=True, assertion_for_inf=True): + nllh = 0 + for maximum, coordinate in zip(maxima_values, coordinate_values): + if maximum_from_obs: + assert len(maximum) == 1, \ + 'So far, only one observation for each coordinate, but code would be easy to change' + maximum = maximum[0] + gev_params = function_from_fit.get_params(coordinate, is_transformed=True) + p = gev_params.density(maximum) + nllh -= np.log(p) + if assertion_for_inf: + assert not np.isinf(nllh), '{} {}'.format(gev_params, coordinate, maximum) + return nllh diff --git a/extreme_trend/ensemble_fit/abstract_ensemble_fit.py b/extreme_trend/ensemble_fit/abstract_ensemble_fit.py index 6e5d9e035dadef3707fd6db4306933d69697aafd..1a209af2e285f5a14c21a5730f1f704f7330acc6 100644 --- a/extreme_trend/ensemble_fit/abstract_ensemble_fit.py +++ b/extreme_trend/ensemble_fit/abstract_ensemble_fit.py @@ -23,3 +23,6 @@ class AbstractEnsembleFit(object): self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method + @property + def altitudes(self): + raise NotImplementedError \ No newline at end of file diff --git a/extreme_trend/ensemble_fit/independent_ensemble_fit/independent_ensemble_fit.py b/extreme_trend/ensemble_fit/independent_ensemble_fit/independent_ensemble_fit.py index 9fb05952e55029dafbad8b4e3322f10e9a94e0a9..c4c3aaca4989d7c45504ac6c3e63039d86bf7cae 100644 --- a/extreme_trend/ensemble_fit/independent_ensemble_fit/independent_ensemble_fit.py +++ b/extreme_trend/ensemble_fit/independent_ensemble_fit/independent_ensemble_fit.py @@ -37,12 +37,19 @@ class IndependentEnsembleFit(AbstractEnsembleFit): self.Mean_merge: np.mean } self.merge_function_name_to_visualizer = { - name: VisualizerMerge(visualizers, self.models_classes, False, self.massif_names, - self.fit_method, self.temporal_covariate_for_fit, - self.only_models_that_pass_goodness_of_fit_test, - self.confidence_interval_based_on_delta_method, - self.remove_physically_implausible_models, - merge_function=merge_function) + name: VisualizerMerge(visualizers, self.models_classes, False, self.massif_names, + self.fit_method, self.temporal_covariate_for_fit, + self.only_models_that_pass_goodness_of_fit_test, + self.confidence_interval_based_on_delta_method, + self.remove_physically_implausible_models, + merge_function=merge_function) for name, merge_function in merge_function_name_to_merge_function.items() } + @property + def visualizer(self): + return list(self.gcm_rcm_couple_to_visualizer.values())[0] + + @property + def altitudes(self): + return self.visualizer.studies.altitudes diff --git a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py index 3149bab5293927f4377120794cc3a352bf89fab9..9d441551205d1a26fb879e74daf526ddb7598a82 100644 --- a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py +++ b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py @@ -1,3 +1,4 @@ +from collections import OrderedDict from typing import List from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \ @@ -24,12 +25,13 @@ class VisualizerForProjectionEnsemble(object): remove_physically_implausible_models=False, gcm_to_year_min_and_year_max=None, ): + self.altitudes_list = altitudes_list self.gcm_rcm_couples = gcm_rcm_couples self.massif_names = massif_names self.ensemble_fit_classes = ensemble_fit_classes # Load all studies - altitude_class_to_gcm_couple_to_studies = {} + altitude_class_to_gcm_couple_to_studies = OrderedDict() for altitudes in altitudes_list: altitude_class = get_altitude_class_from_altitudes(altitudes) gcm_rcm_couple_to_studies = {} @@ -53,7 +55,7 @@ class VisualizerForProjectionEnsemble(object): altitude_class_to_gcm_couple_to_studies[altitude_class] = gcm_rcm_couple_to_studies # Load ensemble fit - self.altitude_class_to_ensemble_class_to_ensemble_fit = {} + self.altitude_class_to_ensemble_class_to_ensemble_fit = OrderedDict() for altitude_class, gcm_rcm_couple_to_studies in altitude_class_to_gcm_couple_to_studies.items(): ensemble_class_to_ensemble_fit = {} for ensemble_fit_class in ensemble_fit_classes: @@ -104,6 +106,7 @@ class VisualizerForProjectionEnsemble(object): plot_shoe_plot_changes_against_altitude(self.massif_names, visualizer_list, relative=relative, with_significance=with_significance) def ensemble_fits(self, ensemble_class): + """Return the ordered ensemble fit for a given ensemble class (in the order of the altitudes)""" return [ensemble_class_to_ensemble_fit[ensemble_class] for ensemble_class_to_ensemble_fit in self.altitude_class_to_ensemble_class_to_ensemble_fit.values()] diff --git a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py index 824328c8d61f25a921ba5ce493cd18abdabdeb52..4b9993bf3c7b74bf5b2df3c4bfdce58519d8fc27 100644 --- a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py +++ b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py @@ -35,7 +35,7 @@ class VisualizerForSensivity(object): self.temp_min_to_visualizer = {} # type: Dict[float, VisualizerForProjectionEnsemble] for temp_min, temp_max in zip(self.temp_min, self.temp_max): - print(temp_min, temp_max) + print("temp min and temp max", temp_min, temp_max) # Build gcm_to_year_min_and_year_max = {} gcm_list = list(set([g for g, r in gcm_rcm_couples])) @@ -60,7 +60,7 @@ class VisualizerForSensivity(object): def plot(self): # todo: before reactivating the subplot, i should ensure that we can modify the prefix - # so that we can have all the subplot for the merge visualizer + # todo: so that we can have all the subplot for the merge visualizer # , and not just the plots for the last t_min # for visualizer in self.temp_min_to_visualizer.values(): # visualizer.plot()