[projection swe] add script for weight computation.

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)

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

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

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

extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py

+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()]

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()