From a97330917c2951fe5cae5fefc9beb684c61c3bfe Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 16 Mar 2020 08:19:52 +0100
Subject: [PATCH] [paper 1] add plot for the presentation. increase the size of
the legend for uncertainty curves. add sample function for the GevParams
class.
---
.../scm_models_data/abstract_study.py | 1 +
extreme_fit/distribution/gev/gev_params.py | 6 ++
.../presentation/__init__.py | 0
.../accumulation_in_winter.py | 0
.../presentation/statistical_model.py | 81 +++++++++++++++++++
.../main_result_trends_and_return_levels.py | 8 +-
.../plot_uncertainty_curves.py | 18 +++--
7 files changed, 104 insertions(+), 10 deletions(-)
create mode 100644 papers/exceeding_snow_loads/presentation/__init__.py
rename papers/exceeding_snow_loads/{ => presentation}/accumulation_in_winter.py (100%)
create mode 100644 papers/exceeding_snow_loads/presentation/statistical_model.py
diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py
index 66ca5708..0e548d1a 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/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):
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index f27beb70..db84b878 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/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
diff --git a/papers/exceeding_snow_loads/presentation/__init__.py b/papers/exceeding_snow_loads/presentation/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/papers/exceeding_snow_loads/accumulation_in_winter.py b/papers/exceeding_snow_loads/presentation/accumulation_in_winter.py
similarity index 100%
rename from papers/exceeding_snow_loads/accumulation_in_winter.py
rename to papers/exceeding_snow_loads/presentation/accumulation_in_winter.py
diff --git a/papers/exceeding_snow_loads/presentation/statistical_model.py b/papers/exceeding_snow_loads/presentation/statistical_model.py
new file mode 100644
index 00000000..9e15740f
--- /dev/null
+++ b/papers/exceeding_snow_loads/presentation/statistical_model.py
@@ -0,0 +1,81 @@
+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
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index c656915b..8769a25b 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/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,
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
index c9e04aae..87b8325e 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
+++ b/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([])
--
GitLab