[projection snowfall] modify altitude_class to altitude_group in...

[projection snowfall] modify altitude_class to altitude_group in visualizer_for_projection_ensemble.py. fix index for ensemble dataset.  add plot_relative_change_in_return_level.py

extreme_fit/estimator/margin_estimator/utils.py

@@ -14,12 +14,14 @@ from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima
 
 
-def fitted_linear_margin_estimator_short(model_class, dataset, fit_method, **model_kwargs) -> LinearMarginEstimator:
-    return fitted_linear_margin_estimator(model_class, dataset.coordinates, dataset, None, fit_method, **model_kwargs)
+def fitted_linear_margin_estimator_short(model_class, dataset, fit_method, drop_duplicates=None, **model_kwargs) -> LinearMarginEstimator:
+    return fitted_linear_margin_estimator(model_class, dataset.coordinates, dataset, None, fit_method, drop_duplicates, **model_kwargs)
 
 
-def fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year, fit_method, **model_kwargs):
+def fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year, fit_method, drop_duplicates=None, **model_kwargs):
     model = model_class(coordinates, starting_point=starting_year, fit_method=fit_method, **model_kwargs)
+    if drop_duplicates is not None:
+        model.drop_duplicates = drop_duplicates
     estimator = LinearMarginEstimator(dataset, model)
     estimator.fit()
     return estimator

extreme_trend/ensemble_fit/independent_ensemble_fit/one_fold_fit_merge.py

@@ -22,14 +22,16 @@ class OneFoldFitMerge(OneFoldFit):
     def get_moment(self, altitude, temporal_covariate, order=1):
         return self.merge_function([o.get_moment(altitude, temporal_covariate, order) for o in self.one_fold_fit_list])
 
-    def changes_of_moment(self, altitudes, order=1):
-        all_changes = [o.changes_of_moment(altitudes, order) for o in self.one_fold_fit_list]
+    def changes_of_moment(self, altitudes, order=1, covariate_before=None, covariate_after=None):
+        all_changes = [o.changes_of_moment(altitudes, order, covariate_before, covariate_after)
+                       for o in self.one_fold_fit_list]
         merged_changes = list(self.merge_function(np.array(all_changes), axis=0))
         assert len(all_changes[0]) == len(merged_changes)
         return merged_changes
 
-    def relative_changes_of_moment(self, altitudes, order=1):
-        all_relative_changes = [o.relative_changes_of_moment(altitudes, order) for o in self.one_fold_fit_list]
+    def relative_changes_of_moment(self, altitudes, order=1, covariate_before=None, covariate_after=None):
+        all_relative_changes = [o.relative_changes_of_moment(altitudes, order, covariate_before, covariate_after)
+                                for o in self.one_fold_fit_list]
         merged_relative_changes = list(self.merge_function(np.array(all_relative_changes), axis=0))
         assert len(all_relative_changes[0]) == len(merged_relative_changes)
         return merged_relative_changes

extreme_trend/ensemble_fit/together_ensemble_fit/visualizer_non_stationary_ensemble.py

@@ -39,8 +39,14 @@ class VisualizerNonStationaryEnsemble(AltitudesStudiesVisualizerForNonStationary
             dataset = studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
             df_coordinates_list.append(dataset.coordinates.df_coordinates(add_climate_informations=True))
             df_maxima_gev_list.append(dataset.observations.df_maxima_gev)
-        observations = AbstractSpatioTemporalObservations(df_maxima_gev=pd.concat(df_maxima_gev_list, axis=0))
-        coordinates = AbstractCoordinates(df=pd.concat(df_coordinates_list, axis=0),
+
+        index = pd.RangeIndex(0, sum([len(df) for df in df_maxima_gev_list]))
+        df_maxima_gev = pd.concat(df_maxima_gev_list, axis=0)
+        df_maxima_gev.index = index
+        observations = AbstractSpatioTemporalObservations(df_maxima_gev=df_maxima_gev)
+        df = pd.concat(df_coordinates_list, axis=0)
+        df.index = index
+        coordinates = AbstractCoordinates(df=df,
                                           slicer_class=type(dataset.slicer))
         dataset = AbstractDataset(observations=observations, coordinates=coordinates)
         return dataset

extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py

@@ -13,10 +13,13 @@ from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import Alt
 from extreme_trend.ensemble_fit.abstract_ensemble_fit import AbstractEnsembleFit
 from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
 from extreme_trend.ensemble_fit.together_ensemble_fit.together_ensemble_fit import TogetherEnsembleFit
-from extreme_trend.one_fold_fit.altitude_group import get_altitude_class_from_altitudes
+from extreme_trend.one_fold_fit.altitude_group import get_altitude_class_from_altitudes, \
+    get_altitude_group_from_altitudes
 from extreme_trend.one_fold_fit.plots.plot_histogram_altitude_studies import \
     plot_histogram_all_trends_against_altitudes, plot_shoe_plot_changes_against_altitude
 from extreme_trend.one_fold_fit.utils_altitude_studies_visualizer import compute_and_assign_max_abs
+from projects.projected_extreme_snowfall.results.plot_relative_change_in_return_level import \
+    plot_relative_dynamic_in_return_level
 
 
 class VisualizerForProjectionEnsemble(object):
@@ -48,9 +51,9 @@ class VisualizerForProjectionEnsemble(object):
                 assert len(years) == 2, years
 
         # Load all studies
-        altitude_class_to_gcm_couple_to_studies = OrderedDict()
+        altitude_group_to_gcm_couple_to_studies = OrderedDict()
         for altitudes in altitudes_list:
-            altitude_class = get_altitude_class_from_altitudes(altitudes)
+            altitude_group = get_altitude_group_from_altitudes(altitudes)
             gcm_rcm_couple_to_studies = {}
             for gcm_rcm_couple in gcm_rcm_couples:
                 if gcm_to_year_min_and_year_max is None:
@@ -69,11 +72,11 @@ class VisualizerForProjectionEnsemble(object):
                 gcm_rcm_couple_to_studies[gcm_rcm_couple] = studies
             if len(gcm_rcm_couple_to_studies) == 0:
                 print('No valid studies for the following couples:', self.gcm_rcm_couples)
-            altitude_class_to_gcm_couple_to_studies[altitude_class] = gcm_rcm_couple_to_studies
+            altitude_group_to_gcm_couple_to_studies[altitude_group] = gcm_rcm_couple_to_studies
 
         # Load 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():
+        self.altitude_group_to_ensemble_class_to_ensemble_fit = OrderedDict()
+        for altitude_group, gcm_rcm_couple_to_studies in altitude_group_to_gcm_couple_to_studies.items():
             ensemble_class_to_ensemble_fit = {}
             for ensemble_fit_class in ensemble_fit_classes:
                 ensemble_fit = ensemble_fit_class(massif_names, gcm_rcm_couple_to_studies, model_classes,
@@ -82,7 +85,7 @@ class VisualizerForProjectionEnsemble(object):
                                                   confidence_interval_based_on_delta_method,
                                                   remove_physically_implausible_models)
                 ensemble_class_to_ensemble_fit[ensemble_fit_class] = ensemble_fit
-            self.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class] = ensemble_class_to_ensemble_fit
+            self.altitude_group_to_ensemble_class_to_ensemble_fit[altitude_group] = ensemble_class_to_ensemble_fit
 
     @property
     def has_elevation_non_stationarity(self):
@@ -99,7 +102,7 @@ class VisualizerForProjectionEnsemble(object):
                 plot_shoe_plot_changes_against_altitude(self.massif_names, visualizer_list, relative=relative,
                                                         with_significance=with_significance)
         else:
-            print('here plot for one altitude visualizer')
+            plot_relative_dynamic_in_return_level(self.massif_names, visualizer_list)
 
     def plot(self):
         # Set limit for the plot
@@ -142,7 +145,7 @@ class VisualizerForProjectionEnsemble(object):
         """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()]
+                in self.altitude_group_to_ensemble_class_to_ensemble_fit.values()]
 
     def plot_preliminary_first_part(self):
         if self.massif_names is None:

extreme_trend/ensemble_fit/visualizer_for_sensitivity.py

@@ -235,12 +235,12 @@ class VisualizerForSensivity(object):
                              merge_visualizer_str):
         if merge_visualizer_str in [AbstractEnsembleFit.Median_merge, AbstractEnsembleFit.Mean_merge]:
             independent_ensemble_fit = \
-                visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][
+                visualizer_projection.altitude_group_to_ensemble_class_to_ensemble_fit[altitude_class][
                     IndependentEnsembleFit]
             merge_visualizer = independent_ensemble_fit.merge_function_name_to_visualizer[merge_visualizer_str]
         else:
             together_ensemble_fit = \
-                visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][
+                visualizer_projection.altitude_group_to_ensemble_class_to_ensemble_fit[altitude_class][
                     TogetherEnsembleFit]
             merge_visualizer = together_ensemble_fit.visualizer
         merge_visualizer.studies.study.gcm_rcm_couple = (merge_visualizer_str, "merge")

extreme_trend/one_fold_fit/one_fold_fit.py

@@ -77,7 +77,8 @@ class OneFoldFit(object):
         return fitted_linear_margin_estimator_short(model_class=model_class,
                                                     dataset=dataset,
                                                     fit_method=self.fit_method,
-                                                    temporal_covariate_for_fit=self.temporal_covariate_for_fit)
+                                                    temporal_covariate_for_fit=self.temporal_covariate_for_fit,
+                                                    drop_duplicates=False)
 
     @classmethod
     def get_moment_str(cls, order):
@@ -115,11 +116,15 @@ class OneFoldFit(object):
     def relative_change_in_return_level_for_reference_altitude(self) -> float:
         return self.relative_changes_of_moment(altitudes=[self.altitude_plot], order=None)[0]
 
-    def changes_of_moment(self, altitudes, order=1):
+    def changes_of_moment(self, altitudes, order=1, covariate_before=None, covariate_after=None):
+        if covariate_before is None:
+            covariate_before = self.covariate_before
+        if covariate_after is None:
+            covariate_after = self.covariate_after
         changes = []
         for altitude in altitudes:
-            mean_after = self.get_moment(altitude, self.covariate_after, order)
-            mean_before = self.get_moment(altitude, self.covariate_before, order)
+            mean_after = self.get_moment(altitude, covariate_after, order)
+            mean_before = self.get_moment(altitude, covariate_before, order)
             change = mean_after - mean_before
             changes.append(change)
         return changes
@@ -151,11 +156,11 @@ class OneFoldFit(object):
                      **d_temperature)
         return s
 
-    def relative_changes_of_moment(self, altitudes, order=1):
+    def relative_changes_of_moment(self, altitudes, order=1, covariate_before=None, covariate_after=None):
         relative_changes = []
         for altitude in altitudes:
-            mean_after = self.get_moment(altitude, self.covariate_after, order)
-            mean_before = self.get_moment(altitude, self.covariate_before, order)
+            mean_after = self.get_moment(altitude, covariate_after, order)
+            mean_before = self.get_moment(altitude, covariate_before, order)
             relative_change = 100 * (mean_after - mean_before) / mean_before
             relative_changes.append(relative_change)
         return relative_changes

projects/projected_extreme_snowfall/results/main_projections_ensemble.py

@@ -43,7 +43,7 @@ from extreme_data.meteo_france_data.scm_models_data.utils import Season
 def main():
     start = time.time()
     study_class = AdamontSnowfall
-    ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][:1]
+    ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][1:]
     temporal_covariate_for_fit = [TimeTemporalCovariate,
                                   AnomalyTemperatureWithSplineTemporalCovariate][1]
     set_seed_for_test()

projects/projected_extreme_snowfall/results/plot_relative_change_in_return_level.py

+from typing import List
+import matplotlib.pyplot as plt
+
+import numpy as np
+
+from extreme_trend.one_fold_fit.altitudes_studies_visualizer_for_non_stationary_models import \
+    AltitudesStudiesVisualizerForNonStationaryModels
+from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
+    AnomalyTemperatureWithSplineTemporalCovariate
+
+
+def plot_relative_dynamic_in_return_level(massif_names, visualizer_list: List[
+    AltitudesStudiesVisualizerForNonStationaryModels]):
+    ax = plt.gca()
+    print(len(visualizer_list))
+    visualizer = visualizer_list[0]
+    assert len(massif_names) == 1
+    assert visualizer.temporal_covariate_for_fit is AnomalyTemperatureWithSplineTemporalCovariate
+    massif_name = massif_names[0]
+    for v in visualizer_list:
+        plot_curve(ax, massif_name, v)
+
+    ax.set_xlabel('Anomaly of global temperature w.r.t. pre-industrial levels (K)')
+    ax.set_ylabel('Relative change in return levels (\%)')
+
+    ax.legend()
+    visualizer.plot_name = 'dynamic of return level'
+    visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
+
+    plt.close()
+
+
+def plot_curve(ax, massif_name, visualizer: AltitudesStudiesVisualizerForNonStationaryModels):
+    temperatures_list = np.linspace(1, 4, num=4)
+    one_fold_fit = visualizer.massif_name_to_one_fold_fit[massif_name]
+    return_levels = [one_fold_fit.relative_changes_of_moment([None], order=None,
+                                                             covariate_before=1,
+                                                             covariate_after=t)[0] for t in temperatures_list]
+    label = '{} m'.format(visualizer.altitude_group.reference_altitude)
+    ax.plot(temperatures_list, return_levels, label=label)
+
+