From 843ef56d0ac7ea28131d8bcb37c60b55a90e3468 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 1 Mar 2021 16:38:06 +0100
Subject: [PATCH] [projection swe] add script for weight computation.
---
.../adamont_data/adamont_scenario.py | 4 ++-
.../abstract_margin_estimator.py | 26 ++++++++++++-------
.../ensemble_fit/abstract_ensemble_fit.py | 3 +++
.../independent_ensemble_fit.py | 19 +++++++++-----
.../visualizer_for_projection_ensemble.py | 7 +++--
.../visualizer_for_sensitivity.py | 4 +--
6 files changed, 42 insertions(+), 21 deletions(-)
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 604a3a12..e9003913 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 6c71b8c5..47aa7b4f 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 6e5d9e03..1a209af2 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 9fb05952..c4c3aaca 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 3149bab5..9d441551 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 824328c8..4b9993bf 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()
--
GitLab