diff --git a/extreme_estimator/estimator/full_estimator.py b/extreme_estimator/estimator/full_estimator.py
index bd3e1618185eadd8a3f0983a7911ac42101e412d..8c6fea1d509d273261ccdd9d97b1ea705193cb60 100644
--- a/extreme_estimator/estimator/full_estimator.py
+++ b/extreme_estimator/estimator/full_estimator.py
@@ -1,6 +1,4 @@
from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
-from extreme_estimator.extreme_models.margin_model.margin_function.independent_margin_function import \
- LinearMarginFunction
from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearMarginModel
from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import AbstractMaxStableModel
from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py
index e5a593b45c1f1241c22fd4e1f075a2dd0c05e33f..ef0928c8e15a4efda995fe7a359dac5c839dd536 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py
@@ -7,17 +7,18 @@ from spatio_temporal_dataset.spatial_coordinates.abstract_spatial_coordinates im
class AbstractMarginFunction(object):
- """
- It represents any function mapping points from a space S (could be 2D, 3D,...) to R^3 (the 3 parameters of the GEV)
- """
+ """ Class of function mapping points from a space S (could be 1D, 2D,...) to R^3 (the 3 parameters of the GEV)"""
def __init__(self, spatial_coordinates: AbstractSpatialCoordinates, default_params: GevParams):
self.spatial_coordinates = spatial_coordinates
self.default_params = default_params.to_dict()
def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
+ """Main function that maps each coordinate to its GEV parameters"""
pass
+ # Visualization function
+
def visualize_2D(self, gev_param_name=GevParams.GEV_LOC, ax=None, show=False):
x = self.spatial_coordinates.x_coordinates
y = self.spatial_coordinates.y_coordinates
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py
index a9f416fa3228082c02985e0529dd00fdbed30d18..422a35848fda151da21122e028b5ab4c77a9dfe9 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py
@@ -1,101 +1,106 @@
-from typing import Dict
+from typing import Dict, List, Tuple
import numpy as np
+from extreme_estimator.extreme_models.margin_model.margin_function.param_function import ConstantParamFunction, \
+ LinearOneAxisParamFunction, ParamFunction, LinearParamFunction
from extreme_estimator.gev_params import GevParams
-from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import AbstractMarginFunction
+from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
+ AbstractMarginFunction
from spatio_temporal_dataset.spatial_coordinates.abstract_spatial_coordinates import AbstractSpatialCoordinates
-class ParamFunction(object):
-
- def get_gev_param_value(self, coordinate: np.ndarray) -> float:
- pass
-
-
class IndependentMarginFunction(AbstractMarginFunction):
+ """Margin Function where each parameter of the GEV are modeled independently"""
def __init__(self, spatial_coordinates: AbstractSpatialCoordinates, default_params: GevParams):
+ """Attribute 'gev_param_name_to_param_function' maps each GEV parameter to its corresponding function"""
super().__init__(spatial_coordinates, default_params)
self.gev_param_name_to_param_function = None # type: Dict[str, ParamFunction]
def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
+ """Each GEV parameter is computed independently through its corresponding param_function"""
assert self.gev_param_name_to_param_function is not None
assert len(self.gev_param_name_to_param_function) == 3
gev_params = {}
for gev_param_name in GevParams.GEV_PARAM_NAMES:
param_function = self.gev_param_name_to_param_function[gev_param_name]
- gev_value = param_function.get_gev_param_value(coordinate)
- gev_params[gev_param_name] = gev_value
+ gev_params[gev_param_name] = param_function.get_gev_param_value(coordinate)
return GevParams.from_dict(gev_params)
-class ConstantParamFunction(ParamFunction):
-
- def __init__(self, constant):
- self.constant = constant
-
- def get_gev_param_value(self, coordinate: np.ndarray) -> float:
- return self.constant
-
-
-class LinearOneAxisParamFunction(ParamFunction):
-
- def __init__(self, linear_axis: int, coordinates_axis: np.ndarray, start: float, end: float = 0.0):
- self.linear_axis = linear_axis
- self.t_min = coordinates_axis.min()
- self.t_max = coordinates_axis.max()
- self.start = start
- self.end = end
-
- def get_gev_param_value(self, coordinate: np.ndarray) -> float:
- t = coordinate[self.linear_axis]
- t_between_zero_and_one = (t - self.t_min) / self.t_max
- return self.start + t_between_zero_and_one * (self.end - self.start)
-
-
class LinearMarginFunction(IndependentMarginFunction):
- """
- On the minimal point along all the dimension, the GevParms will equal default params
- Otherwise, it will augment linearly along a single linear axis
- """
-
- def __init__(self, spatial_coordinates: AbstractSpatialCoordinates, default_params: GevParams,
- gev_param_name_to_linear_axis: Dict[str, int]):
+ """ Margin Function, where each parameter can augment linearly as follows:
+ On the extremal point along all the dimension, the GEV parameters will equal the default_param value
+ Then, it will augment linearly along a single linear axis"""
+
+ def __init__(self, spatial_coordinates: AbstractSpatialCoordinates,
+ default_params: GevParams,
+ gev_param_name_to_linear_axes: Dict[str, List[int]],
+ gev_param_name_and_axis_to_start_fit: Dict[Tuple[str, int], float] = None):
+ """
+ -Attribute 'gev_param_name_to_linear_axis' maps each GEV parameter to its corresponding function
+ -Attribute 'gev_param_name_and_axis_to_start_fit' maps each tuple (GEV parameter, axis) to its start value for
+ fitting (by default equal to 1). Also start value for the intercept is equal to 0 by default."""
super().__init__(spatial_coordinates, default_params)
+ self.gev_param_name_and_axis_to_start_fit = gev_param_name_and_axis_to_start_fit
+ self.gev_param_name_to_linear_axes = gev_param_name_to_linear_axes
- # By default, set a constant function for all the params
- self.gev_param_to_str_form = {gev_param_name: '1' for gev_param_name in GevParams.GEV_PARAM_NAMES}
+ # Check the axes are well-defined with respect to the coordinates
+ for axes in self.gev_param_name_to_linear_axes.values():
+ assert all([axis < np.ndim(spatial_coordinates.coordinates) for axis in axes])
- self.param_to_linear_dims = gev_param_name_to_linear_axis
- assert all([axis < np.ndim(spatial_coordinates.coordinates) for axis in gev_param_name_to_linear_axis.values()])
- # Initialize gev_parameter_to_param_function
- self.gev_param_name_to_param_function = {}
+ # Build gev_parameter_to_param_function dictionary
+ self.gev_param_name_to_param_function = {} # type: Dict[str, ParamFunction]
+ # Map each gev_param_name to its corresponding param_function
for gev_param_name in GevParams.GEV_PARAM_NAMES:
- if gev_param_name not in gev_param_name_to_linear_axis.keys():
+ # By default, if gev_param_name linear_axis is not specified, a constantParamFunction is chosen
+ if gev_param_name not in self.gev_param_name_to_linear_axes.keys():
param_function = ConstantParamFunction(constant=self.default_params[gev_param_name])
+ # Otherwise, we fit a LinearParamFunction
else:
- linear_axis = gev_param_name_to_linear_axis.get(gev_param_name, None)
- coordinates_axis = self.spatial_coordinates.coordinates[:, linear_axis]
- param_function = LinearOneAxisParamFunction(linear_axis=linear_axis, coordinates_axis=coordinates_axis,
- start=self.default_params[gev_param_name])
+ param_function = LinearParamFunction(linear_axes=self.gev_param_name_to_linear_axes[gev_param_name],
+ coordinates=self.spatial_coordinates.coordinates,
+ start=self.default_params[gev_param_name])
+ # Some check on the Linear param function
+ if gev_param_name == GevParams.GEV_SCALE:
+ assert param_function.end > param_function.start, 'Impossible linear rate for Scale parameter'
+
+ # Add the param_function to the dictionary
self.gev_param_name_to_param_function[gev_param_name] = param_function
-
@property
- def fit_marge_form_dict(self):
+ def fit_marge_form_dict(self) -> dict:
"""
- Example:
- loc.form = loc ~ E
- scale.form = scale ~ N
- shape.form = shape ~ E+N
+ Example of formula that could be specified:
+ loc.form = loc ~ coord_x
+ scale.form = scale ~ coord_y
+ shape.form = shape ~ coord_x+coord_y
:return:
"""
- # todo specify all functions
- return {gev_param_name + '.form': gev_param_name + ' ~ ' + self.gev_param_to_str_form[gev_param_name]
- for gev_param_name in GevParams.GEV_PARAM_NAMES}
+ fit_marge_form_dict = {}
+ axis_to_name = {i: name for i, name in enumerate(AbstractSpatialCoordinates.COORDINATE_NAMES)}
+ for gev_param_name in GevParams.GEV_PARAM_NAMES:
+ axes = self.gev_param_name_to_linear_axes.get(gev_param_name, None)
+ formula_str = '1' if axes is None else '+'.join([axis_to_name[axis] for axis in axes])
+ fit_marge_form_dict[gev_param_name + '.form'] = gev_param_name + ' ~ ' + formula_str
+ return fit_marge_form_dict
@property
- def margin_start_dict(self):
- """Example for the names: locCoeff1 locCoeff2 """
- return {gev_param_name + 'Coeff1': 0 for gev_param_name in GevParams.GEV_PARAM_NAMES}
+ def margin_start_dict(self) -> dict:
+ # Define default values
+ default_start_fit_coef = 1.0
+ default_start_fit_intercept = 0.0
+ # Build the dictionary containing all the parameters
+ margin_start_dict = {}
+ for gev_param_name in GevParams.GEV_PARAM_NAMES:
+ coef_template_str = gev_param_name + 'Coeff{}'
+ # Constant param must be specified for all the parameters
+ margin_start_dict[coef_template_str.format(1)] = default_start_fit_intercept
+ for j, axis in enumerate(self.gev_param_name_to_linear_axes.get(gev_param_name, []), 2):
+ if self.gev_param_name_and_axis_to_start_fit is None:
+ coef = default_start_fit_coef
+ else:
+ coef = self.gev_param_name_and_axis_to_start_fit.get((gev_param_name, axis), default_start_fit_coef)
+ margin_start_dict[coef_template_str.format(j)] = coef
+ return margin_start_dict
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/param_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/param_function.py
new file mode 100644
index 0000000000000000000000000000000000000000..5503da2c5aca4157377698df3102a411c270668f
--- /dev/null
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/param_function.py
@@ -0,0 +1,49 @@
+from typing import List
+
+import numpy as np
+
+
+class ParamFunction(object):
+
+ def get_gev_param_value(self, coordinate: np.ndarray) -> float:
+ pass
+
+
+class ConstantParamFunction(ParamFunction):
+
+ def __init__(self, constant):
+ self.constant = constant
+
+ def get_gev_param_value(self, coordinate: np.ndarray) -> float:
+ return self.constant
+
+
+class LinearOneAxisParamFunction(ParamFunction):
+
+ def __init__(self, linear_axis: int, coordinates_axis: np.ndarray, start: float, end: float = 2.0):
+ self.linear_axis = linear_axis
+ self.t_min = coordinates_axis.min()
+ self.t_max = coordinates_axis.max()
+ self.start = start
+ self.end = end
+
+ def get_gev_param_value(self, coordinate: np.ndarray) -> float:
+ t = coordinate[self.linear_axis]
+ t_between_zero_and_one = (t - self.t_min) / (self.t_max - self.t_min)
+ assert 0 <= t_between_zero_and_one <= 1, 'Out of bounds'
+ return self.start + t_between_zero_and_one * (self.end - self.start)
+
+
+class LinearParamFunction(ParamFunction):
+
+ def __init__(self, linear_axes: List[int], coordinates: np.ndarray, start: float, end: float = 2.0):
+ self.linear_one_axis_param_functions = [] # type: List[LinearOneAxisParamFunction]
+ self.start = start
+ self.end = end
+ for linear_axis in linear_axes:
+ param_function = LinearOneAxisParamFunction(linear_axis, coordinates[:, linear_axis], start, end)
+ self.linear_one_axis_param_functions.append(param_function)
+
+ def get_gev_param_value(self, coordinate: np.ndarray) -> float:
+ values = [param_funct.get_gev_param_value(coordinate) for param_funct in self.linear_one_axis_param_functions]
+ return float(np.mean(values))
diff --git a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
index b15b8c6dac1c75fc6b4d68e839d13523bb0a30ea..aeec0112747f237c0dd81054479fc06354090d71 100644
--- a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
@@ -1,7 +1,9 @@
import numpy as np
-from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import AbstractMarginFunction
-from extreme_estimator.extreme_models.margin_model.margin_function.independent_margin_function import LinearMarginFunction
+from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
+ AbstractMarginFunction
+from extreme_estimator.extreme_models.margin_model.margin_function.independent_margin_function import \
+ LinearMarginFunction
from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
from extreme_estimator.gev_params import GevParams
@@ -13,10 +15,10 @@ class LinearMarginModel(AbstractMarginModel):
self.default_params_start_fit = GevParams(1.0, 1.0, 1.0).to_dict()
self.margin_function_sample = LinearMarginFunction(spatial_coordinates=self.spatial_coordinates,
default_params=GevParams.from_dict(self.params_sample),
- gev_param_name_to_linear_axis=gev_param_name_to_linear_axis)
+ gev_param_name_to_linear_axes=gev_param_name_to_linear_axis)
self.margin_function_start_fit = LinearMarginFunction(spatial_coordinates=self.spatial_coordinates,
default_params=GevParams.from_dict(self.params_start_fit),
- gev_param_name_to_linear_axis=gev_param_name_to_linear_axis)
+ gev_param_name_to_linear_axes=gev_param_name_to_linear_axis)
def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates: np.ndarray) -> AbstractMarginFunction:
return self.margin_function_start_fit
@@ -27,18 +29,30 @@ class ConstantMarginModel(LinearMarginModel):
def load_margin_functions(self, gev_param_name_to_linear_axis=None):
super().load_margin_functions({})
- def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates: np.ndarray) -> AbstractMarginFunction:
- return self.margin_function_start_fit
-
class LinearShapeAxis0MarginModel(LinearMarginModel):
def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_axis=None):
- super().load_margin_functions({GevParams.GEV_SHAPE: 0})
+ super().load_margin_functions({GevParams.GEV_SHAPE: [0]})
+
+
+class LinearShapeAxis0and1MarginModel(LinearMarginModel):
+
+ def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_axis=None):
+ super().load_margin_functions({GevParams.GEV_SHAPE: [0, 1]})
+
+
+class LinearAllParametersAxis0MarginModel(LinearMarginModel):
+
+ def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_axis=None):
+ super().load_margin_functions({GevParams.GEV_SHAPE: [0],
+ GevParams.GEV_LOC: [0],
+ GevParams.GEV_SCALE: [0]})
- # def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates: np.ndarray) -> AbstractMarginFunction:
- # pass
+class LinearAllParametersAxis0And1MarginModel(LinearMarginModel):
-if __name__ == '__main__':
- pass
+ def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_axis=None):
+ super().load_margin_functions({GevParams.GEV_SHAPE: [0, 1],
+ GevParams.GEV_LOC: [0, 1],
+ GevParams.GEV_SCALE: [0, 1]})
diff --git a/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py b/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
index 1043f5bc51a3c9879a57ce9f26341b3e9e0323fc..44643e1bf62ac35f4e30758f1aa7e903af2c6ebc 100644
--- a/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
+++ b/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
@@ -3,8 +3,7 @@ from enum import Enum
import numpy as np
import rpy2
-import rpy2.robjects as ro
-from rpy2.robjects import ListVector
+import rpy2.robjects as robjects
from extreme_estimator.extreme_models.abstract_model import AbstractModel
@@ -27,24 +26,23 @@ class AbstractMaxStableModel(AbstractModel):
assert fit_marge_form_dict is not None
assert margin_start_dict is not None
- # Add the colnames to df_coordinates DataFrame to enable specification of the margin functions
- df_coordinates = df_coordinates.copy()
- df_coordinates.colnames = ro.StrVector(list(df_coordinates.columns))
- # Transform the formula string representation into robjects.Formula("y ~ x")
- # Specify the fit params
+ # Prepare the data and the coord objects
+ data = np.transpose(maxima_frech)
+ coord = robjects.vectors.Matrix(df_coordinates.values)
+ coord.colnames = robjects.StrVector(list(df_coordinates.columns))
+
+ # Prepare the fit params
fit_params = self.cov_mod_param.copy()
start_dict = self.params_start_fit
if fit_marge:
start_dict.update(margin_start_dict)
- fit_params.update({k: ro.Formula(v) for k, v in fit_marge_form_dict.items()})
+ fit_params.update({k: robjects.Formula(v) for k, v in fit_marge_form_dict.items()})
fit_params['start'] = self.r.list(**start_dict)
fit_params['fit.marge'] = fit_marge
- # Run the fitmaxstab in R
- # todo: find how to specify the optim function to use
- coord = df_coordinates.values
+ # Run the fitmaxstab in R
try:
- res = self.r.fitmaxstab(data=np.transpose(maxima_frech), coord=coord, **fit_params) # type: ListVector
+ res = self.r.fitmaxstab(data=data, coord=coord, **fit_params) # type: robjects.ListVector
except rpy2.rinterface.RRuntimeError as error:
raise Exception('Some R exception have been launched at RunTime: \n {}'.format(error.__repr__()))
# todo: maybe if the convergence was not successful I could try other starting point several times
diff --git a/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R b/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R
index 4daba7f105b904c578c7df2d1f429feced3f5235..24f36d60416382633ffdc8b5d5c8470994e14729 100644
--- a/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R
+++ b/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R
@@ -20,6 +20,8 @@ if (call_main) {
# res = fitmaxstab(data, coord, "gauss", fit.marge=FALSE, )
# res = fitmaxstab(data, coord, "brown")
# res = fitmaxstab(data, coord, "whitmat", start=)
+ print(class(coord))
+ print(colnames(coord))
loc.form = loc ~ 1
scale.form = scale ~ 1
@@ -32,10 +34,10 @@ if (call_main) {
# namedlist = list(cov11 = 1.0, cov12 = 1.2, cov22 = 2.2)
# res = fitmaxstab(data=data, coord=coord, cov.mod="gauss", start=namedlist)
- for (name in names(res)){
- print(name)
- print(res[name])
- }
+ # for (name in names(res)){
+ # print(name)
+ # print(res[name])
+ # }
print(res['fitted.values'])
# print(res['convergence'])
diff --git a/extreme_estimator/gev_params.py b/extreme_estimator/gev_params.py
index 0f7c7cdd85cceb8a3d9b5af0c8a655942018abdc..31cb8247a5d3ee9e1acd7f0d6c9c79787ee9db7c 100644
--- a/extreme_estimator/gev_params.py
+++ b/extreme_estimator/gev_params.py
@@ -11,6 +11,7 @@ class GevParams(object):
self.location = loc
self.scale = scale
self.shape = shape
+ assert self.scale > 0
@classmethod
def from_dict(cls, params: dict):
diff --git a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
index fe77b35ac48b9bcaa19c5f8faced344a5b8b345c..8729ef3f0455d7dd74b8debd5d1296245319c0d9 100644
--- a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
@@ -2,7 +2,8 @@ import unittest
from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
from extreme_estimator.extreme_models.margin_model.smooth_margin_model import ConstantMarginModel, \
- LinearShapeAxis0MarginModel
+ LinearShapeAxis0MarginModel, LinearShapeAxis0and1MarginModel, LinearAllParametersAxis0MarginModel, \
+ LinearAllParametersAxis0And1MarginModel
from extreme_estimator.estimator.margin_estimator import SmoothMarginEstimator
from extreme_estimator.return_level_plot.spatial_2D_plot import Spatial2DPlot
from spatio_temporal_dataset.dataset.simulation_dataset import MarginDataset
@@ -11,7 +12,9 @@ from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import Ci
class TestSmoothMarginEstimator(unittest.TestCase):
DISPLAY = False
- MARGIN_TYPES = [ConstantMarginModel, LinearShapeAxis0MarginModel][1:]
+ MARGIN_TYPES = [ConstantMarginModel, LinearShapeAxis0MarginModel,
+ LinearShapeAxis0and1MarginModel, LinearAllParametersAxis0MarginModel,
+ LinearAllParametersAxis0And1MarginModel][:]
SMOOTH_MARGIN_ESTIMATORS = [SmoothMarginEstimator]
def setUp(self):
@@ -33,7 +36,7 @@ class TestSmoothMarginEstimator(unittest.TestCase):
# Map name to their margin functions
name_to_margin_function = {
'Ground truth margin function': dataset.margin_model.margin_function_sample,
- # 'Estimated margin function': estimator.margin_function_fitted,
+ 'Estimated margin function': estimator.margin_function_fitted,
}
# Spatial Plot
if self.DISPLAY: