From 11cc9e8a815e0c7127c84b55ccaa7384e3d2be00 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 13 Nov 2019 16:17:37 +0100
Subject: [PATCH] [Thesis committee drawing] improve GEV non stationary plots.
rename ci_delta to ci_mle. catch some ci runtime exceptions. remove
hatch/stripes drawing temporarily. fix unit for snow load. add new crocus
classes & variable classes for snow load. modify conversion factor only for
the new plots.
---
.../eurocode_data/eurocode_visualizer.py | 40 ++++---
.../eurocode_data/main_eurocode_drawing.py | 111 +++++++++++++-----
experiment/eurocode_data/slide_plot.py | 18 +++
experiment/eurocode_data/utils.py | 4 +-
.../scm_models_data/abstract_study.py | 25 ++--
.../scm_models_data/crocus/crocus.py | 18 ++-
.../crocus/crocus_variables.py | 16 +++
.../main_study_visualizer.py | 15 ++-
.../study_visualization/study_visualizer.py | 19 ++-
.../poster_EVAN2019/main_poster_EVAN2019.py | 6 +-
.../confidence_interval_method.py | 4 +-
.../result_from_mle_extremes.py | 6 +-
.../test_model/test_confidence_interval.py | 2 +-
thesis_report/gev_plot.py | 23 ++--
14 files changed, 226 insertions(+), 81 deletions(-)
create mode 100644 experiment/eurocode_data/slide_plot.py
diff --git a/experiment/eurocode_data/eurocode_visualizer.py b/experiment/eurocode_data/eurocode_visualizer.py
index 847b09af..e4f0ff21 100644
--- a/experiment/eurocode_data/eurocode_visualizer.py
+++ b/experiment/eurocode_data/eurocode_visualizer.py
@@ -1,6 +1,7 @@
from typing import Dict, List, Tuple
import matplotlib.pyplot as plt
+import numpy as np
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
EurocodeConfidenceIntervalFromExtremes
@@ -17,7 +18,7 @@ def get_label_name(model_name, ci_method_name: str):
model_name += model_symbol
model_name += '}}'
model_name += parameter
- model_name += ')_{ \\textrm{' + ci_method_name.upper() + '}} $ '
+ model_name += ')_{ \\textrm{' + ci_method_name.upper().split(' ')[1] + '}} $ '
return model_name
@@ -25,7 +26,7 @@ def get_model_name(model_class):
return get_display_name_from_object_type(model_class).split('Stationary')[0] + 'Stationary'
-def plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(d, nb_massif_names, nb_model_names, show=True):
+def plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(d, nb_massif_names, nb_model_names):
"""
Rows correspond to massif names
Columns correspond to stationary/non stationary model name for a given date
@@ -39,11 +40,7 @@ def plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(d, nb_m
plot_model_name_to_uncertainty_method_to_ordered_dict(model_name_to_uncertainty_level,
massif_name, ax)
- plt.suptitle('50-year return levels of extreme snow loads in France for several confiance interval methods.')
-
- if show:
- plt.show()
-
+ # plt.suptitle('50-year return levels of extreme snow loads in France for several confiance interval methods.')
def plot_model_name_to_uncertainty_method_to_ordered_dict(d, massif_name, axes):
if len(d) == 1:
@@ -71,17 +68,30 @@ def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name
# Plot eurocode standards only for the first loop
if j == 0:
eurocode_region.plot_max_loading(ax, altitudes=altitudes)
- conversion_factor = 100 # We divide by 100 to transform the snow water equivalent into snow load
- mean = [r.mean_estimate / conversion_factor for r in ordered_return_level_uncertaines]
+ conversion_factor = 9.8 / 1000
+ mean = [r.mean_estimate * conversion_factor for r in ordered_return_level_uncertaines]
+ # Filter and keep only non nan values
+ not_nan_index = [not np.isnan(m) for m in mean]
+ mean = list(np.array(mean)[not_nan_index])
+ altitudes = list(np.array(altitudes)[not_nan_index])
+ ordered_return_level_uncertaines = list(np.array(ordered_return_level_uncertaines)[not_nan_index])
+
ci_method_name = str(label).split('.')[1].replace('_', ' ')
- ax.plot(altitudes, mean, '-', color=color, label=get_label_name(model_name, ci_method_name))
- lower_bound = [r.confidence_interval[0] / conversion_factor for r in ordered_return_level_uncertaines]
- upper_bound = [r.confidence_interval[1] / conversion_factor for r in ordered_return_level_uncertaines]
+ ax.plot(altitudes, mean, linestyle='--', marker='o', color=color, label=get_label_name(model_name, ci_method_name))
+ lower_bound = [r.confidence_interval[0] * conversion_factor for r in ordered_return_level_uncertaines]
+ upper_bound = [r.confidence_interval[1] * conversion_factor for r in ordered_return_level_uncertaines]
ax.fill_between(altitudes, lower_bound, upper_bound, color=color, alpha=alpha)
ax.legend(loc=2)
- ax.set_ylim([0.0, 4.0])
- ax.set_title(massif_name + ' ' + model_name)
- ax.set_ylabel('50-year return level (N $m^-2$)')
+ ax.set_ylim([0.0, 16])
+ massif_name_str = massif_name.replace('_', ' ')
+ eurocode_region_str = get_display_name_from_object_type(type(eurocode_region))
+ if 'Non' in model_name:
+ model_name = 'non-stationary'
+ else:
+ model_name = 'stationary'
+ title = '{} ({} Eurocodes area) with a {} model'.format(massif_name_str, eurocode_region_str, model_name)
+ ax.set_title(title)
+ ax.set_ylabel('50-year return level (kN $m^-2$)')
ax.set_xlabel('Altitude (m)')
ax.grid()
diff --git a/experiment/eurocode_data/main_eurocode_drawing.py b/experiment/eurocode_data/main_eurocode_drawing.py
index 8cdaeb16..f6c15f6f 100644
--- a/experiment/eurocode_data/main_eurocode_drawing.py
+++ b/experiment/eurocode_data/main_eurocode_drawing.py
@@ -1,13 +1,18 @@
import time
+import os.path as op
+import matplotlib.pyplot as plt
from collections import OrderedDict
+from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
+ StudyVisualizer
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
ConfidenceIntervalMethodFromExtremes
from experiment.eurocode_data.eurocode_visualizer import \
plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict, get_model_name
from experiment.eurocode_data.massif_name_to_departement import MASSIF_NAMES_ALPS
from experiment.eurocode_data.utils import EUROCODE_ALTITUDES, LAST_YEAR_FOR_EUROCODE
-from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days, CrocusSweTotal
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_hypercube_visualizer import \
AltitudeHypercubeVisualizer
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.utils_hypercube import \
@@ -18,6 +23,8 @@ from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_m
# Model class
import matplotlib as mpl
+from root_utils import VERSION_TIME
+
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
@@ -32,7 +39,7 @@ def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for
only_first_one=False, save_to_file=False,
exact_starting_year=1958,
first_starting_year=None,
- study_classes=[CrocusSwe3Days],
+ study_classes=[CrocusSwe3Days, CrocusSweTotal][1:],
trend_test_class=None) # type: AltitudeHypercubeVisualizer
# Loop on the data
assert isinstance(altitude_visualizer.tuple_to_study_visualizer, OrderedDict)
@@ -51,29 +58,7 @@ def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for
return {model_name: massif_name_to_ordered_eurocode_level_uncertainty}
-def main_drawing():
- fast_plot = [True, False][0]
- temporal_covariate = 2017
- # Select parameters
- massif_names = MASSIF_NAMES_ALPS[:]
- model_class_and_last_year = [
- (StationaryTemporalModel, 2017),
- (NonStationaryLocationAndScaleTemporalModel, 2017),
- # Add the temperature here
- ][1:]
- altitudes = EUROCODE_ALTITUDES[:]
- uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
- ConfidenceIntervalMethodFromExtremes.ci_normal]
- show = True
-
- if fast_plot:
- show = True
- model_class_and_last_year = model_class_and_last_year[:1]
- altitudes = altitudes[2:]
- # altitudes = altitudes[:]
- massif_names = massif_names[:1]
- uncertainty_methods = uncertainty_methods[1:]
-
+def plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, temporal_covariate, uncertainty_methods):
model_name_to_massif_name_to_ordered_return_level = {}
for model_class, last_year_for_the_data in model_class_and_last_year:
start = time.time()
@@ -91,8 +76,80 @@ def main_drawing():
# Plot graph
plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(
massif_name_to_model_name_to_ordered_return_level_uncertainties, nb_massif_names=len(massif_names),
- nb_model_names=len(model_class_and_last_year), show=show)
+ nb_model_names=len(model_class_and_last_year))
+ if show:
+ plt.show()
+ else:
+ massif_names_str = '_'.join(massif_names)
+ model_names_str = '_'.join(
+ [model_name for model_name in model_name_to_massif_name_to_ordered_return_level.keys()])
+ filename = op.join(VERSION_TIME, model_names_str + '_' + massif_names_str)
+ StudyVisualizer.savefig_in_results(filename)
+
+
+def main_drawing():
+ fast_plot = [True, False][0]
+ temporal_covariate = 2017
+ # Select parameters
+ massif_names = MASSIF_NAMES_ALPS[:]
+ model_class_and_last_year = [
+ (StationaryTemporalModel, 2017),
+ (NonStationaryLocationAndScaleTemporalModel, 2017),
+ # Add the temperature here
+ ][:]
+ altitudes = EUROCODE_ALTITUDES[:]
+ uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+ ConfidenceIntervalMethodFromExtremes.ci_mle]
+ show = False
+
+ if fast_plot:
+ show = True
+ model_class_and_last_year = model_class_and_last_year[:]
+ altitudes = altitudes[-2:]
+ # altitudes = altitudes[:]
+ massif_names = massif_names[:1]
+ uncertainty_methods = uncertainty_methods[:]
+
+ plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, temporal_covariate, uncertainty_methods)
+
+
+# Create 5 main plots
+def main_5_drawings():
+ model_class_and_last_year = [
+ (StationaryTemporalModel, 2017),
+ (NonStationaryLocationAndScaleTemporalModel, 2017),
+ # Add the temperature here
+ ][:1]
+ altitudes = EUROCODE_ALTITUDES[:]
+ uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+ ConfidenceIntervalMethodFromExtremes.ci_mle]
+ show = False
+ massif_names = MASSIF_NAMES_ALPS[:]
+ temporal_covariate = 2017
+ m = 4
+ n = (23 // m) + 1
+ for i in list(range(n))[:]:
+ massif_names = MASSIF_NAMES_ALPS[m * i: m * (i+1)]
+ print(massif_names)
+ plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, temporal_covariate, uncertainty_methods)
+
+
+def main_3_massif_of_interest():
+ massif_names = ['Parpaillon', 'Chartreuse', 'Maurienne'][1:]
+ model_class_and_last_year = [
+ (StationaryTemporalModel, 2017),
+ (NonStationaryLocationAndScaleTemporalModel, 2017),
+ # Add the temperature here
+ ][:]
+ altitudes = EUROCODE_ALTITUDES[1:]
+ uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+ ConfidenceIntervalMethodFromExtremes.ci_mle][:]
+ temporal_covariate = 2017
+ show = False
+ plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, temporal_covariate, uncertainty_methods)
if __name__ == '__main__':
- main_drawing()
+ # main_drawing()
+ # main_5_drawings()
+ main_3_massif_of_interest()
\ No newline at end of file
diff --git a/experiment/eurocode_data/slide_plot.py b/experiment/eurocode_data/slide_plot.py
new file mode 100644
index 00000000..e46e0f2f
--- /dev/null
+++ b/experiment/eurocode_data/slide_plot.py
@@ -0,0 +1,18 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+from experiment.eurocode_data.eurocode_region import C2, E
+from root_utils import get_display_name_from_object_type
+
+if __name__ == '__main__':
+ ax = plt.gca()
+ altitudes = np.linspace(200, 2000)
+ for region_class in [C2, E][1:]:
+ region_object = region_class()
+ region_object.plot_max_loading(ax, altitudes)
+ # ax.set_title(get_display_name_from_object_type(region_object) + ' Eurocodes region')
+ ax.set_ylabel('50-year return level (kN $m^-2$)')
+ ax.set_xlabel('Altitude (m)')
+ ax.set_ylim([0.0, 11.0])
+ ax.grid()
+ plt.show()
\ No newline at end of file
diff --git a/experiment/eurocode_data/utils.py b/experiment/eurocode_data/utils.py
index 635b3531..764ddbd0 100644
--- a/experiment/eurocode_data/utils.py
+++ b/experiment/eurocode_data/utils.py
@@ -1,8 +1,8 @@
# Eurocode quantile correspond to a 50 year return period
EUROCODE_QUANTILE = 0.98
-# Altitudes (between low and mid altitudes) < 2000m
-EUROCODE_ALTITUDES = [0, 300, 600, 900, 1200, 1500, 1800]
+# Altitudes (between low and mid altitudes) < 2000m and should be > 200m
+EUROCODE_ALTITUDES = [300, 600, 900, 1200, 1500, 1800]
# Last year taken into account for the Eurocode
# Date of publication was 2014, therefore the winter 2013/2014 could not have been measured
# Therefore, the winter 2012/2013 was the last one. Thus, 2012 is the last year for the Eurocode
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 8c407129..07f7a2be 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -335,15 +335,16 @@ class AbstractStudy(object):
else:
ax.fill(*coords_list, **{'color': default_color_for_missing_massif})
- # Add a hatch to visualize the mean & variance variation sign
- hatch_list = ['//', '\\\\']
- if massif_name_to_hatch_boolean_list is not None:
- if massif_name in massif_name_to_hatch_boolean_list:
- a = np.array(coords_list).transpose()
- hatch_boolean_list = massif_name_to_hatch_boolean_list[massif_name]
- for hatch, is_hatch in zip(hatch_list, hatch_boolean_list):
- if is_hatch:
- ax.add_patch(Polygon(xy=a, fill=False, hatch=hatch))
+ # For the moment we comment all the part of this code
+ # # Add a hatch to visualize the mean & variance variation sign
+ # hatch_list = ['//', '\\\\']
+ # if massif_name_to_hatch_boolean_list is not None:
+ # if massif_name in massif_name_to_hatch_boolean_list:
+ # a = np.array(coords_list).transpose()
+ # hatch_boolean_list = massif_name_to_hatch_boolean_list[massif_name]
+ # for hatch, is_hatch in zip(hatch_list, hatch_boolean_list):
+ # if is_hatch:
+ # ax.add_patch(Polygon(xy=a, fill=False, hatch=hatch))
# Display the center of the massif
masssif_coordinate_for_display = cls.massifs_coordinates_for_display(massif_names)
@@ -404,9 +405,9 @@ class AbstractStudy(object):
def map_full_path(self) -> str:
return op.join(self.full_path, 'map')
- @property
- def result_full_path(self) -> str:
- return op.join(self.full_path, 'results')
+ @classproperty
+ def result_full_path(cls) -> str:
+ return op.join(cls.full_path, 'results')
@property
def study_full_path(self) -> str:
diff --git a/experiment/meteo_france_data/scm_models_data/crocus/crocus.py b/experiment/meteo_france_data/scm_models_data/crocus/crocus.py
index 771e5078..0264b02f 100644
--- a/experiment/meteo_france_data/scm_models_data/crocus/crocus.py
+++ b/experiment/meteo_france_data/scm_models_data/crocus/crocus.py
@@ -3,7 +3,7 @@ import numpy as np
from experiment.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusTotalSweVariable, \
- CrocusDepthVariable, CrocusRecentSweVariable
+ CrocusDepthVariable, CrocusRecentSweVariable, TotalSnowLoadVariable, RecentSnowLoadVariable
class Crocus(AbstractStudy):
@@ -12,7 +12,8 @@ class Crocus(AbstractStudy):
"""
def __init__(self, variable_class, *args, **kwargs):
- assert variable_class in [CrocusTotalSweVariable, CrocusDepthVariable, CrocusRecentSweVariable]
+ assert variable_class in [CrocusTotalSweVariable, CrocusDepthVariable, CrocusRecentSweVariable,
+ RecentSnowLoadVariable, TotalSnowLoadVariable]
super().__init__(variable_class, *args, **kwargs)
self.model_name = 'Crocus'
@@ -44,6 +45,19 @@ class CrocusSweTotal(Crocus):
return self.winter_annual_aggregation(time_serie)
+# Create some class that enables to deal directly with the snow load
+
+
+class CrocusSnowLoadTotal(Crocus):
+ def __init__(self, *args, **kwargs):
+ Crocus.__init__(self, TotalSnowLoadVariable, *args, **kwargs)
+
+
+class CrocusSnowLoad3Days(CrocusSweTotal):
+ def __init__(self, *args, **kwargs):
+ Crocus.__init__(self, RecentSnowLoadVariable, *args, **kwargs)
+
+
class ExtendedCrocusSweTotal(AbstractExtendedStudy, CrocusSweTotal):
pass
diff --git a/experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py b/experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py
index d9f63484..bd2a9b27 100644
--- a/experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py
+++ b/experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py
@@ -27,6 +27,22 @@ class CrocusRecentSweVariable(CrocusTotalSweVariable):
return 'SWE_3DY_ISBA'
+class AbstractSnowLoadVariable(CrocusVariable):
+ UNIT = 'kN $m^{-2}$'
+ snow_load_multiplication_factor = 9.81 / 1000
+
+ @property
+ def daily_time_serie_array(self) -> np.ndarray:
+ return self.snow_load_multiplication_factor * super().daily_time_serie_array
+
+
+class RecentSnowLoadVariable(AbstractSnowLoadVariable, CrocusRecentSweVariable):
+ NAME = 'Snow load last 3 days'
+
+class TotalSnowLoadVariable(AbstractSnowLoadVariable, CrocusTotalSweVariable):
+ NAME = 'Snow load total'
+
+
class CrocusDepthVariable(CrocusVariable):
NAME = 'Snow Depth'
UNIT = 'm'
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py
index aeb8f63e..5a5ef8b3 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py
@@ -6,7 +6,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
from experiment.trend_analysis.abstract_score import MannKendall
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusDepth, CrocusSweTotal, ExtendedCrocusDepth, \
- ExtendedCrocusSweTotal, CrocusDaysWithSnowOnGround, CrocusSwe3Days
+ ExtendedCrocusSweTotal, CrocusDaysWithSnowOnGround, CrocusSwe3Days, CrocusSnowLoad3Days, CrocusSnowLoadTotal
from experiment.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, ExtendedSafranSnowfall, \
SafranRainfall, \
SafranTemperature, SafranTotalPrecip
@@ -23,9 +23,11 @@ SCM_STUDIES_NAMES = [get_display_name_from_object_type(k) for k in SCM_STUDIES]
SCM_STUDY_NAME_TO_SCM_STUDY = dict(zip(SCM_STUDIES_NAMES, SCM_STUDIES))
SCM_STUDY_CLASS_TO_ABBREVIATION = {
SafranSnowfall: 'SF3',
- CrocusSweTotal: 'SWET',
+ CrocusSweTotal: 'SWE',
CrocusSwe3Days: 'SWE3',
CrocusDepth: 'SD',
+ CrocusSnowLoadTotal: 'SL',
+ CrocusSnowLoad3Days: 'SL3',
}
altitude_massif_name_and_study_class_for_poster = [
@@ -40,6 +42,12 @@ altitude_massif_name_and_study_class_for_poster_evan = [
(2700, 'Parpaillon', CrocusSwe3Days),
]
+altitude_massif_name_and_study_class_for_committee= [
+ (900, 'Chartreuse', CrocusSnowLoad3Days),
+ (1800, 'Vanoise', CrocusSnowLoad3Days),
+ (2700, 'Parpaillon', CrocusSnowLoad3Days),
+]
+
SCM_STUDY_NAME_TO_ABBREVIATION = {get_display_name_from_object_type(k): v for k, v in
SCM_STUDY_CLASS_TO_ABBREVIATION.items()}
@@ -250,7 +258,8 @@ def max_graph_annual_maxima_poster():
choice_tuple = [
altitude_massif_name_and_study_class_for_poster,
altitude_massif_name_and_study_class_for_poster_evan,
- ][1]
+ altitude_massif_name_and_study_class_for_committee,
+ ][2]
for altitude, massif_name, study_class in choice_tuple:
for study in study_iterator_global([study_class], altitudes=[altitude]):
study_visualizer = StudyVisualizer(study, save_to_file=save_to_file,
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
index 727af31e..505a39d1 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
@@ -625,7 +625,10 @@ class StudyVisualizer(VisualizationParameters):
ax = plt.gca()
x, y = self.smooth_maxima_x_y(massif_names.index(massif_name))
ax.plot(x, y, color=color, linewidth=5)
- ax.set_ylabel('{} (in {})'.format(snow_abbreviation, self.study.variable_unit), color=color, fontsize=15)
+ # ax.set_ylabel('{} (in {})'.format(snow_abbreviation, self.study.variable_unit), color=color, fontsize=15)
+ ax.set_ylabel('{} (in {})'.format(snow_abbreviation, self.study.variable_unit), fontsize=15)
+ ax.set_xlabel('years', fontsize=15)
+ ax.set_title('{} at {} m'.format(massif_name, altitude))
ax.xaxis.set_ticks(x[2::10])
ax.tick_params(axis='both', which='major', labelsize=13)
@@ -824,11 +827,15 @@ class StudyVisualizer(VisualizationParameters):
if not self.only_one_graph:
filename += "{}".format('_'.join(self.plot_name.split())) + '_'
filename += specific_title
- filepath = op.join(self.study.result_full_path, filename + '.png')
- dirname = op.dirname(filepath)
- if not op.exists(dirname):
- os.makedirs(dirname, exist_ok=True)
- plt.savefig(filepath)
+ self.savefig_in_results(filename)
+
+ @classmethod
+ def savefig_in_results(cls, filename):
+ filepath = op.join(AbstractStudy.result_full_path, filename + '.png')
+ dirname = op.dirname(filepath)
+ if not op.exists(dirname):
+ os.makedirs(dirname, exist_ok=True)
+ plt.savefig(filepath)
def visualize_independent_margin_fits(self, threshold=None, axes=None, show=True):
# Fit either a GEV or a GPD
diff --git a/experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py b/experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py
index ac9ddb19..2d81708a 100644
--- a/experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py
+++ b/experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py
@@ -17,7 +17,7 @@ mpl.rcParams['hatch.linewidth'] = 0.3
def main_poster_A_non_stationary_model_choice():
nb = 1
- for altitude in POSTER_ALTITUDES[:nb]:
+ for altitude in POSTER_ALTITUDES[:]:
for trend_test_class in [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest][-nb:]:
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
exact_starting_year=1958, reduce_strength_array=False,
@@ -133,7 +133,7 @@ def main_poster_D_other_quantities_analysis():
if __name__ == '__main__':
- # main_poster_A_non_stationary_model_choice()
- main_poster_B_starting_years_analysis()
+ main_poster_A_non_stationary_model_choice()
+ # main_poster_B_starting_years_analysis()
# main_poster_C_orientation_analysis()
# main_poster_D_other_quantities_analysis()
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
index e151009e..958dae9f 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
@@ -4,7 +4,7 @@ from enum import Enum
class ConfidenceIntervalMethodFromExtremes(Enum):
# Confidence interval from the ci function
ci_bayes = 0
- ci_normal = 1
+ ci_mle = 1
ci_boot = 2
ci_proflik = 3
# Confidence interval from my functions
@@ -12,7 +12,7 @@ class ConfidenceIntervalMethodFromExtremes(Enum):
ci_method_to_method_name = {
- ConfidenceIntervalMethodFromExtremes.ci_normal: 'normal',
+ ConfidenceIntervalMethodFromExtremes.ci_mle: 'normal',
ConfidenceIntervalMethodFromExtremes.ci_boot: 'boot',
ConfidenceIntervalMethodFromExtremes.ci_proflik: 'proflik',
}
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
index 08bbd976..68708932 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
@@ -1,4 +1,5 @@
import numpy as np
+import rpy2
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_result_from_extremes import \
AbstractResultFromExtremes
@@ -23,7 +24,10 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
'method': method_name,
# xrange = NULL, nint = 20
}
- res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
+ try:
+ res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
+ except rpy2.rinterface.RRuntimeError:
+ return np.nan, (np.nan, np.nan)
if self.is_non_stationary:
a = np.array(res)[0]
lower, mean_estimate, upper, _ = a
diff --git a/test/test_extreme_fit/test_model/test_confidence_interval.py b/test/test_extreme_fit/test_model/test_confidence_interval.py
index 9b6a5b9d..7ba6baf0 100644
--- a/test/test_extreme_fit/test_model/test_confidence_interval.py
+++ b/test/test_extreme_fit/test_model/test_confidence_interval.py
@@ -74,7 +74,7 @@ class TestConfidenceInterval(unittest.TestCase):
def test_ci_normal(self):
self.fit_method = TemporalMarginFitMethod.extremes_fevd_mle
- self.ci_method = ConfidenceIntervalMethodFromExtremes.ci_normal
+ self.ci_method = ConfidenceIntervalMethodFromExtremes.ci_mle
self.model_class_to_triplet= {
StationaryTemporalModel: (-4.703945484843988, 30.482318639674023, 65.66858276419204),
NonStationaryLocationTemporalModel: (-4.223086740397132, 30.29842988666537, 64.81994651372787),
diff --git a/thesis_report/gev_plot.py b/thesis_report/gev_plot.py
index 144ac914..91a9d71a 100644
--- a/thesis_report/gev_plot.py
+++ b/thesis_report/gev_plot.py
@@ -38,14 +38,23 @@ def gev_plot_big():
def gev_plot_big_non_stationary_location():
lim = 5
x = np.linspace(-lim, lim, 100)
- scale, shape = 1, 1
- locs = [-1, 1]
- colors = ['red', 'green']
- y = r.dgev(x, 0.0, scale, shape)
- plt.plot(x, y, label='distribution of $Y(1968)$', linewidth=5)
+ scale, shape = 1, 0.0
+ locs = [1, 2, 3]
+ inverse_loc_with_scale = True
+ colors = ['red','k', 'green']
+ greek_leeter = ' $\{}_1'.format('mu' if not inverse_loc_with_scale else 'sigma')
+ plt.title('Density for the distribution of Y(t) with different{}$'.format(greek_leeter))
+ template = greek_leeter + '{} 0$'
for loc, color in zip(locs, colors):
- sign_str = '>' if loc > 0 else '<'
- label = 'distribution of $Y(1969)$ with $\mu_1 {} 0$'.format(sign_str)
+ if loc == locs[1]:
+ sign_str = '='
+ elif loc == locs[2]:
+ sign_str = '>'
+ else:
+ sign_str = '<'
+ label = template.format(sign_str)
+ if inverse_loc_with_scale:
+ loc, scale = scale, loc
y = r.dgev(x, loc, scale, shape)
plt.plot(x, y, label=label, linewidth=5, color=color)
plt.legend(prop={'size': 20})
--
GitLab