diff --git a/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py b/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
index d4ca048ff05cee419c60c8397101cac85855465e..c052dc184085846fef40f06d503bdeb5ff6e4e0d 100644
--- a/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
@@ -1,17 +1,21 @@
+from abc import ABC
+
 import numpy as np
 import pandas as pd
 
 from extreme_estimator.extreme_models.abstract_model import AbstractModel
 from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function \
     import AbstractMarginFunction
+from extreme_estimator.extreme_models.result_from_fit import ResultFromFit
 from extreme_estimator.extreme_models.utils import r
 from extreme_estimator.margin_fits.gev.gev_params import GevParams
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
-class AbstractMarginModel(AbstractModel):
+class AbstractMarginModel(AbstractModel, ABC):
 
-    def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None, params_sample=None):
+    def __init__(self, coordinates: AbstractCoordinates, use_start_value=False,
+                 params_start_fit=None, params_sample=None):
         super().__init__(use_start_value, params_start_fit, params_sample)
         assert isinstance(coordinates, AbstractCoordinates), type(coordinates)
         self.coordinates = coordinates
@@ -20,7 +24,7 @@ class AbstractMarginModel(AbstractModel):
         self.load_margin_functions()
 
     def load_margin_functions(self):
-        pass
+        raise NotImplementedError
 
     def default_load_margin_functions(self, margin_function_class):
         self.margin_function_sample = margin_function_class(coordinates=self.coordinates,
@@ -70,7 +74,7 @@ class AbstractMarginModel(AbstractModel):
     # Fitting methods needs to be defined in child classes
 
     def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, df_coordinates_spatial: pd.DataFrame,
-                                  df_coordinates_temporal: pd.DataFrame) -> AbstractMarginFunction:
+                                  df_coordinates_temporal: pd.DataFrame) -> ResultFromFit:
         raise NotImplementedError
 
 
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 1e77a4ea42be0516f77e79603284081dadf3f30c..25eeba991a0211c7ff43640ea9fe31d1b63d2ec9 100644
--- a/extreme_estimator/extreme_models/margin_model/linear_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/linear_margin_model.py
@@ -1,16 +1,11 @@
-import numpy as np
-import pandas as pd
-
-from extreme_estimator.extreme_models.result_from_fit import ResultFromFit
 from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
 from extreme_estimator.extreme_models.margin_model.margin_function.linear_margin_function import LinearMarginFunction
 from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
-from extreme_estimator.extreme_models.utils import safe_run_r_estimator, r, get_coord, \
-    get_margin_formula
+from extreme_estimator.extreme_models.margin_model.parametric_margin_model import ParametricMarginModel
 from extreme_estimator.margin_fits.gev.gev_params import GevParams
 
 
-class LinearMarginModel(AbstractMarginModel):
+class LinearMarginModel(ParametricMarginModel):
 
     def load_margin_functions(self, gev_param_name_to_linear_dims=None):
         assert gev_param_name_to_linear_dims is not None, 'LinearMarginModel cannot be used for sampling/fitting \n' \
@@ -57,26 +52,6 @@ class LinearMarginModel(AbstractMarginModel):
                 params[(gev_param_name, dim)] = coef
         return cls(coordinates, params_sample=params, params_start_fit=params)
 
-    def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, df_coordinates_spat: pd.DataFrame,
-                                  df_coordinates_temp: pd.DataFrame) -> ResultFromFit:
-        # The reshaping on the line below is only valid if we have a single observation per spatio-temporal point
-        if maxima_gev.shape[1] == 1:
-            maxima_gev = maxima_gev.reshape([len(df_coordinates_temp), len(df_coordinates_spat)])
-        data = np.transpose(maxima_gev)
-
-        fit_params = get_margin_formula(self.margin_function_start_fit.form_dict)
-
-        # Covariables
-        covariables = get_coord(df_coordinates=df_coordinates_spat)
-        fit_params['temp.cov'] = get_coord(df_coordinates=df_coordinates_temp)
-
-        # Start parameters
-        coef_dict = self.margin_function_start_fit.coef_dict
-        fit_params['start'] = r.list(**coef_dict)
-
-        return safe_run_r_estimator(function=r.fitspatgev, use_start=self.use_start_value, data=data,
-                                    covariables=covariables, **fit_params)
-
 
 class ConstantMarginModel(LinearMarginModel):
 
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 6ae93485d6c9fc24d03f40fc4101b61c34c60eeb..44f704459e1b13374a0b144b6ea3b95531b89c05 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
@@ -46,7 +46,7 @@ class AbstractMarginFunction(object):
 
     def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
         """Main method that maps each coordinate to its GEV parameters"""
-        pass
+        raise NotImplementedError
 
     @property
     def gev_value_name_to_serie(self) -> Dict[str, pd.Series]:
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 4196d16e51725d6a2c29e74628607772d07c69e4..34134e31cedbf5fe1d22fbb3f969276334051733 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,7 +1,7 @@
 from typing import Dict, List
 
-from extreme_estimator.extreme_models.margin_model.margin_function.independent_margin_function import \
-    IndependentMarginFunction
+from extreme_estimator.extreme_models.margin_model.margin_function.parametric_margin_function import \
+    ParametricMarginFunction
 from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
 from extreme_estimator.extreme_models.margin_model.param_function.param_function import ConstantParamFunction, \
     ParamFunction, LinearParamFunction
@@ -9,7 +9,7 @@ from extreme_estimator.margin_fits.gev.gev_params import GevParams
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
-class LinearMarginFunction(IndependentMarginFunction):
+class LinearMarginFunction(ParametricMarginFunction):
     """ Margin Function, where each parameter can augment linearly along any dimension.
 
         dim = 0 correspond to the intercept
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/parametric_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/parametric_margin_function.py
new file mode 100644
index 0000000000000000000000000000000000000000..0305bcfbfad0f7e69ab458fa508ebe830c59d654
--- /dev/null
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/parametric_margin_function.py
@@ -0,0 +1,15 @@
+from typing import Dict
+
+from extreme_estimator.extreme_models.margin_model.margin_function.independent_margin_function import \
+    IndependentMarginFunction
+
+
+class ParametricMarginFunction(IndependentMarginFunction):
+
+    @property
+    def form_dict(self) -> Dict[str, str]:
+        raise NotImplementedError
+
+    @property
+    def coef_dict(self) -> Dict[str, float]:
+        raise NotImplementedError
\ No newline at end of file
diff --git a/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py b/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py
new file mode 100644
index 0000000000000000000000000000000000000000..b26d7efb60d5564d2e5787a9a43fa8925639b659
--- /dev/null
+++ b/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py
@@ -0,0 +1,41 @@
+from abc import ABC
+
+import numpy as np
+import pandas as pd
+
+from extreme_estimator.extreme_models.margin_model.margin_function.parametric_margin_function import \
+    ParametricMarginFunction
+from extreme_estimator.extreme_models.result_from_fit import ResultFromFit
+from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
+from extreme_estimator.extreme_models.utils import safe_run_r_estimator, r, get_coord, \
+    get_margin_formula
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+
+
+class ParametricMarginModel(AbstractMarginModel, ABC):
+
+    def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
+                 params_sample=None):
+        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)
+
+    def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, df_coordinates_spat: pd.DataFrame,
+                                  df_coordinates_temp: pd.DataFrame) -> ResultFromFit:
+        # The reshaping on the line below is only valid if we have a single observation per spatio-temporal point
+        if maxima_gev.shape[1] == 1:
+            maxima_gev = maxima_gev.reshape([len(df_coordinates_temp), len(df_coordinates_spat)])
+        data = np.transpose(maxima_gev)
+
+        fit_params = get_margin_formula(self.margin_function_start_fit.form_dict)
+
+        # Covariables
+        covariables = get_coord(df_coordinates=df_coordinates_spat)
+        fit_params['temp.cov'] = get_coord(df_coordinates=df_coordinates_temp)
+
+        # Start parameters
+        coef_dict = self.margin_function_start_fit.coef_dict
+        fit_params['start'] = r.list(**coef_dict)
+
+        return safe_run_r_estimator(function=r.fitspatgev, use_start=self.use_start_value, data=data,
+                                    covariables=covariables, **fit_params)
diff --git a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
deleted file mode 100644
index 1e77a4ea42be0516f77e79603284081dadf3f30c..0000000000000000000000000000000000000000
--- a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
+++ /dev/null
@@ -1,127 +0,0 @@
-import numpy as np
-import pandas as pd
-
-from extreme_estimator.extreme_models.result_from_fit import ResultFromFit
-from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
-from extreme_estimator.extreme_models.margin_model.margin_function.linear_margin_function import LinearMarginFunction
-from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
-from extreme_estimator.extreme_models.utils import safe_run_r_estimator, r, get_coord, \
-    get_margin_formula
-from extreme_estimator.margin_fits.gev.gev_params import GevParams
-
-
-class LinearMarginModel(AbstractMarginModel):
-
-    def load_margin_functions(self, gev_param_name_to_linear_dims=None):
-        assert gev_param_name_to_linear_dims is not None, 'LinearMarginModel cannot be used for sampling/fitting \n' \
-                                                          'load_margin_functions needs to be implemented in child class'
-        # Load sample coef
-        self.default_params_sample = self.default_param_name_and_dim_to_coef()
-        linear_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_linear_coef=linear_coef_sample,
-                                                           gev_param_name_to_linear_dims=gev_param_name_to_linear_dims)
-
-        # Load start fit coef
-        self.default_params_start_fit = self.default_param_name_and_dim_to_coef()
-        linear_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_linear_coef=linear_coef_start_fit,
-                                                              gev_param_name_to_linear_dims=gev_param_name_to_linear_dims)
-
-    @staticmethod
-    def default_param_name_and_dim_to_coef() -> dict:
-        default_intercept = 1
-        default_slope = 0.01
-        gev_param_name_and_dim_to_coef = {}
-        for gev_param_name in GevParams.PARAM_NAMES:
-            gev_param_name_and_dim_to_coef[(gev_param_name, 0)] = default_intercept
-            for dim in [1, 2, 3]:
-                gev_param_name_and_dim_to_coef[(gev_param_name, dim)] = default_slope
-        return gev_param_name_and_dim_to_coef
-
-    @staticmethod
-    def gev_param_name_to_linear_coef(param_name_and_dim_to_coef):
-        gev_param_name_to_linear_coef = {}
-        for gev_param_name in GevParams.PARAM_NAMES:
-            dim_to_coef = {dim: param_name_and_dim_to_coef[(gev_param_name, dim)] for dim in [0, 1, 2, 3]}
-            linear_coef = LinearCoef(gev_param_name=gev_param_name, dim_to_coef=dim_to_coef)
-            gev_param_name_to_linear_coef[gev_param_name] = linear_coef
-        return gev_param_name_to_linear_coef
-
-    @classmethod
-    def from_coef_list(cls, coordinates, gev_param_name_to_coef_list):
-        params = {}
-        for gev_param_name in GevParams.PARAM_NAMES:
-            for dim, coef in enumerate(gev_param_name_to_coef_list[gev_param_name]):
-                params[(gev_param_name, dim)] = coef
-        return cls(coordinates, params_sample=params, params_start_fit=params)
-
-    def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, df_coordinates_spat: pd.DataFrame,
-                                  df_coordinates_temp: pd.DataFrame) -> ResultFromFit:
-        # The reshaping on the line below is only valid if we have a single observation per spatio-temporal point
-        if maxima_gev.shape[1] == 1:
-            maxima_gev = maxima_gev.reshape([len(df_coordinates_temp), len(df_coordinates_spat)])
-        data = np.transpose(maxima_gev)
-
-        fit_params = get_margin_formula(self.margin_function_start_fit.form_dict)
-
-        # Covariables
-        covariables = get_coord(df_coordinates=df_coordinates_spat)
-        fit_params['temp.cov'] = get_coord(df_coordinates=df_coordinates_temp)
-
-        # Start parameters
-        coef_dict = self.margin_function_start_fit.coef_dict
-        fit_params['start'] = r.list(**coef_dict)
-
-        return safe_run_r_estimator(function=r.fitspatgev, use_start=self.use_start_value, data=data,
-                                    covariables=covariables, **fit_params)
-
-
-class ConstantMarginModel(LinearMarginModel):
-
-    def load_margin_functions(self, gev_param_name_to_linear_dims=None):
-        super().load_margin_functions({})
-
-
-class LinearShapeDim1MarginModel(LinearMarginModel):
-
-    def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None):
-        super().load_margin_functions({GevParams.SHAPE: [1]})
-
-
-class LinearScaleDim1MarginModel(LinearMarginModel):
-
-    def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None):
-        super().load_margin_functions({GevParams.SCALE: [1]})
-
-
-class LinearShapeDim1and2MarginModel(LinearMarginModel):
-
-    def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None):
-        super().load_margin_functions({GevParams.SHAPE: [1, 2]})
-
-
-class LinearAllParametersDim1MarginModel(LinearMarginModel):
-
-    def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None):
-        super().load_margin_functions({GevParams.SHAPE: [1],
-                                       GevParams.LOC: [1],
-                                       GevParams.SCALE: [1]})
-
-
-class LinearMarginModelExample(LinearMarginModel):
-
-    def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None):
-        super().load_margin_functions({GevParams.SHAPE: [1],
-                                       GevParams.LOC: [2],
-                                       GevParams.SCALE: [1]})
-
-
-class LinearAllParametersAllDimsMarginModel(LinearMarginModel):
-
-    def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None):
-        all_dims = list(range(1, self.coordinates.nb_coordinates + 1))
-        super().load_margin_functions({GevParams.SHAPE: all_dims.copy(),
-                                       GevParams.LOC: all_dims.copy(),
-                                       GevParams.SCALE: all_dims.copy()})
diff --git a/test/test_extreme_estimator/test_extreme_models/test_margin_model.py b/test/test_extreme_estimator/test_extreme_models/test_margin_model.py
index 552a4b2419dcfc3b44ef4940f6a946f27748893d..5980064c76c6fa338443e8aec27de97b8f006e5b 100644
--- a/test/test_extreme_estimator/test_extreme_models/test_margin_model.py
+++ b/test/test_extreme_estimator/test_extreme_models/test_margin_model.py
@@ -4,6 +4,7 @@ import unittest
 from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
     AbstractMarginFunction
 from extreme_estimator.extreme_models.margin_model.margin_function.linear_margin_function import LinearMarginFunction
+from extreme_estimator.extreme_models.margin_model.spline_margin_model import SplineMarginModel
 from extreme_estimator.margin_fits.gev.gev_params import GevParams
 from extreme_estimator.extreme_models.margin_model.linear_margin_model import LinearShapeDim1MarginModel, \
     LinearAllParametersAllDimsMarginModel
@@ -14,15 +15,17 @@ from spatio_temporal_dataset.coordinates.spatial_coordinates.coordinates_1D impo
 from test.test_utils import load_test_spatiotemporal_coordinates
 
 
-class TestVisualizationMarginModel(unittest.TestCase):
+class TestVisualizationLinearMarginModel(unittest.TestCase):
     DISPLAY = False
     nb_points = 2
-    margin_model_class = [LinearShapeDim1MarginModel, LinearAllParametersAllDimsMarginModel][-1]
+    margin_model_class = LinearAllParametersAllDimsMarginModel
 
     def tearDown(self) -> None:
         self.margin_model.margin_function_sample.visualize_function(show=self.DISPLAY)
         self.assertTrue(True)
 
+
+
     def test_example_visualization_1D(self):
         coordinates = LinSpaceSpatialCoordinates.from_nb_points(nb_points=self.nb_points)
         self.margin_model = self.margin_model_class(coordinates=coordinates, params_sample={(GevParams.SHAPE, 1): 0.02})
@@ -56,6 +59,27 @@ class TestVisualizationMarginModel(unittest.TestCase):
     #     # Load
 
 
+# class TestVisualizationSplineMarginModel(TestVisualizationMarginModel):
+#     margin_model_class = SplineMarginModel
+#
+#     def tearDown(self) -> None:
+#         self.margin_model.margin_function_sample.visualize_function(show=self.DISPLAY)
+#         self.assertTrue(True)
+
+    # def test_example_visualization_1D_spline(self):
+    #     coordinates = LinSpaceSpatialCoordinates.from_nb_points(nb_points=self.nb_points)
+    #     self.margin_model = self.margin_model_class(coordinates=coordinates, params_sample={(GevParams.SHAPE, 1): 0.02})
+
+    # def test_example_visualization_2D_spatial_spline(self):
+    #     spatial_coordinates = LinSpaceSpatial2DCoordinates.from_nb_points(nb_points=self.nb_points)
+    #     self.margin_model = self.margin_model_class(coordinates=spatial_coordinates)
+    #     # Assert that the grid correspond to what we expect in a simple case
+    #     AbstractMarginFunction.VISUALIZATION_RESOLUTION = 2
+    #     grid = self.margin_model.margin_function_sample.grid_2D['loc']
+    #     true_grid = np.array([[0.98, 1.0], [1.0, 1.02]])
+    #     self.assertTrue((grid == true_grid).all(), msg="\nexpected:\n{}, \nfound:\n{}".format(true_grid, grid))
+
+
 if __name__ == '__main__':
     unittest.main()
     # v = TestVisualizationMarginModel()