diff --git a/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
new file mode 100644
index 0000000000000000000000000000000000000000..c911fa5e7c3ed9aa23787afbd21b900681e84906
--- /dev/null
+++ b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
@@ -0,0 +1,88 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.lines import Line2D
+
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_to_color
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, rcp_scenarios, get_gcm_list, \
+ get_linestyle_from_scenario, scenario_to_str, adamont_scenarios_real
+from extreme_data.meteo_france_data.adamont_data.cmip5.climate_explorer_cimp5 import years_and_global_mean_temps
+
+
+def temperature_minmax_to_year_minmax(gcm, scenario, temperature_min, temperature_max):
+ years, global_mean_temps = years_and_global_mean_temps(gcm, scenario, year_min=2005, year_max=2100,
+ rolling=30, anomaly=True)
+ years, global_mean_temps = np.array(years), np.array(global_mean_temps)
+ ind = temperature_min < global_mean_temps
+ ind &= global_mean_temps < temperature_max
+ years_to_select = years[ind]
+ ind2 = years_to_select[:-1] == years_to_select[1:] - 1
+ if not all(ind2):
+ i = list(ind2).index(False)
+ years_to_select = years_to_select[:i + 1]
+ # A minimum of 30 years of data is needed to find a trend
+ if len(years_to_select) >= 30:
+ year_min, year_max = years_to_select[0], years_to_select[-1]
+ assert (year_max - year_min + 1) == len(years_to_select)
+ return year_min, year_max
+ else:
+ return None, None
+
+
+def get_nb_data(gcm, scenario, temperature_min, temperature_max):
+ year_min, year_max = temperature_minmax_to_year_minmax(gcm, scenario, temperature_min, temperature_max)
+ if year_min is None:
+ return 0
+ else:
+ return year_max - year_min + 1
+
+
+def plot_nb_data_one_line(ax, gcm, scenario, temp_min, temp_max):
+
+ nb_data = [get_nb_data(gcm, scenario, mi, ma) for mi, ma in zip(temp_min, temp_max)]
+ color = gcm_to_color[gcm]
+ linestyle = get_linestyle_from_scenario(scenario)
+
+ # Filter out the zero value
+ nb_data, temp_min = np.array(nb_data), np.array(temp_min)
+ ind = np.array(nb_data) > 0
+ nb_data, temp_min = nb_data[ind], temp_min[ind]
+
+ # For the legend
+ if scenario is AdamontScenario.rcp26:
+ ax.plot(temp_min[0], nb_data[0], color=color, linestyle='solid', label=gcm)
+
+ ax.plot(temp_min, nb_data, linestyle=linestyle, color=color)
+
+
+def plot_nb_data():
+ temp_max, temp_min = get_temp_min_and_temp_max()
+
+ ax = plt.gca()
+ for gcm in get_gcm_list(adamont_version=2)[:]:
+ for scenario in rcp_scenarios[:]:
+ plot_nb_data_one_line(ax, gcm, scenario, temp_min, temp_max)
+
+ ax.legend()
+ ticks_labels = ['{}-{}'.format(mi, ma) for mi, ma in zip(temp_min, temp_max)]
+ ax.set_xticks(temp_min)
+ ax.set_xticklabels(ticks_labels)
+ ax.set_xlabel('Temperature interval')
+ ax.set_ylabel('Nb of Data')
+ ax2 = ax.twinx()
+ legend_elements = [
+ Line2D([0], [0], color='k', lw=1, label=scenario_to_str(s),
+ linestyle=get_linestyle_from_scenario(s)) for s in adamont_scenarios_real
+ ]
+ ax2.legend(handles=legend_elements, loc='upper center')
+ ax2.set_yticks([])
+ plt.show()
+
+
+def get_temp_min_and_temp_max():
+ temp_min = np.arange(0.5, 2.5, 0.5)
+ temp_max = temp_min + 1.5
+ return temp_min, temp_max
+
+
+if __name__ == '__main__':
+ plot_nb_data()
diff --git a/extreme_fit/distribution/abstract_extreme_params.py b/extreme_fit/distribution/abstract_extreme_params.py
index 8c98f1b99b78cba3071fd99680a442601a544116..e892342d09b545849dac25f5ac3780c41dae3558 100644
--- a/extreme_fit/distribution/abstract_extreme_params.py
+++ b/extreme_fit/distribution/abstract_extreme_params.py
@@ -1,7 +1,10 @@
+import numpy as np
+
from extreme_fit.distribution.abstract_params import AbstractParams
class AbstractExtremeParams(AbstractParams):
+ SMALL_SCALE_PARAMETERS_ARE_UNDEFINED = True
def __init__(self, loc: float, scale: float, shape: float):
self.location = loc
@@ -11,4 +14,6 @@ class AbstractExtremeParams(AbstractParams):
# (sometimes it happens, when we want to find a quantile for every point of a 2D map
# then it can happen that a corner point that was not used for fitting correspond to a negative scale,
# in the case we set all the parameters as equal to np.nan, and we will not display those points)
- self.has_undefined_parameters = self.scale <= 0
+ self.has_undefined_parameters = (self.scale <= 0)
+ if self.SMALL_SCALE_PARAMETERS_ARE_UNDEFINED:
+ self.has_undefined_parameters = self.has_undefined_parameters or (np.isclose(self.scale, 0))
\ No newline at end of file
diff --git a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
index f8afbdcda4071b873fd009c8e8c7c24a395364f6..6c71b8c53a0bef12871d86bf9c9da2cb31e8de24 100644
--- a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
+++ b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
@@ -61,11 +61,13 @@ class LinearMarginEstimator(AbstractMarginEstimator):
def function_from_fit(self) -> LinearMarginFunction:
return load_margin_function(self, self.margin_model)
+ def coordinates_for_nllh(self, split=Split.all):
+ return pd.concat([self.df_coordinates_spat(split=split), self.df_coordinates_temp(split=split)], axis=1).values
+
def nllh(self, split=Split.all):
nllh = 0
maxima_values = self.dataset.maxima_gev(split=split)
- df = pd.concat([self.df_coordinates_spat(split=split), self.df_coordinates_temp(split=split)], axis=1)
- coordinate_values = df.values
+ coordinate_values = self.coordinates_for_nllh(split=split)
for maximum, coordinate in zip(maxima_values, coordinate_values):
assert len(maximum) == 1, \
'So far, only one observation for each coordinate, but code would be easy to change'
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 30d0eb09f3f3cdd1ff73115b26ba87b1090ae918..5a3f8a16d28e890697027c540b03544c7dde0004 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
@@ -17,6 +17,10 @@ class AbstractAltitudeGroup(object):
def name(self):
raise NotImplementedError
+ @property
+ def altitudes(self):
+ return altitudes_for_groups[self.group_id-1]
+
@property
def reference_altitude(self):
raise NotImplementedError
@@ -141,6 +145,10 @@ class DefaultAltitudeGroup(AbstractAltitudeGroup):
return 500
+def get_altitude_class_from_altitudes(altitudes):
+ return type(get_altitude_group_from_altitudes(altitudes))
+
+
def get_altitude_group_from_altitudes(altitudes):
s = set(altitudes)
if s == set(altitudes_for_groups[0]):
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 65f900e2eb466fd98d9120bb0b21288ff724ad30..a0451c44757d65d4c15c898937e433a5cbc436ed 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
@@ -180,12 +180,31 @@ class OneFoldFit(object):
# Remove wrong shape
estimators = [e for e in estimators if -0.5 < self._compute_shape_for_reference_altitude(e) < 0.5]
# Remove models with undefined parameters for the coordinate of interest
- coordinate = np.array([self.altitude_group.reference_altitude, self.last_year])
- estimators = [e for e in estimators if not e.function_from_fit.get_params(coordinate).has_undefined_parameters]
+ well_defined_estimators = []
+ for e in estimators:
+ coordinate_values_for_the_fit = e.coordinates_for_nllh(Split.all)
+ coordinate_values_for_the_result = [np.array([self.altitude_group.reference_altitude, c])
+ for c in self._covariate_before_and_after]
+ coordinate_values_to_check = list(coordinate_values_for_the_fit) + coordinate_values_for_the_result
+ has_undefined_parameters = False
+ for coordinate in coordinate_values_to_check:
+ gev_params = e.function_from_fit.get_params(coordinate)
+ if gev_params.has_undefined_parameters:
+ has_undefined_parameters = True
+ break
+ if not has_undefined_parameters:
+ well_defined_estimators.append(e)
+ estimators = well_defined_estimators
+
if len(estimators) == 0:
print(self.massif_name, " has only implausible models")
- sorted_estimators = sorted([estimator for estimator in estimators], key=lambda e: e.aic())
+ try:
+ sorted_estimators = sorted([estimator for estimator in estimators], key=lambda e: e.aic())
+ except AssertionError as e:
+ print('Error for')
+ print(self.massif_name, self.altitude_group)
+ raise
return sorted_estimators
def _compute_shape_for_reference_altitude(self, estimator):
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
index e0a556db286c0b1c90d3833323dfebd10e1a2081..484634526ba5032595c4e7a83636ba9c63c5e2bd 100644
--- a/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
+++ b/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
@@ -34,6 +34,9 @@ def compute_and_assign_max_abs(visualizer_list):
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])
+ if max_abs == 0:
+ epsilon = 0.1
+ max_abs = epsilon
method_name_and_order_to_max_abs[c] = max_abs
# Assign the max abs dictionary
for v in visualizer_list:
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py
index 268b38fc16a94ee5972e6f9b8af50a959e93a4b3..23f51598fbd3da1818e43859218098e9b8395dcb 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py
@@ -4,13 +4,16 @@ from typing import List
import matplotlib as mpl
+from projects.projected_snowfall.elevation_temporal_model_for_projections.visualizer_for_sensitivity import \
+ VisualizerForSensivity
+
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
from projects.projected_snowfall.elevation_temporal_model_for_projections.visualizer_for_projection_ensemble import \
- MetaVisualizerForProjectionEnsemble
+ VisualizerForProjectionEnsemble
import matplotlib
from extreme_fit.model.utils import set_seed_for_test
@@ -34,62 +37,72 @@ from extreme_data.meteo_france_data.scm_models_data.utils import Season
def main():
start = time.time()
- study_classes = [AdamontSnowfall][:1]
+ study_class = AdamontSnowfall
ensemble_fit_class = [IndependentEnsembleFit]
temporal_covariate_for_fit = [TimeTemporalCovariate, AnomalyTemperatureTemporalCovariate][1]
set_seed_for_test()
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
fast = None
- scenarios = rcp_scenarios[:1] if fast is False else [AdamontScenario.rcp26]
+ sensitivity_plot = True
+ scenarios = rcp_scenarios if fast is False else [AdamontScenario.rcp85]
for scenario in scenarios:
gcm_rcm_couples = get_gcm_rcm_couples(scenario)
if fast is None:
massif_names = None
- gcm_rcm_couples = gcm_rcm_couples[:]
+ gcm_rcm_couples = None
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
- altitudes_list = altitudes_for_groups[:1]
+ altitudes_list = altitudes_for_groups[3:]
elif fast:
+ massif_names = ['Vanoise', 'Haute-Maurienne']
+ gcm_rcm_couples = gcm_rcm_couples[4:6]
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
- massif_names = None
- gcm_rcm_couples = [('EC-EARTH', 'RACMO22E')]
- altitudes_list = altitudes_for_groups[1:2]
+ altitudes_list = altitudes_for_groups[2:]
else:
massif_names = None
altitudes_list = altitudes_for_groups[:]
assert isinstance(gcm_rcm_couples, list)
- main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_classes, ensemble_fit_class, scenario,
- temporal_covariate_for_fit)
+ main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_class, ensemble_fit_class, scenario,
+ temporal_covariate_for_fit, sensitivity_plot=sensitivity_plot)
end = time.time()
duration = str(datetime.timedelta(seconds=end - start))
print('Total duration', duration)
-def main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_classes, ensemble_fit_classes, scenario,
- temporal_covariate_for_fit):
+def main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_class, ensemble_fit_classes, scenario,
+ temporal_covariate_for_fit, sensitivity_plot=False):
assert isinstance(altitudes_list, List)
assert isinstance(altitudes_list[0], List)
print('Covariate is {}'.format(temporal_covariate_for_fit))
- for study_class in study_classes:
- print('Inner loop', study_class)
- model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
- assert scenario in rcp_scenarios
- visualizer = MetaVisualizerForProjectionEnsemble(
+ model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+ assert scenario in rcp_scenarios
+ remove_physically_implausible_models = True
+
+ if sensitivity_plot:
+ visualizer = VisualizerForSensivity(
+ altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
+ model_classes=model_classes,
+ ensemble_fit_classes=ensemble_fit_classes,
+ massif_names=massif_names,
+ temporal_covariate_for_fit=temporal_covariate_for_fit,
+ remove_physically_implausible_models=remove_physically_implausible_models,
+ )
+ else:
+ visualizer = VisualizerForProjectionEnsemble(
altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
model_classes=model_classes,
ensemble_fit_classes=ensemble_fit_classes,
massif_names=massif_names,
temporal_covariate_for_fit=temporal_covariate_for_fit,
- remove_physically_implausible_models=True,
+ remove_physically_implausible_models=remove_physically_implausible_models,
+ gcm_to_year_min_and_year_max=None,
)
- visualizer.plot()
- del visualizer
- time.sleep(2)
+ visualizer.plot()
if __name__ == '__main__':
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py
index 9d51602b169ccfac2e44ed2262711dfebdfff9ad..e7d9e38ef1aa0b93d6e85395fa257be7123d788a 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py
@@ -5,7 +5,7 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
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
+ get_altitude_group_from_altitudes, get_altitude_class_from_altitudes
from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
plot_histogram_all_trends_against_altitudes, plot_shoe_plot_changes_against_altitude
from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import compute_and_assign_max_abs
@@ -13,7 +13,7 @@ from projects.projected_snowfall.elevation_temporal_model_for_projections.indepe
IndependentEnsembleFit
-class MetaVisualizerForProjectionEnsemble(object):
+class VisualizerForProjectionEnsemble(object):
def __init__(self, altitudes_list, gcm_rcm_couples, study_class, season, scenario,
model_classes: List[AbstractSpatioTemporalPolynomialModel],
@@ -24,25 +24,39 @@ class MetaVisualizerForProjectionEnsemble(object):
display_only_model_that_pass_gof_test=False,
confidence_interval_based_on_delta_method=False,
remove_physically_implausible_models=False,
+ gcm_to_year_min_and_year_max=None,
):
self.gcm_rcm_couples = gcm_rcm_couples
self.massif_names = massif_names
self.ensemble_fit_classes = ensemble_fit_classes
# Load all studies
- altitude_group_to_gcm_couple_to_studies = {}
+ altitude_class_to_gcm_couple_to_studies = {}
for altitudes in altitudes_list:
- altitude_group = get_altitude_group_from_altitudes(altitudes)
+ altitude_class = get_altitude_class_from_altitudes(altitudes)
gcm_rcm_couple_to_studies = {}
for gcm_rcm_couple in gcm_rcm_couples:
+ if gcm_to_year_min_and_year_max is None:
+ kwargs_study = {}
+ else:
+ gcm = gcm_rcm_couple[0]
+ if gcm not in gcm_to_year_min_and_year_max:
+ # It means that for this gcm and scenario,
+ # there is not enough data (less than 30 years) for the fit
+ continue
+ year_min, year_max = gcm_to_year_min_and_year_max[gcm]
+ kwargs_study = {'year_min': year_min, 'year_max': year_max}
studies = AltitudesStudies(study_class, altitudes, season=season,
- scenario=scenario, gcm_rcm_couple=gcm_rcm_couple)
+ scenario=scenario, gcm_rcm_couple=gcm_rcm_couple,
+ **kwargs_study)
gcm_rcm_couple_to_studies[gcm_rcm_couple] = studies
- altitude_group_to_gcm_couple_to_studies[altitude_group] = gcm_rcm_couple_to_studies
+ if len(gcm_rcm_couple_to_studies) == 0:
+ print('No valid studies for the following couples:', self.gcm_rcm_couples)
+ altitude_class_to_gcm_couple_to_studies[altitude_class] = gcm_rcm_couple_to_studies
# Load ensemble fit
- self.altitude_group_to_ensemble_class_to_ensemble_fit = {}
- for altitude_group, gcm_rcm_couple_to_studies in altitude_group_to_gcm_couple_to_studies.items():
+ self.altitude_class_to_ensemble_class_to_ensemble_fit = {}
+ for altitude_class, gcm_rcm_couple_to_studies in altitude_class_to_gcm_couple_to_studies.items():
ensemble_class_to_ensemble_fit = {}
for ensemble_fit_class in ensemble_fit_classes:
ensemble_fit = ensemble_fit_class(massif_names, gcm_rcm_couple_to_studies, model_classes,
@@ -51,7 +65,7 @@ class MetaVisualizerForProjectionEnsemble(object):
confidence_interval_based_on_delta_method,
remove_physically_implausible_models)
ensemble_class_to_ensemble_fit[ensemble_fit_class] = ensemble_fit
- self.altitude_group_to_ensemble_class_to_ensemble_fit[altitude_group] = ensemble_class_to_ensemble_fit
+ self.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class] = ensemble_class_to_ensemble_fit
def plot(self):
if IndependentEnsembleFit in self.ensemble_fit_classes:
@@ -74,6 +88,8 @@ class MetaVisualizerForProjectionEnsemble(object):
# Aggregated at gcm_rcm_level plots
merge_keys = [IndependentEnsembleFit.Median_merge, IndependentEnsembleFit.Mean_merge]
keys = self.gcm_rcm_couples + merge_keys
+ # Only plot Mean for speed
+ keys = [IndependentEnsembleFit.Mean_merge]
for key in keys:
visualizer_list = [independent_ensemble_fit.gcm_rcm_couple_to_visualizer[key]
if key in self.gcm_rcm_couples
@@ -92,7 +108,7 @@ class MetaVisualizerForProjectionEnsemble(object):
def ensemble_fits(self, ensemble_class):
return [ensemble_class_to_ensemble_fit[ensemble_class]
for ensemble_class_to_ensemble_fit
- in self.altitude_group_to_ensemble_class_to_ensemble_fit.values()]
+ in self.altitude_class_to_ensemble_class_to_ensemble_fit.values()]
def plot_together(self):
pass
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_sensitivity.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_sensitivity.py
new file mode 100644
index 0000000000000000000000000000000000000000..f3e40b69e2b800bb894c03b6ef234044a91dfffe
--- /dev/null
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_sensitivity.py
@@ -0,0 +1,115 @@
+from collections import OrderedDict
+import matplotlib.pyplot as plt
+from typing import List, Dict
+
+from extreme_data.meteo_france_data.adamont_data.cmip5.temperature_to_year import get_temp_min_and_temp_max, \
+ temperature_minmax_to_year_minmax
+from extreme_data.meteo_france_data.scm_models_data.utils import Season
+from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
+ 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, get_altitude_class_from_altitudes
+from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
+ plot_histogram_all_trends_against_altitudes, plot_shoe_plot_changes_against_altitude
+from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import compute_and_assign_max_abs
+from projects.projected_snowfall.elevation_temporal_model_for_projections.independent_ensemble_fit.independent_ensemble_fit import \
+ IndependentEnsembleFit
+from projects.projected_snowfall.elevation_temporal_model_for_projections.visualizer_for_projection_ensemble import \
+ VisualizerForProjectionEnsemble
+
+
+class VisualizerForSensivity(object):
+
+ def __init__(self, altitudes_list, gcm_rcm_couples, study_class, season, scenario,
+ model_classes: List[AbstractSpatioTemporalPolynomialModel],
+ ensemble_fit_classes=None,
+ massif_names=None,
+ fit_method=MarginFitMethod.extremes_fevd_mle,
+ temporal_covariate_for_fit=None,
+ display_only_model_that_pass_gof_test=False,
+ confidence_interval_based_on_delta_method=False,
+ remove_physically_implausible_models=False,
+ ):
+ self.altitudes_list = altitudes_list
+ self.massif_names = massif_names
+ self.temp_min, self.temp_max = get_temp_min_and_temp_max()
+ self.temp_min_to_temp_max = OrderedDict(zip(self.temp_min, self.temp_max))
+ self.temp_min_to_visualizer = {} # type: Dict[float, VisualizerForProjectionEnsemble]
+
+ for temp_min, temp_max in zip(self.temp_min, self.temp_max):
+ print(temp_min, temp_max)
+ # Build
+ gcm_to_year_min_and_year_max = {}
+ gcm_list = list(set([g for g, r in gcm_rcm_couples]))
+ for gcm in gcm_list:
+ year_min_and_year_max = temperature_minmax_to_year_minmax(gcm, scenario, temp_min, temp_max)
+ if year_min_and_year_max[0] is not None:
+ gcm_to_year_min_and_year_max[gcm] = year_min_and_year_max
+
+ visualizer = VisualizerForProjectionEnsemble(
+ altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
+ model_classes=model_classes,
+ fit_method=fit_method,
+ ensemble_fit_classes=ensemble_fit_classes,
+ display_only_model_that_pass_gof_test=display_only_model_that_pass_gof_test,
+ confidence_interval_based_on_delta_method=confidence_interval_based_on_delta_method,
+ massif_names=massif_names,
+ temporal_covariate_for_fit=temporal_covariate_for_fit,
+ remove_physically_implausible_models=remove_physically_implausible_models,
+ gcm_to_year_min_and_year_max=gcm_to_year_min_and_year_max
+ )
+ self.temp_min_to_visualizer[temp_min] = visualizer
+
+ def plot(self):
+ # todo: before reactivating the subplot, i should ensure that we can modify the prefix
+ # so that we can have all the subplot, and not juste for the last t_min
+ # for visualizer in self.temp_min_to_visualizer.values():
+ # visualizer.plot()
+ self.sensitivity_plot()
+
+ def sensitivity_plot(self):
+ ax = plt.gca()
+ for altitudes in self.altitudes_list:
+ altitude_class = get_altitude_class_from_altitudes(altitudes)
+ self.temperature_interval_plot(ax, altitude_class)
+
+ ticks_labels = ['+{} to +{}'.format(mi, ma) for mi, ma in zip(self.temp_min, self.temp_max)]
+ ax.set_ylabel('Percentages of massifs (%)')
+ ax.set_xlabel('Range of temperatures used to compute the trends ')
+ ax.set_xticks(self.temp_min)
+ ax.set_xticklabels(ticks_labels)
+ ax.legend()
+ mean_visualizer = self.first_mean_visualizer
+ mean_visualizer.plot_name = 'Sensitivity plot'
+ mean_visualizer.show_or_save_to_file(no_title=True)
+
+ @property
+ def first_mean_visualizer(self):
+ altitude_class = get_altitude_class_from_altitudes(self.altitudes_list[0])
+ visualizer_projection = list(self.temp_min_to_visualizer.values())[0]
+ return self.get_mean_visualizer(altitude_class, visualizer_projection)
+
+ def get_mean_visualizer(self, altitude_class, visualizer_projection: VisualizerForProjectionEnsemble):
+ independent_ensemble_fit = visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][
+ IndependentEnsembleFit]
+ mean_visualizer = independent_ensemble_fit.merge_function_name_to_visualizer[IndependentEnsembleFit.Mean_merge]
+ mean_visualizer.studies.study.gcm_rcm_couple = (IndependentEnsembleFit.Mean_merge, "merge")
+ return mean_visualizer
+
+ def temperature_interval_plot(self, ax, altitude_class):
+ label_to_l = {
+ 'increasing': [],
+ 'decreasing': []
+ }
+ for v in self.temp_min_to_visualizer.values():
+ mean_visualizer = self.get_mean_visualizer(altitude_class, v)
+ _, *trends = mean_visualizer.all_trends(self.massif_names, with_significance=False)
+ label_to_l['decreasing'].append(trends[0])
+ label_to_l['increasing'].append(trends[2])
+ altitude_str = altitude_class().formula
+ for label, l in label_to_l.items():
+ label_improved = 'with {} trends {}'.format(label, altitude_str)
+ ax.plot(self.temp_min, l, label=label_improved)
+
diff --git a/test/test_extreme_fit/test_estimator/test_full_estimators.py b/test/test_extreme_fit/test_estimator/test_full_estimators.py
index 9322f03459ba33915e3d0099492d256eb987e419..abbb240b144f7fabcb45edd8e012d7bd4ea143e9 100644
--- a/test/test_extreme_fit/test_estimator/test_full_estimators.py
+++ b/test/test_extreme_fit/test_estimator/test_full_estimators.py
@@ -1,6 +1,7 @@
import unittest
from itertools import product
+from extreme_fit.distribution.abstract_extreme_params import AbstractExtremeParams
from spatio_temporal_dataset.dataset.simulation_dataset import FullSimulatedDataset
from test.test_utils import load_test_max_stable_models, load_smooth_margin_models, load_test_1D_and_2D_spatial_coordinates, \
load_test_full_estimators
@@ -17,6 +18,7 @@ class TestFullEstimators(unittest.TestCase):
self.max_stable_models = load_test_max_stable_models()
def test_full_estimators(self):
+ AbstractExtremeParams.SMALL_SCALE_PARAMETERS_ARE_UNDEFINED = False
for coordinates in self.spatial_coordinates:
smooth_margin_models = load_smooth_margin_models(coordinates=coordinates)
for margin_model, max_stable_model in product(smooth_margin_models, self.max_stable_models):
diff --git a/test/test_projects/test_altitude_spatial/test_one_fold_fit.py b/test/test_projects/test_altitude_spatial/test_one_fold_fit.py
index 33daf73f40f49f089d600a4bf3f659bc41aad89d..f939768679cbbabbcd27a4173ea1fb0f6662ac93 100644
--- a/test/test_projects/test_altitude_spatial/test_one_fold_fit.py
+++ b/test/test_projects/test_altitude_spatial/test_one_fold_fit.py
@@ -2,6 +2,7 @@ import unittest
from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario
+from extreme_data.meteo_france_data.adamont_data.cmip5.temperature_to_year import temperature_minmax_to_year_minmax
from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import StationaryAltitudinal
@@ -11,6 +12,7 @@ from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pari
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.altitude_group import VeyHighAltitudeGroup
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
TimeTemporalCovariate, AnomalyTemperatureTemporalCovariate
@@ -75,6 +77,24 @@ class TestOneFoldFit(unittest.TestCase):
remove_physically_implausible_models=True)
self.assertFalse(one_fold_fit.has_at_least_one_valid_model)
+ def test_assertion_error_for_a_specific_case(self):
+ self.massif_name = "Thabor"
+ self.model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS[:]
+ self.altitudes = [3000, 3300, 3600]
+ gcm_rcm_couple = ('HadGEM2-ES', 'RegCM4-6')
+ scenario = AdamontScenario.rcp85
+ year_min, year_max = temperature_minmax_to_year_minmax(gcm_rcm_couple[0], scenario, temperature_min=1.0,
+ temperature_max=2.5)
+ dataset = self.load_dataset(AdamontSnowfall,
+ scenario=scenario, gcm_rcm_couple=gcm_rcm_couple,
+ year_min=year_min, year_max=year_max)
+ one_fold_fit = OneFoldFit(self.massif_name, dataset,
+ models_classes=self.model_classes,
+ temporal_covariate_for_fit=AnomalyTemperatureTemporalCovariate,
+ altitude_class=VeyHighAltitudeGroup,
+ only_models_that_pass_goodness_of_fit_test=False,
+ remove_physically_implausible_models=True)
+ self.assertTrue(one_fold_fit.has_at_least_one_valid_model)
if __name__ == '__main__':
unittest.main()