[PAPER 1] create relative change plot

experiment/meteo_france_data/plot/create_shifted_cmap.py

@@ -32,7 +32,7 @@ def create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=None):
     if ticks_values_and_labels is not None:
         cb.ax.set_yticklabels([str(t) for t in ticks_values_and_labels[1]])
     if isinstance(label, str):
-        cb.set_label(label)
+        cb.set_label(label, fontsize=15)
     return norm
 
 

experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -17,14 +17,14 @@ from root_utils import NB_CORES
 
 mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
-
+import matplotlib.pyplot as plt
 
 def minor_result(altitude):
     """Plot trends for a single altitude to be fast"""
     visualizer = StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude), multiprocessing=True,
                                                        )
     visualizer.plot_trends()
-    # plt.show()
+    plt.show()
 
 def compute_minimized_aic(visualizer):
     _ = visualizer.massif_name_to_minimized_aic_non_stationary_trend_test
@@ -50,7 +50,7 @@ def intermediate_result(altitudes, massif_names=None,
     # Load variable object efficiently
     for v in altitude_to_visualizer.values():
         _ = v.study.year_to_variable_object
-    # Plot trends
+    # Compute minimized value efficiently
     visualizers = list(altitude_to_visualizer.values())
     if multiprocessing:
         with Pool(NB_CORES) as p:
@@ -59,9 +59,13 @@ def intermediate_result(altitudes, massif_names=None,
         for visualizer in visualizers:
             _ = compute_minimized_aic(visualizer)
     # Compute common max value for the colorbar
-    max_abs_tdrl = max([visualizer.max_abs_tdrl for visualizer in visualizers])
-    for visualizer in altitude_to_visualizer.values():
-        visualizer.plot_trends(max_abs_tdrl)
+    altitudes_for_plot_trend = [900, 1800, 2700]
+    visualizers_for_altitudes = [visualizer
+                                 for altitude, visualizer in altitude_to_visualizer.items()
+                                 if altitude in altitudes_for_plot_trend]
+    max_abs_tdrl = max([visualizer.max_abs_change for visualizer in visualizers_for_altitudes])
+    for visualizer in visualizers_for_altitudes:
+        visualizer.plot_trends(max_abs_tdrl, add_colorbar=visualizer.study.altitude==2700)
 
     # Plot graph
     plot_uncertainty_massifs(altitude_to_visualizer)
@@ -73,7 +77,7 @@ def major_result():
     uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
                            ConfidenceIntervalMethodFromExtremes.ci_mle][:]
     massif_names = None
-    for study_class in paper_study_classes[1:2]:
+    for study_class in paper_study_classes[:1]:
         if study_class == CrocusSnowLoadEurocode:
             non_stationary_uncertainty = [False]
         else:
@@ -82,16 +86,21 @@ def major_result():
 
 
 if __name__ == '__main__':
-    major_result()
+    # major_result()
     # intermediate_result(altitudes=paper_altitudes, massif_names=['Maurienne'],
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
     #                        ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
     #                     non_stationary_uncertainty=[False, True][1:],
     #                     multiprocessing=True)
+    intermediate_result(altitudes=paper_altitudes, massif_names=['Maurienne'],
+                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+                           ConfidenceIntervalMethodFromExtremes.ci_mle][:],
+                        non_stationary_uncertainty=[False, True][:],
+                        multiprocessing=True)
     # intermediate_result(altitudes=[900, 1200], massif_names=None)
     # intermediate_result(ALL_ALTITUDES_WITHOUT_NAN)
     # intermediate_result(paper_altitudes)
-    # minor_result(altitude=600)
+    # minor_result(altitude=900)
     # intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
     #                                          ConfidenceIntervalMethodFromExtremes.ci_bayes],

experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py

@@ -138,7 +138,9 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, leg
     if len(visualizers) > 0:
         tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in visualizers]
         # Plot bars
-        colors = [v.massif_name_to_tdrl_color[massif_name] for v in visualizers]
+        # colors = [v.massif_name_to_color[massif_name] for v in visualizers]
+        # From snow on, we set a black color for the bars
+        colors = ['white' for v in visualizers]
         ax.bar(valid_altitudes, tdrl_values, width=150, color=colors, label=visualizers[0].label_tdrl_bar,
                edgecolor='black', hatch='//')
         # Plot markers

experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -38,12 +38,14 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                  uncertainty_methods=None,
                  non_stationary_contexts=None,
                  uncertainty_massif_names=None,
-                 effective_temporal_covariate=2017):
+                 effective_temporal_covariate=2017,
+                 relative_change_trend_plot=True):
         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, score_class)
         # Add some attributes
+        self.relative_change_trend_plot = relative_change_trend_plot
         self.effective_temporal_covariate = effective_temporal_covariate
         self.non_stationary_contexts = non_stationary_contexts
         self.uncertainty_methods = uncertainty_methods
@@ -60,7 +62,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         self.non_stationary_trend_test = [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest]
         self.non_stationary_trend_test_to_marker = dict(zip(self.non_stationary_trend_test, ["s", "^", "D"]))
 
-        self.global_max_abs_tdrl = None
+        self.global_max_abs_change = None
 
     # Utils
 
@@ -105,53 +107,64 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
             massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
         return massif_name_to_trend_test_that_minimized_aic
 
-
     # Part 1 - Trends
 
     @property
-    def max_abs_tdrl(self):
-        return max(abs(min(self.massif_name_to_tdrl_value.values())), max(self.massif_name_to_tdrl_value.values()))
+    def max_abs_change(self):
+        return max(abs(min(self.massif_name_to_change_value.values())), max(self.massif_name_to_change_value.values()))
 
     @cached_property
-    def _max_abs_tdrl(self):
-        return self.global_max_abs_tdrl if self.global_max_abs_tdrl is not None else self.max_abs_tdrl
+    def _max_abs_change(self):
+        return self.global_max_abs_change if self.global_max_abs_change is not None else self.max_abs_change
 
-    def plot_trends(self, max_abs_tdrl=None):
+    def plot_trends(self, max_abs_tdrl=None,  add_colorbar=True):
         if max_abs_tdrl is not None:
-            self.global_max_abs_tdrl = max_abs_tdrl
-        ax = self.study.visualize_study(massif_name_to_value=self.massif_name_to_tdrl_value,
+            self.global_max_abs_change = max_abs_tdrl
+        ax = self.study.visualize_study(massif_name_to_value=self.massif_name_to_change_value,
                                         replace_blue_by_white=False,
                                         axis_off=False, show_label=False,
-                                        add_colorbar=True,
+                                        add_colorbar=add_colorbar,
                                         massif_name_to_marker_style=self.massif_name_to_marker_style,
-                                        massif_name_to_color=self.massif_name_to_tdrl_color,
+                                        massif_name_to_color=self.massif_name_to_color,
                                         cmap=self.cmap,
                                         show=False,
                                         ticks_values_and_labels=self.ticks_values_and_labels,
-                                        label=self.label_tdrl_bar + ' for {}'.format(EUROCODE_RETURN_LEVEL_STR))
+                                        label=self.label)
         ax.get_xaxis().set_visible(True)
         ax.set_xticks([])
-        ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=12)
+        ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=15)
 
         self.plot_name = 'tdlr_trends'
         self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
-                                  dpi=dpi_paper1_figure)
+                                  dpi=500)
         plt.close()
 
+    @property
+    def label(self):
+        if self.relative_change_trend_plot:
+            label_tdlr_bar = 'Relative change between {} and {}'.format(self.initial_year, self.final_year)
+        else:
+            label_tdlr_bar = self.label_tdrl_bar
+        label = label_tdlr_bar + '\nfor {}'.format(EUROCODE_RETURN_LEVEL_STR)
+        if self.relative_change_trend_plot:
+            # change units
+            label = label.split('(')[0] + '(\%)'
+        return label
+
     @property
     def label_tdrl_bar(self):
         return 'Change in {} years'.format(AbstractGevTrendTest.nb_years_for_quantile_evolution)
 
     @property
     def ticks_values_and_labels(self):
-        graduation = 0.2
+        graduation = 10 if self.relative_change_trend_plot else 0.2
         positive_ticks = []
         tick = graduation
-        while tick < self._max_abs_tdrl:
+        while tick < self._max_abs_change:
             positive_ticks.append(round(tick, 1))
             tick += graduation
         all_ticks_labels = [-t for t in positive_ticks] + [0] + positive_ticks
-        ticks_values = [((t / self._max_abs_tdrl) + 1) / 2 for t in all_ticks_labels]
+        ticks_values = [((t / self._max_abs_change) + 1) / 2 for t in all_ticks_labels]
         return ticks_values, all_ticks_labels
 
     @cached_property
@@ -159,14 +172,34 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         return {m: t.time_derivative_of_return_level for m, t in
                 self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
 
+    @cached_property
+    def massif_name_to_relative_change_value(self):
+        return {m: t.relative_change_in_return_level(initial_year=self.initial_year, final_year=self.final_year)
+                for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+
+    @property
+    def initial_year(self):
+        return self.final_year - 50
+
+    @property
+    def final_year(self):
+        return 2010
+
+    @cached_property
+    def massif_name_to_change_value(self):
+        if self.relative_change_trend_plot:
+            return self.massif_name_to_relative_change_value
+        else:
+            return self.massif_name_to_tdrl_value
+
     @cached_property
     def cmap(self):
-        return get_shifted_map(-self._max_abs_tdrl, self._max_abs_tdrl)
+        return get_shifted_map(-self._max_abs_change, self._max_abs_change)
 
     @cached_property
-    def massif_name_to_tdrl_color(self):
-        return {m: get_colors([v], self.cmap, -self._max_abs_tdrl, self._max_abs_tdrl)[0]
-                for m, v in self.massif_name_to_tdrl_value.items()}
+    def massif_name_to_color(self):
+        return {m: get_colors([v], self.cmap, -self._max_abs_change, self._max_abs_change)[0]
+                for m, v in self.massif_name_to_change_value.items()}
 
     @cached_property
     def massif_name_to_marker_style(self):

experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py

@@ -143,17 +143,31 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
         return self.constrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([1958]),
                                                                                   is_transformed=False)
 
+    def time_derivative_times_years(self, nb_years):
+        # Compute the slope strength
+        slope = self._slope_strength()
+        # Delta T must in the same unit as were the parameter of slope mu1 and sigma1
+        slope *= nb_years * self.coordinates.transformed_distance_between_two_successive_years[0]
+        return slope
+
     @property
     def time_derivative_of_return_level(self):
         if self.crashed:
             return 0.0
         else:
-            # Compute the slope strength
-            slope = self._slope_strength()
-            # Delta T must in the same unit as were the parameter of slope mu1 and sigma1
-            slope *= self.nb_years_for_quantile_evolution * \
-                     self.coordinates.transformed_distance_between_two_successive_years[0]
-            return slope
+            return self.time_derivative_times_years(self.nb_years_for_quantile_evolution)
+
+    def relative_change_in_return_level(self, initial_year, final_year):
+        return_level_values = []
+        for year in [initial_year, final_year]:
+            gev_params = self.unconstrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([year]),
+                                                                                            is_transformed=False)
+            return_level_values.append(gev_params.quantile(self.quantile_level))
+        change_until_final_year = self.time_derivative_times_years(nb_years=final_year - initial_year)
+        change_in_between = (return_level_values[1] - return_level_values[0])
+        np.testing.assert_almost_equal(change_until_final_year, change_in_between, decimal=5)
+        initial_return_level = return_level_values[0]
+        return 100 * change_until_final_year / initial_return_level
 
     def _slope_strength(self):
         raise NotImplementedError