Commit 51e32a62 authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[projection snowfall] add main_projections_ensemble.py where we try together_ensemble_fit.py

parent aafb602c
No related merge requests found
Showing with 234 additions and 99 deletions
+234 -99
......@@ -86,7 +86,7 @@ def plot_nb_data(is_temperature_interval, is_shift_interval):
]
ax2.legend(handles=legend_elements, loc='upper center')
ax2.set_yticks([])
# plt.show()
plt.show()
def get_interval_limits(is_temperature_interval, is_shift_interval):
......@@ -111,13 +111,13 @@ def get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval):
ticks_labels = [' +${}^o\mathrm{C}$ and +${}^o\mathrm{C}$'.format(left_limit, right_limit, **{'C': '{C}'})
if is_temperature_interval else '{} and {}'.format(left_limit, right_limit)
for left_limit, right_limit in zip(left_limits, right_limits)]
prefix = 'Maxima occured between \n'
prefix = 'Maxima between \n'
ticks_labels = [prefix + l for l in ticks_labels]
return ticks_labels
if __name__ == '__main__':
for shift_interval in [False, True]:
for temp_interval in [True, False]:
for temp_interval in [False, True]:
print("shift = {}, temp_inteval = {}".format(shift_interval, temp_interval))
plot_nb_data(is_temperature_interval=temp_interval, is_shift_interval=shift_interval)
......@@ -2,9 +2,12 @@ from typing import Dict, Tuple
from extreme_fit.model.margin_model.utils import MarginFitMethod
from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies
from extreme_trend.one_fold_fit.one_fold_fit import OneFoldFit
class AbstractEnsembleFit(object):
Median_merge = 'Median'
Mean_merge = 'Mean'
def __init__(self, massif_names, gcm_rcm_couple_to_altitude_studies: Dict[Tuple[str, str], AltitudesStudies],
models_classes,
......@@ -23,6 +26,14 @@ class AbstractEnsembleFit(object):
self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test
self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
# Set appropriate setting
OneFoldFit.last_year = 2100
OneFoldFit.nb_years = 95
@property
def altitudes(self):
raise self.visualizer_list.studies.altitudes
@property
def visualizer_list(self):
raise NotImplementedError
\ No newline at end of file
"""Instead of creating one big group, with all the ensemble together,
and assuming a common temporal trend to all the group.
We could create smaller groups, with an importance proportional to the number of GCM/RCM couples considered
in the group.
For instance, we could group them by GCM, or group them by RCM.
Or we could try to find a metric to group them together.
This is the idea of finding of sweet spot between:
-only independent fits with few assumptions
-one common fit with too much assumption
it links with the idea of "climate model subset".
Generally people try to find one model subset,
the idea here, would be to find group of model subsets
"""
\ No newline at end of file
......@@ -9,15 +9,10 @@ from extreme_trend.one_fold_fit.one_fold_fit import OneFoldFit
class IndependentEnsembleFit(AbstractEnsembleFit):
"""For each gcm_rcm_couple, we create a OneFoldFit"""
Median_merge = 'Median'
Mean_merge = 'Mean'
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
# Set appropriate setting
OneFoldFit.last_year = 2100
OneFoldFit.nb_years = 95
# Load a classical visualizer
self.gcm_rcm_couple_to_visualizer = {}
for gcm_rcm_couple, studies in self.gcm_rcm_couple_to_altitude_studies.items():
......@@ -47,9 +42,7 @@ class IndependentEnsembleFit(AbstractEnsembleFit):
}
@property
def visualizer(self):
return list(self.gcm_rcm_couple_to_visualizer.values())[0]
def visualizer_list(self):
return list(self.gcm_rcm_couple_to_visualizer.values()) \
+ list(self.merge_function_name_to_visualizer.values())
@property
def altitudes(self):
return self.visualizer.studies.altitudes
from extreme_trend.ensemble_fit.abstract_ensemble_fit import AbstractEnsembleFit
from extreme_trend.ensemble_fit.together_ensemble_fit.visualizer_non_stationary_ensemble import \
VisualizerNonStationaryEnsemble
class TogetherEnsembleFit(AbstractEnsembleFit):
"""We create a single OneFoldFit for all gcm_rcm_couples"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self.visualizer = VisualizerNonStationaryEnsemble(self.gcm_rcm_couple_to_altitude_studies,
self.models_classes,
False,
self.massif_names, self.fit_method,
self.temporal_covariate_for_fit,
self.only_models_that_pass_goodness_of_fit_test,
self.confidence_interval_based_on_delta_method,
self.remove_physically_implausible_models
)
@property
def visualizer_list(self):
return [self.visualizer]
from typing import List, Dict, Tuple
from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies
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 extreme_trend.one_fold_fit.altitudes_studies_visualizer_for_non_stationary_models import \
AltitudesStudiesVisualizerForNonStationaryModels
from extreme_trend.trend_test.visualizers.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
class VisualizerNonStationaryEnsemble(AltitudesStudiesVisualizerForNonStationaryModels):
def __init__(self, gcm_rcm_couple_to_studies: Dict[Tuple[str, str], AltitudesStudies], *args, **kwargs):
self.gcm_rcm_couple_to_studies = gcm_rcm_couple_to_studies
studies = list(self.gcm_rcm_couple_to_studies.values())[0]
super().__init__(studies, *args, **kwargs)
def get_massif_altitudes(self, massif_name):
# return self._get_massif_altitudes(massif_name, self.studies)
raise NotImplementedError
def get_dataset(self, massif_altitudes, massif_name):
raise NotImplementedError
# dataset = self.studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
# return dataset
......@@ -5,7 +5,9 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
AbstractSpatioTemporalPolynomialModel
from extreme_fit.model.margin_model.utils import MarginFitMethod
from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies
from extreme_trend.ensemble_fit.abstract_ensemble_fit import AbstractEnsembleFit
from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
from extreme_trend.ensemble_fit.together_ensemble_fit.together_ensemble_fit import TogetherEnsembleFit
from extreme_trend.one_fold_fit.altitude_group import get_altitude_class_from_altitudes
from extreme_trend.one_fold_fit.plots.plot_histogram_altitude_studies import \
plot_histogram_all_trends_against_altitudes, plot_shoe_plot_changes_against_altitude
......@@ -75,29 +77,35 @@ class VisualizerForProjectionEnsemble(object):
ensemble_class_to_ensemble_fit[ensemble_fit_class] = ensemble_fit
self.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class] = ensemble_class_to_ensemble_fit
def plot_for_visualizer_list(self, visualizer_list):
with_significance = False
for v in visualizer_list:
v.plot_moments()
plot_histogram_all_trends_against_altitudes(self.massif_names, visualizer_list,
with_significance=with_significance)
for relative in [True, False]:
plot_shoe_plot_changes_against_altitude(self.massif_names, visualizer_list, relative=relative,
with_significance=with_significance)
def plot(self):
# Set limit for the plot
visualizer_list = []
for ensemble_fit_class in self.ensemble_fit_classes:
for ensemble_fit in self.ensemble_fits(ensemble_fit_class):
visualizer_list.extend(ensemble_fit.visualizer_list)
compute_and_assign_max_abs(visualizer_list)
# Plot
if IndependentEnsembleFit in self.ensemble_fit_classes:
# Set max abs
visualizer_list = []
for ensemble_fit in self.ensemble_fits(IndependentEnsembleFit):
visualizer_list.extend(list(ensemble_fit.gcm_rcm_couple_to_visualizer.values()))
# Potentially I could add more visualizer here...
method_name_and_order_to_max_abs, max_abs_for_shape = compute_and_assign_max_abs(visualizer_list)
# Assign the same max abs for the
for ensemble_fit in self.ensemble_fits(IndependentEnsembleFit):
for v in ensemble_fit.merge_function_name_to_visualizer.values():
v._max_abs_for_shape = max_abs_for_shape
v._method_name_and_order_to_max_abs = method_name_and_order_to_max_abs
# Plot
self.plot_independent()
if TogetherEnsembleFit in self.ensemble_fit_classes:
self.plot_together()
def plot_independent(self):
with_significance = False
# Aggregated at gcm_rcm_level plots
merge_keys = [IndependentEnsembleFit.Median_merge, IndependentEnsembleFit.Mean_merge]
merge_keys = [AbstractEnsembleFit.Median_merge, AbstractEnsembleFit.Mean_merge]
keys = self.gcm_rcm_couples + merge_keys
# Only plot Mean for speed
keys = [IndependentEnsembleFit.Mean_merge]
keys = [AbstractEnsembleFit.Mean_merge]
for key in keys:
visualizer_list = [independent_ensemble_fit.gcm_rcm_couple_to_visualizer[key]
if key in self.gcm_rcm_couples
......@@ -107,17 +115,15 @@ class VisualizerForProjectionEnsemble(object):
if key in merge_keys:
for v in visualizer_list:
v.studies.study.gcm_rcm_couple = (key, "merge")
for v in visualizer_list:
v.plot_moments()
plot_histogram_all_trends_against_altitudes(self.massif_names, visualizer_list, with_significance=with_significance)
for relative in [True, False]:
plot_shoe_plot_changes_against_altitude(self.massif_names, visualizer_list, relative=relative, with_significance=with_significance)
self.plot_for_visualizer_list(visualizer_list)
def plot_together(self):
visualizer_list = [together_ensemble_fit.visualizer
for together_ensemble_fit in self.ensemble_fits(TogetherEnsembleFit)]
self.plot_for_visualizer_list(visualizer_list)
def ensemble_fits(self, ensemble_class):
"""Return the ordered ensemble fit for a given ensemble class (in the order of the altitudes)"""
return [ensemble_class_to_ensemble_fit[ensemble_class]
for ensemble_class_to_ensemble_fit
in self.altitude_class_to_ensemble_class_to_ensemble_fit.values()]
def plot_together(self):
pass
......@@ -8,6 +8,7 @@ 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 extreme_trend.ensemble_fit.abstract_ensemble_fit import AbstractEnsembleFit
from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
from extreme_trend.one_fold_fit.altitude_group import get_altitude_class_from_altitudes, \
......@@ -25,7 +26,7 @@ class VisualizerForSensivity(object):
display_only_model_that_pass_gof_test=False,
confidence_interval_based_on_delta_method=False,
remove_physically_implausible_models=False,
merge_visualizer_str=IndependentEnsembleFit.Median_merge, # if we choose the Mean merge, then it is almost impossible to obtain stationary trends
merge_visualizer_str=AbstractEnsembleFit.Median_merge, # if we choose the Mean merge, then it is almost impossible to obtain stationary trends
is_temperature_interval=False,
is_shift_interval=False,
):
......
......@@ -79,7 +79,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
# Save the massif altitudes only for those who pass the condition
self.massif_name_to_massif_altitudes[massif_name] = massif_altitudes
# Load dataset
dataset = self.studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
dataset = self.get_dataset(massif_altitudes, massif_name)
old_fold_fit = OneFoldFit(massif_name, dataset, self.model_classes, self.fit_method,
self.temporal_covariate_for_fit,
type(self.altitude_group),
......@@ -90,15 +90,22 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
else:
return None
def get_dataset(self, massif_altitudes, massif_name):
dataset = self.studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
return dataset
moment_names = ['moment', 'changes_of_moment', 'relative_changes_of_moment'][:]
orders = [1, 2, None][2:]
def get_massif_altitudes(self, massif_name):
valid_altitudes = [altitude for altitude, study in self.studies.altitude_to_study.items()
return self._get_massif_altitudes(massif_name, self.studies)
def _get_massif_altitudes(self, massif_name, studies):
valid_altitudes = [altitude for altitude, study in studies.altitude_to_study.items()
if massif_name in study.study_massif_names]
massif_altitudes = []
for altitude in valid_altitudes:
study = self.studies.altitude_to_study[altitude]
study = studies.altitude_to_study[altitude]
annual_maxima = study.massif_name_to_annual_maxima[massif_name]
percentage_of_non_zeros = 100 * np.count_nonzero(annual_maxima) / len(annual_maxima)
if percentage_of_non_zeros > 90:
......
......@@ -135,7 +135,7 @@ def main_comparaison_plot():
ax.set_ylabel('Altitude (m)', fontsize=10)
massif_str = 'all massifs' if massif_names is None else 'the {} massif'.format(massif_names[0])
unit = '%' if relative_bias else study.variable_unit
bias_name = 'Relative bias' if relative_bias else 'Bias'
bias_name = 'Relative difference' if relative_bias else 'Difference'
mean_str = 'mean' if mean else 'std'
title = '{} in the {} annual maxima of {} of {}\n' \
'for ADAMONT v{}' \
......@@ -144,7 +144,7 @@ def main_comparaison_plot():
adamont_version,
comparaison_study_class,
unit)
folder = 'relative bias' if relative_bias else 'absolute bias'
folder = 'relative difference' if relative_bias else 'difference'
plot_name = op.join(folder, title)
ax.set_xlabel(title, fontsize=10)
reanalysis_altitude_studies.show_or_save_to_file(plot_name=plot_name, no_title=True)
......
import datetime
import time
from typing import List
import matplotlib
from extreme_trend.ensemble_fit.together_ensemble_fit.together_ensemble_fit import TogetherEnsembleFit
matplotlib.use('Agg')
import matplotlib as mpl
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
......@@ -11,20 +17,16 @@ from extreme_trend.ensemble_fit.visualizer_for_sensitivity import VisualizerForS
from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
AnomalyTemperatureWithSplineTemporalCovariate
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
import matplotlib
from extreme_fit.model.utils import set_seed_for_test
from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
rcp_scenarios
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import TimeTemporalCovariate
matplotlib.use('Agg')
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
TimeTemporalCovariate
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
AbstractExtractEurocodeReturnLevel
......@@ -37,21 +39,20 @@ from extreme_data.meteo_france_data.scm_models_data.utils import Season
def main():
start = time.time()
study_class = AdamontSnowfall
ensemble_fit_class = [IndependentEnsembleFit]
ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][1:]
temporal_covariate_for_fit = [TimeTemporalCovariate,
AnomalyTemperatureWithSplineTemporalCovariate][0]
set_seed_for_test()
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
fast = True
sensitivity_plot = True
scenarios = rcp_scenarios[::-1] 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 = None
gcm_rcm_couples = gcm_rcm_couples
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
altitudes_list = altitudes_for_groups[3:]
elif fast:
......@@ -65,35 +66,15 @@ def main():
assert isinstance(gcm_rcm_couples, list)
main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_class, ensemble_fit_class, scenario,
temporal_covariate_for_fit, sensitivity_plot=sensitivity_plot)
assert isinstance(altitudes_list, List)
assert isinstance(altitudes_list[0], List)
print('Scenario is', scenario)
print('Covariate is {}'.format(temporal_covariate_for_fit))
end = time.time()
duration = str(datetime.timedelta(seconds=end - start))
print('Total duration', duration)
model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
assert scenario in rcp_scenarios
remove_physically_implausible_models = True
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('Scenario is', scenario)
print('Covariate is {}'.format(temporal_covariate_for_fit))
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,
......@@ -103,7 +84,11 @@ def main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_class, ensemb
remove_physically_implausible_models=remove_physically_implausible_models,
gcm_to_year_min_and_year_max=None,
)
visualizer.plot()
visualizer.plot()
end = time.time()
duration = str(datetime.timedelta(seconds=end - start))
print('Total duration', duration)
if __name__ == '__main__':
......
import datetime
import time
from typing import List
import matplotlib
from extreme_trend.ensemble_fit.abstract_ensemble_fit import AbstractEnsembleFit
matplotlib.use('Agg')
import matplotlib as mpl
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
from extreme_trend.ensemble_fit.visualizer_for_sensitivity import VisualizerForSensivity
from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
AnomalyTemperatureWithSplineTemporalCovariate
from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
from extreme_fit.model.utils import set_seed_for_test
from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
rcp_scenarios
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import TimeTemporalCovariate
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
AbstractExtractEurocodeReturnLevel
from extreme_trend.one_fold_fit.altitude_group import altitudes_for_groups
from extreme_data.meteo_france_data.scm_models_data.utils import Season
def main():
start = time.time()
study_class = AdamontSnowfall
ensemble_fit_classes = [IndependentEnsembleFit]
temporal_covariate_for_fit = [TimeTemporalCovariate,
AnomalyTemperatureWithSplineTemporalCovariate][0]
set_seed_for_test()
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
fast = False
scenarios = rcp_scenarios[-1:] 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
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
altitudes_list = altitudes_for_groups[1:2]
elif fast:
massif_names = ['Vanoise', 'Haute-Maurienne']
gcm_rcm_couples = gcm_rcm_couples[4:6]
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
altitudes_list = altitudes_for_groups[:1]
else:
massif_names = None
altitudes_list = altitudes_for_groups[:]
assert isinstance(gcm_rcm_couples, list)
assert isinstance(altitudes_list, List)
assert isinstance(altitudes_list[0], List)
print('Scenario is', scenario)
print('Covariate is {}'.format(temporal_covariate_for_fit))
model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
assert scenario in rcp_scenarios
remove_physically_implausible_models = True
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,
merge_visualizer_str=AbstractEnsembleFit.Median_merge,
is_temperature_interval=False,
is_shift_interval=True,
)
visualizer.plot()
end = time.time()
duration = str(datetime.timedelta(seconds=end - start))
print('Total duration', duration)
if __name__ == '__main__':
main()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment