[contrasting] add histogram for model repartition & evolution of trends....

[contrasting] add histogram for model repartition & evolution of trends. improve preliminary_analysis.py

projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py

-
 import time
 from typing import List
 
 import matplotlib as mpl
 
+from projects.altitude_spatial_model.altitudes_fit.plot_histogram_altitude_studies import \
+    plot_histogram_all_trends_against_altitudes, plot_histogram_all_models_against_altitudes
+
 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 extreme_fit.model.margin_model.polynomial_margin_model.utils import \
     ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups
@@ -22,29 +27,12 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_s
 
 
 def main():
-    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][4:6]
-    # todo: l ecart  pour les saisons de l automne ou de sprint
-    #  vient probablement de certains zéros pas pris en compte pour le fit,
-    # mais pris en compte pour le calcul de mon aic
-    # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
-    # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900][:]
-    # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
-    # altitudes = [600, 900, 1200, 1500, 1800, 2100][:]
-    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200, 4500, 4800][:]
-    # altitudes = [1500, 1800][:]
-    study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
-    study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
-                     SafranPrecipitation7Days][:]
-    study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranPrecipitation1Day
-                        , SafranPrecipitation3Days][:1]
-    altitudes = [1800, 2100, 2400]
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:1]
     # study_classes = [SafranPrecipitation1Day][:1]
 
-    # Common parameters
-    # altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]
     massif_names = None
     # massif_names = ['Mercantour', 'Vercors', 'Ubaye']
+
     seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
 
     main_loop(altitudes_for_groups[:], massif_names, seasons, study_classes)
@@ -55,69 +43,29 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes):
     assert isinstance(altitudes_list[0], List)
     for season in seasons:
         for study_class in study_classes:
-            # if issubclass(study_class, SimulationStudy):
-            #     for ensemble_idx in list(range(14))[:1]:
-            #         studies = AltitudesStudies(study_class, altitudes, season=season,
-            #                                    ensemble_idx=ensemble_idx)
-            #         plot_studies(massif_names, season, studies, study_class)
-            # else:
+            print('Inner loop', season, study_class)
             visualizer_list = load_visualizer_list(season, study_class, altitudes_list, massif_names)
+            plot_visualizers(massif_names, visualizer_list)
             for visualizer in visualizer_list:
-                plots(massif_names, season, visualizer)
+                plot_visualizer(massif_names, visualizer)
             del visualizer_list
             time.sleep(2)
 
 
+def plot_visualizers(massif_names, visualizer_list):
+    plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
+    plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
 
-def load_visualizer_list(season, study_class, altitudes_list, massif_names):
-    model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
-    visualizer_list = []
-    # Load all studies
-    for altitudes in altitudes_list:
-        studies = AltitudesStudies(study_class, altitudes, season=season)
-        visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
-                                                                      model_classes=model_classes,
-                                                                      massif_names=massif_names,
-                                                                      show=False,
-                                                                      temporal_covariate_for_fit=None,
-                                                                      # temporal_covariate_for_fit=MeanAlpsTemperatureCovariate,
-                                                                      )
-        visualizer_list.append(visualizer)
-    # 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
-
-    return visualizer_list
-
-
-def plots(massif_names, season, visualizer):
-    studies = visualizer.studies
-    print('inner loop', season, type(studies.study))
 
+def plot_visualizer(massif_names, visualizer):
     # Plot time series
-    studies.plot_maxima_time_series(massif_names=massif_names)
-
-    # Plot moments
+    # visualizer.studies.plot_maxima_time_series(massif_names=massif_names)
+    # Plot moments against altitude
     # for std in [True, False][:]:
     #     for change in [True, False, None]:
     #         studies.plot_mean_maxima_against_altitude(massif_names=massif_names, std=std, change=change)
-
     # Plot the results for the model that minimizes the individual aic
     plot_individual_aic(visualizer)
-
     # Plot the results for the model that minimizes the total aic
     # plot_total_aic(model_classes, visualizer)
 

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py

+from collections import Counter
 from typing import List, Dict
 
 import matplotlib
@@ -116,23 +117,39 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         moment = ' '.join(method_name.split('_'))
         moment = moment.replace('moment', '{} in 2019'.format(OneFoldFit.get_moment_str(order=order)))
         plot_name = '{}{} for annual maxima of {}'.format(OneFoldFit.folder_for_plots, moment,
-                                                      SCM_STUDY_CLASS_TO_ABBREVIATION[
-                                                          self.studies.study_class])
+                                                          SCM_STUDY_CLASS_TO_ABBREVIATION[
+                                                              self.studies.study_class])
+
+        massif_name_to_text = self.massif_name_to_best_name
         if 'change' in method_name:
-            plot_name += '\nbetween {} and {}'.format(2019-50, 2019)
+            plot_name += '\nbetween {} and {}'.format(2019 - 50, 2019)
+            if 'relative' not in method_name:
+                # Put the relative score as text on the plot for the change.
+                massif_name_to_text = {m: ('+' if v > 0 else '') + str(int(v)) + '\%' for m, v in
+                                       self.method_name_and_order_to_d(self.moment_names[2], order).items()}
+
         parenthesis = self.study.variable_unit if 'relative' not in method_name else '\%'
         ylabel = '{} ({})'.format(plot_name, parenthesis)
 
+        max_abs_change = self.method_name_and_order_to_max_abs(method_name, order)
+        add_colorbar = self.add_colorbar
+
+        is_return_level_plot = (self.moment_names.index(method_name) == 0) and (order is None)
+        if is_return_level_plot:
+            add_colorbar = True
+            max_abs_change = None
+            massif_name_to_text = {m: '+- 0' for m in massif_name_to_text.keys()}
+
+        negative_and_positive_values = self.moment_names.index(method_name) > 0
         # Plot the map
-        a_change_is_displayed = self.moment_names.index(method_name) > 0
         self.plot_map(cmap=plt.cm.coolwarm, fit_method=self.fit_method, graduation=10,
                       label=ylabel, massif_name_to_value=massif_name_to_value,
                       plot_name=plot_name, add_x_label=True,
-                      negative_and_positive_values=a_change_is_displayed,
+                      negative_and_positive_values=negative_and_positive_values,
                       altitude=self.altitude_group.reference_altitude,
-                      add_colorbar=self.add_colorbar,
-                      max_abs_change=self.method_name_and_order_to_max_abs(method_name, order),
-                      massif_name_to_text=self.massif_name_to_name,
+                      add_colorbar=add_colorbar,
+                      max_abs_change=max_abs_change,
+                      massif_name_to_text=massif_name_to_text,
                       )
 
     @property
@@ -175,7 +192,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                 for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items()}
 
     @property
-    def massif_name_to_name(self):
+    def massif_name_to_best_name(self):
         return {massif_name: one_fold_fit.best_name
                 for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items()}
 
@@ -211,23 +228,23 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         negative_and_positive_values = (max(values) > 0) and (min(values) < 0)
         raise NotImplementedError
         self.plot_map(massif_name_to_value=massif_name_to_best_coef,
-                                label='{}Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots,
-                                                                        coef_name,
-                                                                        SCM_STUDY_CLASS_TO_ABBREVIATION[
-                                                                            type(self.study)],
-                                                                        self.study.variable_unit),
-                                add_x_label=False, graduation=graduation, massif_name_to_text=self.massif_name_to_name,
-                                negative_and_positive_values=negative_and_positive_values)
+                      label='{}Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots,
+                                                              coef_name,
+                                                              SCM_STUDY_CLASS_TO_ABBREVIATION[
+                                                                  type(self.study)],
+                                                              self.study.variable_unit),
+                      add_x_label=False, graduation=graduation, massif_name_to_text=self.massif_name_to_best_name,
+                      negative_and_positive_values=negative_and_positive_values)
 
     def plot_shape_map(self):
 
         label = 'Shape parameter for {} maxima of {} in 2019'.format(self.study.season_name,
-                                                                          SCM_STUDY_CLASS_TO_ABBREVIATION[
-                                                                              type(self.study)])
+                                                                     SCM_STUDY_CLASS_TO_ABBREVIATION[
+                                                                         type(self.study)])
         self.plot_map(massif_name_to_value=self.massif_name_to_shape,
                       label=label,
                       plot_name=label,
-                      add_x_label=True, graduation=0.1, massif_name_to_text=self.massif_name_to_name,
+                      add_x_label=True, graduation=0.1, massif_name_to_text=self.massif_name_to_best_name,
                       cmap=matplotlib.cm.get_cmap('BrBG_r'),
                       altitude=self.altitude_group.reference_altitude,
                       add_colorbar=self.add_colorbar,
@@ -350,3 +367,53 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         plot_name = 'Switch altitude'
         self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
         plt.close()
+
+    def all_trends(self, massif_names):
+        """return percents which contain decrease, significant decrease, increase, significant increase percentages"""
+        valid_massif_names = self.get_valid_names(massif_names)
+
+        nb_valid_massif_names = len(valid_massif_names)
+        nbs = np.zeros(4)
+        relative_changes = [[], []]
+        for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
+                                   if m in valid_massif_names]:
+            if one_fold.change_for_reference_altitude == 0:
+                continue
+            # Compute relative changes
+            relative_changes[0].append(one_fold.relative_change_for_reference_altitude)
+            if one_fold.is_significant:
+                relative_changes[1].append(one_fold.relative_change_for_reference_altitude)
+            # Compute nbs
+            idx = 0 if one_fold.change_for_reference_altitude < 0 else 2
+            nbs[idx] += 1
+            if one_fold.is_significant:
+                nbs[idx + 1] += 1
+
+        percents = 100 * nbs / nb_valid_massif_names
+        mean_relative_changes = [np.mean(l) for l in relative_changes]
+        return [nb_valid_massif_names] + mean_relative_changes + list(percents)
+
+    def get_valid_names(self, massif_names):
+        valid_massif_names = set(self.massif_name_to_one_fold_fit.keys())
+        if massif_names is not None:
+            valid_massif_names = valid_massif_names.intersection(set(massif_names))
+        return valid_massif_names
+
+    def model_name_to_percentages(self, massif_names):
+        valid_massif_names = self.get_valid_names(massif_names)
+        nb_valid_massif_names = len(valid_massif_names)
+        best_names = [one_fold_fit.best_estimator.margin_model.name_str
+                      for m, one_fold_fit in self.massif_name_to_one_fold_fit.items()
+                      if m in valid_massif_names]
+        counter = Counter(best_names)
+        d = {name: 100 * c / nb_valid_massif_names for name, c in counter.items()}
+        # Add 0 for the name not present
+        for name in self.model_names:
+            if name not in d:
+                d[name] = 0
+        return d
+
+    @property
+    def model_names(self):
+        massif_name = list(self.massif_name_to_one_fold_fit.keys())[0]
+        return self.massif_name_to_one_fold_fit[massif_name].model_names

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py

@@ -88,6 +88,14 @@ class OneFoldFit(object):
     def moment(self, altitudes, year=2019, order=1):
         return [self.get_moment(altitude, year, order) for altitude in altitudes]
 
+    @property
+    def change_for_reference_altitude(self) -> float:
+        return self.changes_for_moment(altitudes=[self.altitude_plot])[0]
+
+    @property
+    def relative_change_for_reference_altitude(self) -> float:
+        return self.relative_changes_for_moment(altitudes=[self.altitude_plot])[0]
+
     def changes_for_moment(self, altitudes, year=2019, nb_years=50, order=1):
         changes = []
         for altitude in altitudes:
@@ -185,6 +193,10 @@ class OneFoldFit(object):
         except (TypeError, KeyError):
             return None
 
+    @property
+    def model_names(self):
+        return [e.margin_model.name_str for e in self.sorted_estimators]
+
     @property
     def best_name(self):
         name = self.best_estimator.margin_model.name_str

projects/altitude_spatial_model/altitudes_fit/plot_histogram_altitude_studies.py

+from typing import List
+
+import numpy as np
+
+import matplotlib.pyplot as plt
+
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
+    AltitudesStudiesVisualizerForNonStationaryModels
+
+
+def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: List[
+    AltitudesStudiesVisualizerForNonStationaryModels]):
+    visualizer = visualizer_list[0]
+    model_names = visualizer.model_names
+    model_name_to_percentages = {model_name: [] for model_name in model_names}
+    for v in visualizer_list:
+        for model_name, percentage in v.model_name_to_percentages(massif_names).items():
+            model_name_to_percentages[model_name].append(percentage)
+    model_name_to_mean_percentage = {m: np.mean(a) for m, a in model_name_to_percentages.items()}
+    sorted_model_names = sorted(model_names, key=lambda m: model_name_to_mean_percentage[m], reverse=True)
+
+    # Plot part
+    ax = plt.gca()
+    width = 5
+    size = 8
+    legend_fontsize = 10
+    labelsize = 10
+    linewidth = 1
+    tick_list = np.array([((len(visualizer_list) + 2) * i + (1 + len(visualizer_list) / 2)) * width
+                          for i in range(len(sorted_model_names))])
+    for tick_middle, model_name in zip(tick_list, sorted_model_names):
+        x_shifted = [tick_middle + width * shift / 2 for shift in range(-3, 5, 2)]
+        percentages = model_name_to_percentages[model_name]
+        colors = ['white', 'yellow', 'orange', 'red']
+        labels = ['{} m (\% out of {} massifs)'.format(v.altitude_group.reference_altitude,
+                                                      len(v.get_valid_names(massif_names))) for v in visualizer_list]
+        for x, color, percentage, label in zip(x_shifted, colors, percentages, labels):
+            ax.bar([x], [percentage], width=width, label=label,
+                   linewidth=linewidth, edgecolor='black', color=color)
+
+    handles, labels = ax.get_legend_handles_labels()
+    ax.legend(handles[:len(visualizer_list)], labels[:len(visualizer_list)], prop={'size': size})
+    ax.set_xticklabels(sorted_model_names)
+    ax.set_xticks(tick_list)
+    ax.set_ylabel('Percentage of massifs (\%) ', fontsize=legend_fontsize)
+    ax.set_xlabel('Models', fontsize=legend_fontsize)
+    ax.set_ylim(bottom=0)
+    ax.yaxis.grid()
+    ax.tick_params(axis='both', which='major', labelsize=labelsize)
+
+    visualizer.plot_name = 'All models'
+    visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
+    plt.close()
+
+
+def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list):
+    all_trends = [v.all_trends(massif_names) for v in visualizer_list]
+    nb_massifs, mean_relative_change, mean_significant_relative_change, *all_l = zip(*all_trends)
+
+    plt.close()
+    ax = plt.gca()
+    width = 5
+    size = 8
+    legend_fontsize = 10
+    labelsize = 10
+    linewidth = 3
+    x = np.array([3 * width * (i + 1) for i in range(len(nb_massifs))])
+
+    colors = ['blue', 'darkblue', 'red', 'darkred']
+    labels = []
+    for suffix in ['Decrease', 'Increase']:
+        for prefix in ['Non significant', 'Significant']:
+            labels.append('{} {}'.format(prefix, suffix))
+    for l, color, label in zip(all_l, colors, labels):
+        x_shifted = x - width / 2 if 'blue' in color else x + width / 2
+        ax.bar(x_shifted, l, width=width, color=color, edgecolor=color, label=label,
+               linewidth=linewidth)
+    ax.legend(loc='upper left', prop={'size': size})
+    ax.set_ylabel('Percentage of massifs (\%) ', fontsize=legend_fontsize)
+    ax.set_xlabel('Altitudes', fontsize=legend_fontsize)
+    ax.tick_params(axis='both', which='major', labelsize=labelsize)
+    ax.set_xticks(x)
+    ax.yaxis.grid()
+    ax.set_xticklabels([str(v.altitude_group.reference_altitude) for v in visualizer_list])
+
+    plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x)
+
+    # PLot mean relative change against altitude
+    ax_twinx = ax.twinx()
+    colors = ['grey', 'black']
+    relative_changes = [mean_relative_change, mean_significant_relative_change]
+    labels = ['All massifs with a trend', 'Only massifs with a significant trend']
+    for color, y, label in zip(colors, relative_changes, labels):
+        print(x, y)
+        ax_twinx.plot(x, y, label=label, color=color)
+    ax_twinx.legend(loc='upper right', prop={'size': size})
+
+    ax_twinx.set_ylabel('Mean relative change average on massif with trends', fontsize=legend_fontsize)
+
+    visualizer = visualizer_list[0]
+    visualizer.plot_name = 'All trends'
+    visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
+
+    plt.close()
+
+
+def plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x):
+    # Plot number of massifs on the upper axis
+    ax_twiny = ax.twiny()
+    ax_twiny.plot(x, [0 for _ in x], linewidth=0)
+    ax_twiny.set_xticks(x)
+    ax_twiny.tick_params(labelsize=labelsize)
+    ax_twiny.set_xticklabels(nb_massifs)
+    ax_twiny.set_xlim(ax.get_xlim())
+    ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=legend_fontsize)

projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py

+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
+
+
+def load_visualizer_list(season, study_class, altitudes_list, massif_names):
+    model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+    visualizer_list = []
+    # Load all studies
+    for altitudes in altitudes_list:
+        # if issubclass(study_class, SimulationStudy):
+        #     for ensemble_idx in list(range(14))[:1]:
+        #         studies = AltitudesStudies(study_class, altitudes, season=season,
+        #                                    ensemble_idx=ensemble_idx)
+        #         plot_studies(massif_names, season, studies, study_class)
+        # else:
+        studies = AltitudesStudies(study_class, altitudes, season=season)
+        visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
+                                                                      model_classes=model_classes,
+                                                                      massif_names=massif_names,
+                                                                      show=False,
+                                                                      temporal_covariate_for_fit=None,
+                                                                      # temporal_covariate_for_fit=MeanAlpsTemperatureCovariate,
+                                                                      )
+        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/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_study_visualizer.py → projects/altitude_spatial_model/preliminary_analysis.py

@@ -29,25 +29,29 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
             massif_name_to_r2_score = {}
             for massif_name in self.study.all_massif_names()[:]:
                 linear_coef, _, r2 = self._plot_gev_params_against_altitude_one_massif(ax, massif_name, param_name)
-                massif_name_to_linear_coef[massif_name] = linear_coef[0]
+                massif_name_to_linear_coef[massif_name] = 1000 * linear_coef[0]
                 massif_name_to_r2_score[massif_name] = str(round(r2, 2))
             print(massif_name_to_linear_coef, massif_name_to_r2_score)
             # Plot change against altitude
             ax.legend(prop={'size': 7}, ncol=3)
             ax.set_xlabel('Altitude')
             ax.set_ylabel(GevParams.full_name_from_param_name(param_name) + ' parameter for a stationary GEV distribution')
-            plot_name = '{} change /with altitude'.format(param_name)
+            plot_name = '{} change with altitude'.format(param_name)
             self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
             ax.clear()
             # Plot map of slope for each massif
             visualizer = StudyVisualizer(study=self.study, show=False, save_to_file=True)
-            ylabel = 'Linear slope for the {} parameter against the altitude'.format(param_name)
+            ylabel = 'Linear slope for the {} parameter (change every 1000 m of altitude)'.format(param_name)
+            gev_param_name_to_graduation = {
+                GevParams.LOC: 5,
+                GevParams.SCALE: 1,
+                GevParams.SHAPE: 0.1,
+            }
             visualizer.plot_map(cmap=plt.cm.coolwarm, fit_method=self.study.fit_method,
-                                graduation=1,
+                                graduation=gev_param_name_to_graduation[param_name],
                                 label=ylabel, massif_name_to_value=massif_name_to_linear_coef,
                                 plot_name=ylabel, add_x_label=False,
-                                # negative_and_positive_values=param_name == GevParams.SHAPE,
-                                negative_and_positive_values=True,
+                                negative_and_positive_values=param_name == GevParams.SHAPE,
                                 add_colorbar=True,
                                 massif_name_to_text=massif_name_to_r2_score,
                                 )
@@ -190,7 +194,7 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
 
 if __name__ == '__main__':
     altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300]
-    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]
+    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
     # altitudes = paper_altitudes
     # altitudes = [1800, 2100]
     visualizer = PointwiseGevStudyVisualizer(SafranSnowfall1Day, altitudes=altitudes)

test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py

@@ -65,7 +65,7 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
         self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
         self.fit_method = MarginFitMethod.extremes_fevd_mle
 
-    def function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(self, model_class):
+    def function_test_gev_spatio_temporal_margin_fit_non_stationary(self, model_class):
         # Create estimator
         estimator = fitted_linear_margin_estimator_short(model_class=model_class,
                                                          dataset=self.dataset,
@@ -76,32 +76,36 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
         self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
         self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
         # Assert we can compute the return level
-        covariate_for_return_level = np.array([400, 25])[::1]
+        covariate_for_return_level = np.array([400, 20])[::1]
         confidence_interval = EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator,
                                                                                              ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
                                                                                              temporal_covariate=covariate_for_return_level)
         return_level = estimator.function_from_fit.get_params(covariate_for_return_level).return_level(50)
+        print("my return level", return_level)
         self.assertAlmostEqual(return_level, confidence_interval.mean_estimate)
         self.assertFalse(np.isnan(confidence_interval.confidence_interval[0]))
         self.assertFalse(np.isnan(confidence_interval.confidence_interval[1]))
 
     def test_gev_spatio_temporal_margin_fit_1(self):
-        self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(StationaryAltitudinal)
+        self.function_test_gev_spatio_temporal_margin_fit_non_stationary(StationaryAltitudinal)
 
+    #
     def test_gev_spatio_temporal_margin_fit_1_bis(self):
-        self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(AltitudinalShapeLinearTimeStationary)
+        self.function_test_gev_spatio_temporal_margin_fit_non_stationary(AltitudinalShapeLinearTimeStationary)
 
     # def test_gev_spatio_temporal_margin_fit_2(self):
-    #     self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(
+    #     # first model with both a time and altitude non stationarity
+    #     self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
     #         AltitudinalShapeConstantTimeLocationLinear)
-    #
+
     # def test_gev_spatio_temporal_margin_fit_3(self):
-    #     self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(
+    #     self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
     #         AltitudinalShapeLinearTimeLocationLinear)
     #
     # def test_gev_spatio_temporal_margin_fit_4(self):
-    #     self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(
+    #     self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
     #         AltitudinalShapeLinearTimeLocScaleLinear)
 
+
 if __name__ == '__main__':
     unittest.main()

test/test_projects/test_exceeding_snow_loads/test_results.py

@@ -8,11 +8,12 @@ from projects.exceeding_snow_loads.section_results.plot_selection_curves import
 from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_curves, plot_trend_map
 from projects.exceeding_snow_loads.section_results.plot_uncertainty_curves import plot_uncertainty_massifs
 from projects.exceeding_snow_loads.section_results.plot_uncertainty_histogram import plot_uncertainty_histogram
-
+import matplotlib.pyplot as plt
 
 class TestResults(unittest.TestCase):
 
     def test_run_intermediate_results(self):
+        plt.close()
         # Load data
         altitude_to_visualizer = load_altitude_to_visualizer(altitudes=[300, 600], massif_names=None,
                                                              study_class=CrocusSnowLoadTotal,