diff --git a/extreme_data/meteo_france_data/scm_models_data/utils.py b/extreme_data/meteo_france_data/scm_models_data/utils.py
index dec1d7ea6ceb5cb37cb9311f14329973493a7ef1..0502267efe9fb99f7d435f740cd2c100b2e97b41 100644
--- a/extreme_data/meteo_france_data/scm_models_data/utils.py
+++ b/extreme_data/meteo_france_data/scm_models_data/utils.py
@@ -22,6 +22,7 @@ class Season(Enum):
winter = 2
spring = 3
automn = 4
+ summer = 5
def season_to_str(season):
@@ -35,14 +36,18 @@ class FrenchRegion(Enum):
pyrenees = 2
+season_to_start_day_and_last_day = {
+ Season.annual: ('08-01', '07-31'),
+ Season.winter_extended: ('11-01', '05-31'),
+ Season.winter: ('12-01', '02-28'),
+ Season.spring: ('03-01', '05-31'),
+ Season.summer: ('06-01', '08-31'),
+ Season.automn: ('09-01', '11-30'),
+
+}
+
+
def first_day_and_last_day(season):
- season_to_start_day_and_last_day = {
- Season.annual: ('08-01', '07-31'),
- Season.winter_extended: ('11-01', '05-31'),
- Season.winter: ('12-01', '02-28'),
- Season.spring: ('03-01', '05-31'),
- Season.automn: ('09-01', '11-30'),
- }
return season_to_start_day_and_last_day[season]
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
index 20885487073e75fd0bd448b9d324a40924fd7d6c..73881938a0dbdd1543ab9eff38f7db5567fd0b7a 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
@@ -716,14 +716,16 @@ class StudyVisualizer(VisualizationParameters):
# PLot functions that should be common
def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label=True,
- negative_and_positive_values=True, massif_name_to_text=None):
+ negative_and_positive_values=True, massif_name_to_text=None, altitude=None):
+ if altitude is None:
+ altitude = self.study.altitude
if len(massif_name_to_value) > 0:
- load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method),
+ load_plot(cmap, graduation, label, massif_name_to_value, altitude, fitmethod_to_str(fit_method),
add_x_label=add_x_label, negative_and_positive_values=negative_and_positive_values,
massif_name_to_text=massif_name_to_text)
self.plot_name = plot_name
# self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True, dpi=500)
- self.show_or_save_to_file(add_classic_title=False, no_title=True, dpi=500)
+ self.show_or_save_to_file(add_classic_title=False, no_title=True, dpi=500, tight_layout=True)
plt.close()
def plot_abstract(self, massif_name_to_value, label, plot_name, fit_method='', graduation=10.0, cmap=plt.cm.bwr,
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 0afdd4dd8222daebf955a3bf37bfae26cf21a25d..29e3646923eb4df6ef32951dd0b8b73b987dce91 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
@@ -97,20 +97,20 @@ ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS = [
StationaryAltitudinal,
AltitudinalShapeConstantTimeLocationLinear,
AltitudinalShapeConstantTimeScaleLinear,
- AltitudinalShapeConstantTimeShapeLinear,
- AltitudinalShapeConstantTimeLocShapeLinear,
+ # AltitudinalShapeConstantTimeShapeLinear,
+ # AltitudinalShapeConstantTimeLocShapeLinear,
AltitudinalShapeConstantTimeLocScaleLinear,
- AltitudinalShapeConstantTimeScaleShapeLinear,
- AltitudinalShapeConstantTimeLocScaleShapeLinear,
+ # AltitudinalShapeConstantTimeScaleShapeLinear,
+ # AltitudinalShapeConstantTimeLocScaleShapeLinear,
# With a linear shape for the altitude (8 models)
AltitudinalShapeLinearTimeStationary,
AltitudinalShapeLinearTimeLocationLinear,
AltitudinalShapeLinearTimeScaleLinear,
- AltitudinalShapeLinearTimeShapeLinear,
- AltitudinalShapeLinearTimeLocShapeLinear,
+ # AltitudinalShapeLinearTimeShapeLinear,
+ # AltitudinalShapeLinearTimeLocShapeLinear,
AltitudinalShapeLinearTimeLocScaleLinear,
- AltitudinalShapeLinearTimeScaleShapeLinear,
- AltitudinalShapeLinearTimeLocScaleShapeLinear,
+ # AltitudinalShapeLinearTimeScaleShapeLinear,
+ # AltitudinalShapeLinearTimeLocScaleShapeLinear,
]
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 ee31f00096ef1cc0384e8ddbac0a4e285c504ccd..9ea4741ba097a80ca63e0e7d168a9beef92dd217 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -7,6 +7,7 @@ from extreme_data.meteo_france_data.adamont_data.snowfall_simulation import Safr
from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \
ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM, ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES, \
ALTITUDINAL_GEV_MODELS_LOCATION, ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic, \
plot_individual_aic
from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import MeanAlpsTemperatureCovariate
@@ -74,15 +75,20 @@ def main():
SafranPrecipitation7Days][:]
study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranPrecipitation1Day
, SafranPrecipitation3Days][:1]
-
altitudes = [1800, 2100, 2400]
- study_classes = [SafranSnowfall1Day][:]
+ study_classes = [SafranSnowfall1Day, SafranSnowfall3Days][:1]
# Common parameters
# altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]
massif_names = None
- seasons = [Season.annual]
+ seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
+
+ all_groups = altitudes_for_groups[:]
+ for altitudes in all_groups:
+ main_loop(altitudes, massif_names, seasons, study_classes)
+
+def main_loop(altitudes, massif_names, seasons, study_classes):
for season in seasons:
for study_class in study_classes:
if issubclass(study_class, SimulationStudy):
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
index 804abddfac77f7926a775d0b6b9bd271c9e2c542..07e569441e6f95a78674818f9d74a8cfb60ad9ea 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
@@ -1,28 +1,79 @@
from enum import Enum
+# The order is important
+altitudes_for_groups = [
+ [900, 1200, 1500],
+ [1800, 2100, 2400],
+ [2700, 3000, 3300]
+]
+# altitudes_for_groups = [
+# [600, 900, 1200, 1500, 1800],
+# [1500, 1800, 2100, 2400, 2700],
+# [2400, 2700, 3000, 3300, 3600]
+# ]
-class AltitudeGroup(Enum):
- low = 0
- mid = 1
- high = 2
- unspecfied = 3
+class AbstractAltitudeGroup(object):
-def altitude_group_to_reference_altitude():
- return {
- AltitudeGroup.low: 1000,
- AltitudeGroup.mid: 2000,
- AltitudeGroup.high: 3000,
- AltitudeGroup.unspecfied: 1000,
- }
+ @property
+ def name(self):
+ raise NotImplementedError
+
+ @property
+ def reference_altitude(self):
+ raise NotImplementedError
+
+
+class LowAltitudeGroup(AbstractAltitudeGroup):
+
+ @property
+ def name(self):
+ return 'low'
+
+ @property
+ def reference_altitude(self):
+ return 1000
+
+
+class MidAltitudeGroup(AbstractAltitudeGroup):
+
+ @property
+ def name(self):
+ return 'mid'
+
+ @property
+ def reference_altitude(self):
+ return 2000
+
+
+class HighAltitudeGroup(AbstractAltitudeGroup):
+
+ @property
+ def name(self):
+ return 'high'
+
+ @property
+ def reference_altitude(self):
+ return 3000
+
+
+class DefaultAltitudeGroup(AbstractAltitudeGroup):
+
+ @property
+ def name(self):
+ return 'default'
+
+ @property
+ def reference_altitude(self):
+ return 500
def get_altitude_group_from_altitudes(altitudes):
s = set(altitudes)
- if s == {900, 1200, 1500}:
- return AltitudeGroup.low
- elif s == {1800, 2100, 2400}:
- return AltitudeGroup.mid
- elif s == {2700, 3000, 3300}:
- return AltitudeGroup.high
+ if s == set(altitudes_for_groups[0]):
+ return LowAltitudeGroup()
+ elif s == set(altitudes_for_groups[1]):
+ return MidAltitudeGroup()
+ elif s == set(altitudes_for_groups[2]):
+ return HighAltitudeGroup()
else:
- return AltitudeGroup.unspecfied
+ return DefaultAltitudeGroup()
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 3edeca7bea2b9f8d7e644c4066770d48b5920a43..4240af2bd44c51d816fc6966f2361b6c6f6b1625 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,10 +1,14 @@
from typing import List, Dict
+
+import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from cached_property import cached_property
from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import get_shifted_map, \
+ get_colors, ticks_values_and_labels_for_percentages, get_half_colormap, ticks_values_and_labels_for_positive_value
from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
SCM_STUDY_CLASS_TO_ABBREVIATION, ALL_ALTITUDES_WITHOUT_NAN
from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import plot_against_altitude
@@ -16,6 +20,8 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
AbstractSpatioTemporalPolynomialModel
from extreme_fit.model.margin_model.utils import MarginFitMethod
from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import \
+ get_altitude_group_from_altitudes
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import \
OneFoldFit
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -39,27 +45,95 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
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.all_massif_names()
self.massif_name_to_massif_id = {m: i for i, m in enumerate(self.massif_names)}
+ self.altitude_group = get_altitude_group_from_altitudes(self.studies.altitudes)
# Load one fold fit
self._massif_name_to_one_fold_fit = {}
for massif_name in self.massif_names:
- assert top_n_values_to_remove is None
- dataset = studies.spatio_temporal_dataset(massif_name=massif_name,
- top_n_values_to_remove=top_n_values_to_remove)
- old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method,
- self.temporal_covariate_for_fit)
- self._massif_name_to_one_fold_fit[massif_name] = old_fold_fit
+ if any([massif_name in study.study_massif_names for study in self.studies.altitude_to_study.values()]):
+ assert top_n_values_to_remove is None
+ dataset = studies.spatio_temporal_dataset(massif_name=massif_name,
+ top_n_values_to_remove=top_n_values_to_remove)
+ old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method,
+ self.temporal_covariate_for_fit,
+ type(self.altitude_group),
+ self.display_only_model_that_pass_anderson_test)
+ self._massif_name_to_one_fold_fit[massif_name] = old_fold_fit
@property
def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
return {massif_name: old_fold_fit for massif_name, old_fold_fit in self._massif_name_to_one_fold_fit.items()
- if not self.display_only_model_that_pass_anderson_test or old_fold_fit.goodness_of_fit_anderson_test}
+ if not self.display_only_model_that_pass_anderson_test
+ or old_fold_fit.has_at_least_one_valid_non_stationary_model}
def plot_moments(self):
for method_name in ['moment', 'changes_in_the_moment', 'relative_changes_in_the_moment']:
for order in [1, 2, None]:
- self.plot_general(method_name, order)
+ # self.plot_against_years(method_name, order)
+ self.plot_map_moment(method_name, order)
+
+ def plot_map_moment(self, method_name, order):
+ # Compute values
+ massif_name_to_value = {}
+ for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
+ value = one_fold_fit.__getattribute__(method_name)([self.altitude_group.reference_altitude], order=order)[0]
+ massif_name_to_value[massif_name] = value
+
+ # Common plot settings
+ moment = ' '.join(method_name.split('_'))
+ moment = moment.replace('moment', '{} in 2019'.format(OneFoldFit.get_moment_str(order=order)))
+ plot_name = '{}{} annual maxima of {}'.format(OneFoldFit.folder_for_plots, moment,
+ SCM_STUDY_CLASS_TO_ABBREVIATION[
+ self.studies.study_class])
+ ylabel = '{} ({})'.format(plot_name, self.study.variable_unit)
+
+ # Plot the map
+ if any([np.isinf(v) for v in massif_name_to_value.values()]):
+ print("shape to large > 0.5, thus removing std that are infinite")
+ massif_name_to_value = {m: v for m, v in massif_name_to_value.items()
+ if not np.isinf(v)}
+
+ print(massif_name_to_value)
+ negative_and_positive_values = min(massif_name_to_value.values()) < 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=negative_and_positive_values,
+ massif_name_to_text=None, altitude=self.altitude_group.reference_altitude)
- def plot_general(self, method_name, order):
+
+
+ # ax.get_xaxis().set_visible(True)
+ # ax.set_xticks([])
+ # ax.set_xlabel('Altitude = {}m'.format(self.altitude_group.reference_altitude), fontsize=15)
+
+
+ # cmap = get_shifted_map(min_ratio, max_ratio)
+ # print(massif_name_to_value)
+ # massif_name_to_color = {m: get_colors([v], cmap, -max_abs_change, max_abs_change)[0]
+ # for m, v in massif_name_to_value.items()}
+ #
+ #
+ # ticks_values_and_labels = ticks_values_and_labels_for_percentages(graduation, max_abs_change)
+ #
+ # ax = self.study.visualize_study(massif_name_to_value=massif_name_to_value,
+ # replace_blue_by_white=False,
+ # axis_off=False, show_label=False,
+ # add_colorbar=add_colorbar,
+ # # massif_name_to_marker_style=self.massif_name_to_marker_style,
+ # # marker_style_to_label_name=self.selected_marker_style_to_label_name,
+ # massif_name_to_color=massif_name_to_color,
+ # cmap=cmap,
+ # show=False,
+ # ticks_values_and_labels=ticks_values_and_labels,
+ # label=ylabel,
+ # add_legend=False,
+ # )
+
+ # self.plot_name = plot_name
+ # self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
+ # dpi=500)
+ # ax.clear()
+
+ def plot_against_years(self, method_name, order):
ax = plt.gca()
min_altitude, *_, max_altitude = self.studies.altitudes
altitudes_plot = np.linspace(min_altitude, max_altitude, num=50)
@@ -74,20 +148,17 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
ax.legend(prop={'size': 7}, ncol=3)
moment = ' '.join(method_name.split('_'))
moment = moment.replace('moment', '{} in 2019'.format(OneFoldFit.get_moment_str(order=order)))
- plot_name = '{}/Model {} annual maxima of {}'.format(OneFoldFit.folder_for_plots, moment,
+ plot_name = '{}Model {} annual maxima of {}'.format(OneFoldFit.folder_for_plots, moment,
SCM_STUDY_CLASS_TO_ABBREVIATION[self.studies.study_class])
ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
ax.set_xlabel('altitudes', fontsize=15)
- # lim_down, lim_up = ax.get_ylim()
- # lim_up += (lim_up - lim_down) / 3
- # ax.set_ylim([lim_down, lim_up])
ax.tick_params(axis='both', which='major', labelsize=13)
- self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
+ self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True)
ax.clear()
def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm, add_x_label=True,
negative_and_positive_values=True, massif_name_to_text=None):
- plot_name = '{}/{}'.format(OneFoldFit.folder_for_plots, label)
+ plot_name = '{}{}'.format(OneFoldFit.folder_for_plots, label)
self.plot_map(cmap, self.fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label,
negative_and_positive_values,
massif_name_to_text)
@@ -133,7 +204,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
print(coef_name, graduation, max(values), min(values))
negative_and_positive_values = (max(values) > 0) and (min(values) < 0)
self.plot_abstract_fast(massif_name_to_best_coef,
- label='{}/Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots,
+ label='{}Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots,
coef_name,
SCM_STUDY_CLASS_TO_ABBREVIATION[
type(self.study)],
@@ -143,10 +214,13 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
def plot_shape_map(self):
self.plot_abstract_fast(self.massif_name_to_shape,
- label='Shape parameter for {} maxima of {}'.format(self.study.season_name,
- SCM_STUDY_CLASS_TO_ABBREVIATION[
- type(self.study)]),
- add_x_label=False, graduation=0.1, massif_name_to_text=self.massif_name_to_name)
+ label='Shape parameter for {} maxima of {} in 2019 at {}m'.format(
+ self.study.season_name,
+ SCM_STUDY_CLASS_TO_ABBREVIATION[
+ type(self.study)],
+ self.altitude_group.reference_altitude),
+ add_x_label=False, graduation=0.1, massif_name_to_text=self.massif_name_to_name,
+ cmap=matplotlib.cm.get_cmap('BrBG_r'))
def plot_altitude_for_the_peak(self):
pass
@@ -189,7 +263,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
ax.set_ylabel('{} of {} in {} ({})'.format(ylabel, SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)],
massif_name.replace('_', ' '), self.study.variable_unit))
peak_year_folder = 'Peak year ' + ylabel
- plot_name = '{}/{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, peak_year_folder,
+ plot_name = '{}{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, peak_year_folder,
massif_name.replace('_', ''))
self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True, tight_layout=True)
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 62704f80166f2537db7df9453297128b84458e57..a91a23f3cbe683bcbcdc16b1602d8194abcc479c 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
@@ -15,7 +15,8 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_m
from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pariwise_analysis.gev_with_linear_shape_wrt_altitude import \
AltitudinalShapeLinearTimeStationary
from extreme_fit.model.margin_model.utils import MarginFitMethod
-from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import AltitudeGroup
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import AbstractAltitudeGroup, \
+ DefaultAltitudeGroup
from root_utils import classproperty
from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
@@ -26,8 +27,10 @@ class OneFoldFit(object):
def __init__(self, massif_name: str, dataset: AbstractDataset, models_classes,
fit_method=MarginFitMethod.extremes_fevd_mle, temporal_covariate_for_fit=None,
- altitude_group=AltitudeGroup.unspecfied):
- self.altitude_group = altitude_group
+ altitude_class=DefaultAltitudeGroup,
+ only_models_that_pass_anderson_test=True):
+ self.only_models_that_pass_anderson_test = only_models_that_pass_anderson_test
+ self.altitude_group = altitude_class()
self.massif_name = massif_name
self.dataset = dataset
self.models_classes = models_classes
@@ -49,7 +52,7 @@ class OneFoldFit(object):
@classproperty
def folder_for_plots(cls):
- return 'Total aic' if cls.best_estimator_minimizes_total_aic else 'Individual aic'
+ return 'Total aic/' if cls.best_estimator_minimizes_total_aic else ''
@classmethod
def get_moment_str(cls, order):
@@ -67,7 +70,7 @@ class OneFoldFit(object):
elif order == 2:
return gev_params.std
elif order is None:
- return gev_params.return_level(return_period=2019)
+ return gev_params.return_level(return_period=50)
else:
raise NotImplementedError
@@ -106,9 +109,20 @@ class OneFoldFit(object):
return sorted_estimators
@cached_property
- def sorted_estimators_without_stationary(self):
+ def _sorted_estimators_without_stationary(self):
return [e for e in self.sorted_estimators if not isinstance(e.margin_model, StationaryAltitudinal)]
+ @cached_property
+ def sorted_estimators_without_stationary(self):
+ if self.only_models_that_pass_anderson_test:
+ return [e for e in self._sorted_estimators_without_stationary if self.goodness_of_fit_test(e)]
+ else:
+ return self._sorted_estimators_without_stationary
+
+ @property
+ def has_at_least_one_valid_non_stationary_model(self):
+ return len(self.sorted_estimators_without_stationary) > 0
+
@property
def model_class_to_estimator_with_finite_aic(self):
return {type(estimator.margin_model): estimator for estimator in self.sorted_estimators}
@@ -118,8 +132,11 @@ class OneFoldFit(object):
if self.best_estimator_minimizes_total_aic and self.best_estimator_class_for_total_aic is not None:
return self.model_class_to_estimator[self.best_estimator_class_for_total_aic]
else:
- best_estimator = self.sorted_estimators_without_stationary[0]
- return best_estimator
+ if self.has_at_least_one_valid_non_stationary_model:
+ best_estimator = self.sorted_estimators_without_stationary[0]
+ return best_estimator
+ else:
+ raise ValueError('This should not happen')
@property
def best_margin_model(self):
@@ -131,10 +148,11 @@ class OneFoldFit(object):
@property
def best_shape(self):
- # 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=self.altitude_plot, year=2019).shape
- return self.get_gev_params(altitude=1000, year=2019).shape
+ @property
+ def altitude_plot(self):
+ return self.altitude_group.reference_altitude
def best_coef(self, param_name, dim, degree):
try:
@@ -181,19 +199,24 @@ class OneFoldFit(object):
else:
return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=self.degree_freedom_chi2)
- @property
- def goodness_of_fit_anderson_test(self):
- if self.folder_for_plots in self._folder_to_goodness:
- return self._folder_to_goodness[self.folder_for_plots]
- else:
- quantiles = self.compute_empirical_quantiles(estimator=self.best_estimator)
- goodness_of_fit_anderson_test = goodness_of_fit_anderson(quantiles, self.SIGNIFICANCE_LEVEL)
- if not goodness_of_fit_anderson_test:
- print('{} with {} does not pass the anderson test for model {}'.format(self.massif_name,
- self.folder_for_plots,
- type(self.best_margin_model)))
- self._folder_to_goodness[self.folder_for_plots] = goodness_of_fit_anderson_test
- return goodness_of_fit_anderson_test
+ # @property
+ # def goodness_of_fit_anderson_test(self):
+ # if self.folder_for_plots in self._folder_to_goodness:
+ # return self._folder_to_goodness[self.folder_for_plots]
+ # else:
+ # estimator = self.best_estimator
+ # goodness_of_fit_anderson_test = self.goodness_of_fit_test(estimator)
+ # if not goodness_of_fit_anderson_test:
+ # print('{} with {} does not pass the anderson test for model {}'.format(self.massif_name,
+ # self.folder_for_plots,
+ # type(self.best_margin_model)))
+ # self._folder_to_goodness[self.folder_for_plots] = goodness_of_fit_anderson_test
+ # return goodness_of_fit_anderson_test
+
+ def goodness_of_fit_test(self, estimator):
+ quantiles = self.compute_empirical_quantiles(estimator=estimator)
+ goodness_of_fit_anderson_test = goodness_of_fit_anderson(quantiles, self.SIGNIFICANCE_LEVEL)
+ return goodness_of_fit_anderson_test
def compute_empirical_quantiles(self, estimator):
empirical_quantiles = []
diff --git a/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_analysis_cluster_level.py b/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_analysis_cluster_level.py
index 5a7372fa88cb9ed1fb77396a044973b2a66977ef..63457fb441b536dd809905057688d35a2bbbd521 100644
--- a/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_analysis_cluster_level.py
+++ b/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_analysis_cluster_level.py
@@ -91,8 +91,6 @@ class PointWIseGevAnalysisForCluster(AltitudesStudies):
for cluster_id, annual_maxima_list in d.items():
a = np.array(annual_maxima_list)
a = a.flatten()
- print("here")
- print(a.shape)
cluster_id_to_time_annual_maxima_list[cluster_id].append(list(a))
return cluster_id_to_time_annual_maxima_list