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 0000000000000000000000000000000000000000..86688952ddf483ed519009e185cc9203c2372109
--- /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 0000000000000000000000000000000000000000..4eb1267d22ded38da900f5c392708612a41be4fc
--- /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 fdde441a20cf4d7be1325aaa29db33301c7b7496..688444111acbae1f15bc0d5862e82c13be2ecad8 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 0000000000000000000000000000000000000000..d2cd76df5290abec7182e582d7e405355b93d35a
--- /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 1e086777b2fa1bac1bf31e8ea74fc2ece19cef3b..cc38a3eaab72c736e99664e8a848cfd3e7ca0f69 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 4f1a230a107bbcd5eec7742ca6cdd08b0b335a5f..9a6ae22590ad7d375943b627c414cd5dd2d48b46 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 0000000000000000000000000000000000000000..e6af3be6645fc684c6c1bc33dde3ec213748dc53
--- /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 6b10b921cdb052b7841de0233ec1a8581bb08bf1..95f0a251381df66fb5a86f7846c2c401b3593e27 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))