diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index af7bb0a0fdfc2c17bb0ffc03fa5db447886a1a1e..c99cd54855bc127b22326b2e8e8a1b5f2569efcc 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -1,12 +1,17 @@
-
import time
from typing import List
import matplotlib as mpl
+from projects.altitude_spatial_model.altitudes_fit.plot_histogram_altitude_studies import \
+ plot_histogram_all_trends_against_altitudes, plot_histogram_all_models_against_altitudes
+
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
+from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import load_visualizer_list
+
from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups
@@ -22,29 +27,12 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_s
def main():
- altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][4:6]
- # todo: l ecart pour les saisons de l automne ou de sprint
- # vient probablement de certains zéros pas pris en compte pour le fit,
- # mais pris en compte pour le calcul de mon aic
- # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
- # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900][:]
- # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
- # altitudes = [600, 900, 1200, 1500, 1800, 2100][:]
- altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200, 4500, 4800][:]
- # altitudes = [1500, 1800][:]
- study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
- study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
- SafranPrecipitation7Days][:]
- study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranPrecipitation1Day
- , SafranPrecipitation3Days][:1]
- altitudes = [1800, 2100, 2400]
study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:1]
# study_classes = [SafranPrecipitation1Day][:1]
- # Common parameters
- # altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]
massif_names = None
# massif_names = ['Mercantour', 'Vercors', 'Ubaye']
+
seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
main_loop(altitudes_for_groups[:], massif_names, seasons, study_classes)
@@ -55,69 +43,29 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes):
assert isinstance(altitudes_list[0], List)
for season in seasons:
for study_class in study_classes:
- # if issubclass(study_class, SimulationStudy):
- # for ensemble_idx in list(range(14))[:1]:
- # studies = AltitudesStudies(study_class, altitudes, season=season,
- # ensemble_idx=ensemble_idx)
- # plot_studies(massif_names, season, studies, study_class)
- # else:
+ print('Inner loop', season, study_class)
visualizer_list = load_visualizer_list(season, study_class, altitudes_list, massif_names)
+ plot_visualizers(massif_names, visualizer_list)
for visualizer in visualizer_list:
- plots(massif_names, season, visualizer)
+ plot_visualizer(massif_names, visualizer)
del visualizer_list
time.sleep(2)
+def plot_visualizers(massif_names, visualizer_list):
+ plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
+ plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
-def load_visualizer_list(season, study_class, altitudes_list, massif_names):
- model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
- visualizer_list = []
- # Load all studies
- for altitudes in altitudes_list:
- studies = AltitudesStudies(study_class, altitudes, season=season)
- visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
- model_classes=model_classes,
- massif_names=massif_names,
- show=False,
- temporal_covariate_for_fit=None,
- # temporal_covariate_for_fit=MeanAlpsTemperatureCovariate,
- )
- visualizer_list.append(visualizer)
- # Compute the max abs for all metrics
- d = {}
- for method_name in AltitudesStudiesVisualizerForNonStationaryModels.moment_names:
- for order in AltitudesStudiesVisualizerForNonStationaryModels.orders:
- c = (method_name, order)
- max_abs = max([
- max([abs(e) for e in v.method_name_and_order_to_d(method_name, order).values()
- ]) for v in visualizer_list])
- d[c] = max_abs
- # Assign the max abs dictionary
- for v in visualizer_list:
- v._method_name_and_order_to_max_abs = d
- # Compute the max abs for the shape parameter
- max_abs_for_shape = max([max([abs(e) for e in v.massif_name_to_shape.values()]) for v in visualizer_list])
- for v in visualizer_list:
- v._max_abs_for_shape = max_abs_for_shape
-
- return visualizer_list
-
-
-def plots(massif_names, season, visualizer):
- studies = visualizer.studies
- print('inner loop', season, type(studies.study))
+def plot_visualizer(massif_names, visualizer):
# Plot time series
- studies.plot_maxima_time_series(massif_names=massif_names)
-
- # Plot moments
+ # visualizer.studies.plot_maxima_time_series(massif_names=massif_names)
+ # Plot moments against altitude
# for std in [True, False][:]:
# for change in [True, False, None]:
# studies.plot_mean_maxima_against_altitude(massif_names=massif_names, std=std, change=change)
-
# Plot the results for the model that minimizes the individual aic
plot_individual_aic(visualizer)
-
# Plot the results for the model that minimizes the total aic
# plot_total_aic(model_classes, visualizer)
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index a963c82d82ee661e6d05ab1eab49715828cd6f16..ae27f5c5cb09d924143ed5429531229270ed6b94 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -1,3 +1,4 @@
+from collections import Counter
from typing import List, Dict
import matplotlib
@@ -116,23 +117,39 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
moment = ' '.join(method_name.split('_'))
moment = moment.replace('moment', '{} in 2019'.format(OneFoldFit.get_moment_str(order=order)))
plot_name = '{}{} for annual maxima of {}'.format(OneFoldFit.folder_for_plots, moment,
- SCM_STUDY_CLASS_TO_ABBREVIATION[
- self.studies.study_class])
+ SCM_STUDY_CLASS_TO_ABBREVIATION[
+ self.studies.study_class])
+
+ massif_name_to_text = self.massif_name_to_best_name
if 'change' in method_name:
- plot_name += '\nbetween {} and {}'.format(2019-50, 2019)
+ plot_name += '\nbetween {} and {}'.format(2019 - 50, 2019)
+ if 'relative' not in method_name:
+ # Put the relative score as text on the plot for the change.
+ massif_name_to_text = {m: ('+' if v > 0 else '') + str(int(v)) + '\%' for m, v in
+ self.method_name_and_order_to_d(self.moment_names[2], order).items()}
+
parenthesis = self.study.variable_unit if 'relative' not in method_name else '\%'
ylabel = '{} ({})'.format(plot_name, parenthesis)
+ max_abs_change = self.method_name_and_order_to_max_abs(method_name, order)
+ add_colorbar = self.add_colorbar
+
+ is_return_level_plot = (self.moment_names.index(method_name) == 0) and (order is None)
+ if is_return_level_plot:
+ add_colorbar = True
+ max_abs_change = None
+ massif_name_to_text = {m: '+- 0' for m in massif_name_to_text.keys()}
+
+ negative_and_positive_values = self.moment_names.index(method_name) > 0
# Plot the map
- a_change_is_displayed = self.moment_names.index(method_name) > 0
self.plot_map(cmap=plt.cm.coolwarm, fit_method=self.fit_method, graduation=10,
label=ylabel, massif_name_to_value=massif_name_to_value,
plot_name=plot_name, add_x_label=True,
- negative_and_positive_values=a_change_is_displayed,
+ negative_and_positive_values=negative_and_positive_values,
altitude=self.altitude_group.reference_altitude,
- add_colorbar=self.add_colorbar,
- max_abs_change=self.method_name_and_order_to_max_abs(method_name, order),
- massif_name_to_text=self.massif_name_to_name,
+ add_colorbar=add_colorbar,
+ max_abs_change=max_abs_change,
+ massif_name_to_text=massif_name_to_text,
)
@property
@@ -175,7 +192,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items()}
@property
- def massif_name_to_name(self):
+ def massif_name_to_best_name(self):
return {massif_name: one_fold_fit.best_name
for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items()}
@@ -211,23 +228,23 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
negative_and_positive_values = (max(values) > 0) and (min(values) < 0)
raise NotImplementedError
self.plot_map(massif_name_to_value=massif_name_to_best_coef,
- label='{}Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots,
- coef_name,
- SCM_STUDY_CLASS_TO_ABBREVIATION[
- type(self.study)],
- self.study.variable_unit),
- add_x_label=False, graduation=graduation, massif_name_to_text=self.massif_name_to_name,
- negative_and_positive_values=negative_and_positive_values)
+ label='{}Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots,
+ coef_name,
+ SCM_STUDY_CLASS_TO_ABBREVIATION[
+ type(self.study)],
+ self.study.variable_unit),
+ add_x_label=False, graduation=graduation, massif_name_to_text=self.massif_name_to_best_name,
+ negative_and_positive_values=negative_and_positive_values)
def plot_shape_map(self):
label = 'Shape parameter for {} maxima of {} in 2019'.format(self.study.season_name,
- SCM_STUDY_CLASS_TO_ABBREVIATION[
- type(self.study)])
+ SCM_STUDY_CLASS_TO_ABBREVIATION[
+ type(self.study)])
self.plot_map(massif_name_to_value=self.massif_name_to_shape,
label=label,
plot_name=label,
- add_x_label=True, graduation=0.1, massif_name_to_text=self.massif_name_to_name,
+ add_x_label=True, graduation=0.1, massif_name_to_text=self.massif_name_to_best_name,
cmap=matplotlib.cm.get_cmap('BrBG_r'),
altitude=self.altitude_group.reference_altitude,
add_colorbar=self.add_colorbar,
@@ -350,3 +367,53 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
plot_name = 'Switch altitude'
self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
plt.close()
+
+ def all_trends(self, massif_names):
+ """return percents which contain decrease, significant decrease, increase, significant increase percentages"""
+ valid_massif_names = self.get_valid_names(massif_names)
+
+ nb_valid_massif_names = len(valid_massif_names)
+ nbs = np.zeros(4)
+ relative_changes = [[], []]
+ for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
+ if m in valid_massif_names]:
+ if one_fold.change_for_reference_altitude == 0:
+ continue
+ # Compute relative changes
+ relative_changes[0].append(one_fold.relative_change_for_reference_altitude)
+ if one_fold.is_significant:
+ relative_changes[1].append(one_fold.relative_change_for_reference_altitude)
+ # Compute nbs
+ idx = 0 if one_fold.change_for_reference_altitude < 0 else 2
+ nbs[idx] += 1
+ if one_fold.is_significant:
+ nbs[idx + 1] += 1
+
+ percents = 100 * nbs / nb_valid_massif_names
+ mean_relative_changes = [np.mean(l) for l in relative_changes]
+ return [nb_valid_massif_names] + mean_relative_changes + list(percents)
+
+ def get_valid_names(self, massif_names):
+ valid_massif_names = set(self.massif_name_to_one_fold_fit.keys())
+ if massif_names is not None:
+ valid_massif_names = valid_massif_names.intersection(set(massif_names))
+ return valid_massif_names
+
+ def model_name_to_percentages(self, massif_names):
+ valid_massif_names = self.get_valid_names(massif_names)
+ nb_valid_massif_names = len(valid_massif_names)
+ best_names = [one_fold_fit.best_estimator.margin_model.name_str
+ for m, one_fold_fit in self.massif_name_to_one_fold_fit.items()
+ if m in valid_massif_names]
+ counter = Counter(best_names)
+ d = {name: 100 * c / nb_valid_massif_names for name, c in counter.items()}
+ # Add 0 for the name not present
+ for name in self.model_names:
+ if name not in d:
+ d[name] = 0
+ return d
+
+ @property
+ def model_names(self):
+ massif_name = list(self.massif_name_to_one_fold_fit.keys())[0]
+ return self.massif_name_to_one_fold_fit[massif_name].model_names
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index cb2d665122942023f1029fca96ef092f9732e7e8..badae11f99d6f44c2668d1142ab4065c5c65c897 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -88,6 +88,14 @@ class OneFoldFit(object):
def moment(self, altitudes, year=2019, order=1):
return [self.get_moment(altitude, year, order) for altitude in altitudes]
+ @property
+ def change_for_reference_altitude(self) -> float:
+ return self.changes_for_moment(altitudes=[self.altitude_plot])[0]
+
+ @property
+ def relative_change_for_reference_altitude(self) -> float:
+ return self.relative_changes_for_moment(altitudes=[self.altitude_plot])[0]
+
def changes_for_moment(self, altitudes, year=2019, nb_years=50, order=1):
changes = []
for altitude in altitudes:
@@ -185,6 +193,10 @@ class OneFoldFit(object):
except (TypeError, KeyError):
return None
+ @property
+ def model_names(self):
+ return [e.margin_model.name_str for e in self.sorted_estimators]
+
@property
def best_name(self):
name = self.best_estimator.margin_model.name_str
diff --git a/projects/altitude_spatial_model/altitudes_fit/plot_histogram_altitude_studies.py b/projects/altitude_spatial_model/altitudes_fit/plot_histogram_altitude_studies.py
new file mode 100644
index 0000000000000000000000000000000000000000..0fe005e6fde945d1bcd7f87d68f28d78c6a9f5df
--- /dev/null
+++ b/projects/altitude_spatial_model/altitudes_fit/plot_histogram_altitude_studies.py
@@ -0,0 +1,115 @@
+from typing import List
+
+import numpy as np
+
+import matplotlib.pyplot as plt
+
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
+ AltitudesStudiesVisualizerForNonStationaryModels
+
+
+def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: List[
+ AltitudesStudiesVisualizerForNonStationaryModels]):
+ visualizer = visualizer_list[0]
+ model_names = visualizer.model_names
+ model_name_to_percentages = {model_name: [] for model_name in model_names}
+ for v in visualizer_list:
+ for model_name, percentage in v.model_name_to_percentages(massif_names).items():
+ model_name_to_percentages[model_name].append(percentage)
+ model_name_to_mean_percentage = {m: np.mean(a) for m, a in model_name_to_percentages.items()}
+ sorted_model_names = sorted(model_names, key=lambda m: model_name_to_mean_percentage[m], reverse=True)
+
+ # Plot part
+ ax = plt.gca()
+ width = 5
+ size = 8
+ legend_fontsize = 10
+ labelsize = 10
+ linewidth = 1
+ tick_list = np.array([((len(visualizer_list) + 2) * i + (1 + len(visualizer_list) / 2)) * width
+ for i in range(len(sorted_model_names))])
+ for tick_middle, model_name in zip(tick_list, sorted_model_names):
+ x_shifted = [tick_middle + width * shift / 2 for shift in range(-3, 5, 2)]
+ percentages = model_name_to_percentages[model_name]
+ colors = ['white', 'yellow', 'orange', 'red']
+ labels = ['{} m (\% out of {} massifs)'.format(v.altitude_group.reference_altitude,
+ len(v.get_valid_names(massif_names))) for v in visualizer_list]
+ for x, color, percentage, label in zip(x_shifted, colors, percentages, labels):
+ ax.bar([x], [percentage], width=width, label=label,
+ linewidth=linewidth, edgecolor='black', color=color)
+
+ handles, labels = ax.get_legend_handles_labels()
+ ax.legend(handles[:len(visualizer_list)], labels[:len(visualizer_list)], prop={'size': size})
+ ax.set_xticklabels(sorted_model_names)
+ ax.set_xticks(tick_list)
+ ax.set_ylabel('Percentage of massifs (\%) ', fontsize=legend_fontsize)
+ ax.set_xlabel('Models', fontsize=legend_fontsize)
+ ax.set_ylim(bottom=0)
+ ax.yaxis.grid()
+ ax.tick_params(axis='both', which='major', labelsize=labelsize)
+
+ visualizer.plot_name = 'All models'
+ visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
+ plt.close()
+
+
+def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list):
+ all_trends = [v.all_trends(massif_names) for v in visualizer_list]
+ nb_massifs, mean_relative_change, mean_significant_relative_change, *all_l = zip(*all_trends)
+
+ plt.close()
+ ax = plt.gca()
+ width = 5
+ size = 8
+ legend_fontsize = 10
+ labelsize = 10
+ linewidth = 3
+ x = np.array([3 * width * (i + 1) for i in range(len(nb_massifs))])
+
+ colors = ['blue', 'darkblue', 'red', 'darkred']
+ labels = []
+ for suffix in ['Decrease', 'Increase']:
+ for prefix in ['Non significant', 'Significant']:
+ labels.append('{} {}'.format(prefix, suffix))
+ for l, color, label in zip(all_l, colors, labels):
+ x_shifted = x - width / 2 if 'blue' in color else x + width / 2
+ ax.bar(x_shifted, l, width=width, color=color, edgecolor=color, label=label,
+ linewidth=linewidth)
+ ax.legend(loc='upper left', prop={'size': size})
+ ax.set_ylabel('Percentage of massifs (\%) ', fontsize=legend_fontsize)
+ ax.set_xlabel('Altitudes', fontsize=legend_fontsize)
+ ax.tick_params(axis='both', which='major', labelsize=labelsize)
+ ax.set_xticks(x)
+ ax.yaxis.grid()
+ ax.set_xticklabels([str(v.altitude_group.reference_altitude) for v in visualizer_list])
+
+ plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x)
+
+ # PLot mean relative change against altitude
+ ax_twinx = ax.twinx()
+ colors = ['grey', 'black']
+ relative_changes = [mean_relative_change, mean_significant_relative_change]
+ labels = ['All massifs with a trend', 'Only massifs with a significant trend']
+ for color, y, label in zip(colors, relative_changes, labels):
+ print(x, y)
+ ax_twinx.plot(x, y, label=label, color=color)
+ ax_twinx.legend(loc='upper right', prop={'size': size})
+
+ ax_twinx.set_ylabel('Mean relative change average on massif with trends', fontsize=legend_fontsize)
+
+ visualizer = visualizer_list[0]
+ visualizer.plot_name = 'All trends'
+ visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
+
+ plt.close()
+
+
+def plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x):
+ # Plot number of massifs on the upper axis
+ ax_twiny = ax.twiny()
+ ax_twiny.plot(x, [0 for _ in x], linewidth=0)
+ ax_twiny.set_xticks(x)
+ ax_twiny.tick_params(labelsize=labelsize)
+ ax_twiny.set_xticklabels(nb_massifs)
+ ax_twiny.set_xlim(ax.get_xlim())
+ ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=legend_fontsize)
diff --git a/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py b/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
new file mode 100644
index 0000000000000000000000000000000000000000..343de5879e157f1c3124c477e9b2281b15aefa10
--- /dev/null
+++ b/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
@@ -0,0 +1,49 @@
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
+ ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
+ AltitudesStudiesVisualizerForNonStationaryModels
+
+
+def load_visualizer_list(season, study_class, altitudes_list, massif_names):
+ model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+ visualizer_list = []
+ # Load all studies
+ for altitudes in altitudes_list:
+ # if issubclass(study_class, SimulationStudy):
+ # for ensemble_idx in list(range(14))[:1]:
+ # studies = AltitudesStudies(study_class, altitudes, season=season,
+ # ensemble_idx=ensemble_idx)
+ # plot_studies(massif_names, season, studies, study_class)
+ # else:
+ studies = AltitudesStudies(study_class, altitudes, season=season)
+ visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
+ model_classes=model_classes,
+ massif_names=massif_names,
+ show=False,
+ temporal_covariate_for_fit=None,
+ # temporal_covariate_for_fit=MeanAlpsTemperatureCovariate,
+ )
+ visualizer_list.append(visualizer)
+ compute_and_assign_max_abs(visualizer_list)
+
+ return visualizer_list
+
+
+def compute_and_assign_max_abs(visualizer_list):
+ # Compute the max abs for all metrics
+ d = {}
+ for method_name in AltitudesStudiesVisualizerForNonStationaryModels.moment_names:
+ for order in AltitudesStudiesVisualizerForNonStationaryModels.orders:
+ c = (method_name, order)
+ max_abs = max([
+ max([abs(e) for e in v.method_name_and_order_to_d(method_name, order).values()
+ ]) for v in visualizer_list])
+ d[c] = max_abs
+ # Assign the max abs dictionary
+ for v in visualizer_list:
+ v._method_name_and_order_to_max_abs = d
+ # Compute the max abs for the shape parameter
+ max_abs_for_shape = max([max([abs(e) for e in v.massif_name_to_shape.values()]) for v in visualizer_list])
+ for v in visualizer_list:
+ v._max_abs_for_shape = max_abs_for_shape
diff --git a/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_study_visualizer.py b/projects/altitude_spatial_model/preliminary_analysis.py
similarity index 94%
rename from projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_study_visualizer.py
rename to projects/altitude_spatial_model/preliminary_analysis.py
index e771f9476d30786a3cf33250f8072790d7994857..bcc736d53b50388fb6b856913567cc16e1710c95 100644
--- a/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_study_visualizer.py
+++ b/projects/altitude_spatial_model/preliminary_analysis.py
@@ -29,25 +29,29 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
massif_name_to_r2_score = {}
for massif_name in self.study.all_massif_names()[:]:
linear_coef, _, r2 = self._plot_gev_params_against_altitude_one_massif(ax, massif_name, param_name)
- massif_name_to_linear_coef[massif_name] = linear_coef[0]
+ massif_name_to_linear_coef[massif_name] = 1000 * linear_coef[0]
massif_name_to_r2_score[massif_name] = str(round(r2, 2))
print(massif_name_to_linear_coef, massif_name_to_r2_score)
# Plot change against altitude
ax.legend(prop={'size': 7}, ncol=3)
ax.set_xlabel('Altitude')
ax.set_ylabel(GevParams.full_name_from_param_name(param_name) + ' parameter for a stationary GEV distribution')
- plot_name = '{} change /with altitude'.format(param_name)
+ plot_name = '{} change with altitude'.format(param_name)
self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
ax.clear()
# Plot map of slope for each massif
visualizer = StudyVisualizer(study=self.study, show=False, save_to_file=True)
- ylabel = 'Linear slope for the {} parameter against the altitude'.format(param_name)
+ ylabel = 'Linear slope for the {} parameter (change every 1000 m of altitude)'.format(param_name)
+ gev_param_name_to_graduation = {
+ GevParams.LOC: 5,
+ GevParams.SCALE: 1,
+ GevParams.SHAPE: 0.1,
+ }
visualizer.plot_map(cmap=plt.cm.coolwarm, fit_method=self.study.fit_method,
- graduation=1,
+ graduation=gev_param_name_to_graduation[param_name],
label=ylabel, massif_name_to_value=massif_name_to_linear_coef,
plot_name=ylabel, add_x_label=False,
- # negative_and_positive_values=param_name == GevParams.SHAPE,
- negative_and_positive_values=True,
+ negative_and_positive_values=param_name == GevParams.SHAPE,
add_colorbar=True,
massif_name_to_text=massif_name_to_r2_score,
)
@@ -190,7 +194,7 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
if __name__ == '__main__':
altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300]
- altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]
+ altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
# altitudes = paper_altitudes
# altitudes = [1800, 2100]
visualizer = PointwiseGevStudyVisualizer(SafranSnowfall1Day, altitudes=altitudes)
diff --git a/projects/seasonal_analysis/__init__.py b/projects/seasonal_analysis/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/projects/altitude_spatial_model/seasonal_analysis/main_season_repartition_of_maxima.py b/projects/seasonal_analysis/main_season_repartition_of_maxima.py
similarity index 100%
rename from projects/altitude_spatial_model/seasonal_analysis/main_season_repartition_of_maxima.py
rename to projects/seasonal_analysis/main_season_repartition_of_maxima.py
diff --git a/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py b/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py
index 3db108e7c76cfd5ac00c636a79b07ea53ce28cdb..67b6ea1546a725c478013bb369d69a5d139d6298 100644
--- a/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py
+++ b/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py
@@ -65,7 +65,7 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
self.fit_method = MarginFitMethod.extremes_fevd_mle
- def function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(self, model_class):
+ def function_test_gev_spatio_temporal_margin_fit_non_stationary(self, model_class):
# Create estimator
estimator = fitted_linear_margin_estimator_short(model_class=model_class,
dataset=self.dataset,
@@ -76,32 +76,36 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
# Assert we can compute the return level
- covariate_for_return_level = np.array([400, 25])[::1]
+ covariate_for_return_level = np.array([400, 20])[::1]
confidence_interval = EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator,
ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
temporal_covariate=covariate_for_return_level)
return_level = estimator.function_from_fit.get_params(covariate_for_return_level).return_level(50)
+ print("my return level", return_level)
self.assertAlmostEqual(return_level, confidence_interval.mean_estimate)
self.assertFalse(np.isnan(confidence_interval.confidence_interval[0]))
self.assertFalse(np.isnan(confidence_interval.confidence_interval[1]))
def test_gev_spatio_temporal_margin_fit_1(self):
- self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(StationaryAltitudinal)
+ self.function_test_gev_spatio_temporal_margin_fit_non_stationary(StationaryAltitudinal)
+ #
def test_gev_spatio_temporal_margin_fit_1_bis(self):
- self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(AltitudinalShapeLinearTimeStationary)
+ self.function_test_gev_spatio_temporal_margin_fit_non_stationary(AltitudinalShapeLinearTimeStationary)
# def test_gev_spatio_temporal_margin_fit_2(self):
- # self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(
+ # # first model with both a time and altitude non stationarity
+ # self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
# AltitudinalShapeConstantTimeLocationLinear)
- #
+
# def test_gev_spatio_temporal_margin_fit_3(self):
- # self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(
+ # self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
# AltitudinalShapeLinearTimeLocationLinear)
#
# def test_gev_spatio_temporal_margin_fit_4(self):
- # self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(
+ # self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
# AltitudinalShapeLinearTimeLocScaleLinear)
+
if __name__ == '__main__':
unittest.main()
diff --git a/test/test_projects/test_exceeding_snow_loads/test_results.py b/test/test_projects/test_exceeding_snow_loads/test_results.py
index b7434aa278bc4658ecd9c8fdf4afda3d7c8c318d..e5f110727ba6e6b0171fd6a7b09e3ccb1bf3f201 100644
--- a/test/test_projects/test_exceeding_snow_loads/test_results.py
+++ b/test/test_projects/test_exceeding_snow_loads/test_results.py
@@ -8,11 +8,12 @@ from projects.exceeding_snow_loads.section_results.plot_selection_curves import
from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_curves, plot_trend_map
from projects.exceeding_snow_loads.section_results.plot_uncertainty_curves import plot_uncertainty_massifs
from projects.exceeding_snow_loads.section_results.plot_uncertainty_histogram import plot_uncertainty_histogram
-
+import matplotlib.pyplot as plt
class TestResults(unittest.TestCase):
def test_run_intermediate_results(self):
+ plt.close()
# Load data
altitude_to_visualizer = load_altitude_to_visualizer(altitudes=[300, 600], massif_names=None,
study_class=CrocusSnowLoadTotal,