[projections] first version for abstract_ensemble_fit.py and independent_ensemble_fit.py

Try to reuse the previous code as much as possible.

extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py

@@ -13,6 +13,8 @@ import numpy as np
 import pandas as pd
 import seaborn as sns
 
+from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import AbstractAdamontStudy
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_str
 from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import load_plot
 from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
 from extreme_fit.model.margin_model.utils import fitmethod_to_str
@@ -541,6 +543,10 @@ class StudyVisualizer(VisualizationParameters):
 
     def show_or_save_to_file(self, add_classic_title=True, no_title=False, tight_layout=False, tight_pad=None,
                              dpi=None, folder_for_variable=True):
+        if isinstance(self.study, AbstractAdamontStudy):
+            prefix = gcm_rcm_couple_to_str(self.study.gcm_rcm_couple)
+            self.plot_name = prefix + ' ' + self.plot_name
+
         assert self.plot_name is not None
         if add_classic_title:
             title = self.study.title

projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py

@@ -5,7 +5,7 @@ from collections import OrderedDict
 from cached_property import cached_property
 
 from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import AbstractAdamontStudy
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import scenario_to_str
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import scenario_to_str, gcm_rcm_couple_to_str
 from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
     SCM_STUDY_CLASS_TO_ABBREVIATION, STUDY_CLASS_TO_ABBREVIATION

projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py

+from typing import Dict, Tuple, List
+
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import DefaultAltitudeGroup
 
 
 class AbstractEnsembleFit(object):
-    pass
+
\ No newline at end of file
+    def __init__(self, massif_names, gcm_rcm_couple_to_altitude_studies: Dict[Tuple[str, str], AltitudesStudies],
+                 models_classes,
+                 fit_method=MarginFitMethod.extremes_fevd_mle,
+                 temporal_covariate_for_fit=None,
+                 only_models_that_pass_goodness_of_fit_test=True,
+                 confidence_interval_based_on_delta_method=False,
+                 ):
+        self.massif_names = massif_names
+        self.models_classes = models_classes
+        self.gcm_rcm_couple_to_altitude_studies = gcm_rcm_couple_to_altitude_studies
+        self.fit_method = fit_method
+        self.temporal_covariate_for_fit = temporal_covariate_for_fit
+        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
+
+
+    def plot(self):
+        raise NotImplementedError
\ No newline at end of file

projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/independent_ensemble_fit.py

+from typing import Dict, Tuple, List
+
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import DefaultAltitudeGroup
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
+    AltitudesStudiesVisualizerForNonStationaryModels
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
+from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import compute_and_assign_max_abs
+from projects.projected_snowfall.elevation_temporal_model_for_projections.ensemble_fit.abstract_ensemble_fit import \
+    AbstractEnsembleFit
+
+
+class IndependentEnsembleFit(AbstractEnsembleFit):
+    """For each gcm_rcm_couple, we create a OneFoldFit"""
+
+    def __init__(self, massif_names, gcm_rcm_couple_to_altitude_studies: Dict[Tuple[str, str], AltitudesStudies], models_classes,
+                 fit_method=MarginFitMethod.extremes_fevd_mle, temporal_covariate_for_fit=None, only_models_that_pass_goodness_of_fit_test=True,
+                 confidence_interval_based_on_delta_method=False):
+        super().__init__(massif_names, gcm_rcm_couple_to_altitude_studies, models_classes, fit_method, temporal_covariate_for_fit, only_models_that_pass_goodness_of_fit_test,
+                         confidence_interval_based_on_delta_method)
+
+        # Load a classical visualizer
+        self.gcm_rcm_couple_to_visualizer = {}
+        for gcm_rcm_couple, studies in gcm_rcm_couple_to_altitude_studies.items():
+            visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies, 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.gcm_rcm_couple_to_visualizer[gcm_rcm_couple] = visualizer
+
+        # Assign max
+        visualizer_list = list(self.gcm_rcm_couple_to_visualizer.values())
+        compute_and_assign_max_abs(visualizer_list)
+
+    def plot(self):
+        for v in self.gcm_rcm_couple_to_visualizer.values():
+            self.plot_one_visualizer(v)
+
+    @staticmethod
+    def plot_one_visualizer(visualizer):
+        # visualizer.studies.plot_maxima_time_series(['Vanoise'])
+        visualizer.plot_shape_map()
+        visualizer.plot_moments()
+        # visualizer.plot_qqplots()
+

projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py

@@ -2,12 +2,25 @@ import datetime
 import time
 from typing import List
 
+import matplotlib as mpl
+
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
 import matplotlib
+from extreme_fit.model.utils import set_seed_for_test
 
+from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples
+from projects.projected_snowfall.elevation_temporal_model_for_projections.ensemble_fit.independent_ensemble_fit import \
+    IndependentEnsembleFit
 from projects.projected_snowfall.elevation_temporal_model_for_projections.utils_projected_visualizer import \
     load_projected_visualizer_list
 from projects.projected_snowfall.elevation_temporal_model_for_projections.visualizer_for_projection_ensemble import \
     VisualizerForProjectionEnsemble
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
+    AnomalyTemperatureTemporalCovariate
 
 matplotlib.use('Agg')
 
@@ -17,96 +30,60 @@ from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
     AbstractExtractEurocodeReturnLevel
 
-import matplotlib as mpl
 
-from extreme_fit.model.utils import set_seed_for_test
-
-mpl.rcParams['text.usetex'] = True
-mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups
 
-from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
-    SafranSnowfall5Days, SafranSnowfall7Days
 from extreme_data.meteo_france_data.scm_models_data.utils import Season
 
 
 def main():
-    study_classes = [SafranSnowfall1Day
-                        , SafranSnowfall3Days,
-                     SafranSnowfall5Days, SafranSnowfall7Days][:1]
-    seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
-
+    study_classes = [AdamontSnowfall][:1]
+    scenario = AdamontScenario.rcp85
+    gcm_rcm_couples = get_gcm_rcm_couples(scenario)[:2]
+    ensemble_fit_class = [IndependentEnsembleFit]
+    temporal_covariate_for_fit = [None, AnomalyTemperatureTemporalCovariate][0]
     set_seed_for_test()
-    model_must_pass_the_test = False
     AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
 
-    fast = False
+    fast = True
     if fast is None:
         massif_names = None
         AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
         altitudes_list = altitudes_for_groups[2:3]
     elif fast:
         AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
-        massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
+        massif_names = ['Vanoise'][:]
         altitudes_list = altitudes_for_groups[1:2]
     else:
         massif_names = None
         altitudes_list = altitudes_for_groups[:]
 
     start = time.time()
-    main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_pass_the_test)
+    main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_classes, ensemble_fit_class, scenario, temporal_covariate_for_fit)
     end = time.time()
     duration = str(datetime.timedelta(seconds=end - start))
     print('Total duration', duration)
 
 
-def main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_pass_the_test):
+def main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_classes, ensemble_fit_classes, scenario, temporal_covariate_for_fit):
     assert isinstance(altitudes_list, List)
     assert isinstance(altitudes_list[0], List)
-    for season in seasons:
-        for study_class in study_classes:
-            print('Inner loop', season, study_class)
-            visualizer_list = load_projected_visualizer_list(season, study_class, altitudes_list, massif_names,
-                                                             model_must_pass_the_test
-                                                             )
-            plot_visualizers(massif_names, visualizer_list)
-            for visualizer in visualizer_list:
-                plot_visualizer(massif_names, visualizer)
-            del visualizer_list
-            time.sleep(2)
-
-
-def plot_visualizers(massif_names, visualizer_list):
-    # plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
-    plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
-    # plot_shoe_plot_ratio_interval_size_against_altitude(massif_names, visualizer_list)
-    for relative in [True, False]:
-        plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative)
-    # plot_coherence_curves(massif_names, visualizer_list)
-    # plot_coherence_curves(['Vanoise'], visualizer_list)
-    pass
-
-
-def plot_visualizer(massif_names, visualizer):
-    # Plot time series
-    # visualizer.studies.plot_maxima_time_series(massif_names)
-    # visualizer.studies.plot_maxima_time_series(['Vanoise'])
-
-    # Plot the results for the model that minimizes the individual aic
-    plots(visualizer)
-
-    # Plot the results for the model that minimizes the total aic
-    # plot_total_aic(model_classes, visualizer)
-    pass
-
-
-def plots(visualizer: VisualizerForProjectionEnsemble):
-    # visualizer.plot_shape_map()
-    visualizer.plot_moments()
-    # visualizer.plot_qqplots()
-    # for std in [True, False]:
-    #     visualizer.studies.plot_mean_maxima_against_altitude(std=std)
+    for study_class in study_classes:
+        print('Inner loop', study_class)
+        visualizer_list = load_projected_visualizer_list(gcm_rcm_couples=gcm_rcm_couples, ensemble_fit_classes=ensemble_fit_classes,
+                                                         season=Season.annual, study_class=study_class,
+                                                         altitudes_list=altitudes_list, massif_names=massif_names,
+                                                         scenario=scenario,
+                                                         temporal_covariate_for_fit=temporal_covariate_for_fit)
+        for v in visualizer_list:
+            v.plot()
+
+        del visualizer_list
+        time.sleep(2)
+
+
+
 
 
 if __name__ == '__main__':

projects/projected_snowfall/elevation_temporal_model_for_projections/utils_projected_visualizer.py

+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, rcp_scenarios
 from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
     ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
 from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
     AltitudesStudiesVisualizerForNonStationaryModels
+from projects.projected_snowfall.elevation_temporal_model_for_projections.visualizer_for_projection_ensemble import \
+    VisualizerForProjectionEnsemble
 
 
-def load_projected_visualizer_list(season, study_class, altitudes_list, massif_names, model_must_pass_the_test=True, **kwargs_study):
+def load_projected_visualizer_list(gcm_rcm_couples, ensemble_fit_classes,
+                                   season, study_class, altitudes_list, massif_names, model_must_pass_the_test=False,
+                                   scenario=AdamontScenario.rcp85,
+                                   temporal_covariate_for_fit=None, **kwargs_study):
     model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+    assert scenario in rcp_scenarios
     visualizer_list = []
     # Load all studies
     for altitudes in altitudes_list:
-        studies = AltitudesStudies(study_class, altitudes, season=season, **kwargs_study)
-        visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
-                                                                      model_classes=model_classes,
-                                                                      massif_names=massif_names,
-                                                                      show=False,
-                                                                      temporal_covariate_for_fit=None,
-                                                                      confidence_interval_based_on_delta_method=False,
-                                                                      display_only_model_that_pass_anderson_test=model_must_pass_the_test
-                                                                      )
+        gcm_rcm_couple_to_altitude_studies = {}
+        for gcm_rcm_couple in gcm_rcm_couples:
+            studies = AltitudesStudies(study_class, altitudes, season=season,
+                                       scenario=scenario, gcm_rcm_couple=gcm_rcm_couple, **kwargs_study)
+            gcm_rcm_couple_to_altitude_studies[gcm_rcm_couple] = studies
+        visualizer = VisualizerForProjectionEnsemble(gcm_rcm_couple_to_altitude_studies=gcm_rcm_couple_to_altitude_studies,
+                                                     model_classes=model_classes,
+                                                     ensemble_fit_classes=ensemble_fit_classes,
+                                                     massif_names=massif_names,
+                                                     show=False,
+                                                     temporal_covariate_for_fit=temporal_covariate_for_fit,
+                                                     confidence_interval_based_on_delta_method=False,
+                                                     display_only_model_that_pass_gof_test=model_must_pass_the_test
+                                                     )
         visualizer_list.append(visualizer)
-    compute_and_assign_max_abs(visualizer_list)
 
     return visualizer_list
 
 
-def compute_and_assign_max_abs(visualizer_list):
-    # Compute the max abs for all metrics
-    d = {}
-    for method_name in AltitudesStudiesVisualizerForNonStationaryModels.moment_names:
-        for order in AltitudesStudiesVisualizerForNonStationaryModels.orders:
-            c = (method_name, order)
-            max_abs = max([
-                max([abs(e) for e in v.method_name_and_order_to_d(method_name, order).values()
-                     ]) for v in visualizer_list])
-            d[c] = max_abs
-    # Assign the max abs dictionary
-    for v in visualizer_list:
-        v._method_name_and_order_to_max_abs = d
-    # Compute the max abs for the shape parameter
-    max_abs_for_shape = max([max([abs(e) for e in v.massif_name_to_shape.values()]) for v in visualizer_list])
-    for v in visualizer_list:
-        v._max_abs_for_shape = max_abs_for_shape

projects/projected_snowfall/trends/main_projected_trends_separately.py

-import time
-from typing import List
-
-import matplotlib as mpl
-import numpy as np
-
-from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import WrongYearMinOrYearMax
-from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario
-from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days
-from projects.altitude_spatial_model.altitudes_fit.plots.plot_coherence_curves import plot_coherence_curves
-from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
-    plot_histogram_all_models_against_altitudes, plot_histogram_all_trends_against_altitudes, \
-    plot_shoe_plot_changes_against_altitude
-
-mpl.rcParams['text.usetex'] = True
-mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
-
-from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import load_visualizer_list
-
-from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups
-from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_individual_aic
-
-from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
-    SafranSnowfall5Days, SafranSnowfall7Days, SafranDateFirstSnowfall, SafranPrecipitation1Day, SafranPrecipitation3Days
-from extreme_data.meteo_france_data.scm_models_data.utils import Season
-
-
-def main():
-    study_class = AdamontSnowfall
-    scenario = AdamontScenario.rcp85_extended
-    gcm_rcm_couples = [('CNRM-CM5', 'CCLM4-8-17'), ('CNRM-CM5', 'RCA4')][1:]
-
-    fast = True
-    if fast is None:
-        massif_names = None
-        altitudes_list = altitudes_for_groups[1:2]
-    elif fast:
-        massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
-        altitudes_list = altitudes_for_groups[1:3]
-    else:
-        massif_names = None
-        altitudes_list = altitudes_for_groups[:]
-
-    # One plot for each massif name
-    for massif_name in massif_names:
-        print(massif_name)
-        main_loop(altitudes_list, [massif_name], gcm_rcm_couples, study_class, scenario)
-
-
-def main_loop(altitudes_list, massif_names, gcm_rcm_couples, study_class, scenario):
-    assert isinstance(altitudes_list, List)
-    assert isinstance(altitudes_list[0], List)
-    altitudes_for_return_levels = [altitudes[-2] for altitudes in altitudes_list]
-    print(altitudes_for_return_levels)
-    season = Season.annual
-    first_year_min, first_year_max = 1951, 2010
-    nb_years = first_year_max - first_year_min + 1
-    temporal_windows = [(first_year_min + i, first_year_max + i) for i in [30 * j for j in range(4)]]
-    all_changes_in_return_levels = []
-    for gcm_rcm_couple in gcm_rcm_couples:
-        print('Inner', gcm_rcm_couple)
-        changes_in_return_levels = []
-        for temporal_window in temporal_windows:
-            year_min, year_max = temporal_window
-            try:
-                visualizer_list = load_visualizer_list(season, study_class, altitudes_list, massif_names,
-                                                       scenario=scenario, gcm_rcm_couple=gcm_rcm_couple,
-                                                       year_min=year_min, year_max=year_max)
-
-                for visualizer in visualizer_list:
-                    visualizer.studies.plot_maxima_time_series(massif_names)
-                massif_name = massif_names[0]
-                changes_for_temporal_window = [
-                    v.massif_name_to_one_fold_fit[massif_name].changes_of_moment(altitudes=[a],
-                                                                                 year=year_max,
-                                                                                 nb_years=nb_years,
-                                                                                 order=None
-                                                                                 )[0]
-                    for a, v in zip(altitudes_for_return_levels, visualizer_list)]
-
-            except WrongYearMinOrYearMax:
-                changes_for_temporal_window = [np.nan for _ in altitudes_for_return_levels]
-
-            print(changes_for_temporal_window)
-            changes_in_return_levels.append(changes_for_temporal_window)
-        changes_in_return_levels = np.array(changes_in_return_levels)
-        all_changes_in_return_levels.append(changes_in_return_levels)
-    all_changes_in_return_levels = np.array(all_changes_in_return_levels)
-
-    return all_changes_in_return_levels, temporal_windows, altitudes_for_return_levels, nb_years
-
-
-if __name__ == '__main__':
-    main()