From ec53cf8b59c9803ab864d0222762c6d54f7f50a2 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 22 Oct 2019 11:49:35 +0200
Subject: [PATCH] [TEMPORAL LINEAR MARGIN][BAYESIAN ESTIMATION] modify
fevd_fixed prior to ensure that prior is always applied on the shape
parameter. add multiprocessing for eurocode return level uncertainty. add
option to specify a last year for eurocode uncertainty computation.
---
.../eurocode_return_level_uncertainties.py | 63 ++++++++++---------
.../eurocode_data/main_eurocode_drawing.py | 6 +-
experiment/eurocode_data/utils.py | 8 ++-
.../study_visualization/study_visualizer.py | 34 +++++-----
extreme_fit/distribution/gev/fevd_fixed.R | 7 ++-
5 files changed, 63 insertions(+), 55 deletions(-)
diff --git a/experiment/eurocode_data/eurocode_return_level_uncertainties.py b/experiment/eurocode_data/eurocode_return_level_uncertainties.py
index 750adf64..701da19a 100644
--- a/experiment/eurocode_data/eurocode_return_level_uncertainties.py
+++ b/experiment/eurocode_data/eurocode_return_level_uncertainties.py
@@ -1,10 +1,9 @@
from typing import List
-
import numpy as np
from cached_property import cached_property
-from experiment.eurocode_data.utils import EUROCODE_QUANTILE
+from experiment.eurocode_data.utils import EUROCODE_QUANTILE, YEAR_OF_INTEREST_FOR_RETURN_LEVEL
from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
fitted_linear_margin_estimator
from extreme_fit.distribution.gev.gev_params import GevParams
@@ -17,11 +16,42 @@ from extreme_fit.model.margin_model.margin_function.linear_margin_function impor
from extreme_fit.model.result_from_model_fit.result_from_extremes import ResultFromExtremes
+def compute_eurocode_level_uncertainty(last_year_for_the_data, smooth_maxima_x_y, model_class):
+ years, smooth_maxima = smooth_maxima_x_y
+ idx = years.index(last_year_for_the_data) + 1
+ years, smooth_maxima = years[:idx], smooth_maxima[:idx]
+ return EurocodeLevelUncertaintyFromExtremes.from_maxima_years_model_class(smooth_maxima, years, model_class)
+
+
+class EurocodeLevelUncertaintyFromExtremes(object):
+ YEAR_OF_INTEREST = 2017
+
+ def __init__(self, posterior_mean, poster_uncertainty_interval):
+ self.posterior_mean = posterior_mean
+ self.poster_uncertainty_interval = poster_uncertainty_interval
+
+ @classmethod
+ def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator):
+ extractor = ExtractEurocodeReturnLevelFromExtremes(estimator_extremes, cls.YEAR_OF_INTEREST)
+ return cls(extractor.posterior_mean_eurocode_return_level_for_the_year_of_interest,
+ extractor.posterior_eurocode_return_level_uncertainty_interval_for_the_year_of_interest)
+
+ @classmethod
+ def from_maxima_years_model_class(cls, maxima, years, model_class):
+ # Load coordinates and dataset
+ coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
+ # Fitted estimator
+ fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=1958,
+ fit_method=AbstractTemporalLinearMarginModel.EXTREMES_FEVD_BAYESIAN_FIT_METHOD_STR)
+ # Load object from result from extremes
+ return cls.from_estimator_extremes(fitted_estimator)
+
+
class ExtractEurocodeReturnLevelFromExtremes(object):
- def __init__(self, estimator: LinearMarginEstimator, year_of_interest: int = 2017):
+ def __init__(self, estimator: LinearMarginEstimator, year_of_interest: int = YEAR_OF_INTEREST_FOR_RETURN_LEVEL):
self.estimator = estimator
- self.result_from_fit = self.estimator.result_from_model_fit # type: ResultFromExtremes
+ self.result_from_fit = self.estimator.result_from_model_fit # type: ResultFromExtremes
self.year_of_interest = year_of_interest
@property
@@ -53,28 +83,3 @@ class ExtractEurocodeReturnLevelFromExtremes(object):
bottom_and_upper_quantile = (0.250, 0.975)
return [np.quantile(self.posterior_eurocode_return_level_samples_for_year_of_interest, q=q)
for q in bottom_and_upper_quantile]
-
-
-class EurocodeLevelUncertaintyFromExtremes(object):
-
- YEAR_OF_INTEREST = 2017
-
- def __init__(self, posterior_mean, poster_uncertainty_interval):
- self.posterior_mean = posterior_mean
- self.poster_uncertainty_interval = poster_uncertainty_interval
-
- @classmethod
- def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator):
- extractor = ExtractEurocodeReturnLevelFromExtremes(estimator_extremes, cls.YEAR_OF_INTEREST)
- return cls(extractor.posterior_mean_eurocode_return_level_for_the_year_of_interest,
- extractor.posterior_eurocode_return_level_uncertainty_interval_for_the_year_of_interest)
-
- @classmethod
- def from_maxima_years_model_class(cls, maxima, years, model_class):
- # Load coordinates and dataset
- coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
- # Fitted estimator
- fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=1958,
- fit_method=AbstractTemporalLinearMarginModel.EXTREMES_FEVD_BAYESIAN_FIT_METHOD_STR)
- # Load object from result from extremes
- return cls.from_estimator_extremes(fitted_estimator)
\ No newline at end of file
diff --git a/experiment/eurocode_data/main_eurocode_drawing.py b/experiment/eurocode_data/main_eurocode_drawing.py
index 569cb9a8..f5561017 100644
--- a/experiment/eurocode_data/main_eurocode_drawing.py
+++ b/experiment/eurocode_data/main_eurocode_drawing.py
@@ -2,7 +2,7 @@ from collections import OrderedDict
from experiment.eurocode_data.eurocode_visualizer import plot_model_name_to_dep_to_ordered_return_level_uncertainties
from experiment.eurocode_data.massif_name_to_departement import DEPARTEMENT_TYPES
-from experiment.eurocode_data.utils import EUROCODE_ALTITUDES
+from experiment.eurocode_data.utils import EUROCODE_ALTITUDES, LAST_YEAR_FOR_EUROCODE
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_hypercube_visualizer import \
AltitudeHypercubeVisualizer
@@ -31,7 +31,7 @@ def dep_to_ordered_return_level_uncertainties(model_class, last_year_for_the_dat
assert isinstance(altitude_visualizer.tuple_to_study_visualizer, OrderedDict)
dep_to_ordered_return_level_uncertainty = {dep: [] for dep in DEPARTEMENT_TYPES}
for visualizer in altitude_visualizer.tuple_to_study_visualizer.values():
- dep_to_return_level_uncertainty = visualizer.dep_class_to_eurocode_level_uncertainty(model_class)
+ dep_to_return_level_uncertainty = visualizer.dep_class_to_eurocode_level_uncertainty(model_class, last_year_for_the_data)
for dep, return_level_uncertainty in dep_to_return_level_uncertainty.items():
dep_to_ordered_return_level_uncertainty[dep].append(return_level_uncertainty)
@@ -40,7 +40,7 @@ def dep_to_ordered_return_level_uncertainties(model_class, last_year_for_the_dat
def main_drawing():
model_class_and_last_year = [
- (StationaryStationModel, 1991),
+ (StationaryStationModel, LAST_YEAR_FOR_EUROCODE),
(StationaryStationModel, 2017),
(NonStationaryLocationAndScaleModel, 2017),
][:1]
diff --git a/experiment/eurocode_data/utils.py b/experiment/eurocode_data/utils.py
index ef3b88eb..f293fdcd 100644
--- a/experiment/eurocode_data/utils.py
+++ b/experiment/eurocode_data/utils.py
@@ -2,4 +2,10 @@
# Eurocode quantile correspond to a 50 year return period
EUROCODE_QUANTILE = 0.98
# Altitudes (between low and mid altitudes) < 2000m
-EUROCODE_ALTITUDES = [900, 1200, 1500, 1800][:2]
\ No newline at end of file
+EUROCODE_ALTITUDES = [900, 1200, 1500, 1800][:2]
+# Last year taken into account for the Eurocode
+# Date of publication was 2014, therefore the winter 2013/2014 could not have been measured
+# Therefore, the winter 2012/2013 was the last one. Thus, 2012 is the last year for the Eurocode
+LAST_YEAR_FOR_EUROCODE = 2012
+# Year of interest for the EUROCODE
+YEAR_OF_INTEREST_FOR_RETURN_LEVEL = 2017
\ No newline at end of file
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
index 4b361e43..7f624490 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
@@ -1,6 +1,7 @@
import os
import os.path as op
from collections import OrderedDict
+from multiprocessing.pool import Pool
from random import sample, seed
from typing import Dict
@@ -11,7 +12,7 @@ import pandas as pd
import seaborn as sns
from experiment.eurocode_data.eurocode_return_level_uncertainties import ExtractEurocodeReturnLevelFromExtremes, \
- EurocodeLevelUncertaintyFromExtremes
+ EurocodeLevelUncertaintyFromExtremes, compute_eurocode_level_uncertainty
from experiment.eurocode_data.massif_name_to_departement import dep_class_to_massif_names
from experiment.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
from experiment.trend_analysis.abstract_score import MeanScore, AbstractTrendScore
@@ -47,7 +48,7 @@ from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima
from test.test_utils import load_test_max_stable_models
from root_utils import get_display_name_from_object_type, VERSION_TIME, float_to_str_with_only_some_significant_digits, \
- cached_property
+ cached_property, NB_CORES
BLOCK_MAXIMA_DISPLAY_NAME = 'block maxima '
@@ -354,36 +355,29 @@ class StudyVisualizer(VisualizationParameters):
start_year, stop_year = self.study.start_year_and_stop_year
return list(range(start_year, stop_year))
- def massif_name_to_eurocode_level_uncertainty(self, model_class) -> Dict[str, EurocodeLevelUncertaintyFromExtremes]:
- # todo: add multiprocessing
- massif_name_to_eurocode_return_level_uncertainty = {}
- for massif_id, massif_name in enumerate(self.study.study_massif_names):
- print(massif_name)
- years, smooth_maxima = self.smooth_maxima_x_y(massif_id)
- res = EurocodeLevelUncertaintyFromExtremes.from_maxima_years_model_class(smooth_maxima, years, model_class)
- massif_name_to_eurocode_return_level_uncertainty[massif_name] = res
+ def massif_name_to_eurocode_level_uncertainty(self, model_class, last_year_for_the_data) -> Dict[str, EurocodeLevelUncertaintyFromExtremes]:
+ arguments = [[last_year_for_the_data, self.smooth_maxima_x_y(massif_id), model_class] for massif_id, _ in enumerate(self.study.study_massif_names)]
+ if self.multiprocessing:
+ with Pool(NB_CORES) as p:
+ res = p.starmap(compute_eurocode_level_uncertainty, arguments)
+ else:
+ res = [compute_eurocode_level_uncertainty(*argument) for argument in arguments]
+ massif_name_to_eurocode_return_level_uncertainty = OrderedDict(zip(self.study.study_massif_names, res))
return massif_name_to_eurocode_return_level_uncertainty
- def dep_class_to_eurocode_level_uncertainty(self, model_class):
+ def dep_class_to_eurocode_level_uncertainty(self, model_class, last_year_for_the_data):
""" Take the max with respect to the posterior mean for each departement """
- massif_name_to_eurocode_level_uncertainty = self.massif_name_to_eurocode_level_uncertainty(model_class)
- print(massif_name_to_eurocode_level_uncertainty)
+ massif_name_to_eurocode_level_uncertainty = self.massif_name_to_eurocode_level_uncertainty(model_class,
+ last_year_for_the_data)
dep_class_to_eurocode_level_uncertainty = {}
for dep_class, massif_names in dep_class_to_massif_names.items():
# Filter the couple of interest
couples = [(k, v) for k, v in massif_name_to_eurocode_level_uncertainty.items() if k in massif_names]
assert len(couples) > 0
- # Find the massif name with the maximum posterior mean
- for c in couples:
- print(c)
- print(c[1].posterior_mean)
massif_name, eurocode_level_uncertainty = max(couples, key=lambda c: c[1].posterior_mean)
dep_class_to_eurocode_level_uncertainty[dep_class] = eurocode_level_uncertainty
return dep_class_to_eurocode_level_uncertainty
-
-
-
@cached_property
def massif_name_to_detailed_scores(self):
"""
diff --git a/extreme_fit/distribution/gev/fevd_fixed.R b/extreme_fit/distribution/gev/fevd_fixed.R
index 30ee10bc..db0c6742 100644
--- a/extreme_fit/distribution/gev/fevd_fixed.R
+++ b/extreme_fit/distribution/gev/fevd_fixed.R
@@ -5,8 +5,11 @@ library(SpatialExtremes)
# Define a custom Prior function
fevdPriorCustom <- function (theta, q, p, log = FALSE){
- # Select the shape parameter (which is the last parameter)
- shape = theta[length(theta)]
+ # Select the shape parameter (which is always the last parameter either for the stationary or the non-stationary case)
+ print(attributes(theta))
+ theta_names = attr(theta, 'names')
+ shape_idx = match('shape', theta_names)
+ shape = theta[shape_idx]
# + 0.5 enables to shift the Beta law in the interval [-0.5, 0.5]
res = dbeta(shape + 0.5, q, p, log = TRUE)
return(res)
--
GitLab