[projection snowfall] add aicc. add potential framework for model as truth experiments.

extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py

@@ -124,6 +124,9 @@ class LinearMarginEstimator(AbstractMarginEstimator):
     def bic(self, split=Split.all):
         return np.log(self.n(split=split)) * self.nb_params + 2 * self.nllh(split=split)
 
+    def aicc(self, split=Split.all):
+        additional_term = 2 * self.nb_params * (self.nb_params + 1) / (self.n() - self.nb_params - 1)
+        return self.aic(split) + additional_term
 
 def compute_nllh(coordinate_values, maxima_values, function_from_fit, maximum_from_obs=True, assertion_for_inf=True):
     nllh = 0

extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py

@@ -518,6 +518,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
 
             for altitude in self.studies.altitudes:
                 coordinate_for_filter = (altitude, None)
+                # We filter on the transformed gumbel quantiles for the altitude of interest
                 unconstrained_empirical_quantiles = one_fold_fit.best_estimator.sorted_empirical_standard_gumbel_quantiles(
                     coordinate_for_filter=coordinate_for_filter)
                 n = len(unconstrained_empirical_quantiles)

projects/projected_extreme_snowfall/results/part_1/main_model_as_truth_experiment.py

+import datetime
+import time
+
+from extreme_data.meteo_france_data.scm_models_data.utils import Season
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from projects.projected_extreme_snowfall.results.part_1.model_as_truth_experiment import ModelAsTruthExperiment
+from projects.projected_extreme_snowfall.results.part_3.main_projections_ensemble import set_up_and_load
+
+
+def main_model_as_truth_experiment():
+    start = time.time()
+    fast = True
+
+    altitudes_list, climate_coordinates_with_effects_list, gcm_rcm_couples, massif_names, model_classes, scenario, \
+    study_class, temporal_covariate_for_fit = set_up_and_load(fast)
+
+    for altitudes in altitudes_list[:1]:
+        for climate_coordinates_with_effects in climate_coordinates_with_effects_list:
+            xp = ModelAsTruthExperiment(altitudes, gcm_rcm_couples, study_class, Season.annual,
+                                        scenario=scenario,
+                                        model_classes=model_classes,
+                                        massif_names=massif_names,
+                                        fit_method=MarginFitMethod.evgam,
+                                        temporal_covariate_for_fit=temporal_covariate_for_fit,
+                                        remove_physically_implausible_models=True,
+                                        display_only_model_that_pass_gof_test=True,
+                                        climate_coordinates_with_effects=climate_coordinates_with_effects)
+
+    end = time.time()
+    duration = str(datetime.timedelta(seconds=end - start))
+    print('Total duration', duration)
+
+
+if __name__ == '__main__':
+    main_model_as_truth_experiment()

projects/projected_extreme_snowfall/results/part_1/model_as_truth_experiment.py

+from typing import List
+
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+    AbstractTemporalLinearMarginModel
+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.ensemble_fit.together_ensemble_fit.visualizer_non_stationary_ensemble import \
+    VisualizerNonStationaryEnsemble
+
+
+class ModelAsTruthExperiment(object):
+
+    def __init__(self, altitudes, gcm_rcm_couples, study_class, season, scenario,
+                 model_classes: List[AbstractTemporalLinearMarginModel],
+                 massif_names=None,
+                 fit_method=MarginFitMethod.extremes_fevd_mle,
+                 temporal_covariate_for_fit=None,
+                 display_only_model_that_pass_gof_test=False,
+                 remove_physically_implausible_models=False,
+                 climate_coordinates_with_effects=None,
+                 ):
+        self.fit_method = fit_method
+        self.massif_names = massif_names
+        self.temporal_covariate_for_fit = temporal_covariate_for_fit
+        self.display_only_model_that_pass_gof_test = display_only_model_that_pass_gof_test
+        self.remove_physically_implausible_models = remove_physically_implausible_models
+        self.climate_coordinates_with_effects = climate_coordinates_with_effects
+        self.model_classes = model_classes
+        self.scenario = scenario
+        self.season = season
+        self.study_class = study_class
+        self.gcm_rcm_couples = gcm_rcm_couples
+        self.altitudes = altitudes
+
+    def run_all_experiments(self):
+        pass
+
+    def run_one_experiment(self, gcm_rcm_couple_set_as_truth):
+        # Load gcm_rcm_couple_to_studies
+        gcm_rcm_couple_to_studies = self.load_gcm_rcm_couple_to_studies(gcm_rcm_couple_set_as_truth)
+        # Load ensemble fit
+        visualizer = VisualizerNonStationaryEnsemble(gcm_rcm_couple_to_studies=gcm_rcm_couple_to_studies,
+                                                   model_classes=self.model_classes,
+                                                       show=False,
+                                                       massif_names=self.massif_names,
+                                                  fit_method=self.fit_method,
+                                                   temporal_covariate_for_fit=self.temporal_covariate_for_fit,
+                                                  display_only_model_that_pass_gof_test=self.display_only_model_that_pass_gof_test,
+                                                  confidence_interval_based_on_delta_method=False,
+                                                  remove_physically_implausible_models=self.remove_physically_implausible_models,
+                                                  climate_coordinates_with_effects=self.climate_coordinates_with_effects)
+        # Compute the average nllh for the test data
+        studies = self.load_gcm_rcm_couple_to_studies(gcm_rcm_couple_set_as_truth)
+
+    def load_studies_for_test(self, gcm_rcm_couple_set_as_truth):
+        """For gcm_rcm_couple_set_as_truth, load the data from 2020 to 2100"""
+        pass
+
+    def load_gcm_rcm_couple_to_studies(self, gcm_rcm_couple_set_as_truth):
+        """For the gcm_rcm_couple_set_as_truth load only the data from 1959 to 2019"""
+        pass
+        # kwargs_study = {'year_min': year_min, 'year_max': year_max}
+        #
+        # studies = AltitudesStudies(study_class, altitudes, season=season,
+        #                            scenario=scenario, gcm_rcm_couple=gcm_rcm_couple,
+        #                            **kwargs_study)
+        # gcm_rcm_couple_to_studies[gcm_rcm_couple] = studies
+        # # Potentially add the observations
+        #
+        #
+        # if self.safran_study_class is not None:
+        #     studies = AltitudesStudies(self.safran_study_class, altitudes, season=season)
+        #     gcm_rcm_couple_to_studies[(None, None)] = studies
\ No newline at end of file

projects/projected_extreme_snowfall/results/part_3/main_projections_ensemble.py

@@ -47,81 +47,85 @@ from extreme_trend.one_fold_fit.altitude_group import altitudes_for_groups
 from extreme_data.meteo_france_data.scm_models_data.utils import Season
 
 
-def main():
-    start = time.time()
+def set_up_and_load(fast):
     study_class = AdamontSnowfall
-    safran_study_class = [None, SafranSnowfall2019][1]  # None means we do not account for the observations
     climate_coordinates_with_effects_list = [None,
                                              [AbstractCoordinates.COORDINATE_GCM, AbstractCoordinates.COORDINATE_RCM],
                                              [AbstractCoordinates.COORDINATE_GCM],
                                              [AbstractCoordinates.COORDINATE_RCM],
-                                        ][1:2]  # None means we do not create any effect
-    ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][1:]
+                                             ][1:2]  # None means we do not create any effect
     temporal_covariate_for_fit = [TimeTemporalCovariate,
                                   AnomalyTemperatureWithSplineTemporalCovariate][1]
     set_seed_for_test()
-    AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
-    scenarios = [AdamontScenario.rcp85_extended]
+    scenario = AdamontScenario.rcp85_extended
     model_classes = SPLINE_MODELS_FOR_PROJECTION_ONE_ALTITUDE
+    AltitudesStudiesVisualizerForNonStationaryModels.consider_at_least_two_altitudes = False
+    gcm_rcm_couples = get_gcm_rcm_couples(scenario)
+    print('Scenario is', scenario)
+    print('Covariate is {}'.format(temporal_covariate_for_fit))
+    if fast is None:
+        gcm_rcm_couples = gcm_rcm_couples[:]
+        AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
+        altitudes_list = [600, 2100, 3600]
+    elif fast:
+        gcm_rcm_couples = gcm_rcm_couples[:3]
+        AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
+        altitudes_list = [2700, 3000]
+        model_classes = model_classes[:4]
+    else:
+        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
+    massif_names = ['Vanoise']
+    if fast in [None, False]:
+        assert len(set(model_classes)) == 27
+    print('number of models', len(model_classes))
+    return altitudes_list, climate_coordinates_with_effects_list, gcm_rcm_couples, massif_names, model_classes, scenario, study_class, temporal_covariate_for_fit
+
 
+def main():
+    start = time.time()
     fast = None
-    for scenario in scenarios:
-        gcm_rcm_couples = get_gcm_rcm_couples(scenario)
-        if fast is None:
-            gcm_rcm_couples = gcm_rcm_couples[:]
-            AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
-            altitudes_list = [600, 2100, 3600]
-        elif fast:
-            gcm_rcm_couples = gcm_rcm_couples[:3]
-            AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
-            altitudes_list = [2700, 3000]
-            model_classes = model_classes[:4]
-        else:
-            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))
-        print('Take into account the observations: {}'.format(safran_study_class is not None))
-        print('observation class:', get_display_name_from_object_type(safran_study_class))
-
-        # Default parameters
-        gcm_to_year_min_and_year_max = None
-        massif_names = ['Vanoise']
-        if fast in [None, False]:
-            assert len(set(model_classes)) == 27
-        print('number of models', len(model_classes))
-
-        for climate_coordinates_with_effects in climate_coordinates_with_effects_list:
-            print('climate coordinates with effects ', climate_coordinates_with_effects)
-
-            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,
-                safran_study_class=safran_study_class,
-                display_only_model_that_pass_gof_test=True,
-                climate_coordinates_with_effects=climate_coordinates_with_effects,
-            )
-            visualizer.plot()
+
+    altitudes_list, climate_coordinates_with_effects_list, gcm_rcm_couples, massif_names, model_classes, scenario, \
+    study_class, temporal_covariate_for_fit = set_up_and_load(fast)
+
+    # Default parameters
+    gcm_to_year_min_and_year_max = None
+    safran_study_class = [None, SafranSnowfall2019][1]  # None means we do not account for the observations
+    print('Take into account the observations: {}'.format(safran_study_class is not None))
+    print('observation class:', get_display_name_from_object_type(safran_study_class))
+    ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][1:]
+    AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
+
+    for climate_coordinates_with_effects in climate_coordinates_with_effects_list:
+        print('climate coordinates with effects ', climate_coordinates_with_effects)
+
+        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,
+            safran_study_class=safran_study_class,
+            display_only_model_that_pass_gof_test=True,
+            climate_coordinates_with_effects=climate_coordinates_with_effects,
+        )
+        visualizer.plot()
 
     end = time.time()
     duration = str(datetime.timedelta(seconds=end - start))
     print('Total duration', duration)
 
 
+
+
 if __name__ == '__main__':
     main()