[projections] add a constraint to have physically plausible models for the fitted models.

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py

 from collections import Counter
 from math import ceil, floor
+from multiprocessing import Pool
 from typing import List, Dict
 
 import matplotlib
@@ -27,6 +28,7 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_gr
     get_altitude_group_from_altitudes, HighAltitudeGroup, VeyHighAltitudeGroup, MidAltitudeGroup
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import \
     OneFoldFit
+from root_utils import NB_CORES
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 
@@ -40,11 +42,12 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                  fit_method=MarginFitMethod.extremes_fevd_mle,
                  temporal_covariate_for_fit=None,
                  display_only_model_that_pass_anderson_test=True,
-                 confidence_interval_based_on_delta_method=False
+                 confidence_interval_based_on_delta_method=False,
+                 remove_physically_implausible_models=False,
                  ):
         super().__init__(studies.study, show=show, save_to_file=not show)
         self.studies = studies
-        self.non_stationary_models = model_classes
+        self.model_classes = model_classes
         self.fit_method = fit_method
         self.temporal_covariate_for_fit = temporal_covariate_for_fit
         self.display_only_model_that_pass_test = display_only_model_that_pass_anderson_test
@@ -52,23 +55,20 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         self.massif_name_to_massif_id = {m: i for i, m in enumerate(self.massif_names)}
         self.altitude_group = get_altitude_group_from_altitudes(self.studies.altitudes)
         self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
+        self.remove_physically_implausible_models = remove_physically_implausible_models
         # Load one fold fit
         self.massif_name_to_massif_altitudes = {}
-        self._massif_name_to_one_fold_fit = {}
-        for massif_name in self.massif_names:
-            # Load valid massif altitudes
-            massif_altitudes = self.get_massif_altitudes(massif_name)
-            if self.load_condition(massif_altitudes):
-                # Save the massif altitudes only for those who pass the condition
-                self.massif_name_to_massif_altitudes[massif_name] = massif_altitudes
-                # Load dataset
-                dataset = studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
-                old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method,
-                                          self.temporal_covariate_for_fit,
-                                          type(self.altitude_group),
-                                          self.display_only_model_that_pass_test,
-                                          self.confidence_interval_based_on_delta_method)
-                self._massif_name_to_one_fold_fit[massif_name] = old_fold_fit
+
+        # Multiprocess
+        multiprocessing = False
+        if multiprocessing:
+            with Pool(NB_CORES) as p:
+                one_fold_fit_list = p.map(self.fit_one_fold, self.massif_names)
+        else:
+            one_fold_fit_list = [self.fit_one_fold(massif_name) for massif_name in self.massif_names]
+        self._massif_name_to_one_fold_fit = {m: o for m, o in zip(self.massif_names, one_fold_fit_list) if
+                                             o is not None}
+
         # Print number of massif without any validated fit
         massifs_without_any_validated_fit = [massif_name
                                              for massif_name, old_fold_fit in self._massif_name_to_one_fold_fit.items()
@@ -79,6 +79,24 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         self._method_name_and_order_to_max_abs = {}
         self._max_abs_for_shape = None
 
+    def fit_one_fold(self, massif_name):
+        # Load valid massif altitudes
+        massif_altitudes = self.get_massif_altitudes(massif_name)
+        if self.load_condition(massif_altitudes):
+            # Save the massif altitudes only for those who pass the condition
+            self.massif_name_to_massif_altitudes[massif_name] = massif_altitudes
+            # Load dataset
+            dataset = self.studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
+            old_fold_fit = OneFoldFit(massif_name, dataset, self.model_classes, self.fit_method,
+                                      self.temporal_covariate_for_fit,
+                                      type(self.altitude_group),
+                                      self.display_only_model_that_pass_test,
+                                      self.confidence_interval_based_on_delta_method,
+                                      self.remove_physically_implausible_models)
+            return old_fold_fit
+        else:
+            return None
+
     moment_names = ['moment', 'changes_of_moment', 'relative_changes_of_moment'][:]
     orders = [1, 2, None][2:]
 
@@ -106,8 +124,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
     @property
     def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
         return {massif_name: old_fold_fit for massif_name, old_fold_fit in self._massif_name_to_one_fold_fit.items()
-                if not self.display_only_model_that_pass_test
-                or old_fold_fit.has_at_least_one_valid_model}
+                if old_fold_fit.has_at_least_one_valid_model}
 
     def plot_moments(self):
         for method_name in self.moment_names[:2]:
@@ -452,7 +469,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
             idx = 0 if one_fold.change_in_return_level_for_reference_altitude < 0 else 2
             nbs[idx] += 1
             if with_significance and one_fold.is_significant:
-                    nbs[idx + 1] += 1
+                nbs[idx + 1] += 1
 
         percents = 100 * nbs / nb_valid_massif_names
         return [nb_valid_massif_names] + list(percents)

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py

@@ -53,7 +53,9 @@ class OneFoldFit(object):
                  altitude_class=DefaultAltitudeGroup,
                  only_models_that_pass_goodness_of_fit_test=True,
                  confidence_interval_based_on_delta_method=False,
+                 remove_physically_implausible_models=False,
                  ):
+        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()
@@ -142,14 +144,21 @@ class OneFoldFit(object):
     @cached_property
     def sorted_estimators(self):
         estimators = list(self.model_class_to_estimator.values())
-        # print(self.massif_name)
-        # print(self.altitude_group)
-        # for estimator in estimators:
-        #     print(estimator.margin_model)
-        #     print(estimator.aic())
+        print(self.massif_name, self.altitude_group)
+        if self.remove_physically_implausible_models:
+            estimators = [e for e in estimators if -0.5 < self._compute_shape_for_reference_altitude(e) < 0.5]
+            if len(estimators) == 0:
+                print(self.massif_name, " has only implausible models")
+
         sorted_estimators = sorted([estimator for estimator in estimators], key=lambda e: e.aic())
         return sorted_estimators
 
+    def _compute_shape_for_reference_altitude(self, estimator):
+        coordinate = np.array([self.altitude_plot, self.last_year])
+        gev_params = estimator.function_from_fit.get_params(coordinate, is_transformed=False)
+        shape = gev_params.shape
+        return shape
+
     @cached_property
     def sorted_estimators_with_stationary(self):
         if self.only_models_that_pass_goodness_of_fit_test:
@@ -159,7 +168,8 @@ class OneFoldFit(object):
                     # and self.sensitivity_of_fit_test_last_years(e)
                     ]
         else:
-            assert len(self.sorted_estimators) == len(self.models_classes)
+            if not self.remove_physically_implausible_models:
+                assert len(self.sorted_estimators) == len(self.models_classes)
             return self.sorted_estimators
 
     @property
@@ -180,7 +190,8 @@ class OneFoldFit(object):
                 best_estimator = self.sorted_estimators_with_stationary[0]
                 return best_estimator
             else:
-                raise ValueError('This should not happen')
+                raise ValueError('This object should not have been called because '
+                                 'has_at_least_one_valid_model={}'.format(self.has_at_least_one_valid_model))
 
     @property
     def best_margin_model(self):

projects/altitude_spatial_model/preliminary_analysis.py

@@ -27,8 +27,6 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
                  **kwargs_study):
         super().__init__(study_class, altitudes, spatial_transformation_class, temporal_transformation_class,
                          **kwargs_study)
-        # self.altitudes_for_temporal_hypothesis = [min(self.altitudes), 2100, max(self.altitudes)]
-        self.altitudes_for_temporal_hypothesis = [600, 1500, 2400, 3300]
 
     def plot_gev_params_against_altitude(self):
         legend = False
@@ -332,8 +330,10 @@ def main_paper2():
 def main_paper3():
     altitudes = list(chain.from_iterable(altitudes_for_groups))
     # altitudes = [1200, 1500, 1800]
-    for scenario in rcp_scenarios[:]:
-        for gcm_rcm_couple in get_gcm_rcm_couples(scenario):
+    for scenario in rcp_scenarios[2:]:
+        gcm_rcm_couples = get_gcm_rcm_couples(scenario)
+        gcm_rcm_couples =[('CNRM-CM5', 'CCLM4-8-17')]
+        for gcm_rcm_couple in gcm_rcm_couples:
             visualizer = PointwiseGevStudyVisualizer(AdamontSnowfall, altitudes=altitudes, scenario=scenario,
                                                      gcm_rcm_couple=gcm_rcm_couple)
             visualizer.plot_gev_params_against_altitude()

projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py

@@ -13,7 +13,9 @@ class AbstractEnsembleFit(object):
                  temporal_covariate_for_fit=None,
                  only_models_that_pass_goodness_of_fit_test=True,
                  confidence_interval_based_on_delta_method=False,
+                 remove_physically_implausible_models=False,
                  ):
+        self.remove_physically_implausible_models = remove_physically_implausible_models
         self.massif_names = massif_names
         self.models_classes = models_classes
         self.gcm_rcm_couple_to_altitude_studies = gcm_rcm_couple_to_altitude_studies

projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py

@@ -44,7 +44,8 @@ class MetaVisualizerForProjectionEnsemble(object):
                  fit_method=MarginFitMethod.extremes_fevd_mle,
                  temporal_covariate_for_fit=None,
                  display_only_model_that_pass_gof_test=False,
-                 confidence_interval_based_on_delta_method=False
+                 confidence_interval_based_on_delta_method=False,
+                 remove_physically_implausible_models=False,
                  ):
         self.gcm_rcm_couples = gcm_rcm_couples
         self.massif_names = massif_names
@@ -69,7 +70,8 @@ class MetaVisualizerForProjectionEnsemble(object):
                 ensemble_fit = ensemble_fit_class(massif_names, gcm_rcm_couple_to_studies, model_classes,
                                                   fit_method, temporal_covariate_for_fit,
                                                   display_only_model_that_pass_gof_test,
-                                                  confidence_interval_based_on_delta_method)
+                                                  confidence_interval_based_on_delta_method,
+                                                  remove_physically_implausible_models)
                 ensemble_class_to_ensemble_fit[ensemble_fit_class] = ensemble_fit
             self.altitude_group_to_ensemble_class_to_ensemble_fit[altitude_group] = ensemble_class_to_ensemble_fit