diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index 407fae9ec5daa6cefd79e6518c74614a14325385..258c03251ff26123fa2936339cd648644ac88c1d 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -8,6 +8,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
ALL_ALTITUDES_WITHOUT_NAN
from experiment.paper_past_snow_loads.paper_main_utils import load_altitude_to_visualizer
from experiment.paper_past_snow_loads.paper_utils import paper_study_classes, paper_altitudes, ModelSubsetForUncertainty
+from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_diagnosis_risk import plot_diagnosis_risk
from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_trend_curves import plot_trend_curves, \
plot_trend_map
from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import \
@@ -68,10 +69,11 @@ def intermediate_result(altitudes, massif_names=None,
_ = compute_minimized_aic(visualizer)
# Plots
- plot_trend_map(altitude_to_visualizer)
+ # plot_trend_map(altitude_to_visualizer)
+ plot_diagnosis_risk(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)
def major_result():
@@ -92,7 +94,7 @@ def major_result():
if __name__ == '__main__':
# major_result()
- intermediate_result(altitudes=ALL_ALTITUDES_WITHOUT_NAN[:2], massif_names=None,
+ intermediate_result(altitudes=[1500, 1800], massif_names=None,
uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
multiprocessing=True)
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_diagnosis_risk.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_diagnosis_risk.py
new file mode 100644
index 0000000000000000000000000000000000000000..a6112812dd2349afe9e567e1ab22ff92805b3f2e
--- /dev/null
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_diagnosis_risk.py
@@ -0,0 +1,36 @@
+import matplotlib.pyplot as plt
+
+from experiment.paper_past_snow_loads.paper_utils import ModelSubsetForUncertainty, dpi_paper1_figure
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+ ConfidenceIntervalMethodFromExtremes
+
+
+def plot_diagnosis_risk(altitude_to_visualizer):
+ ax = plt.gca()
+ altitudes = list(altitude_to_visualizer.keys())
+ visualizers = list(altitude_to_visualizer.values())
+ visualizer = visualizers[0]
+ model_subset_for_uncertainty = ModelSubsetForUncertainty.non_stationary_gumbel_and_gev
+ ci_method = ConfidenceIntervalMethodFromExtremes.ci_mle
+
+ plot_mean_exceedance(visualizers, altitudes, ax, model_subset_for_uncertainty, ci_method)
+
+ visualizer.plot_name = 'Diagnosis plot'
+ visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
+ plt.close()
+
+ plt.show()
+
+
+def plot_mean_exceedance(visualizers, altitudes, ax, model_subset_for_uncertainty, ci_method):
+ l = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[3:] for v in visualizers]
+ diff, diff_c, diff_e = zip(*l)
+ ax.plot(altitudes, diff, marker='o', color='red', label='diff')
+ ax.plot(altitudes, diff_c, marker='o', color='yellow', label='diff c')
+ ax.plot(altitudes, diff_e, marker='o', color='orange', label='diff e')
+
+ ax.set_ylabel("Mean exceedance (kN)", fontsize=15)
+ ax.set_xlabel("Altitude", fontsize=15)
+ ax.yaxis.grid()
+
+ ax.legend()
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
index dd5ff29783a363a5767e66cfcc56859ffecd3a6f..680c7e888ec85e8b6b14c077be04b1c7b4f9ba8d 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
import numpy as np
from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_ALTITUDES
-from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure
+from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure, ModelSubsetForUncertainty
from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color, \
@@ -33,6 +33,10 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
ax = plt.gca()
altitudes = np.array(list(altitude_to_visualizer.keys()))
bincenters = altitudes
+
+ fontsize_label = 15
+ legend_size = 15
+ # Plot histogram
for j, ci_method in enumerate(visualizer.uncertainty_methods):
if len(visualizer.uncertainty_methods) == 2:
offset = -50 if j == 0 else 50
@@ -41,13 +45,15 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
else:
width = 200
plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width=width)
- fontsize_label = 15
- legend_size = 15
+
+
ax.set_xticks(altitudes)
ax.tick_params(labelsize=fontsize_label)
if not (len(visualizer.uncertainty_methods) == 1
and visualizer.uncertainty_methods[0] == ConfidenceIntervalMethodFromExtremes.ci_mle):
ax.legend(loc='upper left', prop={'size': legend_size})
+ # ax.set_ylabel('Massifs whose 50-year return level\n'
+ # 'exceeds French standards (\%)', fontsize=fontsize_label)
ax.set_ylabel('Massifs exceeding French standards (\%)', fontsize=fontsize_label)
ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
ax.set_ylim([0, 100])
@@ -59,7 +65,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
def plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width):
- three_percentages_of_excess = [v.three_percentages_of_excess(ci_method, model_subset_for_uncertainty) for v in
+ 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]
diff --git a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
index 0116fe33f7ea201cae951c24a2ac70aa908d00d0..ab2a9c0ee0259a6213804c04e6b49ef5f862af45 100644
--- a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -6,6 +6,7 @@ import matplotlib.pyplot as plt
import numpy as np
from cached_property import cached_property
+from experiment.eurocode_data.eurocode_region import C2, C1, E
from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
from experiment.eurocode_data.utils import EUROCODE_QUANTILE, EUROCODE_RETURN_LEVEL_STR
from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map, get_colors
@@ -359,16 +360,31 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
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, model_subset_for_uncertainty):
- eurocode_and_uncertainties = [(self.massif_name_to_eurocode_values[massif_name],
- self.triplet_to_eurocode_uncertainty[
- (ci_method, model_subset_for_uncertainty, massif_name)])
- for massif_name in self.uncertainty_massif_names]
- a = np.array([(uncertainty.confidence_interval[0] > eurocode,
- uncertainty.mean_estimate > eurocode,
- uncertainty.confidence_interval[1] > eurocode)
- for eurocode, uncertainty in eurocode_and_uncertainties])
- return 100 * np.mean(a, axis=0)
+ def excess_metrics(self, ci_method, model_subset_for_uncertainty):
+ triplet = [(massif_name_to_eurocode_region[massif_name],
+ self.massif_name_to_eurocode_values[massif_name],
+ self.triplet_to_eurocode_uncertainty[(ci_method, model_subset_for_uncertainty, massif_name)])
+ for massif_name in self.uncertainty_massif_names]
+ # First array for histogram
+ a = 100 * np.array([(uncertainty.confidence_interval[0] > eurocode,
+ uncertainty.mean_estimate > eurocode,
+ uncertainty.confidence_interval[1] > eurocode)
+ for _, eurocode, uncertainty in triplet])
+ a = np.mean(a, axis=0)
+ # Second array for curve
+ b = [[] for _ in range(3)]
+ for eurocode_region, eurocode, uncertainty in triplet:
+ diff = uncertainty.mean_estimate - eurocode
+ b[0].append(diff)
+ if eurocode_region in [C1, C2]:
+ b[1].append(diff)
+ if eurocode_region in [E]:
+ b[2].append(diff)
+
+ b = [np.mean(np.array(e)) for e in b]
+ # Return the concatenated results
+ concatenated_result = list(a) + list(b)
+ return concatenated_result
# Part 3 - QQPLOT