[PAPER 1] improve uncertainty plot. fix eurocode value for the histogram....

[PAPER 1] improve uncertainty plot. fix eurocode value for the histogram. modify number of bayesian iterations to 20000 for uncertainty plot

experiment/eurocode_data/eurocode_region.py

@@ -18,9 +18,9 @@ class AbstractEurocodeRegion(object):
             return max(self.sad, valeur_caracteritique)
 
     def valeur_caracteristique(self, altitude):
-        return self.sk0 + self.lois_de_variation_de_la_valeur_caracteristique(altitude)
+        return self.sk0 + self._lois_de_variation_de_la_valeur_caracteristique(altitude)
 
-    def lois_de_variation_de_la_valeur_caracteristique(self, altitude):
+    def _lois_de_variation_de_la_valeur_caracteristique(self, altitude):
         if 200 <= altitude <= 2000:
             if 200 <= altitude <= 500:
                 a, b = self.lois_de_variation_200_and_500

experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py

@@ -634,7 +634,6 @@ class StudyVisualizer(VisualizationParameters):
         ax.plot(x, y, color=color, linewidth=5, label=label, linestyle=linestyle)
         # ax.set_ylabel('{} (in {})'.format(snow_abbreviation, self.study.variable_unit), color=color, fontsize=15)
 
-
         ax.xaxis.set_ticks(x[2::10])
         ax.tick_params(axis='both', which='major', labelsize=13)
         plot_name = 'Annual maxima of {} in {} at {}m'.format(snow_abbreviation, massif_name, altitude)

experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode.py

@@ -26,7 +26,13 @@ def max_graph_annual_maxima_comparison():
                      # CrocusSnowDepthAtMaxofSwe,
                      CrocusSnowDepthDifference,
                      ][:]
+    study_class_to_ylim_and_yticks = {
+        CrocusSnowDensityAtMaxofSwe: ([100, 500], [50*i for i in range(2, 11)]),
+        CrocusDifferenceSnowLoad: ([0, 12], [2*i for i in range(0, 7)]),
+        CrocusSnowDepthDifference: ([0, 1], [0.2*i for i in range(0, 6)]),
+    }
     for study_class in study_classes:
+        ylim, yticks = study_class_to_ylim_and_yticks[study_class]
 
         marker_altitude_massif_name = [
             ('magenta', 900, 'Ubaye'),
@@ -45,13 +51,15 @@ def max_graph_annual_maxima_comparison():
                                                              False, ax)
                 last_plot = altitude == 2700
                 if last_plot:
-                    constant = 150 if study_class == CrocusSnowDensityAtMaxofSwe else 0
-                    label = '{} Eurocode'.format(
-                        snow_density_str) if study_class == CrocusSnowDensityAtMaxofSwe else None
-                    snow_density_eurocode = [constant for _ in study.ordered_years]
-                    ax.plot(study.ordered_years, snow_density_eurocode, color='k', label=label)
+                    if study_class == CrocusSnowDensityAtMaxofSwe:
+                        label = '{} Eurocode'.format(snow_density_str)
+                        snow_density_eurocode = [150 for _ in study.ordered_years]
+                        ax.plot(study.ordered_years, snow_density_eurocode, color='k', label=label)
                     ax.legend()
                     tight_pad = {'h_pad': 0.2}
+                    ax.set_ylim(ylim)
+                    ax.set_xlim([1957, 2018])
+                    ax.yaxis.set_ticks(yticks)
                     study_visualizer.show_or_save_to_file(no_title=True, tight_layout=True, tight_pad=tight_pad)
                     ax.clear()
 

experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -63,12 +63,11 @@ def major_result():
 
 
 if __name__ == '__main__':
-    major_result()
-    # intermediate_result(altitudes=paper_altitudes, massif_names=['Chartreuse', 'Vanoise', 'Mercantour'],
-    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
-    #                        ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
-    #                     non_stationary_uncertainty=[False, True][1:],
-    #                     study_class=CrocusSnowLoadEurocode)
+    # major_result()
+    intermediate_result(altitudes=[600, 900], massif_names=['Maurienne', 'Oisans'],
+                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+                           ConfidenceIntervalMethodFromExtremes.ci_mle][:],
+                        non_stationary_uncertainty=[False, True][:])
     # intermediate_result(altitudes=[900, 1200], massif_names=None)
     # intermediate_result(ALL_ALTITUDES_WITHOUT_NAN)
     # intermediate_result(paper_altitudes)

experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py

@@ -42,29 +42,27 @@ def plot_subgroup_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVis
     visualizer = list(altitude_to_visualizer.values())[0]
     nb_massif_names = len(massif_names)
     assert nb_massif_names <= 5
-    axes = create_adjusted_axes(nb_massif_names, visualizer.nb_contexts)
-    if nb_massif_names == 1:
-        axes = [axes]
-    for ax, massif_name in zip(axes, massif_names):
-        plot_single_uncertainty_massif(altitude_to_visualizer,
-                                       massif_name, ax)
-
-    # Save plot
-    massif_names_str = '_'.join(massif_names)
-    model_names_str = 'NonStationarity=' + '_'.join([str(e) for e in visualizer.non_stationary_contexts])
-    visualizer.plot_name = model_names_str + '_' + massif_names_str
-    visualizer.show_or_save_to_file(no_title=True)
-    plt.close()
+    # axes = create_adjusted_axes(nb_massif_names, visualizer.nb_contexts)
+    # if nb_massif_names == 1:
+    #     axes = [axes]
+    for massif_name in massif_names:
+        plot_single_uncertainty_massif(altitude_to_visualizer, massif_name)
 
 
 def plot_single_uncertainty_massif(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
-                                   massif_name, axes):
+                                   massif_name):
     visualizer = list(altitude_to_visualizer.values())[0]
-    if visualizer.nb_contexts == 1:
-        axes = [axes]
-    for ax, non_stationary_context in zip(axes, visualizer.non_stationary_contexts):
+
+    for non_stationary_context in visualizer.non_stationary_contexts:
+        ax = create_adjusted_axes(1, 1)
         plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, non_stationary_context,
                                                                   altitude_to_visualizer)
+        # Save plot
+        massif_names_str = massif_name
+        model_names_str = 'NonStationarity={}'.format(non_stationary_context)
+        visualizer.plot_name = model_names_str + '_' + massif_names_str
+        visualizer.show_or_save_to_file(no_title=True)
+        plt.close()
 
 
 def get_label_name(non_stationary_context, ci_method_name):
@@ -114,7 +112,7 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
         non_stationary_context = 'selected non-stationary models'
     else:
         non_stationary_context = 'the stationary model'
-    title = '{} massif with {}'.format(massif_name_str,  non_stationary_context)
+    title = '{} massif with {}'.format(massif_name_str, non_stationary_context)
     ax.set_title(title)
     ax.set_xticks(altitudes)
     ylabel = EUROCODE_RETURN_LEVEL_STR.replace('GSL', SCM_STUDY_CLASS_TO_ABBREVIATION[type(visualizer.study)])

experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py

@@ -46,6 +46,7 @@ def plot_histogram(altitude_to_visualizer, non_stationary_context):
     visualizer.plot_name = 'Percentages of exceedance with non_stationary={}'.format(non_stationary_context)
     # visualizer.show = True
     visualizer.show_or_save_to_file(no_title=True)
+    ax.clear()
 
 
 def plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width):

experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -242,7 +242,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     @cached_property
     def massif_name_to_eurocode_values(self):
         """Eurocode values for the altitude"""
-        return {m: r().lois_de_variation_de_la_valeur_caracteristique(altitude=self.study.altitude)
+        return {m: r().valeur_caracteristique(altitude=self.study.altitude)
                 for m, r in massif_name_to_eurocode_region.items() if m in self.uncertainty_massif_names}
 
     def three_percentages_of_excess(self, ci_method, non_stationary_context):

extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py

@@ -17,18 +17,13 @@ ci_method_to_method_name = {
     ConfidenceIntervalMethodFromExtremes.ci_proflik: 'proflik',
 }
 
-# ci_method_to_color = {
-#     ConfidenceIntervalMethodFromExtremes.ci_mle: 'tab:brown',
-#     ConfidenceIntervalMethodFromExtremes.my_bayes: 'tab:green'
-# }
-
 ci_method_to_color = {
-    ConfidenceIntervalMethodFromExtremes.ci_mle: 'yellowgreen',
+    ConfidenceIntervalMethodFromExtremes.ci_mle: 'mediumseagreen',
     ConfidenceIntervalMethodFromExtremes.my_bayes: 'darkgreen'
 }
 
-common_part_and_uncertainty = 'Return level and its uncertainty with the'
+# common_part_and_uncertainty = 'Return level and its uncertainty with the'
 ci_method_to_label = {
-    ConfidenceIntervalMethodFromExtremes.ci_mle: 'Bayesian procedure',
-    ConfidenceIntervalMethodFromExtremes.my_bayes: 'Delta method '
+    ConfidenceIntervalMethodFromExtremes.ci_mle: 'Delta method',
+    ConfidenceIntervalMethodFromExtremes.my_bayes: 'Bayesian procedure'
 }
\ No newline at end of file

extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py

@@ -53,7 +53,7 @@ class EurocodeConfidenceIntervalFromExtremes(object):
             fit_method = TemporalMarginFitMethod.extremes_fevd_mle
         # Fitted estimator
         fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=None,
-                                                          fit_method=fit_method)
+                                                          fit_method=fit_method, nb_iterations_for_bayesian_fit=20000)
         # Load object from result from extremes
         return cls.from_estimator_extremes(fitted_estimator, ci_method, temporal_covariate, quantile_level)