diff --git a/extreme_data/eurocode_data/utils.py b/extreme_data/eurocode_data/utils.py
index 3e5312cadbac0823237cb891c2de27c0b3e2394d..ed431a9931b6bf7f790de5b4200487323dc2d354 100644
--- a/extreme_data/eurocode_data/utils.py
+++ b/extreme_data/eurocode_data/utils.py
@@ -1,6 +1,6 @@
 
 # Eurocode quantile str
-EUROCODE_RETURN_LEVEL_STR = '50-year return level of GSL (kN $m^-2$)'
+EUROCODE_RETURN_LEVEL_STR = '50-year return levels of GSL (kN $m^-2$)'
 # Eurocode quantile correspond to a 50 year return period
 EUROCODE_QUANTILE = 0.98
 # Altitudes (between low and mid altitudes) < 2000m and should be > 200m
diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index 43696208f02528068fd63d3b83d98e40613c9dd9..c6212fef48cc9131262ab220becfe4e448fd770a 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
@@ -589,8 +589,8 @@ class AbstractStudy(object):
             legend_elements = cls.get_legend_for_model_symbol(marker_style_to_label_name, markersize=7)
             ax.legend(handles=legend_elements[:], loc='upper right', prop={'size': 9})
             # ax.legend(handles=legend_elements[4:], bbox_to_anchor=(0.01, 0.03), loc='lower left')
-            ax.annotate(' '.join(filled_marker_legend_list),
-                        xy=(0.05, 0.015), xycoords='axes fraction', fontsize=7)
+            # ax.annotate(' '.join(filled_marker_legend_list),
+            #             xy=(0.05, 0.015), xycoords='axes fraction', fontsize=7)
 
         if show:
             plt.show()
diff --git a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
index 36a4f22c044a3d449c88a682884d94664847395f..b07a5eca8e0167fd5898fdaefc5ebb87005c086f 100644
--- a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
+++ b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
@@ -58,12 +58,14 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                  select_only_acceptable_shape_parameter=True,
                  fit_gev_only_on_non_null_maxima=False,
                  fit_only_time_series_with_ninety_percent_of_non_null_values=True,
+                 select_only_model_that_pass_anderson_test=False,
                  ):
         super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
                          year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
                          transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis,
                          normalization_under_one_observations)
         # Add some attributes
+        self.select_only_model_that_pass_anderson_test = select_only_model_that_pass_anderson_test
         self.fit_only_time_series_with_ninety_percent_of_non_null_values = fit_only_time_series_with_ninety_percent_of_non_null_values
         self.fit_gev_only_on_non_null_maxima = fit_gev_only_on_non_null_maxima
         self.select_only_acceptable_shape_parameter = select_only_acceptable_shape_parameter
@@ -76,6 +78,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         self.uncertainty_massif_names = uncertainty_massif_names
         # Assign some default arguments
         if self.model_subsets_for_uncertainty is None:
+            # We regroup the results either with all the models,
+            # or only the stationary gumbel model which correspond to the building standards
             self.model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
                                                   ModelSubsetForUncertainty.non_stationary_gumbel_and_gev][:]
         if self.uncertainty_methods is None:
@@ -155,6 +159,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     @cached_property
     def massif_name_to_trend_test_tuple(self) -> Tuple[
         Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest]]:
+        print('cached', self.altitude, id(self))
 
         massif_name_to_trend_test_that_minimized_aic = {}
         massif_name_to_stationary_trend_test_that_minimized_aic = {}
@@ -165,7 +170,13 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
             # Extract the stationary or non-stationary model that minimized AIC
             trend_test_that_minimized_aic = sorted_trend_test[0]
-            massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
+            if self.select_only_model_that_pass_anderson_test and \
+                    (not trend_test_that_minimized_aic.goodness_of_fit_anderson_test):
+                    print('not anderson')
+            else:
+                print('ok')
+                massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
+
             # Extract the stationary model that minimized AIC
             stationary_trend_tests_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
                                                          [GumbelVersusGumbel, GevStationaryVersusGumbel]]
diff --git a/extreme_trend/visualizers/utils.py b/extreme_trend/visualizers/utils.py
index f056e22aa07475e882098ed7204488fd2414595f..01c1b76ba2bb56fc7dc4073a96247ca92b59c9e3 100644
--- a/extreme_trend/visualizers/utils.py
+++ b/extreme_trend/visualizers/utils.py
@@ -24,7 +24,8 @@ def load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncer
                                                   model_subsets_for_uncertainty=model_subsets_for_uncertainty,
                                                   fit_method=fit_method, select_only_acceptable_shape_parameter=True,
                                                   fit_gev_only_on_non_null_maxima=False,
-                                                  fit_only_time_series_with_ninety_percent_of_non_null_values=True)
+                                                  fit_only_time_series_with_ninety_percent_of_non_null_values=True,
+                                                  select_only_model_that_pass_anderson_test=True,
+                                                  )
         altitude_to_visualizer[altitude] = study_visualizer
-
     return altitude_to_visualizer
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
index 4fe9ff3ae18cccca9cb44940730fa095d6e19bfe..7d5c0100b544335d80f03ca07d2186fbaefc5261 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
@@ -81,7 +81,7 @@ def intermediate_result(altitudes, massif_names=None,
     visualizers = list(altitude_to_visualizer.values())
     if multiprocessing:
         with Pool(4) as p:
-            _ = p.map(compute_minimized_aic, visualizers)
+            _ = p.imap(compute_minimized_aic, visualizers)
     else:
         for visualizer in visualizers:
             _ = compute_minimized_aic(visualizer)
diff --git a/projects/exceeding_snow_loads/section_data/main_eurocode_plot.py b/projects/exceeding_snow_loads/section_data/main_eurocode_plot.py
index 0cdde26558e08fa5c65965135307e3f7c5aeb56f..850d09150d5a187dbbe969de467a15ee6d4b6932 100644
--- a/projects/exceeding_snow_loads/section_data/main_eurocode_plot.py
+++ b/projects/exceeding_snow_loads/section_data/main_eurocode_plot.py
@@ -30,7 +30,7 @@ def main_eurocode_norms(ax=None, poster_plot=False):
                 ax.set_yticks(list(range(0, 13, 2)))
                 ax.set_xlim([0, 2000])
                 ax.set_xticks(list(range(0, 2001, 500)))
-                prop = {'size': 25} if poster_plot else {}
+                prop = {'size': 20} if poster_plot else {}
                 ax.legend(prop=prop, loc='upper left')
 
                 ax.grid()
diff --git a/projects/exceeding_snow_loads/section_data/main_example_swe_total_plot.py b/projects/exceeding_snow_loads/section_data/main_example_swe_total_plot.py
index 5719e15f7da9022bfd904341d80578f98f9a1428..bb78ef8dfe0f60caa629844f13daeaa043b8cba2 100644
--- a/projects/exceeding_snow_loads/section_data/main_example_swe_total_plot.py
+++ b/projects/exceeding_snow_loads/section_data/main_example_swe_total_plot.py
@@ -55,7 +55,7 @@ def max_graph_annual_maxima_poster():
             last_plot = color == "magenta"
             label = '{} massif at {}m'.format(massif_name, altitude)
             tight_pad = {'h_pad': 0.2}
-            snow_abbreviation = 'max ' + snow_abbreviation
+            snow_abbreviation = 'annual maxima of ' + snow_abbreviation
             study_visualizer.visualize_max_graphs_poster(massif_name, altitude, snow_abbreviation, color, label,
                                                          last_plot, ax, tight_pad=tight_pad,
                                                          dpi=dpi_paper1_figure,
diff --git a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
index 61a2fe2773d054ecc8cfe9bcd95f3eec28675ce0..46f93e5223ea9a8262f0b88bbeee0f63e203478d 100644
--- a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
+++ b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
@@ -5,7 +5,8 @@ import matplotlib as mpl
 from projects.exceeding_snow_loads.checks.qqplot.plot_qqplot import \
     plot_intensity_against_gumbel_quantile_for_3_examples, plot_full_diagnostic
 from projects.exceeding_snow_loads.section_results.plot_selection_curves import plot_selection_curves
-from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map
+from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map, plot_trend_curves
+from projects.exceeding_snow_loads.section_results.plot_uncertainty_histogram import plot_uncertainty_histogram
 
 mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
@@ -33,6 +34,7 @@ def minor_result(altitude):
 
 
 def compute_minimized_aic(visualizer):
+    print(visualizer.altitude)
     _ = visualizer.massif_name_to_trend_test_that_minimized_aic
     return True
 
@@ -58,22 +60,22 @@ def intermediate_result(altitudes, massif_names=None,
     for v in altitude_to_visualizer.values():
         _ = v.study.year_to_variable_object
     # Compute minimized value efficiently
-    visualizers = list(altitude_to_visualizer.values())
+    # visualizers = list()
     if multiprocessing:
         with Pool(NB_CORES) as p:
-            _ = p.map(compute_minimized_aic, visualizers)
+            _ = p.imap(compute_minimized_aic, altitude_to_visualizer.values())
     else:
-        for visualizer in visualizers:
+        for visualizer in altitude_to_visualizer.values():
             _ = compute_minimized_aic(visualizer)
 
     # Plots
-    # plot_trend_map(altitude_to_visualizer)
-    # plot_trend_curves(altitude_to_visualizer={a: v for a, v in altitude_to_visualizer.items() if a >= 900})
-    # plot_uncertainty_massifs(altitude_to_visualizer)
-    # plot_uncertainty_histogram(altitude_to_visualizer)
+    plot_trend_map(altitude_to_visualizer)
+    plot_trend_curves(altitude_to_visualizer={a: v for a, v in altitude_to_visualizer.items() if a >= 900})
+    plot_uncertainty_massifs(altitude_to_visualizer)
+    plot_uncertainty_histogram(altitude_to_visualizer)
     plot_selection_curves(altitude_to_visualizer)
     # uncertainty_interval_size(altitude_to_visualizer)
-    # plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)
+    plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)
     # plot_full_diagnostic(altitude_to_visualizer)
 
 
@@ -89,6 +91,7 @@ def major_result():
     #                                  ModelSubsetForUncertainty.non_stationary_gumbel_and_gev]
     model_subsets_for_uncertainty = None
     # study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
+    altitudes = [300, 600]
     altitudes = paper_altitudes
     # altitudes = [900, 1800, 2700][:1]
     for study_class in study_classes:
diff --git a/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py b/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py
index 04403d6b2d0408eb52924b18113894cab39547a6..61ea9ccc3a3f8ae9669fb159ff3827ba6565f077 100644
--- a/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py
+++ b/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py
@@ -183,9 +183,8 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, leg
         # ax2.legend(handles=legend_elements, bbox_to_anchor=(0.93, 0.7), loc='upper right')
         # ax2.annotate("Filled symbol = significant trend ", xy=(0.85, 0.5), xycoords='axes fraction', fontsize=7)
         upper_legend_y = 0.55
-        ax2.annotate('\n'.join(filled_marker_legend_list2), xy=(0.23, upper_legend_y - 0.2), xycoords='axes fraction', fontsize=fontsize)
-        ax2.annotate('Markers show selected model $\mathcal{M}_N$', xy=(0.02, upper_legend_y), xycoords='axes fraction', fontsize=fontsize)
-        print(legend_size)
+        # ax2.annotate('\n'.join(filled_marker_legend_list2), xy=(0.23, upper_legend_y - 0.2), xycoords='axes fraction', fontsize=fontsize)
+        # ax2.annotate('Markers show selected model $\mathcal{M}_N$', xy=(0.02, upper_legend_y), xycoords='axes fraction', fontsize=fontsize)
         ax2.legend(handles=legend_elements, prop={'size': legend_size},
                    loc='upper left', bbox_to_anchor=(0.00, upper_legend_y))
         # for handle in lgnd.legendHandles: