[contrasting] add assertion that the max_abs_change is different from np.inf...

[contrasting] add assertion that the max_abs_change is different from np.inf otherwise it makes the computer freeze
IMPORTANT MODIFICATION, THAT MIGHT CHANGE PAPER 1: change unconstrained_estimator_gev_params to unconstrained_estimator_gev_params_last_year

extreme_data/meteo_france_data/scm_models_data/visualization/create_shifted_cmap.py

@@ -77,6 +77,7 @@ def imshow_shifted(ax, param_name, values, visualization_extend, mask_2D=None):
 
 
 def ticks_values_and_labels_for_percentages(graduation, max_abs_change):
+    assert max_abs_change != np.inf
     positive_ticks = []
     tick = graduation
     while tick < max_abs_change:
@@ -86,7 +87,9 @@ def ticks_values_and_labels_for_percentages(graduation, max_abs_change):
     ticks_values = [((t / max_abs_change) + 1) / 2 for t in all_ticks_labels]
     return ticks_values, all_ticks_labels
 
+
 def ticks_values_and_labels_for_positive_value(graduation, max_abs_change):
+    assert max_abs_change != np.inf
     positive_ticks = []
     tick = 0
     while tick < max_abs_change:

extreme_trend/abstract_gev_trend_test.py

@@ -139,10 +139,12 @@ class AbstractGevTrendTest(object):
                                                                        AbstractCoordinates.COORDINATE_T)
 
     @cached_property
-    def unconstrained_estimator_gev_params(self) -> GevParams:
-        # Constant parameters correspond to the gev params in 1958
-        return self.unconstrained_estimator.function_from_fit.get_params(coordinate=np.array([1958]),
-                                                                         is_transformed=False)
+    def unconstrained_estimator_gev_params_last_year(self) -> GevParams:
+        return self.get_unconstrained_gev_params(year=self.years[-1])
+
+    @cached_property
+    def unconstrained_estimator_gev_params_first_year(self) -> GevParams:
+        return self.get_unconstrained_gev_params(year=self.years[0])
 
     def time_derivative_times_years(self, nb_years):
         # Compute the slope strength
@@ -186,7 +188,7 @@ class AbstractGevTrendTest(object):
 
     @property
     def test_trend_constant_quantile(self):
-        return self.unconstrained_estimator_gev_params.quantile(p=self.quantile_level)
+        return self.unconstrained_estimator_gev_params_last_year.quantile(p=self.quantile_level)
 
     # Some class properties for display purpose
 
@@ -471,7 +473,8 @@ class AbstractGevTrendTest(object):
     def unconstrained_average_mean_value(self, year_min, year_max) -> float:
         mean_values = []
         for year in range(year_min, year_max + 1):
-            mean = self.get_unconstrained_gev_params(year).mean
+            gev_params = self.get_unconstrained_gev_params(year)
+            mean = gev_params.mean
             mean_values.append(mean)
         average_mean_value = np.mean(mean_values)
         assert isinstance(average_mean_value, float)

extreme_trend/trend_test_one_parameter/gev_trend_test_one_parameter.py

@@ -73,12 +73,12 @@ class GevLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
                          fit_method=fit_method)
 
     def _slope_strength(self):
-        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
-                                                                                       mu1=self.non_stationary_linear_coef)
+        return self.unconstrained_estimator_gev_params_last_year.time_derivative_of_return_level(p=self.quantile_level,
+                                                                                                 mu1=self.non_stationary_linear_coef)
 
     @property
     def mean_difference_same_sign_as_slope_strenght(self) -> bool:
-        zeta0 = self.unconstrained_estimator_gev_params.shape
+        zeta0 = self.unconstrained_estimator_gev_params_last_year.shape
         mean_difference = self.mean_difference(zeta0=zeta0, mu1=self.non_stationary_linear_coef)
         return self.same_sign(mean_difference, self._slope_strength())
 
@@ -100,13 +100,13 @@ class GevScaleTrendTest(GevTrendTestOneParameterAgainstStationary):
                          fit_method=fit_method)
 
     def _slope_strength(self):
-        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(
+        return self.unconstrained_estimator_gev_params_last_year.time_derivative_of_return_level(
             p=self.quantile_level,
             sigma1=self.non_stationary_linear_coef)
 
     @property
     def mean_difference_same_sign_as_slope_strenght(self) -> bool:
-        zeta0 = self.unconstrained_estimator_gev_params.shape
+        zeta0 = self.unconstrained_estimator_gev_params_last_year.shape
         mean_difference = self.mean_difference(zeta0=zeta0, sigma1=self.non_stationary_linear_coef)
         return self.same_sign(mean_difference, self._slope_strength())
 

extreme_trend/trend_test_one_parameter/gumbel_trend_test_one_parameter.py

@@ -79,8 +79,8 @@ class GumbelLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
         return 3
 
     def _slope_strength(self):
-        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
-                                                                                       mu1=self.non_stationary_linear_coef)
+        return self.unconstrained_estimator_gev_params_last_year.time_derivative_of_return_level(p=self.quantile_level,
+                                                                                                 mu1=self.non_stationary_linear_coef)
 
     @classproperty
     def label(self):
@@ -103,7 +103,7 @@ class GumbelScaleTrendTest(GevTrendTestOneParameterAgainstStationary):
                          fit_method=fit_method)
 
     def _slope_strength(self):
-        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(
+        return self.unconstrained_estimator_gev_params_last_year.time_derivative_of_return_level(
             p=self.quantile_level,
             sigma1=self.non_stationary_linear_coef)
 

extreme_trend/trend_test_three_parameters/gev_trend_test_three_parameters.py

@@ -50,9 +50,9 @@ class GevLocationAndScaleTrendTestAgainstGumbel(GevTrendTestThreeParameters):
         return self.get_non_stationary_linear_coef(param_name=GevParams.SCALE)
 
     def _slope_strength(self):
-        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
-                                                                                       mu1=self.mu1,
-                                                                                       sigma1=self.sigma1)
+        return self.unconstrained_estimator_gev_params_last_year.time_derivative_of_return_level(p=self.quantile_level,
+                                                                                                 mu1=self.mu1,
+                                                                                                 sigma1=self.sigma1)
 
     @classproperty
     def label(self):

extreme_trend/trend_test_two_parameters/gev_trend_test_two_parameters.py

@@ -66,13 +66,13 @@ class GevLocationAndScaleTrendTest(GevTrendTestTwoParametersAgainstGev):
         return self.get_non_stationary_linear_coef(param_name=GevParams.SCALE)
 
     def _slope_strength(self):
-        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
-                                                                                       mu1=self.mu1,
-                                                                                       sigma1=self.sigma1)
+        return self.unconstrained_estimator_gev_params_last_year.time_derivative_of_return_level(p=self.quantile_level,
+                                                                                                 mu1=self.mu1,
+                                                                                                 sigma1=self.sigma1)
 
     @property
     def mean_difference_same_sign_as_slope_strenght(self) -> bool:
-        zeta0 = self.unconstrained_estimator_gev_params.shape
+        zeta0 = self.unconstrained_estimator_gev_params_last_year.shape
         mean_difference = self.mean_difference(zeta0=zeta0, mu1=self.mu1, sigma1=self.sigma1)
         return self.same_sign(mean_difference, self._slope_strength())
 

extreme_trend/trend_test_two_parameters/gumbel_test_two_parameters.py

@@ -31,9 +31,9 @@ class GumbelLocationAndScaleTrendTest(GevTrendTestTwoParameters):
         return self.get_non_stationary_linear_coef(param_name=GevParams.SCALE)
 
     def _slope_strength(self):
-        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
-                                                                                       mu1=self.mu1,
-                                                                                       sigma1=self.sigma1)
+        return self.unconstrained_estimator_gev_params_last_year.time_derivative_of_return_level(p=self.quantile_level,
+                                                                                                 mu1=self.mu1,
+                                                                                                 sigma1=self.sigma1)
     @classproperty
     def label(self):
         return super().label % '\\mu_1, \\sigma_1'

extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py

@@ -169,7 +169,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
             if self.select_only_acceptable_shape_parameter:
                 acceptable_shape_parameter = lambda s: -0.5 <= s <= 0.5  # physically acceptable prior
                 all_trend_test = [t for t in all_trend_test
-                                  if acceptable_shape_parameter(t.unconstrained_estimator_gev_params.shape)]
+                                  if (acceptable_shape_parameter(t.unconstrained_estimator_gev_params_last_year.shape)
+                                  and acceptable_shape_parameter(t.unconstrained_estimator_gev_params_first_year.shape))]
             sorted_trend_test = sorted(all_trend_test, key=lambda t: t.aic)
 
             # Extract the stationary or non-stationary model that minimized AIC
@@ -424,14 +425,14 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         psnow = self.massif_name_to_psnow[massif_name]
         trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
         trend_test.intensity_plot_wrt_standard_gumbel_simple(massif_name, self.altitude, color, psnow)
-        self.plot_name = 'intensity_plot_{}_{}_{}_{}'.format(self.altitude, massif_name, psnow, trend_test.unconstrained_estimator_gev_params.shape)
+        self.plot_name = 'intensity_plot_{}_{}_{}_{}'.format(self.altitude, massif_name, psnow, trend_test.unconstrained_estimator_gev_params_last_year.shape)
         self.show_or_save_to_file(add_classic_title=False, no_title=True)
         plt.close()
 
     def qqplot(self, massif_name, color=None):
         trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
         trend_test.qqplot_wrt_standard_gumbel(massif_name, self.altitude)
-        self.plot_name = 'qpplot_plot_{}_{}_{}'.format(self.altitude, massif_name, trend_test.unconstrained_estimator_gev_params.shape)
+        self.plot_name = 'qpplot_plot_{}_{}_{}'.format(self.altitude, massif_name, trend_test.unconstrained_estimator_gev_params_last_year.shape)
         self.show_or_save_to_file(add_classic_title=False, no_title=True)
         plt.close()
 

extreme_trend/visualizers/study_visualizer_for_shape_repartition.py

@@ -14,7 +14,7 @@ class StudyVisualizerForShape(StudyVisualizerForNonStationaryTrends):
 
     @cached_property
     def massif_name_to_unconstrained_shape_parameter(self):
-        return {m: t.unconstrained_estimator_gev_params.shape
+        return {m: t.unconstrained_estimator_gev_params_last_year.shape
                 for m, t in self.massif_name_to_trend_test_that_minimized_aic.items()}
 
     @cached_property

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py

@@ -88,13 +88,14 @@ def major_result():
     model_subsets_for_uncertainty = None
     altitudes = paper_altitudes
     # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
-    altitudes = [900, 1200, 1500, 1800]
-    # draft_altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
+    # altitudes = [900, 1200, 1500, 1800][:2]
+    # altitudes = [2400, 2700][:2]
+    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
     # altitudes = draft_altitudes
 
     for study_class in study_classes:
         intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
-                            uncertainty_methods, study_class)
+                            uncertainty_methods, study_class, multiprocessing=True)
 
 
 if __name__ == '__main__':

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py

@@ -77,19 +77,18 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
         return massif_name_to_empirical_value
 
     @cached_property
-    def massif_name_to_model_mean(self):
-        massif_name_to_model_value = {}
+    def massif_name_to_model_mean_last_year(self):
+        massif_name_to_model_value_last_year = {}
         for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items():
-            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
+            massif_name_to_model_value_last_year[massif_name] = trend_test.unconstrained_estimator_gev_params_last_year.mean
+        return massif_name_to_model_value_last_year
 
     @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]
-            m = self.massif_name_to_model_mean[massif_name]
+            m = self.massif_name_to_model_mean_last_year[massif_name]
             relative_diference = 100 * (m - e) / e
             massif_name_to_relative_difference[massif_name] = relative_diference
         return massif_name_to_relative_difference

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/validation_plot.py

+from multiprocessing import Pool
 from typing import Dict
 
 import matplotlib.pyplot as plt
@@ -8,25 +9,27 @@ from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and
     StudyVisualizerForMeanValues
 from projects.ogorman.gorman_figures.figure1.study_visualizer_for_double_stationary_fit import \
     StudyVisualizerForReturnLevelChange
+from root_utils import NB_CORES
 
 
 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 = {}
+    visualizers = list(altitude_to_visualizer.values())
     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
+    altitude_to_relative_differences = {}
     for altitude, visualizer in altitude_to_visualizer.items():
         altitude_to_relative_differences[altitude] = plot_function(visualizer)
-        study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
+        # 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, study_visualizer,
     #                                             order_derivative)
-    study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
+    # study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
     plt.close()
 
 
@@ -53,11 +56,12 @@ def plot_relative_difference_map_order_zero(visualizer: StudyVisualizerForMeanVa
     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,
     #                               label='Empirical' + label, negative_and_positive_values=False)
-    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_model_mean,
+    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_model_mean_last_year,
                                   label='Model' + label, negative_and_positive_values=False, add_text=True)
     # visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_mean,
     #                               label='Relative difference of the model mean w.r.t. the empirical mean \n'
     #                                     'for the ' + label, graduation=1)
+    visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
     return list(visualizer.massif_name_to_relative_difference_for_mean.values())
 
 

projects/exceeding_snow_loads/checks/qqplot/main_qqplot_for_big_shapes.py

@@ -23,7 +23,7 @@ def get_tuple_ordered_by_shape(fast=False):
     # Extract all the values
     l = []
     for a, v in altitude_to_visualizer.items():
-        l.extend([(a, v, m, t.unconstrained_estimator_gev_params.shape) for m, t in
+        l.extend([(a, v, m, t.unconstrained_estimator_gev_params_last_year.shape) for m, t in
                   v.massif_name_to_trend_test_that_minimized_aic.items()])
     # Sort them and keep the highest examples
     l = sorted(l, key=lambda t: t[-1])

projects/exceeding_snow_loads/checks/without_maximum/study_visualizer_for_fit_witout_maximum.py

@@ -40,7 +40,7 @@ class StudyVisualizerForFitWithoutMaximum(StudyVisualizerForNonStationaryTrends)
         for massif_name, maximum in self.massif_name_to_maximum.items():
             t = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
             msg = '{} {}m'.format(massif_name, self.study.altitude)
-            upper_bound = t.unconstrained_estimator_gev_params.density_upper_bound
+            upper_bound = t.unconstrained_estimator_gev_params_last_year.density_upper_bound
             if not np.isinf(upper_bound):
                 diff.append(upper_bound - maximum)
             assert maximum <= upper_bound, msg

projects/exceeding_snow_loads/section_results/plot_selection_curves.py

@@ -2,7 +2,7 @@ from typing import Dict
 import matplotlib.pyplot as plt
 
 from extreme_data.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
-from projects.exceeding_snow_loads.utils import dpi_paper1_figure
+from projects.exceeding_snow_loads.utils import dpi_paper1_figure, get_trend_test_name
 from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
 
 
@@ -10,7 +10,8 @@ def permute(l, permutation):
     # permutation = [i//2  if i % 2 == 0 else 4 + i //2 for i in range(8)]
     return [l[i] for i in permutation]
 
-def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
+                          paper1=True):
     """
     Plot a single trend curves
     :return:
@@ -20,8 +21,7 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
     ax = create_adjusted_axes(1, 1)
 
     selected_counter = merge_counter([v.selected_trend_test_class_counter for v in altitude_to_visualizer.values()])
-    # selected_and_significative_counter = merge_counter([v.selected_and_significative_trend_test_class_counter for v in altitude_to_visualizer.values()])
-    selected_and_significative_counter = merge_counter([v.selected_and_anderson_goodness_of_fit_trend_test_class_counter for v in altitude_to_visualizer.values()])
+    selected_and_significative_counter = merge_counter([v.selected_and_significative_trend_test_class_counter for v in altitude_to_visualizer.values()])
     # selected_and_significative_counter = merge_counter([v.selected_and_kstest_goodness_of_fit_trend_test_class_counter for v in altitude_to_visualizer.values()])
     total_of_selected_models = sum(selected_counter.values())
     l = sorted(enumerate(visualizer.non_stationary_trend_test), key=lambda e: selected_counter[e[1]])
@@ -29,7 +29,11 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
 
     select_list = get_ordered_list_from_counter(selected_counter, total_of_selected_models, visualizer, permutation)
     selected_and_signifcative_list = get_ordered_list_from_counter(selected_and_significative_counter, total_of_selected_models, visualizer, permutation)
-    labels = permute(['${}$'.format(t.label) for t in visualizer.non_stationary_trend_test], permutation)
+    if paper1:
+        labels = ['${}$'.format(t.label) for t in visualizer.non_stationary_trend_test]
+    else:
+        labels = [get_trend_test_name(t) for t in visualizer.non_stationary_trend_test]
+    labels = permute(labels, permutation)
 
     print(l)
     print(sum(select_list), select_list)
@@ -42,7 +46,10 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
     width = 5
     size = 30
     legend_fontsize = 30
-    labelsize = 25
+    if paper1:
+        labelsize = 25
+    else:
+        labelsize = 15
     linewidth = 3
     x = [10 * (i+1) for i in range(len(select_list))]
     ax.bar(x, select_list, width=width, color='grey', edgecolor='grey', label='Non significant model',
@@ -95,8 +102,18 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
     #     ax_twinx.plot(altitudes, mean_decrease, label=label, linewidth=linewidth, marker='o')
     #     ax_twinx.legend(loc='upper right', prop={'size': size})
 
+    # Compute the ratio of plots that pass the goodness of fit test
+    if paper1:
+        ratio = '100'
+    else:
+        nb_selected_total = 0
+        nb_selected_total_and_goodness_of_fit = 0
+        for v in altitude_to_visualizer.values():
+            nb_selected_total += len(v.super_massif_name_to_trend_test_that_minimized_aic)
+            nb_selected_total_and_goodness_of_fit += len(v.massif_name_to_trend_test_that_minimized_aic)
+        ratio = 100 * nb_selected_total_and_goodness_of_fit / nb_selected_total
     # Save plot
-    visualizer.plot_name = 'Selection curves'
+    visualizer.plot_name = 'Selection curves ({}% of selected models were kept)'.format(ratio)
     visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
     plt.close()