From 2e6731a7d546d715c2a9b4c5d5a76889b7fcd986 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 22 Apr 2021 15:49:01 +0200
Subject: [PATCH] [projection snowfall] refactoring & move folders. adapt the
 relative change formula for negative values.activate the goodness of fit
 test. handle an additional type of climate coordinates. and plot of the value
 & of the GEV parameters.

---
 .../independent_margin_function.py            |  21 ++--
 .../visualizer_for_projection_ensemble.py     |  25 +++--
 extreme_trend/one_fold_fit/one_fold_fit.py    |  12 ++-
 .../__init__.py                               |   0
 .../comparison_window_for_the_max.py          |   0
 .../checks}/__init__.py                       |   0
 .../checks/main_inquiry_weird_pattern.py      |   0
 .../{ => data}/evaluation/__init__.py         |   0
 .../comparison_historical_visualizer.py       |   0
 .../{ => data}/evaluation/comparison_plot.py  |   6 +-
 .../main_comparison_on_quantile_period.py     |   0
 .../evaluation/main_comparison_reanalysis.py  |   0
 .../data/main_data.py                         |   5 +-
 .../part_1}/__init__.py                       |   0
 .../part_1}/main_preliminary_works.py         |   0
 .../results/part_2/__init__.py                |   0
 .../results/part_3/__init__.py                |   0
 .../{ => part_3}/main_projections_ensemble.py |  15 +--
 .../{ => part_3}/plot_gcm_rcm_effects.py      |  43 +++++---
 .../plot_relative_change_in_return_level.py   | 100 ++++++++++++++++++
 .../plot_relative_change_in_return_level.py   |  55 ----------
 .../coordinates/abstract_coordinates.py       |  15 ++-
 22 files changed, 180 insertions(+), 117 deletions(-)
 rename projects/{projected_extreme_snowfall/checks => archive/comparison_window_for_the_max}/__init__.py (100%)
 rename projects/{projected_extreme_snowfall => archive}/comparison_window_for_the_max/comparison_window_for_the_max.py (100%)
 rename projects/projected_extreme_snowfall/{comparison_window_for_the_max => data/checks}/__init__.py (100%)
 rename projects/projected_extreme_snowfall/{ => data}/checks/main_inquiry_weird_pattern.py (100%)
 rename projects/projected_extreme_snowfall/{ => data}/evaluation/__init__.py (100%)
 rename projects/projected_extreme_snowfall/{ => data}/evaluation/comparison_historical_visualizer.py (100%)
 rename projects/projected_extreme_snowfall/{ => data}/evaluation/comparison_plot.py (96%)
 rename projects/projected_extreme_snowfall/{ => data}/evaluation/main_comparison_on_quantile_period.py (100%)
 rename projects/projected_extreme_snowfall/{ => data}/evaluation/main_comparison_reanalysis.py (100%)
 rename projects/projected_extreme_snowfall/{preliminary_works => results/part_1}/__init__.py (100%)
 rename projects/projected_extreme_snowfall/{preliminary_works => results/part_1}/main_preliminary_works.py (100%)
 create mode 100644 projects/projected_extreme_snowfall/results/part_2/__init__.py
 create mode 100644 projects/projected_extreme_snowfall/results/part_3/__init__.py
 rename projects/projected_extreme_snowfall/results/{ => part_3}/main_projections_ensemble.py (93%)
 rename projects/projected_extreme_snowfall/results/{ => part_3}/plot_gcm_rcm_effects.py (65%)
 create mode 100644 projects/projected_extreme_snowfall/results/part_3/plot_relative_change_in_return_level.py
 delete mode 100644 projects/projected_extreme_snowfall/results/plot_relative_change_in_return_level.py

diff --git a/extreme_fit/function/margin_function/independent_margin_function.py b/extreme_fit/function/margin_function/independent_margin_function.py
index a381debd..d5001e69 100644
--- a/extreme_fit/function/margin_function/independent_margin_function.py
+++ b/extreme_fit/function/margin_function/independent_margin_function.py
@@ -36,14 +36,21 @@ class IndependentMarginFunction(AbstractMarginFunction):
         assert len(self.param_name_to_param_function) == len(self.params_class.PARAM_NAMES)
 
         # Potentially separate the coordinate into two groups: the spatio temporal coordnate & the climatic coordinate
+        # The climatic coordinate can be of two types either 1 and 0 vectors,
+        # or a vector with several information such as the GCM str, RCM str and the climate coordinates with effects
         if len(coordinate) > self.coordinates.nb_coordinates:
             assert self.param_name_to_ordered_climate_effects is not None
             assert AbstractCoordinates.COORDINATE_X not in self.coordinates.coordinates_names, \
                 'check the order of coordinates that everything is ok'
             climate_coordinate = coordinate[self.coordinates.nb_coordinates:].copy()
+            # Transform the climate coordinate if they are represent with strings
+            if not isinstance(climate_coordinate[0], float):
+                climate_coordinates_with_effects, gcm_rcm_couple = climate_coordinate
+                climate_coordinate = self.coordinates.get_climate_coordinate(climate_coordinates_with_effects, gcm_rcm_couple)
+            # Then build the param_name_to_total_effect dictionary
             param_name_to_total_effect = {param_name: np.dot(effects, climate_coordinate)
                                           for param_name, effects in self.param_name_to_ordered_climate_effects.items()}
-            coordinate = coordinate[:self.coordinates.nb_coordinates]
+            coordinate = np.array(coordinate[:self.coordinates.nb_coordinates])
         else:
             param_name_to_total_effect = None
 
@@ -59,18 +66,6 @@ class IndependentMarginFunction(AbstractMarginFunction):
             params[GevParams.SCALE] = np.exp(params[GevParams.SCALE])
         return self.params_class.from_dict(params)
 
-    def get_params_from_climate_model(self, coordinate: np.ndarray, climate_coordinate) -> GevParams:
-        """Main method to retrieve the additional parameters that take into account the effect of
-        a given rcp/gcm/rcm triplet"""
-        assert self.param_name_to_ordered_climate_effects is not None
-        assert all([c in AbstractCoordinates.COORDINATE_CLIMATE_MODEL_NAMES for c in climate_coordinate])
-        # the method should transform the rcp/gcm/rcm triplet into the appropriate coordinates then call get_params
-        # todo: or we could integrate this mode, into gev_params, so that we can call gev_params with three manners:
-        # with just spatio temporal corodinate
-        # spatio temporal coordinate +  the extended climate coordinate with 0 and 1
-        # spatio temproal coordinates + the climate corodinates with three columns
-        raise NotImplementedError
-
     def get_first_derivative_param(self, coordinate: np.ndarray, is_transformed: bool, dim: int = 0):
         transformed_coordinate = coordinate if is_transformed else self.transform(coordinate)
         return {
diff --git a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
index 6d1603c9..1c3f9052 100644
--- a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
+++ b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
@@ -8,21 +8,17 @@ from extreme_data.meteo_france_data.scm_models_data.abstract_study import Abstra
 from extreme_fit.distribution.gev.gev_params import GevParams
 from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
     AbstractTemporalLinearMarginModel
-from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
-    AbstractSpatioTemporalPolynomialModel
 from extreme_fit.model.margin_model.utils import MarginFitMethod
 from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies
 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, \
-    get_altitude_group_from_altitudes
+from extreme_trend.one_fold_fit.altitude_group import 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_gcm_rcm_effects import plot_gcm_rcm_effects
-from projects.projected_extreme_snowfall.results.plot_relative_change_in_return_level import \
-    plot_relative_dynamic_in_return_level
+from projects.projected_extreme_snowfall.results.part_3.plot_gcm_rcm_effects import plot_gcm_rcm_effects
+from projects.projected_extreme_snowfall.results.part_3.plot_relative_change_in_return_level import \
+    plot_relative_dynamic
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
@@ -115,11 +111,14 @@ class VisualizerForProjectionEnsemble(object):
                 plot_shoe_plot_changes_against_altitude(self.massif_names, visualizer_list, relative=relative,
                                                         with_significance=with_significance)
         else:
-            for relative in [True, False]:
-                plot_relative_dynamic_in_return_level(self.massif_names, visualizer_list,
-                                                      self.climate_coordinates_with_effects,
-                                                      self.safran_study_class,
-                                                      relative)
+            for relative in [None, True, False]:
+                for order in [None] + GevParams.PARAM_NAMES:
+                    plot_relative_dynamic(self.massif_names, visualizer_list,
+                                          self.climate_coordinates_with_effects,
+                                          self.safran_study_class,
+                                          relative,
+                                          order,
+                                          self.gcm_rcm_couples)
             if self.climate_coordinates_with_effects is not None:
                 climate_coordinate_with_effects_to_list = {
                     (AbstractCoordinates.COORDINATE_GCM, AbstractCoordinates.COORDINATE_RCM): self.gcm_rcm_couples,
diff --git a/extreme_trend/one_fold_fit/one_fold_fit.py b/extreme_trend/one_fold_fit/one_fold_fit.py
index a79e4cf9..f98f09e8 100644
--- a/extreme_trend/one_fold_fit/one_fold_fit.py
+++ b/extreme_trend/one_fold_fit/one_fold_fit.py
@@ -93,6 +93,9 @@ class OneFoldFit(object):
         elif order is None:
             return '{}-year return levels'.format(cls.return_period)
 
+    def get_moment_for_plots(self, altitudes, order=1, covariate_before=None, covariate_after=None):
+        return [self.get_moment(altitudes[0], covariate_after, order)]
+
     def get_moment(self, altitude, temporal_covariate, order=1):
         gev_params = self.get_gev_params(altitude, temporal_covariate)
         if order == 1:
@@ -101,6 +104,8 @@ class OneFoldFit(object):
             return gev_params.std
         elif order is None:
             return gev_params.return_level(return_period=self.return_period)
+        elif order in GevParams.PARAM_NAMES:
+            return gev_params.to_dict()[order]
         else:
             raise NotImplementedError
 
@@ -165,7 +170,7 @@ class OneFoldFit(object):
         for altitude in altitudes:
             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_change = 100 * (mean_after - mean_before) / np.abs(mean_before)
             relative_changes.append(relative_change)
         return relative_changes
 
@@ -213,7 +218,10 @@ class OneFoldFit(object):
 
     def get_coordinate(self, altitude, year):
         if isinstance(self.altitude_group, DefaultAltitudeGroup):
-            coordinate = np.array([year])
+            if not isinstance(altitude, list):
+                coordinate = np.array([year])
+            else:
+                coordinate = [year] + altitude
         else:
             coordinate = np.array([altitude, year])
         return coordinate
diff --git a/projects/projected_extreme_snowfall/checks/__init__.py b/projects/archive/comparison_window_for_the_max/__init__.py
similarity index 100%
rename from projects/projected_extreme_snowfall/checks/__init__.py
rename to projects/archive/comparison_window_for_the_max/__init__.py
diff --git a/projects/projected_extreme_snowfall/comparison_window_for_the_max/comparison_window_for_the_max.py b/projects/archive/comparison_window_for_the_max/comparison_window_for_the_max.py
similarity index 100%
rename from projects/projected_extreme_snowfall/comparison_window_for_the_max/comparison_window_for_the_max.py
rename to projects/archive/comparison_window_for_the_max/comparison_window_for_the_max.py
diff --git a/projects/projected_extreme_snowfall/comparison_window_for_the_max/__init__.py b/projects/projected_extreme_snowfall/data/checks/__init__.py
similarity index 100%
rename from projects/projected_extreme_snowfall/comparison_window_for_the_max/__init__.py
rename to projects/projected_extreme_snowfall/data/checks/__init__.py
diff --git a/projects/projected_extreme_snowfall/checks/main_inquiry_weird_pattern.py b/projects/projected_extreme_snowfall/data/checks/main_inquiry_weird_pattern.py
similarity index 100%
rename from projects/projected_extreme_snowfall/checks/main_inquiry_weird_pattern.py
rename to projects/projected_extreme_snowfall/data/checks/main_inquiry_weird_pattern.py
diff --git a/projects/projected_extreme_snowfall/evaluation/__init__.py b/projects/projected_extreme_snowfall/data/evaluation/__init__.py
similarity index 100%
rename from projects/projected_extreme_snowfall/evaluation/__init__.py
rename to projects/projected_extreme_snowfall/data/evaluation/__init__.py
diff --git a/projects/projected_extreme_snowfall/evaluation/comparison_historical_visualizer.py b/projects/projected_extreme_snowfall/data/evaluation/comparison_historical_visualizer.py
similarity index 100%
rename from projects/projected_extreme_snowfall/evaluation/comparison_historical_visualizer.py
rename to projects/projected_extreme_snowfall/data/evaluation/comparison_historical_visualizer.py
diff --git a/projects/projected_extreme_snowfall/evaluation/comparison_plot.py b/projects/projected_extreme_snowfall/data/evaluation/comparison_plot.py
similarity index 96%
rename from projects/projected_extreme_snowfall/evaluation/comparison_plot.py
rename to projects/projected_extreme_snowfall/data/evaluation/comparison_plot.py
index afce2c2e..7768bb69 100644
--- a/projects/projected_extreme_snowfall/evaluation/comparison_plot.py
+++ b/projects/projected_extreme_snowfall/data/evaluation/comparison_plot.py
@@ -4,12 +4,10 @@ import matplotlib.pyplot as plt
 from typing import Dict
 
 import numpy as np
-from matplotlib.lines import Line2D
 
 from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_rcm_couple_to_color
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_str, \
-    scenario_to_str
-from projects.projected_extreme_snowfall.evaluation.comparison_historical_visualizer import \
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import scenario_to_str
+from projects.projected_extreme_snowfall.data.evaluation.comparison_historical_visualizer import \
     ComparisonHistoricalVisualizer
 
 
diff --git a/projects/projected_extreme_snowfall/evaluation/main_comparison_on_quantile_period.py b/projects/projected_extreme_snowfall/data/evaluation/main_comparison_on_quantile_period.py
similarity index 100%
rename from projects/projected_extreme_snowfall/evaluation/main_comparison_on_quantile_period.py
rename to projects/projected_extreme_snowfall/data/evaluation/main_comparison_on_quantile_period.py
diff --git a/projects/projected_extreme_snowfall/evaluation/main_comparison_reanalysis.py b/projects/projected_extreme_snowfall/data/evaluation/main_comparison_reanalysis.py
similarity index 100%
rename from projects/projected_extreme_snowfall/evaluation/main_comparison_reanalysis.py
rename to projects/projected_extreme_snowfall/data/evaluation/main_comparison_reanalysis.py
diff --git a/projects/projected_extreme_snowfall/data/main_data.py b/projects/projected_extreme_snowfall/data/main_data.py
index bbda2b4a..40fc0ea6 100644
--- a/projects/projected_extreme_snowfall/data/main_data.py
+++ b/projects/projected_extreme_snowfall/data/main_data.py
@@ -1,5 +1,3 @@
-from collections import OrderedDict
-
 import matplotlib as mpl
 
 mpl.use('Agg')
@@ -7,11 +5,10 @@ mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
 import matplotlib.pyplot as plt
-from projects.projected_extreme_snowfall.evaluation.comparison_plot import individual_plot
 
 from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
-    scenario_to_real_scenarios, rcp_scenarios, rcm_scenarios_extended
+    scenario_to_real_scenarios
 from extreme_data.meteo_france_data.adamont_data.adamont_studies import AdamontStudies
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
 from extreme_data.meteo_france_data.scm_models_data.utils import Season
diff --git a/projects/projected_extreme_snowfall/preliminary_works/__init__.py b/projects/projected_extreme_snowfall/results/part_1/__init__.py
similarity index 100%
rename from projects/projected_extreme_snowfall/preliminary_works/__init__.py
rename to projects/projected_extreme_snowfall/results/part_1/__init__.py
diff --git a/projects/projected_extreme_snowfall/preliminary_works/main_preliminary_works.py b/projects/projected_extreme_snowfall/results/part_1/main_preliminary_works.py
similarity index 100%
rename from projects/projected_extreme_snowfall/preliminary_works/main_preliminary_works.py
rename to projects/projected_extreme_snowfall/results/part_1/main_preliminary_works.py
diff --git a/projects/projected_extreme_snowfall/results/part_2/__init__.py b/projects/projected_extreme_snowfall/results/part_2/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/projects/projected_extreme_snowfall/results/part_3/__init__.py b/projects/projected_extreme_snowfall/results/part_3/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py b/projects/projected_extreme_snowfall/results/part_3/main_projections_ensemble.py
similarity index 93%
rename from projects/projected_extreme_snowfall/results/main_projections_ensemble.py
rename to projects/projected_extreme_snowfall/results/part_3/main_projections_ensemble.py
index bdd44c08..f7bb1583 100644
--- a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py
+++ b/projects/projected_extreme_snowfall/results/part_3/main_projections_ensemble.py
@@ -55,15 +55,16 @@ def main():
                                              [AbstractCoordinates.COORDINATE_GCM, AbstractCoordinates.COORDINATE_RCM],
                                              [AbstractCoordinates.COORDINATE_GCM],
                                              [AbstractCoordinates.COORDINATE_RCM],
-                                        ][:1]  # None means we do not create any effect
+                                        ][1:2]  # None means we do not create any effect
     ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][1:]
     temporal_covariate_for_fit = [TimeTemporalCovariate,
                                   AnomalyTemperatureWithSplineTemporalCovariate][1]
     set_seed_for_test()
     AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
     scenarios = [AdamontScenario.rcp85_extended]
+    model_classes = SPLINE_MODELS_FOR_PROJECTION_ONE_ALTITUDE
 
-    fast = True
+    fast = None
     for scenario in scenarios:
         gcm_rcm_couples = get_gcm_rcm_couples(scenario)
         if fast is None:
@@ -71,9 +72,10 @@ def main():
             AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
             altitudes_list = [600, 2100, 3600]
         elif fast:
-            gcm_rcm_couples = gcm_rcm_couples[:2]
+            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]
 
@@ -94,8 +96,8 @@ def main():
         # Default parameters
         gcm_to_year_min_and_year_max = None
         massif_names = ['Vanoise']
-        model_classes = SPLINE_MODELS_FOR_PROJECTION_ONE_ALTITUDE
-        assert len(set(model_classes)) == 27
+        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:
@@ -111,7 +113,8 @@ def main():
                 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,
-                climate_coordinates_with_effects=climate_coordinates_with_effects
+                display_only_model_that_pass_gof_test=True,
+                climate_coordinates_with_effects=climate_coordinates_with_effects,
             )
             visualizer.plot()
 
diff --git a/projects/projected_extreme_snowfall/results/plot_gcm_rcm_effects.py b/projects/projected_extreme_snowfall/results/part_3/plot_gcm_rcm_effects.py
similarity index 65%
rename from projects/projected_extreme_snowfall/results/plot_gcm_rcm_effects.py
rename to projects/projected_extreme_snowfall/results/part_3/plot_gcm_rcm_effects.py
index 7deca889..f54934f7 100644
--- a/projects/projected_extreme_snowfall/results/plot_gcm_rcm_effects.py
+++ b/projects/projected_extreme_snowfall/results/part_3/plot_gcm_rcm_effects.py
@@ -1,32 +1,37 @@
-
 from typing import List
 
 import matplotlib.pyplot as plt
 import numpy as np
 
-from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_rcm_couple_to_color
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_rcm_couple_to_color, gcm_to_color
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_str
 from extreme_fit.distribution.gev.gev_params import GevParams
 from extreme_trend.one_fold_fit.altitudes_studies_visualizer_for_non_stationary_models import \
     AltitudesStudiesVisualizerForNonStationaryModels
-from projects.projected_extreme_snowfall.results.plot_relative_change_in_return_level import set_plot_name
-from root_utils import get_display_name_from_object_type
+from projects.projected_extreme_snowfall.results.part_3.plot_relative_change_in_return_level import set_plot_name
 
 
 def plot_gcm_rcm_effects(massif_names, visualizer_list: List[
     AltitudesStudiesVisualizerForNonStationaryModels], climate_coordinates_for_plot, climate_coordinates_with_effects,
-                                          safran_study_class,
-gcm_rcm_couples, param_name
-                                          ):
+                         safran_study_class,
+                         gcm_rcm_couples, param_name,
+                         ):
     ax = plt.gca()
     altitudes = [v.study.altitude for v in visualizer_list]
     visualizer = visualizer_list[0]
     assert len(massif_names) == 1
     assert climate_coordinates_with_effects is not None
     massif_name = massif_names[0]
+    all_effects = []
     for gcm_rcm_couple in gcm_rcm_couples:
-        plot_curve_gcm_rcm_effect(ax, massif_name, visualizer_list,
-                                  climate_coordinates_with_effects, gcm_rcm_couple, param_name)
+        effects = plot_curve_gcm_rcm_effect(ax, massif_name, visualizer_list,
+                                            climate_coordinates_with_effects, gcm_rcm_couple, param_name)
+        all_effects.append(effects)
+    all_effects = np.array(all_effects)
+    mean_effects = np.mean(all_effects, axis=0)
+    assert len(mean_effects) == len(altitudes)
+
+    ax.plot(mean_effects, altitudes, label='Mean effect', color='k', linewidth=4)
 
     effect_name = '-'.join([c.replace('coord_', '').upper() for c in climate_coordinates_for_plot])
     param_name_str = GevParams.full_name_from_param_name(param_name)
@@ -36,6 +41,7 @@ gcm_rcm_couples, param_name
     ax.set_ylabel('Altitude (m)')
     title = '{} parameter'.format(param_name_str)
     ax.set_ylim(top=altitudes[-1] + 1300)
+    ax.yaxis.set_ticks(altitudes)
     size = 7 if len(climate_coordinates_for_plot) == 2 else 10
     ax.legend(ncol=3, prop={'size': size})
     set_plot_name(climate_coordinates_for_plot, safran_study_class, title, visualizer)
@@ -44,19 +50,17 @@ gcm_rcm_couples, param_name
     plt.close()
 
 
-
 def plot_curve_gcm_rcm_effect(ax, massif_name, visualizer_list: List[AltitudesStudiesVisualizerForNonStationaryModels],
-climate_coordinates_with_effects,
-               gcm_rcm_couple, param_name):
+                              climate_coordinates_with_effects,
+                              gcm_rcm_couple, param_name):
     altitudes = [v.study.altitude for v in visualizer_list]
     effects = []
-    print("\n\n",climate_coordinates_with_effects, gcm_rcm_couple)
     for visualizer in visualizer_list[:]:
         one_fold_fit = visualizer.massif_name_to_one_fold_fit[massif_name]
-        indices = one_fold_fit.dataset.coordinates.get_indices_for_effects(climate_coordinates_with_effects, gcm_rcm_couple)
+        indices = one_fold_fit.dataset.coordinates.get_indices_for_effects(climate_coordinates_with_effects,
+                                                                           gcm_rcm_couple)
         assert len(indices) <= 2, indices
         ordered_climate_effects = one_fold_fit.best_function_from_fit.param_name_to_ordered_climate_effects[param_name]
-        print(ordered_climate_effects, indices)
         sum_effects = sum([ordered_climate_effects[i] for i in indices])
         effects.append(sum_effects)
     if len(gcm_rcm_couple) == 2:
@@ -64,4 +68,11 @@ climate_coordinates_with_effects,
         label = gcm_rcm_couple_to_str(gcm_rcm_couple)
         ax.plot(effects, altitudes, label=label, color=color)
     else:
-        ax.plot(effects, altitudes, label=gcm_rcm_couple[0])
+        is_gcm = gcm_rcm_couple[0] in gcm_to_color
+        if is_gcm:
+            gcm = gcm_rcm_couple[0]
+            ax.plot(effects, altitudes, label=gcm, color=gcm_to_color[gcm])
+        else:
+            rcm = gcm_rcm_couple[0]
+            ax.plot(effects, altitudes, label=rcm)
+    return effects
diff --git a/projects/projected_extreme_snowfall/results/part_3/plot_relative_change_in_return_level.py b/projects/projected_extreme_snowfall/results/part_3/plot_relative_change_in_return_level.py
new file mode 100644
index 00000000..4d325246
--- /dev/null
+++ b/projects/projected_extreme_snowfall/results/part_3/plot_relative_change_in_return_level.py
@@ -0,0 +1,100 @@
+from typing import List
+import matplotlib.pyplot as plt
+
+import numpy as np
+
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_trend.one_fold_fit.altitudes_studies_visualizer_for_non_stationary_models import \
+    AltitudesStudiesVisualizerForNonStationaryModels
+from extreme_trend.one_fold_fit.one_fold_fit import OneFoldFit
+from root_utils import get_display_name_from_object_type
+from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
+    AnomalyTemperatureWithSplineTemporalCovariate
+
+
+def plot_relative_dynamic(massif_names, visualizer_list: List[
+    AltitudesStudiesVisualizerForNonStationaryModels], climate_coordinates_with_effects,
+                          safran_study_class, relative, order,
+                          gcm_rcm_couples
+                          ):
+    ax = plt.gca()
+    visualizer = visualizer_list[0]
+    assert len(massif_names) == 1
+    is_temp_covariate = visualizer.temporal_covariate_for_fit is AnomalyTemperatureWithSplineTemporalCovariate
+    massif_name = massif_names[0]
+    for v in visualizer_list:
+        plot_curve(ax, massif_name, v, relative, is_temp_covariate, order, gcm_rcm_couples,
+                   climate_coordinates_with_effects)
+
+    xlabel = 'Anomaly of global temperature w.r.t. pre-industrial levels (K)' if is_temp_covariate else "Years"
+    ax.set_xlabel(xlabel)
+    change = 'Relative change in' if relative is True else ("Value of" if relative is None else "Change in")
+    unit = '\%' if relative is True else (visualizer.study.variable_unit if order != GevParams.SHAPE else 'no unit')
+    name = 'the {} parameter'.format(GevParams.full_name_from_param_name(order)) if order is not None \
+        else '{}-year return levels'.format(OneFoldFit.return_period)
+    ylabel = '{} {} ({})'.format(change, name, unit)
+    ax.set_ylabel(ylabel)
+
+    ax.legend(ncol=2, prop={'size': 9}, loc='upper left')
+    title = ylabel.split('(')[0]
+    set_plot_name(climate_coordinates_with_effects, safran_study_class, title, visualizer)
+    visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
+
+    plt.close()
+
+
+def plot_curve(ax, massif_name, visualizer: AltitudesStudiesVisualizerForNonStationaryModels,
+               relative, is_temp_cov, order, gcm_rcm_couples, climate_coordinates_with_effects
+               ):
+    if is_temp_cov:
+        x_list = np.linspace(1, 4.5, num=400)
+        covariate_before = 1
+    else:
+        x_list = np.linspace(1951, 2100, num=150)
+        covariate_before = 1951
+    one_fold_fit = visualizer.massif_name_to_one_fold_fit[massif_name]
+    print(relative, order)
+    print(get_display_name_from_object_type(type(one_fold_fit.best_margin_model)),
+          "significant={}".format(one_fold_fit.is_significant))
+    if relative is None:
+        f = one_fold_fit.get_moment_for_plots
+    else:
+        f = one_fold_fit.relative_changes_of_moment if relative else one_fold_fit.changes_of_moment
+    color = altitude_to_color[visualizer.study.altitude]
+    # Plot the main trend
+    changes = [f([None], order=order, covariate_before=covariate_before, covariate_after=t)[0] for t in x_list]
+    label = '{} m'.format(visualizer.altitude_group.reference_altitude)
+    ax.plot(x_list, changes, label=label, color=color, linewidth=4)
+    # Plot the sub trend, i.e. for each GCM-RCM couples
+    for gcm_rcm_couple in gcm_rcm_couples[:]:
+        fake_altitude = [climate_coordinates_with_effects, gcm_rcm_couple]
+        changes = [f([fake_altitude], order=order, covariate_before=covariate_before, covariate_after=t)[0] for t in x_list]
+        ax.plot(x_list, changes, color=color, linewidth=1, linestyle='dotted')
+
+
+altitude_to_color = {
+    # Low altitude group
+    600: 'darkred',
+    900: 'red',
+    # Mid altitude gruop
+    1200: 'darkorange',
+    1500: 'orange',
+    1800: 'gold',
+    # High altitude group
+    2100: 'lightskyblue',
+    2400: 'skyblue',
+    2700: 'dodgerblue',
+    # Very high altitude group
+    3000: 'b',
+    3300: 'mediumblue',
+    3600: 'darkblue',
+}
+
+
+def set_plot_name(climate_coordinates_with_effects, safran_study_class, title, visualizer):
+    plot_name = ' %s' % title
+    plot_name += ' with {} effects'.format('no' if climate_coordinates_with_effects is None
+                                           else ' and '.join(
+        [c.replace('coord_', '') for c in climate_coordinates_with_effects]))
+    plot_name += ' with{} observations'.format('out' if safran_study_class is None else '')
+    visualizer.plot_name = plot_name
diff --git a/projects/projected_extreme_snowfall/results/plot_relative_change_in_return_level.py b/projects/projected_extreme_snowfall/results/plot_relative_change_in_return_level.py
deleted file mode 100644
index af61ad06..00000000
--- a/projects/projected_extreme_snowfall/results/plot_relative_change_in_return_level.py
+++ /dev/null
@@ -1,55 +0,0 @@
-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 root_utils import get_display_name_from_object_type
-from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
-    AnomalyTemperatureWithSplineTemporalCovariate
-
-
-def plot_relative_dynamic_in_return_level(massif_names, visualizer_list: List[
-    AltitudesStudiesVisualizerForNonStationaryModels], climate_coordinates_with_effects,
-                                          safran_study_class, relative
-                                          ):
-    ax = plt.gca()
-    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, relative)
-
-    ax.set_xlabel('Anomaly of global temperature w.r.t. pre-industrial levels (K)')
-    change = 'Relative change' if relative else "Change"
-    unit = '\%' if relative else visualizer.study.variable_unit
-    ax.set_ylabel((change + (' in return levels (' + unit + ')')))
-
-    ax.legend(ncol=2, prop={'size': 9}, loc='upper left')
-    title = 'change {}of return level'.format('relative ' if relative else '')
-    set_plot_name(climate_coordinates_with_effects, safran_study_class, title, visualizer)
-    visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
-
-    plt.close()
-
-def set_plot_name(climate_coordinates_with_effects, safran_study_class, title, visualizer):
-    plot_name = ' %s' % title
-    plot_name += ' with {} effects'.format('no' if climate_coordinates_with_effects is None
-                                         else ' and '.join(
-        [c.replace('coord_', '') for c in climate_coordinates_with_effects]))
-    plot_name += ' with{} observations'.format('out' if safran_study_class is None else '')
-    visualizer.plot_name = plot_name
-
-
-def plot_curve(ax, massif_name, visualizer: AltitudesStudiesVisualizerForNonStationaryModels,
-               relative):
-    temperatures_list = np.linspace(1, 4.5, num=400)
-    one_fold_fit = visualizer.massif_name_to_one_fold_fit[massif_name]
-    print(get_display_name_from_object_type(type(one_fold_fit.best_margin_model)),
-          "significant={}".format(one_fold_fit.is_significant))
-    f = one_fold_fit.relative_changes_of_moment if relative else one_fold_fit.changes_of_moment
-    return_levels = [f([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)
diff --git a/spatio_temporal_dataset/coordinates/abstract_coordinates.py b/spatio_temporal_dataset/coordinates/abstract_coordinates.py
index 729236fb..a195fc6a 100644
--- a/spatio_temporal_dataset/coordinates/abstract_coordinates.py
+++ b/spatio_temporal_dataset/coordinates/abstract_coordinates.py
@@ -316,22 +316,29 @@ class AbstractCoordinates(object):
         unique_values_without_nan = [v for v in unique_values if isinstance(v, str)]
         return s, unique_values, unique_values_without_nan
 
-    def load_ordered_columns_names(self, climate_coordinates, split=Split.all):
+    def load_ordered_columns_names(self, climate_coordinates_names_with_effects, split=Split.all):
         column_names = []
-        for climate_coordinate in climate_coordinates:
+        for climate_coordinate in climate_coordinates_names_with_effects:
             _, _, names = self.load_unique_values(climate_coordinate, split)
             column_names.extend(names)
         return column_names
 
-    def get_indices_for_effects(self, climate_coordinates, gcm_rcm_couple):
+    def get_indices_for_effects(self, climate_coordinates_names_with_effects, gcm_rcm_couple):
         indices = []
-        columns_names = self.load_ordered_columns_names(climate_coordinates)
+        columns_names = self.load_ordered_columns_names(climate_coordinates_names_with_effects)
         for name in gcm_rcm_couple:
             name_for_fit = self.climate_model_coordinate_name_to_name_for_fit(name)
             index = columns_names.index(name_for_fit)
             indices.append(index)
         return indices
 
+    def get_climate_coordinate(self, climate_coordinates_names_with_effects, gcm_rcm_couple):
+        columns_names = self.load_ordered_columns_names(climate_coordinates_names_with_effects)
+        climate_coordinates = np.zeros(len(columns_names))
+        for indice in self.get_indices_for_effects(climate_coordinates_names_with_effects, gcm_rcm_couple):
+            climate_coordinates[indice] = 1
+        return climate_coordinates
+
     def df_climate_models(self, split=Split.all):
         return df_sliced(df=self.df_coordinate_climate_model, split=split, slicer=self.slicer)
 
-- 
GitLab