[MARGIN FUNCTION] refactor LinearMarginFunction object

extreme_estimator/estimator/full_estimator.py

 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

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

extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py

-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

extreme_estimator/extreme_models/margin_model/margin_function/param_function.py

+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))

extreme_estimator/extreme_models/margin_model/smooth_margin_model.py

 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]})

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

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'])
 

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):

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: