[contrasting] add a validation for the change in the ratio

extreme_trend/abstract_gev_trend_test.py

@@ -385,10 +385,9 @@ class AbstractGevTrendTest(object):
 
     # Mean value
 
-    @property
-    def unconstrained_average_mean_value(self) -> float:
+    def unconstrained_average_mean_value(self, year_min, year_max) -> float:
         mean_values = []
-        for year in self.years:
+        for year in range(year_min, year_max+1):
             mean = self.get_unconstrained_gev_params(year).mean
             mean_values.append(mean)
         average_mean_value = np.mean(mean_values)

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py

@@ -72,9 +72,10 @@ def intermediate_result(altitudes, massif_names=None,
             _ = compute_minimized_aic(visualizer)
 
     # Plots
-    validation_plot(altitude_to_visualizer)
-    plot_snowfall_mean(altitude_to_visualizer)
-    plot_snowfall_time_derivative_mean(altitude_to_visualizer)
+    validation_plot(altitude_to_visualizer, order_derivative=0)
+    validation_plot(altitude_to_visualizer, order_derivative=1)
+    # plot_snowfall_mean(altitude_to_visualizer)
+    # plot_snowfall_time_derivative_mean(altitude_to_visualizer)
 
 
 def major_result():

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py

@@ -27,6 +27,13 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
                          non_stationary_trend_test_to_marker, fit_method, select_only_acceptable_shape_parameter,
                          fit_gev_only_on_non_null_maxima, fit_only_time_series_with_ninety_percent_of_non_null_values)
 
+
+    def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm, add_x_label=True,
+                           negative_and_positive_values=True):
+        super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap, add_x_label, negative_and_positive_values)
+
+    # Study of the mean
+
     @cached_property
     def massif_name_to_empirical_mean(self):
         massif_name_to_empirical_value = {}
@@ -36,23 +43,48 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
 
     @cached_property
     def massif_name_to_model_mean(self):
-        massif_name_to_parameter_value = {}
+        massif_name_to_model_value = {}
         for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items():
-            parameter_value = trend_test.unconstrained_average_mean_value
-            massif_name_to_parameter_value[massif_name] = parameter_value
-        return massif_name_to_parameter_value
+            parameter_value = trend_test.unconstrained_average_mean_value(self.study.year_min, self.study.year_max)
+            massif_name_to_model_value[massif_name] = parameter_value
+        return massif_name_to_model_value
 
     @cached_property
     def massif_name_to_relative_difference_for_mean(self):
         massif_name_to_relative_difference = {}
         for massif_name in self.massif_name_to_trend_test_that_minimized_aic.keys():
             e = self.massif_name_to_empirical_mean[massif_name]
-            p = self.massif_name_to_model_mean[massif_name]
-            relative_diference = 100 * (p-e) / e
+            m = self.massif_name_to_model_mean[massif_name]
+            relative_diference = 100 * (m-e) / e
             massif_name_to_relative_difference[massif_name] = relative_diference
         return massif_name_to_relative_difference
 
-    def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm, add_x_label=True,
-                           negative_and_positive_values=True):
-        super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap, add_x_label, negative_and_positive_values)
+    # Study of the change in the mean
 
+    @cached_property
+    def massif_name_to_change_ratio_in_empirical_mean(self):
+        massif_name_to_empirical_value = {}
+        for massif_name, maxima in self.study.massif_name_to_annual_maxima.items():
+            index = self.study.ordered_years.index(1989)
+            maxima_before, maxima_after = maxima[:index+1], maxima[index+1:]
+            massif_name_to_empirical_value[massif_name] = np.mean(maxima_after) / np.mean(maxima_before)
+        return massif_name_to_empirical_value
+
+    @cached_property
+    def massif_name_to_change_ratio_in_model_mean(self):
+        massif_name_to_parameter_value = {}
+        for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items():
+            model_mean_before = trend_test.unconstrained_average_mean_value(year_min=self.study.year_min, year_max=1989)
+            model_mean_after = trend_test.unconstrained_average_mean_value(year_min=1990, year_max=self.study.year_max)
+            massif_name_to_parameter_value[massif_name] = model_mean_after / model_mean_before
+        return massif_name_to_parameter_value
+
+    @cached_property
+    def massif_name_to_relative_difference_for_change_ratio_in_mean(self):
+        massif_name_to_relative_difference = {}
+        for massif_name in self.massif_name_to_trend_test_that_minimized_aic.keys():
+            e = self.massif_name_to_change_ratio_in_empirical_mean[massif_name]
+            m = self.massif_name_to_change_ratio_in_model_mean[massif_name]
+            relative_diference = 100 * (m-e) / e
+            massif_name_to_relative_difference[massif_name] = relative_diference
+        return massif_name_to_relative_difference
\ No newline at end of file

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py

@@ -10,34 +10,42 @@ from projects.ogorman.gorman_figures.figure1.study_visualizer_for_double_station
     StudyVisualizerForReturnLevelChange
 
 
-def validation_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
+def validation_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], order_derivative=0):
     # Plot the mean empirical, the mean parametric and the relative difference between the two
     altitudes = list(altitude_to_visualizer.keys())
     study_visualizer = list(altitude_to_visualizer.values())[0]
     altitude_to_relative_differences = {}
+    if order_derivative == 0:
+        plot_function = plot_relative_difference_map_order_zero
+    else:
+        plot_function = plot_relative_difference_map_order_one
     # Plot map for the repartition of the difference
     for altitude, visualizer in altitude_to_visualizer.items():
-        altitude_to_relative_differences[altitude] = plot_relative_difference_map(visualizer)
+        altitude_to_relative_differences[altitude] = plot_function(visualizer)
         study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
     # Shoe plot with respect to the altitude.
-    plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes)
+    plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, study_visualizer, order_derivative)
     study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
     plt.close()
 
 
-def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes):
+def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, visualizer, order_derivative):
+    study = visualizer.study
     ax = plt.gca()
     width = 150
     ax.boxplot([altitude_to_relative_differences[a] for a in altitudes], positions=altitudes, widths=width)
     ax.set_xlim([min(altitudes) - width, max(altitudes) + width])
-    ylabel = 'Relative difference of the model mean w.r.t. the empirical mean (\%)'
+    moment = '' if order_derivative == 0 else 'time derivative of '
+    ylabel = 'Relative difference of the {} model mean \n' \
+             'w.r.t. the {}empirical mean of {} (\%)'.format(moment, moment, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)])
+    visualizer.plot_trends = ylabel
     ax.set_ylabel(ylabel)
     ax.set_xlabel('Altitude (m)')
     ax.legend()
     ax.grid()
 
 
-def plot_relative_difference_map(visualizer: StudyVisualizerForMeanValues):
+def plot_relative_difference_map_order_zero(visualizer: StudyVisualizerForMeanValues):
     study = visualizer.study
     label = ' mean annual maxima of {} ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit)
     visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_empirical_mean,
@@ -48,3 +56,15 @@ def plot_relative_difference_map(visualizer: StudyVisualizerForMeanValues):
                                   label='Relative difference of the model mean w.r.t. the empirical mean \n'
                                         'for the ' + label, graduation=1)
     return list(visualizer.massif_name_to_relative_difference_for_mean.values())
+
+def plot_relative_difference_map_order_one(visualizer: StudyVisualizerForMeanValues):
+    study = visualizer.study
+    label = ' time derivative of mean annual maxima of {} ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit)
+    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_change_ratio_in_empirical_mean,
+                                  label='Empirical' + label, negative_and_positive_values=False, graduation=0.5)
+    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_change_ratio_in_model_mean,
+                                  label='Model' + label, negative_and_positive_values=False, graduation=0.5)
+    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_change_ratio_in_mean,
+                                  label='Relative difference of the model mean w.r.t. the empirical mean \n'
+                                        'for the ' + label, graduation=1)
+    return list(visualizer.massif_name_to_relative_difference_for_mean.values())
\ No newline at end of file