[paper 1] fix main plots for the paper

papers/exceeding_snow_loads/check_mle_convergence_for_trends/shape/main_shape_repartition.py

@@ -27,8 +27,8 @@ def main_shape_repartition(altitudes, massif_names=None,
 
 if __name__ == '__main__':
     # main_shape_repartition([900], save_to_file=False)
-    # main_shape_repartition([900, 1800, 2700])
+    main_shape_repartition([900, 1800, 2700])
     # main_shape_repartition([300, 600, 900, 1200, 1500, 1800, 2700])
-    main_shape_repartition(paper_altitudes, study_visualizer_class=StudyVisualizerForShape, save_to_file=True)
+    # main_shape_repartition(paper_altitudes, study_visualizer_class=StudyVisualizerForShape, save_to_file=True)
     # main_shape_repartition([300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200],
     #                        study_visualizer_class=StudyVisualizerForShape, save_to_file=True)

papers/exceeding_snow_loads/paper_main_utils.py

@@ -11,16 +11,13 @@ def load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncer
                                 study_visualizer_class=StudyVisualizerForNonStationaryTrends,
                                 save_to_file=True):
     fit_method = TemporalMarginFitMethod.extremes_fevd_mle
-    select_only_acceptable_shape_parameter = True
-    print('Fit method: {}, Select only acceptable shape parameter: {}'
-          .format(fit_method, select_only_acceptable_shape_parameter))
     altitude_to_visualizer = OrderedDict()
     for altitude in altitudes:
         altitude_to_visualizer[altitude] = study_visualizer_class(
             study=study_class(altitude=altitude), multiprocessing=True, save_to_file=save_to_file,
             uncertainty_massif_names=massif_names, uncertainty_methods=uncertainty_methods,
             model_subsets_for_uncertainty=model_subsets_for_uncertainty, fit_method=fit_method,
-            select_only_acceptable_shape_parameter=select_only_acceptable_shape_parameter)
+            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)
     return altitude_to_visualizer
-
-

papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -69,7 +69,7 @@ def intermediate_result(altitudes, massif_names=None,
     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_uncertainty_histogram(altitude_to_visualizer)
     # plot_selection_curves(altitude_to_visualizer)
 
 
@@ -77,7 +77,7 @@ def major_result():
     uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
                            ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
     massif_names = None
-    study_classes = paper_study_classes[:1]
+    study_classes = paper_study_classes[:2]
     # model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
     #                                  ModelSubsetForUncertainty.stationary_gumbel_and_gev,
     #                                  ModelSubsetForUncertainty.non_stationary_gumbel,
@@ -90,11 +90,11 @@ def major_result():
 
 
 if __name__ == '__main__':
-    # major_result()
-    intermediate_result(altitudes=[1800], massif_names=None,
-                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
-                                             ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
-                        multiprocessing=True)
+    major_result()
+    # intermediate_result(altitudes=[1800], massif_names=None,
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+    #                                          ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+    #                     multiprocessing=True)
     # intermediate_result(altitudes=[900, 1200], massif_names=['Maurienne'],
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
     #                                          ConfidenceIntervalMethodFromExtremes.ci_mle][1:],

papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py

@@ -31,12 +31,12 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
     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)
 
-    # print(l)
-    # print(select_list)
-    # print(selected_and_signifcative_list)
-    # [(5, <class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters.GevLocationAgainstGumbel'> ), (6, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters.GevScaleAgainstGumbel' > ), (2, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelScaleTrendTest' > ), (1, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelLocationTrendTest' > ), (7, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_three_parameters.gev_trend_test_three_parameters.GevLocationAndScaleTrendTestAgainstGumbel' > ), (3, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gumbel_test_two_parameters.GumbelLocationAndScaleTrendTest' > ), (4, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GevStationaryVersusGumbel' > ), (0, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelVersusGumbel' > )]
-    # [32.520325203252035, 25.609756097560975, 12.195121951219512, 9.34959349593496, 9.34959349593496, 4.471544715447155, 3.252032520325203, 3.252032520325203]
-    # [0, 13.821138211382113, 7.317073170731708, 8.536585365853659, 5.691056910569106, 1.6260162601626016, 1.2195121951219512, 2.4390243902439024]
+    print(l)
+    print(select_list)
+    print(selected_and_signifcative_list)
+    # [(5, <    class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters.GevLocationAgainstGumbel'> ), (6, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters.GevScaleAgainstGumbel' > ), (2, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelScaleTrendTest' > ), (1, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelLocationTrendTest' > ), (7, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_three_parameters.gev_trend_test_three_parameters.GevLocationAndScaleTrendTestAgainstGumbel' > ), (3, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gumbel_test_two_parameters.GumbelLocationAndScaleTrendTest' > ), (4, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GevStationaryVersusGumbel' > ), (0, < class 'experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelVersusGumbel' > )]
+    # [32.64462809917355, 24.380165289256198, 12.396694214876034, 9.50413223140496, 9.090909090909092, 5.785123966942149, 3.71900826446281, 2.479338842975207]
+    # [0, 13.223140495867769, 7.851239669421488, 8.264462809917354, 4.958677685950414, 2.479338842975207, 0.8264462809917356, 2.0661157024793386]
 
     # parameters
     width = 5

papers/exceeding_snow_loads/result_trends_and_return_levels/plot_trend_curves.py

@@ -17,7 +17,7 @@ def plot_trend_map(altitude_to_visualizer):
 
     for altitude, visualizer in altitude_to_visualizer.items():
         if 900 <= altitude <= 4200:
-            add_color = (visualizer.study.altitude - 1500) % 900 == 0
+            add_color = (visualizer.study.altitude - 1500) % 1200 == 0
             visualizer.plot_trends(max_abs_tdrl_above_900, add_colorbar=add_color)
             # Plot 2700 also with a colorbar
             if altitude == 2700:
@@ -52,7 +52,7 @@ def plot_trend_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNonSta
     ax.bar(altitudes, percent_decrease, width=width, color=color, edgecolor='blue', label='decreasing trend',
            linewidth=linewidth)
     ax.bar(altitudes, percent_decrease_signi, width=width, color=color, edgecolor='black',
-           label='significative decreasing trend',
+           label='significant decreasing trend',
            linewidth=linewidth)
     ax.legend(loc='upper left', prop={'size': size})
 

papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py

@@ -151,7 +151,7 @@ def add_title(ax, eurocode_region, massif_name, non_stationary_context):
 
 def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, legend_size, fontsize):
     visualizers = [v for a, v in altitude_to_visualizer.items()
-                   if a in valid_altitudes and massif_name in v.uncertainty_massif_names]
+                   if a in valid_altitudes and massif_name in v.massif_names_fitted]
     if len(visualizers) > 0:
         tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in visualizers]
         # Plot bars

papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py

@@ -46,7 +46,6 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
             width = 200
         plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width=width)
 
-
     ax.set_xticks(altitudes)
     ax.tick_params(labelsize=fontsize_label)
     if not (len(visualizer.uncertainty_methods) == 1
@@ -57,19 +56,36 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
     ax.set_ylabel('Massifs exceeding French standards (\%)', fontsize=fontsize_label)
     ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
     ax.set_ylim([0, 100])
+    ax.yaxis.grid()
+
+    ax_twiny = ax.twiny()
+    ax_twiny.plot(altitudes, [0 for _ in altitudes], linewidth=0)
+    ax_twiny.tick_params(labelsize=fontsize_label)
+    ax_twiny.set_xlim(ax.get_xlim())
+    ax_twiny.set_xticks(altitudes)
+    nb_massif_names = [v.study.nb_study_massif_names for v in altitude_to_visualizer.values()]
+    print(nb_massif_names)
+    ax_twiny.set_xticklabels(nb_massif_names)
+    ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=fontsize_label)
+
     ax.set_yticks([10 * i for i in range(11)])
-    visualizer.plot_name = 'Percentages of exceedance with {}'.format(get_display_name_from_object_type(model_subset_for_uncertainty))
+    visualizer.plot_name = 'Percentages of exceedance with {}'.format(
+        get_display_name_from_object_type(model_subset_for_uncertainty))
     # visualizer.show = True
     visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
     ax.clear()
+    ax_twiny.clear()
+    plt.close()
 
 
 def plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width):
     three_percentages_of_excess = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[:3] for v in
                                    visualizers]
     epsilon = 0.5
-    three_percentages_of_excess = [(a, b, c) if a == b else (max(epsilon, a), b, c) for (a, b, c) in three_percentages_of_excess]
-    three_percentages_of_excess = [(a, b, c) if b == c else (a, b, min(100 - epsilon, c)) for (a, b, c) in three_percentages_of_excess]
+    three_percentages_of_excess = [(a, b, c) if a == b else (max(epsilon, a), b, c) for (a, b, c) in
+                                   three_percentages_of_excess]
+    three_percentages_of_excess = [(a, b, c) if b == c else (a, b, min(100 - epsilon, c)) for (a, b, c) in
+                                   three_percentages_of_excess]
     y = [d[1] for d in three_percentages_of_excess]
     yerr = np.array([[d[1] - d[0], d[2] - d[1]] for d in three_percentages_of_excess]).transpose()
     label = ci_method_to_label[ci_method]

papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -117,6 +117,10 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         else:
             return {m: EUROCODE_QUANTILE for m in self.massif_name_to_psnow.keys()}
 
+    @property
+    def massif_names_fitted(self):
+        return list(self.massif_name_to_years_and_maxima_for_model_fitting.keys())
+
     @cached_property
     def massif_name_to_years_and_maxima_for_model_fitting(self):
         if self.fit_gev_only_on_non_null_maxima:
@@ -162,6 +166,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                 all_trend_test = [t for t in all_trend_test
                                   if acceptable_shape_parameter(t.unconstrained_estimator_gev_params.shape)]
             sorted_trend_test = sorted(all_trend_test, key=lambda t: t.aic)
+
             # 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
@@ -175,6 +180,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                                                     [GumbelVersusGumbel, GumbelLocationTrendTest, GumbelScaleTrendTest,
                                                      GumbelLocationAndScaleTrendTest]][0]
             massif_name_to_gumbel_trend_test_that_minimized_aic[massif_name] = gumbel_trend_test_that_minimized_aic
+
         return massif_name_to_trend_test_that_minimized_aic, massif_name_to_stationary_trend_test_that_minimized_aic, massif_name_to_gumbel_trend_test_that_minimized_aic
 
     # Part 1 - Trends
@@ -342,9 +348,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                 res = p.starmap(compute_eurocode_confidence_interval, arguments)
         else:
             res = [compute_eurocode_confidence_interval(*argument) for argument in arguments]
-        massif_name_to_eurocode_return_level_uncertainty = dict(zip(self.uncertainty_massif_names, res))
+        massif_name_to_eurocode_return_level_uncertainty = dict(zip(massifs_names, res))
         # For the rest of the massif names. Create a Eurocode Return Level Uncertainty as nan
-        for massif_name in set(self.study.all_massif_names) - set(self.uncertainty_massif_names):
+        for massif_name in set(self.study.all_massif_names) - set(massifs_names):
             massif_name_to_eurocode_return_level_uncertainty[massif_name] = self.default_eurocode_uncertainty
         return massif_name_to_eurocode_return_level_uncertainty