From 5bc105b833d78d6f747ce9fe3d3ec7a46aeddcbd Mon Sep 17 00:00:00 2001 From: Le Roux Erwan <erwan.le-roux@irstea.fr> Date: Tue, 16 Mar 2021 17:16:01 +0100 Subject: [PATCH] [projection swe] add model_as_truth.py. add default weight solver and a file for the method with bootstrap --- .../main_model_as_truth.py | 68 ++++++++++ .../model_as_truth.py | 118 ++++++++++++++++++ .../weight_solver/abtract_weight_solver.py | 52 +++++++- .../weight_solver/default_weight_solver.py | 16 +++ .../projected_swe/weight_solver/indicator.py | 12 +- .../weight_solver/knutti_weight_solver.py | 30 ++--- .../knutti_weight_solver_with_bootstrap.py | 32 +++++ .../test_projected_swe/test_model_as_truth.py | 7 +- 8 files changed, 301 insertions(+), 34 deletions(-) create mode 100644 projects/projected_swe/model_as_truth_visualizer/main_model_as_truth.py create mode 100644 projects/projected_swe/model_as_truth_visualizer/model_as_truth.py create mode 100644 projects/projected_swe/weight_solver/default_weight_solver.py create mode 100644 projects/projected_swe/weight_solver/knutti_weight_solver_with_bootstrap.py diff --git a/projects/projected_swe/model_as_truth_visualizer/main_model_as_truth.py b/projects/projected_swe/model_as_truth_visualizer/main_model_as_truth.py new file mode 100644 index 00000000..86688952 --- /dev/null +++ b/projects/projected_swe/model_as_truth_visualizer/main_model_as_truth.py @@ -0,0 +1,68 @@ +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 +from extreme_data.meteo_france_data.scm_models_data.safran.safran_max_snowf import SafranSnowfall2020 +from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \ + AbstractExtractEurocodeReturnLevel +from projects.projected_swe.model_as_truth_visualizer.model_as_truth import ModelAsTruth +from projects.projected_swe.weight_solver.default_weight_solver import EqualWeight +from projects.projected_swe.weight_solver.indicator import AnnualMaximaMeanIndicator, ReturnLevel30YearsIndicator +from projects.projected_swe.weight_solver.knutti_weight_solver import KnuttiWeightSolver +from projects.projected_swe.weight_solver.knutti_weight_solver_with_bootstrap import \ + KnuttiWeightSolverWithBootstrapVersion1, KnuttiWeightSolverWithBootstrapVersion2 + + +def main(): + altitude = 900 + year_min_histo = 1982 + year_max_histo = 2011 + year_min_projected = 2070 + year_max_projected = 2099 + scenario = AdamontScenario.rcp85_extended + fast = None + gcm_rcm_couples = get_gcm_rcm_couples(adamont_scenario=scenario) + + if fast is None: + AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10 + massif_names = None + knutti_weight_solver_classes = [EqualWeight, + KnuttiWeightSolver, + KnuttiWeightSolverWithBootstrapVersion1, + KnuttiWeightSolverWithBootstrapVersion2] + indicator_class = ReturnLevel30YearsIndicator + gcm_rcm_couples = gcm_rcm_couples[:3] + sigma_list = [10, 100, 1000, 10000] + + elif fast: + AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10 + massif_names = ['Vercors'] + + knutti_weight_solver_classes = [EqualWeight, KnuttiWeightSolver, + KnuttiWeightSolverWithBootstrapVersion1] + indicator_class = ReturnLevel30YearsIndicator + gcm_rcm_couples = gcm_rcm_couples[:3] + sigma_list = [100, 1000] + + else: + massif_names = None + indicator_class = AnnualMaximaMeanIndicator + knutti_weight_solver_classes = [KnuttiWeightSolver, + KnuttiWeightSolverWithBootstrapVersion1, + KnuttiWeightSolverWithBootstrapVersion2] + + observation_study = SafranSnowfall2020(altitude=altitude, year_min=year_min_histo, year_max=year_max_histo) + couple_to_historical_study = {c: AdamontSnowfall(altitude=altitude, scenario=scenario, + year_min=year_min_histo, year_max=year_max_histo, + gcm_rcm_couple=c) for c in gcm_rcm_couples} + couple_to_projected_study = {c: AdamontSnowfall(altitude=altitude, scenario=scenario, + year_min=year_min_projected, year_max=year_max_projected, + gcm_rcm_couple=c) for c in gcm_rcm_couples + } + + model_as_truth = ModelAsTruth(observation_study, couple_to_projected_study, couple_to_historical_study, + indicator_class, knutti_weight_solver_classes, massif_names, + add_interdependence_weight=False) + model_as_truth.plot_against_sigma(sigma_list) + + +if __name__ == '__main__': + main() diff --git a/projects/projected_swe/model_as_truth_visualizer/model_as_truth.py b/projects/projected_swe/model_as_truth_visualizer/model_as_truth.py new file mode 100644 index 00000000..4eb1267d --- /dev/null +++ b/projects/projected_swe/model_as_truth_visualizer/model_as_truth.py @@ -0,0 +1,118 @@ +import matplotlib.pyplot as plt +from typing import Dict, Tuple, List + +from matplotlib.lines import Line2D +from scipy.special import softmax +import numpy as np + +from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy +from projects.projected_swe.weight_solver.abtract_weight_solver import AbstractWeightSolver +from projects.projected_swe.weight_solver.default_weight_solver import EqualWeight +from projects.projected_swe.weight_solver.indicator import AbstractIndicator +from projects.projected_swe.weight_solver.knutti_weight_solver import KnuttiWeightSolver +from projects.projected_swe.weight_solver.knutti_weight_solver_with_bootstrap import \ + KnuttiWeightSolverWithBootstrapVersion2, KnuttiWeightSolverWithBootstrapVersion1 +from root_utils import get_display_name_from_object_type + + +class ModelAsTruth(object): + + def __init__(self, observation_study: AbstractStudy, + couple_to_study_projected: Dict[Tuple[str, str], AbstractStudy], + couple_to_study_historical: Dict[Tuple[str, str], AbstractStudy], + indicator_class: type, + knutti_weight_solver_classes, + massif_names=None, + add_interdependence_weight=False): + self.knutti_weight_solver_classes = knutti_weight_solver_classes + self.massif_names = massif_names + self.add_interdependence_weight = add_interdependence_weight + self.indicator_class = indicator_class + self.couple_to_study_historical = couple_to_study_historical + self.couple_to_study_projected = couple_to_study_projected + self.observation_study = observation_study + # Parameters + self.width = 2 + self.solver_class_to_color = { + EqualWeight: 'white', + KnuttiWeightSolver: 'yellow', + KnuttiWeightSolverWithBootstrapVersion1: 'orange', + KnuttiWeightSolverWithBootstrapVersion2: 'red', + } + + # Update massif names + study_list = [observation_study] + list(couple_to_study_historical.values()) + list(couple_to_study_projected.values()) + self.massif_names = self.get_massif_names_subset_from_study_list(study_list) + print('Nb of massifs') + print(len(self.massif_names)) + # Set some parameters to speed up results (by caching some results) + for study in study_list: + study._massif_names_for_cache = self.massif_names + + def plot_against_sigma(self, sigma_list): + ax = plt.gca() + solver_class_to_score_list = self.get_solver_class_to_score_list(sigma_list) + all_x = [] + labels = [] + colors = [] + for j, (solver_class, score_list) in enumerate(list(solver_class_to_score_list.items())): + x_list = self.get_x_list(j, sigma_list) + all_x.extend(x_list) + assert len(x_list) == len(sigma_list) + label = get_display_name_from_object_type(solver_class) + color = self.solver_class_to_color[solver_class] + bplot = ax.boxplot(score_list, positions=x_list, widths=self.width, patch_artist=True, showmeans=True, + labels=[str(sigma) for sigma in sigma_list]) + for patch in bplot['boxes']: + patch.set_facecolor(color) + colors.append(color) + labels.append(label) + + custom_lines = [Line2D([0], [0], color=color, lw=4) for color in colors] + ax.legend(custom_lines, labels, prop={'size': 8}, loc='upper left') + ax.set_xlim(min(all_x) - self.width, max(all_x) + self.width) + _, max_y = ax.get_ylim() + ax.set_ylim((0, max_y * 1.1)) + + plt.show() + + def get_x_list(self, j, sigma_list): + shift = len(self.knutti_weight_solver_classes) + 1 + x_list = [((j + 1) * (self.width * 1.1)) + shift * i * self.width for i in range(len(sigma_list))] + return x_list + + def get_solver_class_to_score_list(self, sigma_list): + return {solver_class: [self.compute_score(solver_class, sigma) for sigma in sigma_list] + for solver_class in self.knutti_weight_solver_classes} + + def compute_score(self, solver_class, sigma): + # return [sigma, sigma*2] + score_list = [] + for gcm_rcm_couple in self.couple_to_study_historical.keys(): + historical_observation_study = self.couple_to_study_historical[gcm_rcm_couple] + projected_observation_study = self.couple_to_study_projected[gcm_rcm_couple] + couple_to_study_historical = {c: s for c, s in self.couple_to_study_historical.items() if + c != gcm_rcm_couple} + couple_to_study_projected = {c: s for c, s in self.couple_to_study_projected.items() if c != gcm_rcm_couple} + + if issubclass(solver_class, KnuttiWeightSolver): + weight_solver = solver_class(sigma, None, historical_observation_study, couple_to_study_historical, + self.indicator_class, self.massif_names, self.add_interdependence_weight, + ) # type: AbstractWeightSolver + else: + weight_solver = solver_class(historical_observation_study, couple_to_study_historical, + self.indicator_class, self.massif_names, self.add_interdependence_weight, + ) # type: AbstractWeightSolver + + print(solver_class, sigma, weight_solver.couple_to_weight.values()) + mean_score = weight_solver.mean_prediction_score(self.massif_names, couple_to_study_projected, + projected_observation_study) + score_list.append(mean_score) + return np.array(score_list) + + def get_massif_names_subset_from_study_list(self, study_list: List[AbstractStudy]): + massifs_names_list = [set(s.study_massif_names) for s in study_list] + massif_names_intersection = massifs_names_list[0].intersection(*massifs_names_list[1:]) + if self.massif_names is not None: + massif_names_intersection = massif_names_intersection.intersection(set(self.massif_names)) + return list(massif_names_intersection) diff --git a/projects/projected_swe/weight_solver/abtract_weight_solver.py b/projects/projected_swe/weight_solver/abtract_weight_solver.py index fdde441a..68844411 100644 --- a/projects/projected_swe/weight_solver/abtract_weight_solver.py +++ b/projects/projected_swe/weight_solver/abtract_weight_solver.py @@ -1,21 +1,22 @@ +import properscoring as ps from typing import Dict, Tuple from scipy.special import softmax import numpy as np from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy -from projects.projected_swe.weight_solver.indicator import AbstractIndicator +from projects.projected_swe.weight_solver.indicator import AbstractIndicator, ReturnLevelComputationException class AbstractWeightSolver(object): def __init__(self, observation_study: AbstractStudy, - couple_to_study: Dict[Tuple[str, str], AbstractStudy], + couple_to_historical_study: Dict[Tuple[str, str], AbstractStudy], indicator_class: type, massif_names=None, add_interdependence_weight=False): self.observation_study = observation_study - self.couple_to_study = couple_to_study + self.couple_to_historical_study = couple_to_historical_study self.indicator_class = indicator_class self.add_interdependence_weight = add_interdependence_weight # Compute intersection massif names @@ -27,9 +28,48 @@ class AbstractWeightSolver(object): assert set(massif_names).issubset(intersection_massif_names) self.massif_names = massif_names + # Prediction on the future period + + def weighted_projected_expected_indicator(self, massif_name, couple_to_projected_study): + couple_to_projected_expected_indicator = self.couple_to_projected_expected_indicator(massif_name, + couple_to_projected_study) + return sum([couple_to_projected_expected_indicator[c] * w for c, w in self.couple_to_weight.items()]) + + def couple_to_projected_expected_indicator(self, massif_name, couple_to_projected_study): + assert issubclass(self.indicator_class, AbstractIndicator) + return {c: self.indicator_class.get_indicator(s, massif_name) for c, s in couple_to_projected_study.items()} + + def prediction_score(self, massif_name, couple_to_projected_study, projected_observation_study): + try: + target = self.target(massif_name, projected_observation_study) + couple_to_projected_indicator = self.couple_to_projected_expected_indicator(massif_name, + couple_to_projected_study) + couples, ensemble = zip(*list(couple_to_projected_indicator.items())) + couple_to_weight = self.couple_to_weight + weights = [couple_to_weight[c] for c in couples] + return ps.crps_ensemble(target, ensemble, weights=weights) + except ReturnLevelComputationException: + return np.nan + + def mean_prediction_score(self, massif_names, couple_to_projected_study, projected_observation_study): + scores = [self.prediction_score(massif_name, couple_to_projected_study, projected_observation_study) for + massif_name in massif_names] + scores_filtered = [s for s in scores if not np.isnan(s)] + assert len(scores_filtered) > 0 + nb_massif_names_removed = len(scores) - len(scores_filtered) + if nb_massif_names_removed > 0: + print('{} massifs removed'.format(nb_massif_names_removed)) + return np.mean(scores_filtered) + + def target(self, massif_name, projected_observation_study): + assert issubclass(self.indicator_class, AbstractIndicator) + return self.indicator_class.get_indicator(projected_observation_study, massif_name, bootstrap=True).mean() + + # Weight computation on the historical period + @property def study_list(self): - return [self.observation_study] + list(self.couple_to_study.values()) + return [self.observation_study] + list(self.couple_to_historical_study.values()) @property def couple_to_weight(self): @@ -48,7 +88,7 @@ class AbstractWeightSolver(object): @property def couple_to_nllh_skill(self): return {couple: self.compute_skill_nllh(couple_study=couple_study) - for couple, couple_study in self.couple_to_study.items()} + for couple, couple_study in self.couple_to_historical_study.items()} def compute_skill_nllh(self, couple_study): raise NotImplementedError @@ -56,7 +96,7 @@ class AbstractWeightSolver(object): @property def couple_to_nllh_interdependence(self): return {couple: self.compute_interdependence_nllh(couple_study=couple_study) - for couple, couple_study in self.couple_to_study.items()} + for couple, couple_study in self.couple_to_historical_study.items()} def compute_interdependence_nllh(self, couple_study): raise NotImplementedError diff --git a/projects/projected_swe/weight_solver/default_weight_solver.py b/projects/projected_swe/weight_solver/default_weight_solver.py new file mode 100644 index 00000000..d2cd76df --- /dev/null +++ b/projects/projected_swe/weight_solver/default_weight_solver.py @@ -0,0 +1,16 @@ +from typing import Dict, Tuple + +from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy +from projects.projected_swe.weight_solver.abtract_weight_solver import AbstractWeightSolver + + +class EqualWeight(AbstractWeightSolver): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + @property + def couple_to_weight(self): + couples = list(self.couple_to_historical_study.keys()) + weight = 1 / len(couples) + return {c: weight for c in couples} diff --git a/projects/projected_swe/weight_solver/indicator.py b/projects/projected_swe/weight_solver/indicator.py index 1e086777..cc38a3ea 100644 --- a/projects/projected_swe/weight_solver/indicator.py +++ b/projects/projected_swe/weight_solver/indicator.py @@ -34,10 +34,10 @@ class ReturnLevel30YearsIndicator(AbstractIndicator): @classmethod def get_indicator(cls, study: AbstractStudy, massif_name, bootstrap=False): - if bootstrap: - return study.massif_name_to_return_level_list_from_bootstrap(return_period=30)[massif_name] - else: - try: + try: + if bootstrap: + return study.massif_name_to_return_level_list_from_bootstrap(return_period=30)[massif_name] + else: return study.massif_name_to_return_level(return_period=30)[massif_name] - except KeyError: - raise ReturnLevelComputationException + except KeyError: + raise ReturnLevelComputationException diff --git a/projects/projected_swe/weight_solver/knutti_weight_solver.py b/projects/projected_swe/weight_solver/knutti_weight_solver.py index 4f1a230a..9a6ae225 100644 --- a/projects/projected_swe/weight_solver/knutti_weight_solver.py +++ b/projects/projected_swe/weight_solver/knutti_weight_solver.py @@ -13,18 +13,22 @@ class KnuttiWeightSolver(AbstractWeightSolver): super().__init__(*args, **kwargs) self.sigma_skill = sigma_skill self.sigma_interdependence = sigma_interdependence + if self.add_interdependence_weight: + assert self.sigma_interdependence is not None # Set some parameters for the bootstrap ReturnLevelBootstrap.only_physically_plausible_fits = True # Set some parameters to speed up results (by caching some results) - self.observation_study._massif_names_for_cache = self.massif_names - for study in self.couple_to_study.values(): - study._massif_names_for_cache = self.massif_names + study_list = [self.observation_study] + list(self.couple_to_historical_study.values()) + for study in study_list: + if study._massif_names_for_cache is None: + study._massif_names_for_cache = self.massif_names # Compute the subset of massif_names used for the computation self.massif_names_for_computation = [] for massif_name in self.massif_names: try: [self.compute_skill_one_massif(couple_study, massif_name) for couple_study in self.study_list] - [self.compute_interdependence_nllh_one_massif(couple_study, massif_name) for couple_study in self.study_list] + if self.add_interdependence_weight: + [self.compute_interdependence_nllh_one_massif(couple_study, massif_name) for couple_study in self.study_list] except WeightComputationException: continue self.massif_names_for_computation.append(massif_name) @@ -47,7 +51,7 @@ class KnuttiWeightSolver(AbstractWeightSolver): def compute_interdependence_nllh_one_massif(self, couple_study, massif_name): sum_proba = 0 - for other_couple_study in self.couple_to_study.values(): + for other_couple_study in self.couple_to_historical_study.values(): if other_couple_study is not couple_study: nllh = self.compute_nllh_from_two_study(couple_study, other_couple_study, self.sigma_interdependence, massif_name) @@ -72,23 +76,11 @@ class KnuttiWeightSolver(AbstractWeightSolver): assert issubclass(self.indicator_class, AbstractIndicator) return np.array([self.indicator_class.get_indicator(study_1, massif_name) - self.indicator_class.get_indicator(study_2, massif_name)]) + + -class KnuttiWeightSolverWithBootstrapVersion1(KnuttiWeightSolver): - def differences(self, study_1, study_2, massif_name): - assert issubclass(self.indicator_class, AbstractIndicator) - bootstrap_study_1 = self.indicator_class.get_indicator(study_1, massif_name, bootstrap=True) - bootstrap_study_2 = self.indicator_class.get_indicator(study_2, massif_name, bootstrap=True) - differences = bootstrap_study_1 - bootstrap_study_2 - return differences -class KnuttiWeightSolverWithBootstrapVersion2(KnuttiWeightSolver): - def differences(self, study_1, study_2, massif_name): - assert issubclass(self.indicator_class, AbstractIndicator) - bootstrap_study_1 = self.indicator_class.get_indicator(study_1, massif_name, bootstrap=True) - bootstrap_study_2 = self.indicator_class.get_indicator(study_2, massif_name, bootstrap=True) - differences = np.subtract.outer(bootstrap_study_1, bootstrap_study_2) - return differences diff --git a/projects/projected_swe/weight_solver/knutti_weight_solver_with_bootstrap.py b/projects/projected_swe/weight_solver/knutti_weight_solver_with_bootstrap.py new file mode 100644 index 00000000..e6af3be6 --- /dev/null +++ b/projects/projected_swe/weight_solver/knutti_weight_solver_with_bootstrap.py @@ -0,0 +1,32 @@ +import numpy as np + +from projects.projected_swe.weight_solver.indicator import AbstractIndicator +from projects.projected_swe.weight_solver.knutti_weight_solver import KnuttiWeightSolver + + +class KnuttiWeightSolverWithBootstrap(KnuttiWeightSolver): + + def couple_to_projected_expected_indicator(self, massif_name, couple_to_projected_study): + assert issubclass(self.indicator_class, AbstractIndicator) + return {c: self.indicator_class.get_indicator(s, massif_name, bootstrap=True).mean() + for c, s in couple_to_projected_study.items()} + + +class KnuttiWeightSolverWithBootstrapVersion1(KnuttiWeightSolverWithBootstrap): + + def differences(self, study_1, study_2, massif_name): + assert issubclass(self.indicator_class, AbstractIndicator) + bootstrap_study_1 = self.indicator_class.get_indicator(study_1, massif_name, bootstrap=True) + bootstrap_study_2 = self.indicator_class.get_indicator(study_2, massif_name, bootstrap=True) + differences = bootstrap_study_1 - bootstrap_study_2 + return differences + + +class KnuttiWeightSolverWithBootstrapVersion2(KnuttiWeightSolverWithBootstrap): + + def differences(self, study_1, study_2, massif_name): + assert issubclass(self.indicator_class, AbstractIndicator) + bootstrap_study_1 = self.indicator_class.get_indicator(study_1, massif_name, bootstrap=True) + bootstrap_study_2 = self.indicator_class.get_indicator(study_2, massif_name, bootstrap=True) + differences = np.subtract.outer(bootstrap_study_1, bootstrap_study_2) + return differences diff --git a/test/test_projects/test_projected_swe/test_model_as_truth.py b/test/test_projects/test_projected_swe/test_model_as_truth.py index 6b10b921..95f0a251 100644 --- a/test/test_projects/test_projected_swe/test_model_as_truth.py +++ b/test/test_projects/test_projected_swe/test_model_as_truth.py @@ -8,7 +8,8 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran_max_snowf impo from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \ AbstractExtractEurocodeReturnLevel from projects.projected_swe.weight_solver.indicator import AnnualMaximaMeanIndicator, ReturnLevel30YearsIndicator -from projects.projected_swe.weight_solver.knutti_weight_solver import KnuttiWeightSolver, \ +from projects.projected_swe.weight_solver.knutti_weight_solver import KnuttiWeightSolver +from projects.projected_swe.weight_solver.knutti_weight_solver_with_bootstrap import \ KnuttiWeightSolverWithBootstrapVersion1, KnuttiWeightSolverWithBootstrapVersion2 @@ -37,11 +38,11 @@ class TestModelAsTruth(unittest.TestCase): knutti_weight = knutti_weight_solver_class(sigma_skill=100.0, sigma_interdependence=100.0, massif_names=massif_names, observation_study=observation_study, - couple_to_study=couple_to_study, + couple_to_historical_study=couple_to_study, indicator_class=indicator_class, add_interdependence_weight=add_interdependence_weight ) - print(knutti_weight.couple_to_weight) + # print(knutti_weight.couple_to_weight) weight = knutti_weight.couple_to_weight[('CNRM-CM5', 'CCLM4-8-17')] self.assertFalse(np.isnan(weight)) -- GitLab