From 249ea6a49b2d8ae17a55c290d092298e7d6c299a Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 29 Jul 2020 16:35:43 +0200
Subject: [PATCH] [contrasting] add cubic models
---
.../gev_altitudinal_models.py | 16 +++++++++++++++-
.../gumbel_altitudinal_models.py | 2 +-
.../polynomial_margin_model.py | 4 ++++
.../spatio_temporal_polynomial_model.py | 2 +-
.../polynomial_margin_model/utils.py | 17 ++++++++++++++++-
.../model/result_from_model_fit/utils.py | 3 ---
.../altitudes_fit/main_altitudes_studies.py | 16 ++++++++++------
...ies_visualizer_for_non_stationary_models.py | 8 ++++----
.../one_fold_analysis/one_fold_fit.py | 3 ++-
.../main_result_trends_and_return_levels.py | 2 +-
.../plot_uncertainty_histogram.py | 18 +++++++++++-------
11 files changed, 65 insertions(+), 26 deletions(-)
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py
index 218b2780..8d8068f7 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py
@@ -36,7 +36,7 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
name += self.dim_to_str_number(GevParams.LOC, self.coordinates.idx_temporal_coordinates)
name += self.dim_to_str_number(GevParams.SCALE, self.coordinates.idx_temporal_coordinates)
if isinstance(self, AbstractAddCrossTermForLocation):
- name += 'l'
+ name += 'x'
if isinstance(self, AbstractAddCrossTermForScale):
name += 's'
return name
@@ -81,6 +81,15 @@ class NonStationaryAltitudinalLocationQuadratic(AbstractAltitudinalModel):
GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
}
+class NonStationaryAltitudinalLocationCubic(AbstractAltitudinalModel):
+
+ @property
+ def param_name_to_list_dim_and_degree(self):
+ return {
+ GevParams.LOC: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 3)],
+ GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1)]
+ }
+
class NonStationaryAltitudinalLocationLinearScaleLinear(AbstractAltitudinalModel):
@@ -141,6 +150,10 @@ class NonStationaryAltitudinalLocationQuadraticCrossTermForLocation(AbstractAddC
NonStationaryAltitudinalLocationQuadratic):
pass
+class NonStationaryAltitudinalLocationCubicCrossTermForLocation(AbstractAddCrossTermForLocation,
+ NonStationaryAltitudinalLocationCubic,
+ ):
+ pass
class NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation,
NonStationaryAltitudinalLocationLinearScaleLinear):
@@ -155,3 +168,4 @@ class NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation(A
class NonStationaryAltitudinalScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation,
NonStationaryAltitudinalScaleLinear):
pass
+
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py b/extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py
index 6ff8e3b6..e2a6e7dc 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py
@@ -20,7 +20,7 @@ class AbstractGumbelAltitudinalModel(AbstractAltitudinalModel):
def __init__(self, coordinates: AbstractCoordinates, params_user=None, starting_point=None,
fit_method=MarginFitMethod.extremes_fevd_mle, nb_iterations_for_bayesian_fit=5000,
- params_initial_fit_bayesian=None, params_class=GevParams, max_degree=2):
+ params_initial_fit_bayesian=None, params_class=GevParams, max_degree=3):
super().__init__(coordinates, params_user, starting_point, fit_method, nb_iterations_for_bayesian_fit,
params_initial_fit_bayesian, "Gumbel", params_class, max_degree)
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/polynomial_margin_model.py b/extreme_fit/model/margin_model/polynomial_margin_model/polynomial_margin_model.py
index 7d08ad61..0f044d8b 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/polynomial_margin_model.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/polynomial_margin_model.py
@@ -40,6 +40,10 @@ class PolynomialMarginModel(AbstractTemporalLinearMarginModel):
assert dims.index(self.coordinates.idx_x_coordinates) == 0
if self.coordinates.has_temporal_coordinates and self.coordinates.idx_temporal_coordinates in dims:
assert dims.index(self.coordinates.idx_temporal_coordinates) == len(dims) - 1
+ # Assert that the degree are inferior to the max degree
+ for list_dim_and_degree in param_name_to_list_dim_and_degree.values():
+ for _, max_degree in list_dim_and_degree:
+ assert max_degree <= self.max_degree, 'Max degree (={}) specified is too high'.format(max_degree)
# Load param_name_to_polynomial_all_coef
param_name_to_polynomial_all_coef = self.param_name_to_polynomial_all_coef(
param_name_to_list_dim_and_degree=param_name_to_list_dim_and_degree,
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py b/extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py
index fefeb13b..4de87b47 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py
@@ -8,7 +8,7 @@ class AbstractSpatioTemporalPolynomialModel(PolynomialMarginModel):
def __init__(self, coordinates: AbstractCoordinates, params_user=None, starting_point=None,
fit_method=MarginFitMethod.extremes_fevd_mle, nb_iterations_for_bayesian_fit=5000,
- params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=2):
+ params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=3):
super().__init__(coordinates, params_user, starting_point, fit_method, nb_iterations_for_bayesian_fit,
params_initial_fit_bayesian, type_for_MLE, params_class, max_degree)
self.drop_duplicates = False
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
index ae5b147d..5906a9a7 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
@@ -5,7 +5,8 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_mode
NonStationaryAltitudinalLocationQuadraticCrossTermForLocation, \
NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation, \
NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, NonStationaryAltitudinalScaleLinear, \
- NonStationaryAltitudinalScaleLinearCrossTermForLocation
+ NonStationaryAltitudinalScaleLinearCrossTermForLocation, NonStationaryAltitudinalLocationCubicCrossTermForLocation, \
+ NonStationaryAltitudinalLocationCubic
from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_cross_term_in_scale import \
NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForScale, \
NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForScale, \
@@ -35,6 +36,20 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
NonStationaryLocationSpatioTemporalLinearityModelAssertError3, NonStationaryLocationSpatioTemporalLinearityModel6
+ALTITUDINAL_GEV_MODELS_LOCATION = [
+ StationaryAltitudinal,
+ NonStationaryAltitudinalLocationLinear,
+ NonStationaryAltitudinalLocationQuadratic,
+ NonStationaryAltitudinalLocationCubic,
+ # First order cross term
+ NonStationaryCrossTermForLocation,
+ NonStationaryAltitudinalLocationLinearCrossTermForLocation,
+ NonStationaryAltitudinalLocationQuadraticCrossTermForLocation,
+ NonStationaryAltitudinalLocationCubicCrossTermForLocation
+
+]
+
+
ALTITUDINAL_GEV_MODELS = [
StationaryAltitudinal,
NonStationaryAltitudinalScaleLinear,
diff --git a/extreme_fit/model/result_from_model_fit/utils.py b/extreme_fit/model/result_from_model_fit/utils.py
index b0ff5edc..8a51ceee 100644
--- a/extreme_fit/model/result_from_model_fit/utils.py
+++ b/extreme_fit/model/result_from_model_fit/utils.py
@@ -33,9 +33,6 @@ def get_margin_coef_ordered_dict(param_name_to_dims, mle_values, type_for_mle="G
coef_dict[coef_name] = mle_values[i]
i += 1
else:
- # We found (thanks to the test) that time was the first parameter when len(param_name_to_dims) == 1
- # otherwise time is the second parameter in the order of the mle parameters
- # inverted_dims = dims[::-1] if len(param_name_to_dims) == 1 else dims
for dim, max_degree in dims[:]:
coefficient_name = LinearCoef.coefficient_name(dim, dim_to_coordinate_name)
coef_template = LinearCoef.coef_template_str(param_name, coefficient_name)
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 ec41adb2..068f7154 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -2,7 +2,8 @@ import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
-from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \
+ ALTITUDINAL_GEV_MODELS_LOCATION
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic, \
plot_individual_aic
@@ -33,7 +34,7 @@ def plot_moments(studies, massif_names=None):
def plot_altitudinal_fit(studies, massif_names=None):
- model_classes = ALTITUDINAL_GEV_MODELS
+ model_classes = ALTITUDINAL_GEV_MODELS_LOCATION
# model_classes = ALTITUDINAL_GEV_MODELS_QUADRATIC
visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
model_classes=model_classes,
@@ -42,7 +43,7 @@ def plot_altitudinal_fit(studies, massif_names=None):
# 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)
+ plot_total_aic(model_classes, visualizer)
def main():
@@ -50,12 +51,15 @@ def main():
# 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 = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900][:]
+ altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200][3:]
+ altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
+ print(altitudes)
+ altitudes = [900, 1200, 1500][:]
study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
SafranPrecipitation7Days][:]
- study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day,
- SafranSnowfall3Days, SafranPrecipitation3Days][:1]
+ study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranPrecipitation1Day
+ , SafranPrecipitation3Days][:2]
# seasons = [Season.automn, Season.winter, Season.spring][::-1]
seasons = [Season.winter]
# seasons = [Season.winter_extended]
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 351f9c7c..b01b505b 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
@@ -32,7 +32,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self.non_stationary_models = model_classes
self.fit_method = fit_method
self.display_only_model_that_pass_anderson_test = display_only_model_that_pass_anderson_test
- self.massif_names = massif_names if massif_names is not None else self.study.study_massif_names
+ self.massif_names = massif_names if massif_names is not None else self.study.all_massif_names()
self.massif_name_to_massif_id = {m: i for i, m in enumerate(self.massif_names)}
# Load one fold fit
self._massif_name_to_one_fold_fit = {}
@@ -143,7 +143,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
pass
def plot_year_for_the_peak(self):
- t_list = 1900 + np.arange(200)
+ t_list = 1800 + np.arange(400)
for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
ax = plt.gca()
# One plot for each altitude
@@ -162,7 +162,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
ax.legend()
ax.set_xlabel('Year')
ax.set_xlabel('Mean annual maxima of ')
- self.plot_name = 'Peak year for {}'.format(massif_name)
- self.show_or_save_to_file(no_title=True)
+ plot_name = '{}/{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, 'Peak year', massif_name.replace('_', ''))
+ self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
plt.close()
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 56cf37f9..898f825f 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
@@ -116,7 +116,8 @@ class OneFoldFit(object):
@property
def best_shape(self):
- # We take any altitude (altitude=1000 for instance) as the shape is constant w.r.t the altitude
+ # We take any altitude (altitude=1000 for instance) and any year
+ # as the shape is constant w.r.t the altitude and the year
return self.get_gev_params(altitude=1000, year=2019).shape
def best_coef(self, param_name, dim, degree):
diff --git a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
index 62055821..c77b8374 100644
--- a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
+++ b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
@@ -74,7 +74,7 @@ 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_massifs(altitude_to_visualizer)
plot_uncertainty_histogram(altitude_to_visualizer)
# plot_selection_curves(altitude_to_visualizer)
# plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)
diff --git a/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py b/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py
index 923e9211..88272fb2 100644
--- a/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py
+++ b/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py
@@ -1,4 +1,4 @@
-from typing import Dict
+from typing import Dict, Tuple
import matplotlib.pyplot as plt
import numpy as np
@@ -39,6 +39,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
fontsize_label = 15
legend_size = 15
# Plot histogram
+ ylim = (0, 110)
for j, ci_method in enumerate(visualizer.uncertainty_methods):
if len(visualizer.uncertainty_methods) == 2:
offset = -50 if j == 0 else 50
@@ -51,7 +52,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width, legend_size)
# Plot percentages of return level excess on the right axis
- plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax.twinx(), fontsize_label, legend_size)
+ ylim = plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax.twinx(), fontsize_label, legend_size)
ax.set_xticks(altitudes)
@@ -63,7 +64,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
# '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, 110])
+ ax.set_ylim(ylim)
ax.yaxis.grid()
ax_twiny = ax.twiny()
@@ -76,7 +77,8 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
ax_twiny.set_xticklabels(nb_massif_names)
ax_twiny.set_xlabel('Number of massifs at each altitude (for the percentage and the mean)', fontsize=fontsize_label)
- ax.set_yticks([20 * i for i in range(6)])
+ nb_ticks = 1 + (ylim[1] - ylim[0]) // 20
+ ax.set_yticks([100 - 20 * i for i in range(nb_ticks)][::-1])
visualizer.plot_name = 'Percentages of exceedance with {}'.format(
get_display_name_from_object_type(model_subset_for_uncertainty))
visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure, tight_layout=True)
@@ -109,7 +111,7 @@ def plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_metho
ax.legend(loc='upper left', prop={'size': legend_size})
-def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax, fontsize_label, legend_size=10):
+def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax, fontsize_label, legend_size=10) -> Tuple[int, int]:
l = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[2] for v in visualizers]
lower_bound, mean, upper_bound = list(zip(*l))
other_mean = [e[1] for e in l]
@@ -121,7 +123,6 @@ def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_meth
'50-year return levels and French standards (\%)'
label_name = 'Mean relative difference'
- print('mean relative difference', mean)
ax.plot(altitudes, mean, linestyle='--', marker='o', color=color,
label=label_name)
@@ -132,7 +133,10 @@ def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_meth
ax.tick_params(labelsize=fontsize_label)
ax.legend(loc='upper right', prop={'size': legend_size})
ax.set_ylabel(full_label_name, fontsize=fontsize_label)
- ax.set_ylim([0, 110])
+
+ ylim = (-85, 110)
+ ax.set_ylim(ylim)
+ return ylim
def get_label_confidence_interval(label_name):
--
GitLab