[paper 1] add plot for the presentation. increase the size of the legend for...

[paper 1] add plot for the presentation. increase the size of the legend for uncertainty curves. add sample function for the GevParams class.

experiment/meteo_france_data/scm_models_data/abstract_study.py

@@ -40,6 +40,7 @@ with redirect_stdout(f):
     from simpledbf import Dbf5
 
 filled_marker_legend_list = ['Filled marker =', 'Selected model is significant', 'w.r.t $\mathcal{M}_0$']
+filled_marker_legend_list2 = ['Filled marker = Selected', 'model is significant', 'w.r.t $\mathcal{M}_0$']
 
 
 class AbstractStudy(object):

extreme_fit/distribution/gev/gev_params.py

@@ -24,6 +24,12 @@ class GevParams(AbstractParams):
         if accept_zero_scale_parameter and scale == 0.0:
             self.has_undefined_parameters = False
 
+    def sample(self, n) -> float:
+        if self.has_undefined_parameters:
+            return np.nan
+        else:
+            return r.rgev(n, self.location, self.scale, self.shape)
+
     def quantile(self, p) -> float:
         if self.has_undefined_parameters:
             return np.nan

papers/exceeding_snow_loads/presentation/statistical_model.py

+import matplotlib.pyplot as plt
+import numpy as np
+
+from extreme_fit.distribution.gev.gevmle_fit import GevMleFit
+from extreme_fit.distribution.gev.ismev_gev_fit import ismev_gev_fit
+
+
+def binomial_observation():
+    marker_altitude_massif_name_for_paper1 = [
+        ('magenta', 'stationary', [0.5, 0.5]),
+        # ('darkmagenta', 'increase', [0.5, 0.8]),
+        # ('mediumpurple', 'decrease', [0.5, 0.3]),
+    ]
+    ax = plt.gca()
+    for color, label, l in marker_altitude_massif_name_for_paper1:
+        before, after = l
+        total_time = 60
+        data_before = np.random.binomial(1, before, int(total_time / 3))
+        data_after = np.random.binomial(1, after, int(2 * total_time / 3))
+        data = np.concatenate([data_before, data_after], axis=0)
+        time = list(range(total_time))
+        ax.tick_params(axis='both', which='major', labelsize=15)
+        fontsize = 20
+        ax.set_xlabel('Time', fontsize=fontsize)
+        ax.set_ylabel('Coin value', fontsize=fontsize)
+        plt.yticks([0.0, 1.0], ['Heads', 'Tails'])
+        ax.plot(time, data, color=color, label=label, linewidth=4)
+
+def histogram_for_gev():
+    import matplotlib.pyplot as plt
+    from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+    from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import AbstractSnowLoadVariable
+    ax = plt.gca()
+    study = CrocusSnowLoadTotal(altitude=1800)
+    s = study.observations_annual_maxima.df_maxima_gev.loc['Vercors']
+    x_gev = s.values
+    gev = GevMleFit(x_gev)
+    gev_params = gev.gev_params
+    samples = gev_params.sample(10000)
+    nb = 10
+    epsilon = 0.0
+    x, bins, p = ax.hist(samples, bins=nb, color='white', edgecolor='grey', density=True, stacked=True,
+                         linewidth=3, bottom=[-epsilon for _ in range(nb)])
+    for item in p:
+        item.set_height((item.get_height() / sum(x)))
+    # print(gev_params)
+    # x = np.linspace(0.0, 10, 1000)
+    # y = gev_params.density(x)
+    # ax.plot(x, y, linewidth=5)
+    ax.set_xlabel('Annual maxima of GSL ({})'.format(AbstractSnowLoadVariable.UNIT), fontsize=15)
+    ax.set_ylabel('Probability', fontsize=15)
+    ax.tick_params(axis='both', which='major', labelsize=15)
+    ax.set_yticks([0, 0.1, 0.2, 0.3])
+    ax.set_xlim([0, 10])
+    ax.set_ylim([0, 0.3])
+
+
+
+
+
+def histogram_for_normal():
+    ax = plt.gca()
+    linewidth = 5
+    absisse = [0.25, 0.75]
+    ax.bar(absisse, [0.6, 0.6], width=0.2, color='white', edgecolor='grey', linewidth=linewidth, bottom=-0.1)
+    plt.xticks(absisse, ['Heads', 'Tail'])
+    ax.set_yticks([0, 0.5])
+    ax.set_xlabel('Coin value', fontsize=15)
+    ax.tick_params(axis='both', which='major', labelsize=15)
+
+    ax.set_ylabel('Probability', fontsize=15)
+    ax.set_xlim([0, 1])
+    ax.set_ylim([0, 0.6])
+
+
+
+if __name__ == '__main__':
+    # binomial_observation()
+    # histogram_for_gev()
+    histogram_for_normal()
+    plt.show()
\ No newline at end of file

papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -71,8 +71,8 @@ def intermediate_result(altitudes, massif_names=None,
     # Plots
     # 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_massifs(altitude_to_visualizer)
+    # plot_uncertainty_histogram(altitude_to_visualizer)
     # plot_selection_curves(altitude_to_visualizer)
     # uncertainty_interval_size(altitude_to_visualizer)
 
@@ -80,8 +80,8 @@ def intermediate_result(altitudes, massif_names=None,
 def major_result():
     uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
                            ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
-    massif_names = None
-    study_classes = paper_study_classes[:2]
+    massif_names = ['Beaufortain', 'Vercors']
+    study_classes = paper_study_classes[:1]
     # model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
     #                                  ModelSubsetForUncertainty.stationary_gumbel_and_gev,
     #                                  ModelSubsetForUncertainty.non_stationary_gumbel,

papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py

@@ -5,7 +5,8 @@ import numpy as np
 
 from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_ALTITUDES, \
     YEAR_OF_INTEREST_FOR_RETURN_LEVEL
-from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy, filled_marker_legend_list
+from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy, filled_marker_legend_list, \
+    filled_marker_legend_list2
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
     SCM_STUDY_CLASS_TO_ABBREVIATION
 from papers.exceeding_snow_loads.paper_utils import dpi_paper1_figure, ModelSubsetForUncertainty
@@ -90,7 +91,7 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, m
     altitudes = list(altitude_to_visualizer.keys())
     visualizer = list(altitude_to_visualizer.values())[0]
     alpha = 0.2
-    legend_size = 25
+    legend_size = 30
     fontsize_label = 35
     # Display the EUROCODE return level
     eurocode_region = massif_name_to_eurocode_region[massif_name]()
@@ -130,7 +131,8 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, m
     # add_title(ax, eurocode_region, massif_name, non_stationary_context)
     ax.set_xticks(altitudes)
     ax.tick_params(labelsize=fontsize_label)
-    ylabel = EUROCODE_RETURN_LEVEL_STR.replace('GSL', SCM_STUDY_CLASS_TO_ABBREVIATION[type(visualizer.study)])
+    ylabel = EUROCODE_RETURN_LEVEL_STR
+    # ylabel = EUROCODE_RETURN_LEVEL_STR.replace('GSL', SCM_STUDY_CLASS_TO_ABBREVIATION[type(visualizer.study)])
     ax.set_ylabel(ylabel, fontsize=fontsize_label)
     ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
     ax.grid()
@@ -179,12 +181,16 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, leg
         ax2 = ax.twinx()
         # ax2.legend(handles=legend_elements, bbox_to_anchor=(0.93, 0.7), loc='upper right')
         # ax2.annotate("Filled symbol = significant trend ", xy=(0.85, 0.5), xycoords='axes fraction', fontsize=7)
-        ax2.legend(handles=legend_elements, loc='center left', prop={'size': legend_size})
+        upper_legend_y = 0.55
+        ax2.annotate('\n'.join(filled_marker_legend_list2), xy=(0.23, upper_legend_y - 0.2), xycoords='axes fraction', fontsize=fontsize)
+        ax2.annotate('Markers show selected model $\mathcal{M}_N$', xy=(0.02, upper_legend_y), xycoords='axes fraction', fontsize=fontsize)
+        print(legend_size)
+        ax2.legend(handles=legend_elements, prop={'size': legend_size},
+                   loc='upper left', bbox_to_anchor=(0.00, upper_legend_y))
         # for handle in lgnd.legendHandles:
         #     handle.set_sizes([6.0])
         # ax2.annotate("Filled symbol =\nsignificant trend  \nw.r.t $\mathcal{M}_0$", xy=(0.6, 0.85), xycoords='axes fraction', fontsize=fontsize)
-        ax2.annotate('\n'.join(filled_marker_legend_list), xy=(0.23, 0.43), xycoords='axes fraction', fontsize=fontsize)
-        ax2.annotate('Markers show selected model $\mathcal{M}_N$', xy=(0.02, 0.605), xycoords='axes fraction', fontsize=fontsize)
+
         ax2.set_yticks([])