diff --git a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
index 4dfb2e52e232d59de897549da28599f98dfe6ff6..561387abc63859aae890198e9124f600c13c5031 100644
--- a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
+++ b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
@@ -26,6 +26,7 @@ class VisualizerForProjectionEnsemble(object):
gcm_to_year_min_and_year_max=None,
):
self.altitudes_list = altitudes_list
+ self.temporal_covariate_for_fit = temporal_covariate_for_fit
self.scenario = scenario
self.gcm_rcm_couples = gcm_rcm_couples
self.massif_names = massif_names
diff --git a/projects/projected_swe/abstract_weight_computer.py b/projects/projected_swe/abstract_weight_computer.py
index dbff2a20edd74a35fc39461a7e234d9925eeee18..669419f1f0b75141fd4ee9790fb3e70002e794b4 100644
--- a/projects/projected_swe/abstract_weight_computer.py
+++ b/projects/projected_swe/abstract_weight_computer.py
@@ -53,7 +53,7 @@ class AbstractWeightComputer(object):
# Save csv
save_to_filepath(df, self.visualizer.gcm_rcm_couples, self.visualizer.altitudes_list,
self.year_min, self.year_max,
- self.visualizer.scenario)
+ self.visualizer.scenario, type(self))
def compute_weight(self, gcm_rcm_couples_subset, gcm_rcm_couple_missing):
# Initial the dictionary
@@ -70,13 +70,11 @@ class AbstractWeightComputer(object):
reference_studies = ensemble_fit.gcm_rcm_couple_to_visualizer[gcm_rcm_couple_missing].studies
# Computer the nllh for each "station", i.e. each altitude and massif_name
for altitude, reference_study in reference_studies.altitude_to_study.items():
- for massif_name, maxima in reference_study.massif_name_to_annual_maxima.items():
+ for massif_name in reference_study.study_massif_names:
assert all([self.year_min <= year <= self.year_max for year in reference_study.ordered_years])
gcm_rcm_couple_to_local_nllh = self.compute_gcm_rcm_couple_to_local_nllh(ensemble_fit,
gcm_rcm_couples_subset,
massif_name,
- maxima,
- altitude,
reference_study)
if gcm_rcm_couple_to_local_nllh is not None:
for c, nllh in gcm_rcm_couple_to_local_nllh.items():
@@ -89,6 +87,22 @@ class AbstractWeightComputer(object):
gcm_rcm_couple_to_weight[gcm_rcm_couple_missing] = np.nan
return gcm_rcm_couple_to_weight
- def compute_gcm_rcm_couple_to_local_nllh(self, ensemble_fit, gcm_rcm_couples_subset, massif_name, maxima,
- altitude, reference_study):
+ def compute_gcm_rcm_couple_to_local_nllh(self, ensemble_fit, gcm_rcm_couples_subset, massif_name, reference_study):
+ # Check that all the gcm_rcm_couple have a model for this massif_name
+ if self.condition_to_compute_nllh(ensemble_fit, massif_name, self.visualizer):
+ print(ensemble_fit.altitudes, massif_name)
+ nllh_list = []
+ for gcm_rcm_couple in gcm_rcm_couples_subset:
+ nllh = self.compute_nllh_local(ensemble_fit, gcm_rcm_couple, massif_name, reference_study)
+ nllh_list.append(nllh)
+
+ if all([not np.isinf(nllh) for nllh in nllh_list]):
+ return dict(zip(gcm_rcm_couples_subset, nllh_list))
+
+ def condition_to_compute_nllh(self, ensemble_fit, massif_name, visualizer):
+ return all(
+ [massif_name in ensemble_fit.gcm_rcm_couple_to_visualizer[c].massif_name_to_one_fold_fit for c in
+ visualizer.gcm_rcm_couples])
+
+ def compute_nllh_local(self, ensemble_fit, gcm_rcm_couple, massif_name, reference_study):
raise NotImplementedError
diff --git a/projects/projected_swe/knutti_weight_computer.py b/projects/projected_swe/knutti_weight_computer.py
new file mode 100644
index 0000000000000000000000000000000000000000..41b43b1331109aa963231273a9fef39ad18cf585
--- /dev/null
+++ b/projects/projected_swe/knutti_weight_computer.py
@@ -0,0 +1,25 @@
+import numpy as np
+
+from projects.projected_swe.abstract_weight_computer import AbstractWeightComputer
+from projects.projected_swe.non_stationary_weight_computer import NllhWeightComputer
+
+
+class KnuttiWeightComputer(AbstractWeightComputer):
+
+ def __init__(self, *args, sigma_D, **kwargs):
+ super().__init__(*args, **kwargs)
+ self.sigma_D = sigma_D
+
+ def compute_nllh_local(self, ensemble_fit, gcm_rcm_couple, massif_name, reference_study):
+ mean_maxima_reference = reference_study.massif_name_to_annual_maxima[massif_name].mean()
+ studies = ensemble_fit.gcm_rcm_couple_to_visualizer[gcm_rcm_couple].studies
+ study = studies.altitude_to_study[reference_study.altitude]
+ mean_maxima_couple = study.massif_name_to_annual_maxima[massif_name].mean()
+ ratio = (mean_maxima_couple - mean_maxima_reference) / self.sigma_D
+ proba = np.exp(-np.power(ratio, 2))
+ return -np.log(proba)
+
+
+
+class MixedWeightComputer(NllhWeightComputer, KnuttiWeightComputer):
+ pass
\ No newline at end of file
diff --git a/projects/projected_swe/main_weight_computation.py b/projects/projected_swe/main_weight_computation.py
index 579346bf97351a44b00bb2a8d4093a64acba8067..8a6b87513be25d40e0a87eeb6f56280c42c9ad08 100644
--- a/projects/projected_swe/main_weight_computation.py
+++ b/projects/projected_swe/main_weight_computation.py
@@ -26,6 +26,7 @@ from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fi
from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
from extreme_trend.one_fold_fit.altitude_group import altitudes_for_groups
from projects.projected_swe.abstract_weight_computer import AbstractWeightComputer
+from projects.projected_swe.knutti_weight_computer import KnuttiWeightComputer
from projects.projected_swe.non_stationary_weight_computer import NllhWeightComputer
from projects.projected_swe.utils import load_gcm_rcm_couple_to_weight
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
@@ -49,7 +50,7 @@ def main_weight_computation():
gcm_rcm_couples = get_gcm_rcm_couples(scenario)
year_min = 1982
year_max = 2019
- weight_computer_class = NllhWeightComputer
+ weight_computer_class = [NllhWeightComputer, KnuttiWeightComputer][1]
fast = True
if fast is None:
@@ -74,7 +75,9 @@ def main_weight_computation():
gcm_to_year_min_and_year_max={c[0]: (year_min, year_max) for c in gcm_rcm_couples},
)
- weight_computer = weight_computer_class(visualizer, scm_study_class, year_min, year_max) # type:AbstractWeightComputer
+
+ weight_computer = weight_computer_class(visualizer, scm_study_class, year_min, year_max,
+ sigma_D=10) # type:AbstractWeightComputer
weight_computer.compute_weights_and_save_them()
end = time.time()
@@ -84,5 +87,6 @@ def main_weight_computation():
if __name__ == '__main__':
main_weight_computation()
- d = load_gcm_rcm_couple_to_weight(['sd', 'sdf'], [23], 1982, 2019, AdamontScenario.rcp85_extended)
- print(d)
+ # d = load_gcm_rcm_couple_to_weight(['sd', 'sdf'], [23], 1982, 2019, AdamontScenario.rcp85_extended,
+ # weight_class=NllhWeightComputer, gcm_rcm_couple_missing=None)
+ # print(d)
diff --git a/projects/projected_swe/non_stationary_weight_computer.py b/projects/projected_swe/non_stationary_weight_computer.py
index c7ffaec6f164ac819460f9947aeaa47197de9916..dd305395470c24252ad5a85cc65d741e2d8643f0 100644
--- a/projects/projected_swe/non_stationary_weight_computer.py
+++ b/projects/projected_swe/non_stationary_weight_computer.py
@@ -1,6 +1,7 @@
import numpy as np
from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import compute_nllh
+from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
from projects.projected_swe.abstract_weight_computer import AbstractWeightComputer
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
TimeTemporalCovariate
@@ -8,32 +9,19 @@ from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_
class NllhWeightComputer(AbstractWeightComputer):
- def compute_gcm_rcm_couple_to_local_nllh(self, ensemble_fit, gcm_rcm_couples_subset, massif_name, maxima,
- altitude, reference_study):
- assert ensemble_fit.temporal_covariate_for_fit is TimeTemporalCovariate, \
+ def __init__(self, *args, **kwargs):
+ super().__init__(*args, **kwargs)
+ assert self.visualizer.temporal_covariate_for_fit is TimeTemporalCovariate, \
'for temperature covariate you should transform the coordinates'
- gcm_rcm_couple_to_local_nllh = {}
- # Check that all the gcm_rcm_couple have a model for this massif_name
- if self.condition_to_compute_nllh(ensemble_fit, massif_name, self.visualizer):
- print(ensemble_fit.altitudes, massif_name)
- coordinates = [np.array([altitude, year]) for year in reference_study.ordered_years]
- nllh_list = []
- for gcm_rcm_couple in gcm_rcm_couples_subset:
- best_function_from_fit = self.get_function_from_fit(ensemble_fit, massif_name, gcm_rcm_couple)
- # It is normal that it could crash, because some models where fitted with data smaller than
- # the data used to compute the nllh
- nllh = compute_nllh(coordinates, maxima, best_function_from_fit,
- maximum_from_obs=False, assertion_for_inf=False)
- nllh_list.append(nllh)
-
- if all([not np.isinf(nllh) for nllh in nllh_list]):
- return dict(zip(gcm_rcm_couples_subset, nllh_list))
-
- def condition_to_compute_nllh(self, ensemble_fit, massif_name, visualizer):
- return all(
- [massif_name in ensemble_fit.gcm_rcm_couple_to_visualizer[c].massif_name_to_one_fold_fit for c in
- visualizer.gcm_rcm_couples])
+ def compute_nllh_local(self, ensemble_fit, gcm_rcm_couple, massif_name, reference_study):
+ coordinates = [np.array([reference_study.altitude, year]) for year in reference_study.ordered_years]
+ maxima = reference_study.massif_name_to_annual_maxima[massif_name]
+ best_function_from_fit = self.get_function_from_fit(ensemble_fit, massif_name, gcm_rcm_couple)
+ # It is normal that it could crash, because some models where fitted with data smaller than
+ # the data used to compute the nllh, therefore we use assertion_for_inf=False below
+ return compute_nllh(coordinates, maxima, best_function_from_fit,
+ maximum_from_obs=False, assertion_for_inf=False)
def get_function_from_fit(self, ensemble_fit, massif_name, gcm_rcm_couple):
visualizer = ensemble_fit.gcm_rcm_couple_to_visualizer[gcm_rcm_couple]
diff --git a/projects/projected_swe/utils.py b/projects/projected_swe/utils.py
index 2c04c17d8065f5373b52c9f35874a61aae86607d..0a7202e0f2fb072dbb6630e5e4f8b0769f572e48 100644
--- a/projects/projected_swe/utils.py
+++ b/projects/projected_swe/utils.py
@@ -23,6 +23,7 @@ from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
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 altitudes_for_groups
+from root_utils import get_display_name_from_object_type
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
TimeTemporalCovariate
from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
@@ -33,23 +34,24 @@ WEIGHT_COLUMN_NAME = "all weights"
WEIGHT_FOLDER = "ensemble weight"
-def get_csv_filepath(gcm_rcm_couples, altitudes_list, year_min, year_max, scenario):
+def get_csv_filepath(gcm_rcm_couples, altitudes_list, year_min, year_max, scenario, weight_computer_class):
nb_gcm_rcm_couples = len(gcm_rcm_couples)
nb_altitudes_list = len(altitudes_list)
ensemble_folder_path = op.join(DATA_PATH, WEIGHT_FOLDER)
if not op.exists(ensemble_folder_path):
os.makedirs(ensemble_folder_path, exist_ok=True)
scenario_str = scenario_to_str(scenario)
- csv_filename = "weights_{}_{}_{}_{}_{}.csv" \
- .format(nb_gcm_rcm_couples, nb_altitudes_list, year_min, year_max, scenario_str)
+ class_name = get_display_name_from_object_type(weight_computer_class)
+ csv_filename = "weights_{}_{}_{}_{}_{}_{}.csv" \
+ .format(nb_gcm_rcm_couples, nb_altitudes_list, year_min, year_max, scenario_str, class_name)
weight_csv_filepath = op.join(ensemble_folder_path, csv_filename)
return weight_csv_filepath
def save_to_filepath(df, gcm_rcm_couples, altitudes_list,
year_min, year_max,
- scenario):
- filepath = get_csv_filepath(gcm_rcm_couples, altitudes_list, year_min, year_max, scenario)
+ scenario, weight_computer_class):
+ filepath = get_csv_filepath(gcm_rcm_couples, altitudes_list, year_min, year_max, scenario, weight_computer_class)
df = df.round(decimals=3)
df.index = [gcm_rcm_couple_to_str(i) for i in df.index]
df.columns = [gcm_rcm_couple_to_str(i) if j > 0 else i for j, i in enumerate(df.columns)]
@@ -60,8 +62,8 @@ def save_to_filepath(df, gcm_rcm_couples, altitudes_list,
def load_gcm_rcm_couple_to_weight(gcm_rcm_couples, altitudes_list, year_min, year_max, scenario,
- gcm_rcm_couple_missing=None):
- filepath = get_csv_filepath(gcm_rcm_couples, altitudes_list, year_min, year_max, scenario)
+ weight_class, gcm_rcm_couple_missing):
+ filepath = get_csv_filepath(gcm_rcm_couples, altitudes_list, year_min, year_max, scenario, weight_class)
df = pd.read_csv(filepath, index_col=0)
df.index = [str_to_gcm_rcm_couple(i) for i in df.index]
df.columns = [str_to_gcm_rcm_couple(i) if j > 0 else i for j, i in enumerate(df.columns)]