diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
index 43e3d2eda4e46801f2205187a48a28b3219ba2f6..19a01c235df96a65fed835df420d4a23d5d43821 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
@@ -69,10 +69,10 @@ def max_stable_process_vizu_compare_gaume_study(altitude=1800, nb_days=1):
def normal_visualization(temporal_non_stationarity=False):
- save_to_file = True
+ save_to_file = False
only_first_one = True
# for study_class in SCM_STUDIES[:1]:
- for study_class in [CrocusDepth, SafranSnowfall, SafranRainfall, SafranTemperature][:]:
+ for study_class in [CrocusDepth, SafranSnowfall, SafranRainfall, SafranTemperature][1:2]:
for study in study_iterator(study_class, only_first_one=only_first_one):
study_visualizer = StudyVisualizer(study, save_to_file=save_to_file, temporal_non_stationarity=temporal_non_stationarity)
# study_visualizer.visualize_independent_margin_fits(threshold=[None, 20, 40, 60][0])
@@ -96,8 +96,8 @@ def complete_analysis(only_first_one=False):
if __name__ == '__main__':
- annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True, altitude=2100)
- # normal_visualization(temporal_non_stationarity=True)
+ # annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True, altitude=2100)
+ normal_visualization(temporal_non_stationarity=True)
# max_stable_process_vizu_compare_gaume_study(altitude=1800, nb_days=1)
# extended_visualization()
# complete_analysis()
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
index d1dad192052b48a3c8ef9f4b057f2b2d6ce948e2..a57c41b2bb04a0032a89bae148a4d64f34576ea1 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
@@ -14,7 +14,8 @@ from experiment.utils import average_smoothing_with_sliding_window
from extreme_estimator.estimator.full_estimator.abstract_full_estimator import \
FullEstimatorInASingleStepWithSmoothMargin
from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import SmoothMarginEstimator
-from extreme_estimator.extreme_models.margin_model.linear_margin_model import LinearAllParametersAllDimsMarginModel
+from extreme_estimator.extreme_models.margin_model.linear_margin_model import LinearAllParametersAllDimsMarginModel, \
+ LinearNonStationaryMarginModel, LinearStationaryMarginModel
from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
AbstractMarginFunction
from extreme_estimator.extreme_models.margin_model.param_function.param_function import AbstractParamFunction
@@ -248,7 +249,7 @@ class StudyVisualizer(object):
def visualize_linear_margin_fit(self, only_first_max_stable=False):
default_covariance_function = CovarianceFunction.powexp
- margin_class = LinearAllParametersAllDimsMarginModel
+ margin_class = LinearNonStationaryMarginModel if self.temporal_non_stationarity else LinearStationaryMarginModel
plot_name = 'Full Likelihood with Linear marginals and max stable dependency structure'
plot_name += '\n(with {} covariance structure when a covariance is needed)'.format(
str(default_covariance_function).split('.')[-1])
@@ -281,13 +282,13 @@ class StudyVisualizer(object):
self.visualize_independent_margin_fits(threshold=None, axes=axes[-1], show=False)
# Plot the smooth margin only
- margin_model = margin_class(coordinates=self.coordinates)
+ margin_model = margin_class(coordinates=self.coordinates, starting_point=None)
estimator = SmoothMarginEstimator(dataset=self.dataset, margin_model=margin_model)
self.fit_and_visualize_estimator(estimator, axes[0], title='without max stable')
# Plot the smooth margin fitted with a max stable
for i, max_stable_model in enumerate(max_stable_models, 1):
- margin_model = margin_class(coordinates=self.coordinates)
+ margin_model = margin_class(coordinates=self.coordinates, starting_point=None)
estimator = FullEstimatorInASingleStepWithSmoothMargin(self.dataset, margin_model, max_stable_model)
title = get_display_name_from_object_type(type(max_stable_model))
self.fit_and_visualize_estimator(estimator, axes[i], title=title)
diff --git a/extreme_estimator/estimator/full_estimator/abstract_full_estimator.py b/extreme_estimator/estimator/full_estimator/abstract_full_estimator.py
index fa9c0e768dad34a20f8cbc5c43e3e4ec743505ec..07ed06c1577951130fb97bd5d0705e56c5e1304e 100644
--- a/extreme_estimator/estimator/full_estimator/abstract_full_estimator.py
+++ b/extreme_estimator/estimator/full_estimator/abstract_full_estimator.py
@@ -50,9 +50,12 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
super().__init__(dataset)
self.max_stable_model = max_stable_model
self.linear_margin_model = margin_model
- self.linear_margin_function_to_fit = self.linear_margin_model.margin_function_start_fit
assert isinstance(self.linear_margin_function_to_fit, LinearMarginFunction)
+ @property
+ def linear_margin_function_to_fit(self):
+ return self.linear_margin_model.margin_function_start_fit
+
def _fit(self):
# Estimate both the margin and the max-stable structure
self._result_from_fit = self.max_stable_model.fitmaxstab(
diff --git a/extreme_estimator/extreme_models/margin_model/fitspatgev.R b/extreme_estimator/extreme_models/margin_model/fitspatgev.R
index f8c7d264259775a5956c144f0103ad391d3ef61a..e7fd380e2f9ccec0a42126dda1a552ee262ac988 100644
--- a/extreme_estimator/extreme_models/margin_model/fitspatgev.R
+++ b/extreme_estimator/extreme_models/margin_model/fitspatgev.R
@@ -74,8 +74,8 @@ fitspatgev_3D_test <- function (n.obs){
# Add the temporal covariates
temp.form.loc = loc ~ T
temp.form.scale = scale ~ T
- temp.form.shape = shape ~ T
- # temp.form.shape = NULL
+ # temp.form.shape = shape ~ T
+ temp.form.shape = NULL
temp.cov = matrix(1:n.obs, ncol=1)
colnames(temp.cov) = c("T")
@@ -89,6 +89,8 @@ fitspatgev_3D_test <- function (n.obs){
start = list(locCoeff1=0.0, locCoeff2=1.0, scaleCoeff1=0.1, shapeCoeff1=1.0, scaleCoeff2=0.1, shapeCoeff2=1.0,
tempCoeffLoc1=0.0, tempCoeffScale1=0.1, tempCoeffShape1=0.1)
+ start = list(locCoeff1=0.0, locCoeff2=1.0, scaleCoeff1=0.1, shapeCoeff1=1.0, scaleCoeff2=0.1, shapeCoeff2=1.0,
+ tempCoeffLoc1=0.0, tempCoeffScale1=0.1)
res = fitspatgev(data=data, covariables=covariables, start=start, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form,
temp.cov=temp.cov, temp.form.loc=temp.form.loc, temp.form.scale=temp.form.scale, temp.form.shape = temp.form.shape)
@@ -105,6 +107,7 @@ fitspatgev_3D_test <- function (n.obs){
loc.form <- update(loc.form, y ~ .)
scale.form <- update(scale.form, y ~ .)
shape.form <- update(shape.form, y ~ .)
+ print('here')
print(use.temp.cov)
if (use.temp.cov[1])
temp.form.loc <- update(temp.form.loc, y ~ . + 0)
diff --git a/extreme_estimator/extreme_models/margin_model/linear_margin_model.py b/extreme_estimator/extreme_models/margin_model/linear_margin_model.py
index 9408fb72316c5365560ad874ad16e2190fe89928..adbb16d926c3a426f3114677042dc766070472d3 100644
--- a/extreme_estimator/extreme_models/margin_model/linear_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/linear_margin_model.py
@@ -17,7 +17,7 @@ class LinearMarginModel(ParametricMarginModel):
def load_margin_functions(self, gev_param_name_to_dims=None):
assert gev_param_name_to_dims is not None, 'LinearMarginModel cannot be used for sampling/fitting \n' \
- 'load_margin_functions needs to be implemented in child class'
+ 'load_margin_functions needs to be implemented in child class'
# Load default params (with a dictionary format to enable quick replacement)
# IMPORTANT: Using a dictionary format enable using the default/user params methodology
self.default_params_sample = self.default_param_name_and_dim_to_coef
@@ -27,13 +27,15 @@ class LinearMarginModel(ParametricMarginModel):
coef_sample = self.gev_param_name_to_linear_coef(param_name_and_dim_to_coef=self.params_sample)
self.margin_function_sample = LinearMarginFunction(coordinates=self.coordinates,
gev_param_name_to_coef=coef_sample,
- gev_param_name_to_dims=gev_param_name_to_dims)
+ gev_param_name_to_dims=gev_param_name_to_dims,
+ starting_point=self.starting_point)
# Load start fit coef
coef_start_fit = self.gev_param_name_to_linear_coef(param_name_and_dim_to_coef=self.params_start_fit)
self.margin_function_start_fit = LinearMarginFunction(coordinates=self.coordinates,
gev_param_name_to_coef=coef_start_fit,
- gev_param_name_to_dims=gev_param_name_to_dims)
+ gev_param_name_to_dims=gev_param_name_to_dims,
+ starting_point=self.starting_point)
@property
def default_param_name_and_dim_to_coef(self) -> dict:
@@ -49,7 +51,8 @@ class LinearMarginModel(ParametricMarginModel):
def gev_param_name_to_linear_coef(self, param_name_and_dim_to_coef):
gev_param_name_to_linear_coef = {}
for gev_param_name in GevParams.PARAM_NAMES:
- idx_to_coef = {idx: param_name_and_dim_to_coef[(gev_param_name, idx)] for idx in [-1] + self.coordinates.coordinates_dims}
+ idx_to_coef = {idx: param_name_and_dim_to_coef[(gev_param_name, idx)] for idx in
+ [-1] + self.coordinates.coordinates_dims}
linear_coef = LinearCoef(gev_param_name=gev_param_name, idx_to_coef=idx_to_coef)
gev_param_name_to_linear_coef[gev_param_name] = linear_coef
return gev_param_name_to_linear_coef
@@ -101,3 +104,28 @@ class LinearAllParametersAllDimsMarginModel(LinearMarginModel):
super().load_margin_functions({GevParams.SHAPE: self.coordinates.coordinates_dims,
GevParams.LOC: self.coordinates.coordinates_dims,
GevParams.SCALE: self.coordinates.coordinates_dims})
+
+
+class LinearAllParametersTwoFirstCoordinatesMarginModel(LinearMarginModel):
+
+ def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_dims=None):
+ super().load_margin_functions({GevParams.SHAPE: [0, 1],
+ GevParams.LOC: [0, 1],
+ GevParams.SCALE: [0, 1]})
+
+
+class LinearAllTwoStatialCoordinatesLocationLinearMarginModel(LinearMarginModel):
+
+ def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_dims=None):
+ super().load_margin_functions({GevParams.SHAPE: [0, 1],
+ GevParams.LOC: [0, 1, 2],
+ GevParams.SCALE: [0, 1]})
+
+
+# Some renaming that defines the stationary and non-stationary models of reference
+class LinearStationaryMarginModel(LinearAllParametersTwoFirstCoordinatesMarginModel):
+ pass
+
+
+class LinearNonStationaryMarginModel(LinearAllTwoStatialCoordinatesLocationLinearMarginModel):
+ pass
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py
index 427d0d38ee175a1a85638a3bc30a34b79dd8ac5c..103d6d382be7fdda5e870c2456b106cc86ae0113 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py
@@ -1,5 +1,7 @@
from typing import Dict, List
+import numpy as np
+
from extreme_estimator.extreme_models.margin_model.margin_function.parametric_margin_function import \
ParametricMarginFunction
from extreme_estimator.extreme_models.margin_model.param_function.abstract_coef import AbstractCoef
@@ -27,7 +29,10 @@ class LinearMarginFunction(ParametricMarginFunction):
COEF_CLASS = LinearCoef
def __init__(self, coordinates: AbstractCoordinates, gev_param_name_to_dims: Dict[str, List[int]],
- gev_param_name_to_coef: Dict[str, AbstractCoef]):
+ gev_param_name_to_coef: Dict[str, AbstractCoef],
+ starting_point=None):
+ # Starting point for the trend is the same for all the parameters
+ self.starting_point = starting_point
self.gev_param_name_to_coef = None # type: Dict[str, LinearCoef]
super().__init__(coordinates, gev_param_name_to_dims, gev_param_name_to_coef)
@@ -36,6 +41,16 @@ class LinearMarginFunction(ParametricMarginFunction):
coordinates=self.coordinates.coordinates_values(),
linear_coef=self.gev_param_name_to_coef[gev_param_name])
+ def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
+ if self.starting_point is not None:
+ # Shift temporal coordinate to enable to model temporal trend with starting point
+ assert self.coordinates.has_temporal_coordinates
+ idx_coordinate_t = self.coordinates.coordinates_names.index(self.coordinates.COORDINATE_T)
+ assert 0 <= idx_coordinate_t < len(coordinate)
+ if coordinate[idx_coordinate_t] < self.starting_point:
+ coordinate[idx_coordinate_t] = self.starting_point
+ return super().get_gev_params(coordinate)
+
@classmethod
def idx_to_coefficient_name(cls, coordinates: AbstractCoordinates) -> Dict[int, str]:
# Intercept correspond to the dimension 0
diff --git a/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py b/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py
index c9946d727a9b97585227afd6001bd1e1ebd9ae3b..2a933e9da11d919de95990dd09ed325658a7351a 100644
--- a/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py
@@ -15,7 +15,11 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
class ParametricMarginModel(AbstractMarginModel, ABC):
def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
- params_sample=None):
+ params_sample=None, starting_point=None):
+ """
+ :param starting_point: starting coordinate for the temporal trend
+ """
+ self.starting_point = starting_point # type: int
self.margin_function_sample = None # type: ParametricMarginFunction
self.margin_function_start_fit = None # type: ParametricMarginFunction
super().__init__(coordinates, use_start_value, params_start_fit, params_sample)
@@ -23,6 +27,15 @@ class ParametricMarginModel(AbstractMarginModel, ABC):
def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,
df_coordinates_temp: pd.DataFrame) -> ResultFromFit:
assert data.shape[1] == len(df_coordinates_spat)
+ # assert data.shape[0] == len(df_coordinates_temp)
+
+ # Enforce a starting point for the temporal trend
+ if self.starting_point is not None:
+ ind_to_modify = df_coordinates_temp.iloc[:, 0] <= self.starting_point # type: pd.Series
+ # Assert that some coordinates are selected but not all (at least 20 data should be left for temporal trend)
+ assert 0 < sum(ind_to_modify) < len(ind_to_modify) - 20
+ # Modify the temporal coordinates to enforce the stationarity
+ df_coordinates_temp.loc[ind_to_modify] = self.starting_point
fit_params = get_margin_formula(self.margin_function_start_fit.form_dict)
diff --git a/extreme_estimator/extreme_models/utils.py b/extreme_estimator/extreme_models/utils.py
index 057dfdfb6663d2973800c3717088ccf556f7d965..327ea0d6d95f42aa830312617de049e8f25983fc 100644
--- a/extreme_estimator/extreme_models/utils.py
+++ b/extreme_estimator/extreme_models/utils.py
@@ -83,7 +83,9 @@ def get_coord(df_coordinates: pd.DataFrame):
def get_margin_formula(fit_marge_form_dict) -> Dict:
- return {k: robjects.Formula(v) for k, v in fit_marge_form_dict.items()}
+ as_null = r['as.null']
+ margin_formula = {k: robjects.Formula(v) if v != 'NULL' else as_null(1.0) for k, v in fit_marge_form_dict.items()}
+ return margin_formula
# def conversion_to_FloatVector(data):
# """Convert DataFrame or numpy array into FloatVector for r"""
@@ -91,3 +93,4 @@ def get_margin_formula(fit_marge_form_dict) -> Dict:
# data = data.values
# assert isinstance(data, np.ndarray)
# return npr.numpy2ri(data)
+