From 83beed61aa474124efc5c022f02967fd07394e61 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 19 Feb 2021 16:38:07 +0100
Subject: [PATCH] [projections] add a constraint to have physically plausible
models for the fitted models.
---
...es_visualizer_for_non_stationary_models.py | 57 ++++++++++++-------
.../one_fold_analysis/one_fold_fit.py | 25 +++++---
.../preliminary_analysis.py | 8 +--
.../ensemble_fit/abstract_ensemble_fit.py | 2 +
.../visualizer_for_projection_ensemble.py | 6 +-
5 files changed, 65 insertions(+), 33 deletions(-)
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index f49d5ac3..d4292774 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -1,5 +1,6 @@
from collections import Counter
from math import ceil, floor
+from multiprocessing import Pool
from typing import List, Dict
import matplotlib
@@ -27,6 +28,7 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_gr
get_altitude_group_from_altitudes, HighAltitudeGroup, VeyHighAltitudeGroup, MidAltitudeGroup
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import \
OneFoldFit
+from root_utils import NB_CORES
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
@@ -40,11 +42,12 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
fit_method=MarginFitMethod.extremes_fevd_mle,
temporal_covariate_for_fit=None,
display_only_model_that_pass_anderson_test=True,
- confidence_interval_based_on_delta_method=False
+ confidence_interval_based_on_delta_method=False,
+ remove_physically_implausible_models=False,
):
super().__init__(studies.study, show=show, save_to_file=not show)
self.studies = studies
- self.non_stationary_models = model_classes
+ self.model_classes = model_classes
self.fit_method = fit_method
self.temporal_covariate_for_fit = temporal_covariate_for_fit
self.display_only_model_that_pass_test = display_only_model_that_pass_anderson_test
@@ -52,23 +55,20 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self.massif_name_to_massif_id = {m: i for i, m in enumerate(self.massif_names)}
self.altitude_group = get_altitude_group_from_altitudes(self.studies.altitudes)
self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
+ self.remove_physically_implausible_models = remove_physically_implausible_models
# Load one fold fit
self.massif_name_to_massif_altitudes = {}
- self._massif_name_to_one_fold_fit = {}
- for massif_name in self.massif_names:
- # Load valid massif altitudes
- massif_altitudes = self.get_massif_altitudes(massif_name)
- if self.load_condition(massif_altitudes):
- # Save the massif altitudes only for those who pass the condition
- self.massif_name_to_massif_altitudes[massif_name] = massif_altitudes
- # Load dataset
- dataset = studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
- old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method,
- self.temporal_covariate_for_fit,
- type(self.altitude_group),
- self.display_only_model_that_pass_test,
- self.confidence_interval_based_on_delta_method)
- self._massif_name_to_one_fold_fit[massif_name] = old_fold_fit
+
+ # Multiprocess
+ multiprocessing = False
+ if multiprocessing:
+ with Pool(NB_CORES) as p:
+ one_fold_fit_list = p.map(self.fit_one_fold, self.massif_names)
+ else:
+ one_fold_fit_list = [self.fit_one_fold(massif_name) for massif_name in self.massif_names]
+ self._massif_name_to_one_fold_fit = {m: o for m, o in zip(self.massif_names, one_fold_fit_list) if
+ o is not None}
+
# Print number of massif without any validated fit
massifs_without_any_validated_fit = [massif_name
for massif_name, old_fold_fit in self._massif_name_to_one_fold_fit.items()
@@ -79,6 +79,24 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self._method_name_and_order_to_max_abs = {}
self._max_abs_for_shape = None
+ def fit_one_fold(self, massif_name):
+ # Load valid massif altitudes
+ massif_altitudes = self.get_massif_altitudes(massif_name)
+ if self.load_condition(massif_altitudes):
+ # 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)
+ old_fold_fit = OneFoldFit(massif_name, dataset, self.model_classes, self.fit_method,
+ self.temporal_covariate_for_fit,
+ type(self.altitude_group),
+ self.display_only_model_that_pass_test,
+ self.confidence_interval_based_on_delta_method,
+ self.remove_physically_implausible_models)
+ return old_fold_fit
+ else:
+ return None
+
moment_names = ['moment', 'changes_of_moment', 'relative_changes_of_moment'][:]
orders = [1, 2, None][2:]
@@ -106,8 +124,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
@property
def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
return {massif_name: old_fold_fit for massif_name, old_fold_fit in self._massif_name_to_one_fold_fit.items()
- if not self.display_only_model_that_pass_test
- or old_fold_fit.has_at_least_one_valid_model}
+ if old_fold_fit.has_at_least_one_valid_model}
def plot_moments(self):
for method_name in self.moment_names[:2]:
@@ -452,7 +469,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
idx = 0 if one_fold.change_in_return_level_for_reference_altitude < 0 else 2
nbs[idx] += 1
if with_significance and one_fold.is_significant:
- nbs[idx + 1] += 1
+ nbs[idx + 1] += 1
percents = 100 * nbs / nb_valid_massif_names
return [nb_valid_massif_names] + list(percents)
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 93ce1c0b..e37365f9 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
@@ -53,7 +53,9 @@ class OneFoldFit(object):
altitude_class=DefaultAltitudeGroup,
only_models_that_pass_goodness_of_fit_test=True,
confidence_interval_based_on_delta_method=False,
+ remove_physically_implausible_models=False,
):
+ self.remove_physically_implausible_models = remove_physically_implausible_models
self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test
self.altitude_group = altitude_class()
@@ -142,14 +144,21 @@ class OneFoldFit(object):
@cached_property
def sorted_estimators(self):
estimators = list(self.model_class_to_estimator.values())
- # print(self.massif_name)
- # print(self.altitude_group)
- # for estimator in estimators:
- # print(estimator.margin_model)
- # print(estimator.aic())
+ print(self.massif_name, self.altitude_group)
+ if self.remove_physically_implausible_models:
+ estimators = [e for e in estimators if -0.5 < self._compute_shape_for_reference_altitude(e) < 0.5]
+ 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())
return sorted_estimators
+ def _compute_shape_for_reference_altitude(self, estimator):
+ coordinate = np.array([self.altitude_plot, self.last_year])
+ gev_params = estimator.function_from_fit.get_params(coordinate, is_transformed=False)
+ shape = gev_params.shape
+ return shape
+
@cached_property
def sorted_estimators_with_stationary(self):
if self.only_models_that_pass_goodness_of_fit_test:
@@ -159,7 +168,8 @@ class OneFoldFit(object):
# and self.sensitivity_of_fit_test_last_years(e)
]
else:
- assert len(self.sorted_estimators) == len(self.models_classes)
+ if not self.remove_physically_implausible_models:
+ assert len(self.sorted_estimators) == len(self.models_classes)
return self.sorted_estimators
@property
@@ -180,7 +190,8 @@ class OneFoldFit(object):
best_estimator = self.sorted_estimators_with_stationary[0]
return best_estimator
else:
- raise ValueError('This should not happen')
+ raise ValueError('This object should not have been called because '
+ 'has_at_least_one_valid_model={}'.format(self.has_at_least_one_valid_model))
@property
def best_margin_model(self):
diff --git a/projects/altitude_spatial_model/preliminary_analysis.py b/projects/altitude_spatial_model/preliminary_analysis.py
index 5148971a..dcc35592 100644
--- a/projects/altitude_spatial_model/preliminary_analysis.py
+++ b/projects/altitude_spatial_model/preliminary_analysis.py
@@ -27,8 +27,6 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
**kwargs_study):
super().__init__(study_class, altitudes, spatial_transformation_class, temporal_transformation_class,
**kwargs_study)
- # self.altitudes_for_temporal_hypothesis = [min(self.altitudes), 2100, max(self.altitudes)]
- self.altitudes_for_temporal_hypothesis = [600, 1500, 2400, 3300]
def plot_gev_params_against_altitude(self):
legend = False
@@ -332,8 +330,10 @@ def main_paper2():
def main_paper3():
altitudes = list(chain.from_iterable(altitudes_for_groups))
# altitudes = [1200, 1500, 1800]
- for scenario in rcp_scenarios[:]:
- for gcm_rcm_couple in get_gcm_rcm_couples(scenario):
+ for scenario in rcp_scenarios[2:]:
+ gcm_rcm_couples = get_gcm_rcm_couples(scenario)
+ gcm_rcm_couples =[('CNRM-CM5', 'CCLM4-8-17')]
+ for gcm_rcm_couple in gcm_rcm_couples:
visualizer = PointwiseGevStudyVisualizer(AdamontSnowfall, altitudes=altitudes, scenario=scenario,
gcm_rcm_couple=gcm_rcm_couple)
visualizer.plot_gev_params_against_altitude()
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py
index 6a1d7b01..ad85be10 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py
@@ -13,7 +13,9 @@ class AbstractEnsembleFit(object):
temporal_covariate_for_fit=None,
only_models_that_pass_goodness_of_fit_test=True,
confidence_interval_based_on_delta_method=False,
+ remove_physically_implausible_models=False,
):
+ self.remove_physically_implausible_models = remove_physically_implausible_models
self.massif_names = massif_names
self.models_classes = models_classes
self.gcm_rcm_couple_to_altitude_studies = gcm_rcm_couple_to_altitude_studies
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 7e9185cd..e4f0c23e 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
@@ -44,7 +44,8 @@ class MetaVisualizerForProjectionEnsemble(object):
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
+ confidence_interval_based_on_delta_method=False,
+ remove_physically_implausible_models=False,
):
self.gcm_rcm_couples = gcm_rcm_couples
self.massif_names = massif_names
@@ -69,7 +70,8 @@ class MetaVisualizerForProjectionEnsemble(object):
ensemble_fit = ensemble_fit_class(massif_names, gcm_rcm_couple_to_studies, model_classes,
fit_method, temporal_covariate_for_fit,
display_only_model_that_pass_gof_test,
- confidence_interval_based_on_delta_method)
+ 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
--
GitLab