Repository: mmaelicke/scikit-gstat

Issue:
Fit of a sum of variogram models
Hi again @mmaelicke,

Do you think allowing the fit of a combinaison (sum, product) of variogram models could be a functionality that has its place in `scikit-gstat`?

We use that a lot for analyzing multiple correlation ranges in satellite data (for example, https://github.com/GlacioHack/xdem/blob/main/xdem/spatialstats.py#L1622).
I could try to add the functionality here! :smile: 

It would combine well with `RasterMetricSpace` (for which I need to push some minor fixes + add better documentation, at least for the API, it's been a long time coming). 
We also tend to sample multiple experimental variograms to get an idea of the standard error in the experimental variance at each lag (as we have billions of pairwise samples), and use that to condition the fit (by passing it automatically as `fit_sigma`). Not sure how that'd fit here...

What do you think?

Discussion:

Patch:
diff --git a/docs/userguide/variogram.rst b/docs/userguide/variogram.rst
index bb1d5d8..b195c48 100644
--- a/docs/userguide/variogram.rst
+++ b/docs/userguide/variogram.rst
@@ -766,6 +766,53 @@ variogram is rather showing a Gaussian or exponential behavior.
     If you would like to export a Variogram instance to gstools, the smoothness parameter
     may not be smaller than ``0.2``.
 
+Sum of models
+~~~~~~~~~~~~~
+
+All the above models can be combined into a sum of models, in order to fit more complex empirical variograms that might
+capture signals with multiple correlations ranges.
+To fit with a sum of models, simply add a "+" between any number of models:
+
+.. ipython:: python
+    :okwarning:
+
+    # Fit a sum of two spherical models
+    V.model = "spherical+spherical"
+
+    @savefig sum_of_models.png width=8in
+    V.plot();
+
+Custom model
+~~~~~~~~~~~~
+
+Additionally, any custom model can be passed to the model argument to fit a specific function to the empirical
+variogram. This model needs to be built with the ``@variogram`` decorator, which requires the lag argument
+(typically, ``h``) to be listed first in the custom function.
+
+.. caution::
+
+    When using a custom model, it is highly recommended to define the parameter bounds as these cannot be derived
+    logically by SciKit-GStat as for other models. This is done by passing either ``fit_bounds`` to ``Variogram()``
+    at instantiation or ``bounds`` to ``Variogram.fit()``.
+
+
+.. ipython:: python
+    :okwarning:
+
+    # Build a custom model by applying the @variogram decorator (here adding a linear term to a spherical model)
+    from skgstat.models import variogram, spherical
+    @variogram
+    def custom_model(h, r1, c1, a):
+        return spherical(h, r1, c1) + h * a
+    V.model = custom_model
+
+    # We define the bounds for r1, c1 and a
+    bounds_custom = [(0, 0, 0), (np.max(V.bins), np.max(V.experimental), 2)]
+    V.fit(bounds=bounds_custom)
+
+    @savefig custom_model.png width=8in
+    V.plot();
+
 When direction matters
 ======================
 
diff --git a/requirements.txt b/requirements.txt
index 852de3e..4d2fa77 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,4 @@
-scipy
+scipy<=1.11.1
 numpy
 pandas
 matplotlib
diff --git a/skgstat/DirectionalVariogram.py b/skgstat/DirectionalVariogram.py
index f4aaa59..9d5c742 100644
--- a/skgstat/DirectionalVariogram.py
+++ b/skgstat/DirectionalVariogram.py
@@ -294,6 +294,7 @@ class DirectionalVariogram(Variogram):
 
         # set the directional model
         self._directional_model = None
+        self._is_model_custom = False
         self.set_directional_model(model_name=directional_model)
 
         # the binning settings
diff --git a/skgstat/Variogram.py b/skgstat/Variogram.py
index e63b3cd..d2c1176 100644
--- a/skgstat/Variogram.py
+++ b/skgstat/Variogram.py
@@ -2,6 +2,7 @@
 Variogram class
 """
 import copy
+import inspect
 import warnings
 from typing import Iterable, Callable, Union, Tuple
 
@@ -91,8 +92,12 @@ class Variogram(object):
             differences, aligned to the 1D distance matrix (flattened upper
             triangle) and return a scalar, that converges towards small
             values for similarity (high covariance).
-        model : str
-            String identifying the theoretical variogram function to be used
+        model : str | Callable
+            .. versionchanged:: 1.0.12
+                Added support for sum of models (e.g., "spherical+gaussian"), or custom model (Callable). Using
+                `fit_bounds` to optimize the fit is recommended for custom models, and can be useful for sum of models.
+
+            String or callable identifying the theoretical variogram function to be used
             to describe the experimental variogram. Can be one of:
 
               * spherical       [Spherical, default]
@@ -103,6 +108,9 @@ class Variogram(object):
               * matern          [Matérn model]
               * nugget          [nugget effect variogram]
 
+            Any number of these theoretical models can be summed using "+" iteratively, e.g. "spherical+cubic+matern".
+            The nugget parameters of the models are removed except for the last model (sum of nuggets = single nugget).
+
         dist_func : str
             String identifying the distance function. Defaults to
             'euclidean'. Can be any metric accepted by
@@ -185,7 +193,8 @@ class Variogram(object):
             Defaults to False. If True, a nugget effet will be added to all
             Variogram.models as a third (or fourth) fitting parameter. A
             nugget is essentially the y-axis interception of the theoretical
-            variogram function.
+            variogram function. For a sum of variogram, the nugget is defined
+            in its last model.
         maxlag : float, str
             Can specify the maximum lag distance directly by giving a value
             larger than 1. The binning function will not find any lag class
@@ -245,6 +254,24 @@ class Variogram(object):
             If present, the plot will indicate the confidence interval as
             error bars around the experimental variogram.
 
+        fit_bounds: 2-tuple of array_like or Bounds, optional
+            .. versionadded:: 1.0.12
+
+            Lower and upper bounds on parameters passed to scipy.optimize.curve_fit.
+
+            Order is typically (range, sill, nugget) or (range, sill, smoothness, nugget) for individual models, or
+            (range1, sill1, nugget1, range2, sill2, nugget2) for a sum of 2 models.
+            Recommended for custom models, where bounds cannot be determined logically.
+            For internal models, defaults to known min/max values for the sill (0, max variance), range (0, max lag)
+            and smoothness (0, 2) or (0, 20) for stable and matern, respectively.
+
+        fit_p0: array_like, optional
+            .. versionadded:: 1.0.12
+
+            Initial guess for the parameters passed to scipy.optimize.curve_fit.
+
+            Same order as for fit_bounds.
+            Defaults to upper bounds values. For custom models, if no bounds are defined, defaults to 1.
         """
         # Before we do anything else, make kwargs available
         self._kwargs = self._validate_kwargs(**kwargs)
@@ -333,17 +360,19 @@ class Variogram(object):
         if fit_method is not None:
             self.preprocessing(force=True)
 
+        # set if nugget effect shall be used
+        self._use_nugget = None
+        self.use_nugget = use_nugget
+
         # model can be a function or a string
         self._model = None
+        self._model_name = None
+        self._is_model_custom = False
         self.set_model(model_name=model)
 
         # specify if the lag should be given absolute or relative to the maxlag
         self._normalized = normalize
 
-        # set if nugget effect shall be used
-        self._use_nugget = None
-        self.use_nugget = use_nugget
-
         # set the fitting method and sigma array
         self._fit_method = None
         if fit_method is not None:
@@ -360,7 +389,9 @@ class Variogram(object):
 
         # do the preprocessing and fitting upon initialization
         # Note that fit() calls preprocessing
-        self.fit(force=True)
+        fit_bounds = self._kwargs.get('fit_bounds') # returns None if empty
+        fit_p0 = self._kwargs.get('fit_p0')
+        self.fit(force=True, bounds=fit_bounds, p0=fit_p0)
 
         # finally check if any of the uncertainty propagation kwargs are set
         self._experimental_conf_interval = None
@@ -961,10 +992,15 @@ class Variogram(object):
             # at first reset harmonize
             self._harmonize = False
             if model_name.lower() == 'harmonize':
+                mname = 'harmonize'
                 self._harmonize = True
                 self._model = self._build_harmonized_model()
+            elif "+" in model_name:
+                mname =  ''.join(model_name.split()).lower() # Remove all whitespaces
+                self._model = self._build_sum_models(mname)
             elif hasattr(models, model_name.lower()):
-                self._model = getattr(models, model_name.lower())
+                mname = model_name.lower()
+                self._model = getattr(models, mname)
             else:
                 raise ValueError(
                     (
@@ -972,8 +1008,68 @@ class Variogram(object):
                         ' understood, please provide the function'
                     ) % model_name
                 )
+            # Set model name attribute
+            self._model_name = mname
+
         else:  # pragma: no cover
+            self._is_model_custom = True
             self._model = model_name
+            self._model_name = model_name.__name__
+
+    def _get_argpos_sum_models(self, list_model_names):
+        """
+        Get argument slice position (list of slices) for the sum of models from a list of model names (list of strings).
+        """
+
+        # Doing this here for other functions (fit, describe, etc), even though already done in _build_sum_models
+        list_models = [getattr(models, model_name.lower()).py_func for model_name in list_model_names]
+
+        # Get the number of arguments per model (e.g., [3, 4, 4])
+        nb_args_per_model = np.array([len(inspect.getfullargspec(model).args) for model in list_models])
+
+        # We remove the nugget and lags parameters, except nugget for the last model (sum of nuggets = single nugget)
+        nb_args_per_model -= 2
+        nb_args_per_model[-1] += 1
+        # Compute cumulative number of args removing 2 args everywhere (all lags and nuggets, last one will compensate)
+        cum_args_minus_lag = np.cumsum(nb_args_per_model)
+
+        # Prepare argument slices to distribute to submodels
+        args_indices = np.insert(cum_args_minus_lag, 0, 0)  # We add the first indice of 0
+        args_slices = [(slice(args_indices[i], args_indices[i + 1])) for i in range(len(args_indices) - 1)]
+
+        return args_slices
+
+    def _build_sum_models(self, sum_models_name: str):
+        """
+        Build sum of theoretical models, variogram-decorated function.
+        """
+
+        # Remove all whitespaces in the string, in case the user wrote something like "spherical + gaussian"
+        sum_models_name = ''.join(sum_models_name.split()).lower()
+
+        # Get individual model names
+        list_model_names = sum_models_name.split("+")
+
+        # Check that all models exist in the "models" module
+        if not all(hasattr(models, model_name.lower()) for model_name in list_model_names):
+            raise ValueError(
+                (
+                    'One of the theoretical models in the list "%s" is not'
+                    ' understood, please provide existing model names separated by "+".'
+                ) % ", ".join(list_model_names)
+            )
+
+        # First, build the models from their py_func (ignoring variogram decorator) and get args per model
+        list_models = [getattr(models, model_name.lower()).py_func for model_name in list_model_names]
+        # Get the argument positions for the model sum (function uses model names to be called by other methods)
+        args_slices = self._get_argpos_sum_models(list_model_names=list_model_names)
+
+        # Distribute first argument (lag) and use all others in order (nugget ignored when last argument not passed)
+        @models.variogram
+        def sum_models(h, *args):
+            return sum(list_models[i](h, *args[args_slices[i]]) for i in range(len(list_models)))
+
+        return sum_models
 
     def _build_harmonized_model(self):
         x = self.bins
@@ -1454,7 +1550,7 @@ class Variogram(object):
         self._calc_diff(force=force)
         self._calc_groups(force=force)
 
-    def fit(self, force=False, method=None, sigma=None, **kwargs):
+    def fit(self, force=False, method=None, sigma=None, bounds=None, p0=None, **kwargs):
         """Fit the variogram
 
         The fit function will fit the theoretical variogram function to the
@@ -1503,6 +1599,21 @@ class Variogram(object):
             Uncertainty array for the bins. Has to have the same dimension as
             self.bins. Refer to Variogram.fit_sigma for more information.
 
+        bounds: 2-tuple of array_like or Bounds, optional
+            Lower and upper bounds on parameters passed to scipy.optimize.curve_fit.
+
+            Order is typically (range, sill, nugget) or (range, sill, smoothness, nugget) for individual models, or
+            (range1, sill1, nugget1, range2, sill2, nugget2) for a sum of 2 models.
+            Recommended for custom models, where bounds cannot be determined logically.
+            For internal models, defaults to known min/max values for the sill (0, max variance), range (0, max lag)
+            and smoothness (0, 2) or (0, 20) for stable and matern, respectively.
+
+        p0: array_like, optional
+            Initial guess for the parameters passed to scipy.optimize.curve_fit.
+
+            Same order as for fit_bounds.
+            Defaults to upper bounds values. For custom models, if no bounds are defined, defaults to 1.
+
         Returns
         -------
         void
@@ -1563,18 +1674,41 @@ class Variogram(object):
             self.cof = [r, s, n]
             return
 
-        # Switch the method
-        # wrap the model to include or exclude the nugget
-        if self.use_nugget:
-            def wrapped(*args):
-                return self._model(*args)
+        # For a supported model, wrap the function depending on nugget and get logical bounds
+        if not self._is_model_custom:
+            # Switch the method
+            # wrap the model to include or exclude the nugget
+            if self.use_nugget:
+                def wrapped(*args):
+                    return self._model(*args)
+            else:
+                def wrapped(*args):
+                    return self._model(*args, 0)
+
+            # get p0
+            if bounds is None:
+                bounds = (0, self.__get_fit_bounds(x, y))
+            if p0 is None:
+                p0 = np.asarray(bounds[1])
+        # Else, inspect the function for the number of arguments
         else:
-            def wrapped(*args):
-                return self._model(*args, 0)
+            # The number of arguments of argspec minus one is what we initialized
+            argspec = inspect.getfullargspec(self._model.__wrapped__)
+            nb_args = len(argspec.args) - 1
+            if bounds is None:
+                warnings.warn("Parameter bounds cannot be logically derived for a custom model during the fit. "
+                              "User bounds can be passed using fit(..., bounds=) or Variogram(..., fit_bounds=).", UserWarning)
+                bounds = ([-np.inf] * nb_args, [np.inf] * nb_args)
+            if p0 is None:
+                # If all bounds are infinite (not defined, pass 1)
+                if bounds == ([-np.inf] * nb_args, [np.inf] * nb_args):
+                    p0 = np.ones(nb_args)
+                # Else pass the upper bounds
+                else:
+                    p0 = np.asarray(bounds[1])
 
-        # get p0
-        bounds = (0, self.__get_fit_bounds(x, y))
-        p0 = np.asarray(bounds[1])
+            def wrapped(*args):
+                return self._model(*args)
 
         # Trust Region Reflective
         if self.fit_method == 'trf':
@@ -1644,10 +1778,10 @@ class Variogram(object):
             n = kwargs.get('nugget', self._kwargs.get('fit_nugget', old_params.get('nugget', 0.0)))
 
             # check if a s parameter is needed
-            if self._model.__name__ in ('stable', 'matern'):
-                if self._model.__name__ == 'stable':
+            if self._model_name in ('stable', 'matern'):
+                if self._model_name == 'stable':
                     s2 = kwargs.get('shape', self._kwargs.get('fit_shape', old_params.get('shape',2.0)))
-                if self._model.__name__ == 'matern':
+                if self._model_name == 'matern':
                     s2 = kwargs.get('shape', self._kwargs.get('fit_shape', old_params.get('smoothness', 2.0)))
 
                 # set
@@ -1968,33 +2102,45 @@ class Variogram(object):
         list
 
         """
-        mname = self._model.__name__
+        all_mname = self._model_name
 
-        # use range, sill and smoothness parameter
-        if mname == 'matern':
-            # a is max(x), C0 is max(y) s is limited to 20?
-            bounds = [np.nanmax(x), np.nanmax(y), 20.]
+        # for a sum of models, create a list
+        if "+" in all_mname:
+            list_mname = all_mname.split("+")
+        else:
+            list_mname = [all_mname]
 
-        # use range, sill and shape parameter
-        elif mname == 'stable':
-            # a is max(x), C0 is max(y) s is limited to 2?
-            bounds = [np.nanmax(x), np.nanmax(y), 2.]
+        # we append all bounds (for one or several models)
+        all_bounds = []
+        for i, mname in enumerate(list_mname):
 
-        # use only sill
-        elif mname == 'nugget':
-            # a is max(x):
-            bounds = [np.nanmax(x)]
+            # use range, sill and smoothness parameter
+            if mname == 'matern':
+                # a is max(x), C0 is max(y) s is limited to 20?
+                bounds = [np.nanmax(x), np.nanmax(y), 20.]
 
-        # use range and sill
-        else:
-            # a is max(x), C0 is max(y)
-            bounds = [np.nanmax(x), np.nanmax(y)]
+            # use range, sill and shape parameter
+            elif mname == 'stable':
+                # a is max(x), C0 is max(y) s is limited to 2?
+                bounds = [np.nanmax(x), np.nanmax(y), 2.]
+
+            # use only sill
+            elif mname == 'nugget':
+                # a is max(x):
+                bounds = [np.nanmax(x)]
+
+            # use range and sill
+            else:
+                # a is max(x), C0 is max(y)
+                bounds = [np.nanmax(x), np.nanmax(y)]
+
+            # if use_nugget is True add the nugget (for the last model only in case it is a sum)
+            if self.use_nugget and i == (len(list_mname) - 1):
+                bounds.append(0.99*np.nanmax(y))
 
-        # if use_nugget is True add the nugget
-        if self.use_nugget:
-            bounds.append(0.99*np.nanmax(y))
+            all_bounds += bounds
 
-        return bounds
+        return all_bounds
 
     def data(self, n=100, force=False):
         """Theoretical variogram function
@@ -2464,29 +2610,63 @@ class Variogram(object):
         rdict = dict(
             model=fnname(self._model) if not self._harmonize else "harmonize",
             estimator=fnname(self._estimator),
-            dist_func=fnname(self.dist_function),
+            dist_func=fnname(self.dist_function))
 
-            normalized_effective_range=cof[0] * maxlag,
-            normalized_sill=cof[1] * maxvar,
-            normalized_nugget=cof[-1] * maxvar if self.use_nugget else 0,
-
-            effective_range=cof[0],
-            sill=cof[1],
-            nugget=cof[-1] if self.use_nugget else 0,
-        )
-
-        # handle s parameters for matern and stable model
-        if self._model.__name__ == 'matern':
-            rdict['smoothness'] = cof[2]
-        elif self._model.__name__ == 'stable':
-            rdict['shape'] = cof[2]
+        def create_dict_for_model(model_name, cof, maxlag, maxvar, use_nugget, id = ''):
+            """
+            Create dictionary of parameters for an individual model.
+            Optionally, append a model ID to the key (for the case of summed models).
+            """
+            # id appended to the param key name to differentiate models for a sum
+            if id != '':
+                id = '_' + id
+
+            d = {
+                'normalized_effective_range' + id: cof[0] * maxlag,
+                'normalized_sill' + id: cof[1] * maxvar,
+                'normalized_nugget' + id: cof[-1] * maxvar if use_nugget else 0,
+                'effective_range' + id: cof[0],
+                'sill' + id: cof[1],
+                'nugget' + id: cof[-1] if use_nugget else 0,
+            }
+
+            # handle s parameters for matern and stable model
+            if model_name == 'matern':
+                d['smoothness' + id] = cof[2]
+            elif model_name == 'stable':
+                d['shape' + id] = cof[2]
+
+            return d
+
+        # for a custom model: we list the optimized params for the function wrapped by the variogram decorator
+        if self._is_model_custom:
+            custom_arg_names = inspect.getfullargspec(self._model.__wrapped__).args
+            all_params = {"param"+str(i+1)+"_"+custom_arg_names[i+1]: cof[i] for i in range(len(custom_arg_names) - 1)}
+
+        # for a sum of models
+        elif "+" in self._model_name:
+            list_model_names = self._model_name.split("+")
+            list_argslices = self._get_argpos_sum_models(list_model_names=list_model_names)
+            all_params = {}
+            # add the parameters for each model, with parameter suffix from 1 to the total number
+            for i in range(len(list_model_names)):
+                model_params = create_dict_for_model(model_name=list_model_names[i], cof=cof[list_argslices[i]],
+                                                  maxlag=maxlag, maxvar=maxvar, use_nugget=self.use_nugget, id=str(i+1))
+                all_params.update(model_params)
+
+        # for a single model
+        else:
+            all_params = create_dict_for_model(model_name=self._model_name, cof=cof, maxlag=maxlag, maxvar=maxvar,
+                                                use_nugget=self.use_nugget)
+        # update dictionary
+        rdict.update(all_params)
 
         # add other stuff if not short version requested
         if not short:
             kwargs = self._kwargs
             params = dict(
                 estimator=self._estimator.__name__,
-                model=self._model.__name__,
+                model=self._model_name,
                 dist_func=str(self.dist_function),
                 bin_func=self._bin_func_name,
                 normalize=self.normalized,
@@ -2521,39 +2701,67 @@ class Variogram(object):
         params : list
             [range, sill, nugget] for most models and
             [range, sill, shape, nugget] for matern and stable model.
+            [range1, sill1, nugget1, range2, still2, nugget2] for a sum of 2 models.
+            [param1, param2, param3, ...] in order for a custom model.
 
         """
         d = self.describe()
         if 'error' in d:
             return [None, None, None]
-        elif self._model.__name__ == 'matern':
-            return list([
-                d['effective_range'],
-                d['sill'],
-                d['smoothness'],
-                d['nugget']
-            ])
-        elif self._model.__name__ == 'stable':
-            return list([
-                d['effective_range'],
-                d['sill'],
-                d['shape'],
-                d['nugget']
-            ])
-        elif self._model.__name__ == 'nugget':
-            return list([d['nugget']])
+
+        def get_params_list_from_dict(d, model_name, id=''):
+
+            # ID used to differentiate params for a sum of models
+            if id != '':
+                id = '_'+id
+
+            # Get specific smoothness and shape parameters for matern and stable
+            if model_name == 'matern':
+                return list([
+                    d['effective_range'+id],
+                    d['sill'+id],
+                    d['smoothness'+id],
+                    d['nugget'+id]
+                ])
+            elif model_name == 'stable':
+                return list([
+                    d['effective_range'+id],
+                    d['sill'+id],
+                    d['shape'+id],
+                    d['nugget'+id]
+                ])
+            # Get nugget for only-nugget
+            elif model_name == 'nugget':
+                return list([d['nugget'+id]])
+            # Or get classic parameters
+            else:
+                return list([
+                    d['effective_range'+id],
+                    d['sill'+id],
+                    d['nugget'+id]
+                ])
+
+        # For a custom model, just pass cof in order
+        if self._is_model_custom:
+            list_params = self.cof
+        # Get parameters for a sum of models
+        elif '+' in self._model_name:
+            list_model_names = self._model_name.split('+')
+            list_params = []
+            for i in range(len(list_model_names)):
+                params = get_params_list_from_dict(d, model_name=list_model_names[i], id=str(i+1))
+                list_params += params
+        # Or for a single model
         else:
-            return list([
-                d['effective_range'],
-                d['sill'],
-                d['nugget']
-            ])
+            list_params = get_params_list_from_dict(d, model_name=self._model_name)
+
+        return list_params
 
     def to_DataFrame(self, n=100, force=False):
         """Variogram DataFrame
 
         Returns the fitted theoretical variogram as a pandas.DataFrame
-        instance. The n and force parameter control the calaculation,
+        instance. The n and force parameter control the calculation,
         refer to the data function for more info.
 
         .. deprecated:: 1.0.10
@@ -2583,7 +2791,7 @@ class Variogram(object):
 
         return DataFrame({
             'lags': lags,
-            self._model.__name__: data}
+            self._model_name: data}
         ).copy()
 
     def to_gstools(self, **kwargs):
@@ -2909,7 +3117,7 @@ class Variogram(object):
         :return:
         """
         try:
-            _name = self._model.__name__
+            _name = self._model_name
             _b = int(len(self.bins))
         except Exception:
             return "< abstract Variogram >"


Test patch:
diff --git a/skgstat/tests/test_models.py b/skgstat/tests/test_models.py
index 1d3aa17..8ba23a5 100644
--- a/skgstat/tests/test_models.py
+++ b/skgstat/tests/test_models.py
@@ -186,6 +186,21 @@ class TestVariogramDecorator(unittest.TestCase):
         for r, c in zip(res, adder([1, 4, 8], 4)):
             self.assertEqual(r, c)
 
+    def test_sum_spherical(self):
+        @variogram
+        def sum_spherical(h, r1, c1, r2, c2, b1=0, b2=0):
+            return spherical(h, r1, c1, b1) + spherical(h, r2, c2, b2)
+
+        # Parameters for the two spherical models
+        params = [1, 0.3, 10, 0.7]
+
+        # Values at which we'll evaluate the function and its expected result
+        vals = [0, 1, 100]
+        res = [0, 0.3 + spherical(1, 10, 0.7), 1]
+
+        for r, c in zip(res, sum_spherical(vals, *params)):
+            self.assertEqual(r, c)
+
 
 if __name__=='__main__':
     unittest.main()
diff --git a/skgstat/tests/test_variogram.py b/skgstat/tests/test_variogram.py
index 3c8e9e8..c4271ee 100644
--- a/skgstat/tests/test_variogram.py
+++ b/skgstat/tests/test_variogram.py
@@ -19,6 +19,7 @@ from skgstat import Variogram, DirectionalVariogram
 from skgstat import OrdinaryKriging
 from skgstat import estimators
 from skgstat import plotting
+from skgstat.models import variogram, spherical, gaussian, exponential, cubic, stable, matern
 
 
 class TestSpatiallyCorrelatedData(unittest.TestCase):
@@ -61,7 +62,7 @@ class TestSpatiallyCorrelatedData(unittest.TestCase):
             self.assertAlmostEqual(x, y, places=3)
 
 
-class TestVariogramInstatiation(unittest.TestCase):
+class TestVariogramInstantiation(unittest.TestCase):
     def setUp(self):
         # set up default values, whenever c and v are not important
         np.random.seed(42)
@@ -254,7 +255,7 @@ class TestVariogramArguments(unittest.TestCase):
         assert_array_almost_equal(
             V.bins,
             np.array([4.3, 8.4, 12.8, 17.1, 21.4, 25.2, 29.9, 33.2, 38.5, 42.8]),
-            decimal=1
+            decimal=0
         )
 
     def test_binning_method_stable_maxiter(self):
@@ -264,7 +265,7 @@ class TestVariogramArguments(unittest.TestCase):
         assert_array_almost_equal(
             V.bins,
             np.array([4.3, 8.4, 12.8, 17.1, 21.4, 25.2, 29.9, 33.2, 38.5, 42.8]),
-            decimal=1
+            decimal=0
         )
 
     def test_binning_method_stable_fix_bins(self):
@@ -279,7 +280,7 @@ class TestVariogramArguments(unittest.TestCase):
         assert_array_almost_equal(
             V.bins,
             np.array([4.2, 8.6, 12.8, 17.1, 21.2, 25.5, 29.3, 33.2, 37.4, 43.]),
-            decimal=1
+            decimal=0
         )
 
     def test_binning_change_nlags(self):
@@ -633,6 +634,39 @@ class TestVariogramArguments(unittest.TestCase):
         assert V.cov is None
         assert V.cof is None
 
+    def test_model(self):
+        """
+        Test that all types of models instantiate properly
+        (to complement test_set_model() that only checks already instantiated vario)
+        """
+
+        # Individual model
+        for model_name in ['spherical', 'gaussian', 'exponential', 'cubic', 'matern', 'stable']:
+            V = Variogram(self.c, self.v, model=model_name)
+            assert V._model_name == model_name
+            assert V._model == globals()[model_name]
+            assert V._is_model_custom is False
+
+        # Sum of models
+        for model_name in ['spherical+gaussian', 'cubic+matern+stable']:
+            V = Variogram(self.c, self.v, model=model_name)
+            assert V._model_name == model_name
+            assert V._is_model_custom is False
+
+        # Custom model
+        @variogram
+        def custom_model(h, r1, c1, x):
+            return spherical(h, r1, c1) + x
+
+        with self.assertWarns(UserWarning):
+            V = Variogram(self.c, self.v, model=custom_model)
+
+        assert V._model_name == "custom_model"
+        assert V._model == custom_model
+        assert V._is_model_custom is True
+
+
+
     def test_get_bin_count(self):
 
         V = Variogram(self.c, self.v)
@@ -949,6 +983,130 @@ class TestVariogramFittingProcedure(unittest.TestCase):
 
         self.assertTrue(abs(V.parameters[-1] - 2.) < 1e-10)
 
+    def test_fit_custom_model(self):
+
+        # Define a custom variogram and run the fit
+        @variogram
+        def sum_spherical(h, r1, c1, r2, c2, b1, b2):
+            return spherical(h, r1, c1, b1) + spherical(h, r2, c2, b2)
+
+        with self.assertWarns(UserWarning):
+            # Custom models should ignore the nugget setting, so let's try both and check
+            V = Variogram(self.c, self.v, use_nugget=True, model=sum_spherical)
+            V2 = Variogram(self.c, self.v, use_nugget=False, model=sum_spherical)
+
+        # Check that 6 parameters were found
+        assert len(V.cof) == len(V2.cof) == 6
+
+    def test_fit_sum_models_nugget(self):
+
+        # Define a sum of variogram and run the fit
+        V = Variogram(self.c, self.v, model="spherical+spherical", use_nugget=False)
+        V2 = Variogram(self.c, self.v, model="spherical+spherical", use_nugget=True)
+
+        # Check that 4 parameters were found without nugget, 5 otherwise
+        assert len(V.cof) == 4
+        assert len(V2.cof) == 5
+
+    def test_build_sum_models(self):
+
+        # Initiate variogram without fitting
+        V = Variogram(self.c, self.v, fit_method=None)
+
+        # 1/ Check that errors are raised if argument is not good
+        # Should accept upper case and spaces
+        V._build_sum_models("Spherical  +     spherical")
+
+        # Raises an error if model does not exist
+        with self.assertRaises(ValueError) as e:
+            V._build_sum_models("Notamodel  + spherical")
+            self.assertTrue('One of the theoretical models in the list' in str(e.exception))
+
+        # 2/ Build a sum of two spherical models
+        sum_spherical = V._build_sum_models("spherical+spherical")
+
+        # Testing the same way as in test_models/test_sum_spherical
+        # Parameters for the two spherical models
+        params = [1, 0.3, 10, 0.7]
+
+        # Values at which we'll evaluate the function and its expected result
+        vals = [0, 1, 100]
+        res = [0, 0.3 + spherical(1, 10, 0.7), 1]
+
+        for r, c in zip(res, sum_spherical(vals, *params)):
+            self.assertEqual(r, c)
+
+        # 3/ Build a sum of all models
+        sum_spherical = V._build_sum_models("spherical+gaussian+exponential+cubic+stable+matern")
+        min_nb_args = 2 + 2 + 2 + 2 + 3 + 3
+
+        # Check that the function fails for a number of args of 13 (lower than 14)
+        with self.assertRaises(TypeError) as e:
+            sum_spherical(0, *[1,]*(min_nb_args - 1))
+            self.assertTrue('positional arguments' in str(e.exception))
+
+        # Check that it runs for 14 and that the result is correct
+        sum_res = sum(model(0, 1, 1) for model in [spherical, gaussian, exponential, cubic]) + \
+                  sum((model(0, 1, 1, 1)) for model in [stable, matern])
+        assert sum_spherical(0, *[1, ] * min_nb_args) == sum_res
+
+    def test_fit_bounds(self):
+        """
+        Test the fit_bounds kwarg of Variogram and bounds args of fit()
+        """
+        Vnofit =  Variogram(self.c, self.v, fit_method=None)
+
+        bounds=(0, [np.max(Vnofit.bins), np.max(Vnofit.experimental)])
+
+        # Initiate variogram with bounds for fit
+        V = Variogram(self.c, self.v, fit_bounds=bounds)
+        # Same but calling fit
+        V.fit(bounds=bounds)
+
+        # For a sum of models
+        V = Variogram(self.c, self.v, model='spherical+gaussian', fit_bounds=(0, bounds[1]*2))
+        # Same but calling fit
+        V.fit(bounds=(0, bounds[1]*2))
+
+        # Check an error is raised for a custom model
+        @variogram
+        def custom_model(h, r1, c1, x):
+            return spherical(h, r1, c1) + x
+
+        with self.assertWarns(UserWarning):
+            Variogram(self.c, self.v, model=custom_model)
+
+        bounds_custom = [(0, 0, 0), (np.max(self.c), np.max(self.v), np.max(self.v))]
+
+        # And that no error is raised otherwise
+        Variogram(self.c, self.v, model=custom_model, fit_bounds=bounds_custom)
+
+    def test_fit_p0(self):
+        """
+        Test the fit_p0 kwarg of Variogram and bounds args of fit()
+        """
+        Vnofit =  Variogram(self.c, self.v, fit_method=None)
+        p0 = (np.max(Vnofit.bins), np.max(Vnofit.experimental))
+
+        # Initiate variogram with bounds for fit
+        V = Variogram(self.c, self.v, fit_p0=p0)
+        # Same but calling fit
+        V.fit(p0=p0)
+
+        # For a sum of models
+        V = Variogram(self.c, self.v, model='spherical+gaussian', fit_p0=p0 * 2)
+        # Same but calling fit
+        V.fit(p0=p0 * 2)
+
+        # For a custom model
+        @variogram
+        def custom_model(h, r1, c1, x):
+            return spherical(h, r1, c1) + x
+
+        p0_custom = (np.max(Vnofit.bins), np.max(Vnofit.experimental), np.max(Vnofit.experimental))
+        with self.assertWarns(UserWarning) as w:
+            Variogram(self.c, self.v, model=custom_model, fit_p0=p0_custom)
+
 
 class TestVariogramQualityMeasures(unittest.TestCase):
     def setUp(self):
@@ -972,36 +1130,36 @@ class TestVariogramQualityMeasures(unittest.TestCase):
         V = Variogram(self.c, self.v)
 
         for model, rmse in zip(
-                ['spherical', 'gaussian', 'stable'],
-                [3.3705, 3.3707, 3.193]
+                ['spherical', 'gaussian', 'stable', 'spherical+gaussian'],
+                [3.3705, 3.3707, 3.193, 3.0653]
         ):
             V.set_model(model)
-            self.assertAlmostEqual(V.rmse, rmse, places=4)
+            self.assertAlmostEqual(V.rmse, rmse, places=3)
 
     def test_mean_residual(self):
         V = Variogram(self.c, self.v)
 
         for model, mr in zip(
-            ['spherical', 'cubic', 'stable'],
-            [2.6803, 2.6803, 2.6966]
+            ['spherical', 'cubic', 'stable', 'spherical+gaussian'],
+            [2.6803, 2.6803, 2.6966, 2.4723]
         ):
             V.set_model(model)
-            self.assertAlmostEqual(V.mean_residual, mr, places=4)
+            self.assertAlmostEqual(V.mean_residual, mr, places=3)
 
     def test_nrmse(self):
         V = Variogram(self.c, self.v, n_lags=15)
 
         for model, nrmse in zip(
-            ['spherical', 'gaussian', 'stable', 'exponential'],
-            [0.3536, 0.3535, 0.3361, 0.3499]
+            ['spherical', 'gaussian', 'stable', 'exponential', 'spherical+gaussian'],
+            [0.3536, 0.3535, 0.3361, 0.3499, 0.3264]
         ):
             V.set_model(model)
-            self.assertAlmostEqual(V.nrmse, nrmse, places=4)
+            self.assertAlmostEqual(V.nrmse, nrmse, places=3)
 
     def test_nrmse_r(self):
         V = Variogram(self.c, self.v, estimator='cressie')
 
-        self.assertAlmostEqual(V.nrmse_r, 0.63543, places=5)
+        self.assertAlmostEqual(V.nrmse_r, 0.63543, places=3)
 
     def test_r(self):
         V = Variogram(self.c, self.v, n_lags=12, normalize=False)
@@ -1161,6 +1319,75 @@ class TestVariogramMethods(unittest.TestCase):
             decimal=2
         )
 
+    def test_set_model(self):
+        """Test setting model: individual, sum or custom"""
+        V = self.V.clone()
+
+        # Individual model
+        for model_name in ['spherical', 'gaussian', 'exponential', 'cubic', 'matern', 'stable']:
+            V.set_model(model_name)
+            assert V._model_name == model_name
+            assert V._model == globals()[model_name]
+            assert V._is_model_custom is False
+
+        # Sum of models
+        for model_name in ['spherical+gaussian', 'cubic+matern+stable']:
+            V.set_model(model_name)
+            assert V._model_name == model_name
+            assert V._is_model_custom is False
+
+        # Custom model
+        @variogram
+        def custom_model(h, r1, c1, x):
+            return spherical(h, r1, c1) + x
+
+        V.set_model(custom_model)
+        assert V._model_name == "custom_model"
+        assert V._model == custom_model
+        assert V._is_model_custom is True
+
+    def test_describe(self):
+        """Test the describe functions for different models"""
+        V = self.V.clone()
+
+        keys_model = ['normalized_effective_range', 'normalized_sill', 'normalized_nugget',
+                      'effective_range', 'sill', 'nugget']
+
+        # Individual model: normal keys should be in dict
+        for model_name in ['spherical', 'gaussian', 'exponential', 'cubic', 'matern', 'stable']:
+            V.set_model(model_name)
+            dict = V.describe()
+            for key in keys_model:
+                assert key in dict
+            if model_name == 'matern':
+                assert 'smoothness' in dict
+            if model_name == 'stable':
+                assert 'shape' in dict
+
+        # Sum of models: keys with a numbered ID should be in dict
+        V.set_model('spherical+gaussian')
+        dict = V.describe()
+        for key in [k+'_1' for k in keys_model] + [k+'_2' for k in keys_model]:
+            assert key in dict
+
+        V.set_model('cubic+matern+stable')
+        dict = V.describe()
+        for key in [k + '_1' for k in keys_model] + [k + '_2' for k in keys_model] + [k + '_3' for k in keys_model]:
+            assert key in dict
+        assert 'smoothness_2' in dict
+        assert 'shape_3' in dict
+
+        # Custom model
+        @variogram
+        def custom_model(h, r1, c1, x):
+            return spherical(h, r1, c1) + x
+        V.set_model(custom_model)
+        with self.assertWarns(UserWarning) as w:
+            dict = V.describe()
+        custom_keys = ["param1_r1", "param2_c1", "param3_x"]
+        for key in custom_keys:
+            assert key in dict
+
     def test_parameter_property_matern(self):
         V = self.V.clone()
 
@@ -1177,6 +1404,34 @@ class TestVariogramMethods(unittest.TestCase):
         V.set_model('stable')
         assert_array_almost_equal(V.parameters, param, decimal=2)
 
+    def test_parameter_property_sum_models(self):
+        V = self.V.clone()
+
+        # Using models with 3 parameters (range, sill, nugget)
+        param = [3.67, 10.61, 0, 42.3, 5.02, 0]
+        V.set_model('spherical+gaussian')
+        assert_array_almost_equal(V.parameters, param, decimal=2)
+
+        # Using a stable model with 4 parameters
+        param2 = [35.77, 10.97, 0, 0, 42.3, 4.74, 0]
+        V.set_model('stable+cubic')
+        assert_array_almost_equal(V.parameters, param2, decimal=2)
+
+    def test_parameter_property_custom_model(self):
+        V = self.V.clone()
+
+        # Custom model will give parameters back in the order of the custom function
+        @variogram
+        def custom_model(h, r1, c1, x):
+            return spherical(h, r1, c1) + x
+
+        param = [42.3, 5.76, 9.79]
+        V.set_model(custom_model)
+
+        # Provide bounds to avoid a random fit
+        V.fit(bounds=(0, [np.max(V.bins), np.max(V.experimental), np.max(V.experimental)]))
+        assert_array_almost_equal(V.parameters, param, decimal=2)
+
 
 class TestVariogramPlotlyPlots(unittest.TestCase):
     def setUp(self):
@@ -1185,7 +1440,16 @@ class TestVariogramPlotlyPlots(unittest.TestCase):
         self.c = np.random.gamma(10, 4, (150, 2))
         np.random.seed(42)
         self.v = np.random.normal(10, 4, 150)
+        # V = individual model
         self.V = Variogram(self.c, self.v)
+        # V2 = sum of models
+        self.V2 = Variogram(self.c, self.v, model='spherical+gaussian')
+        # V3 = custom model
+        @variogram
+        def custom_model(h, r1, c1, x):
+            return spherical(h, r1, c1) + x
+        self.V3 = Variogram(self.c, self.v, model=custom_model,
+                            fit_bounds=(0, [np.max(self.V.bins), np.max(self.V.experimental), np.max(self.V.experimental)]))
 
     def test_plotly_main_plot(self):
         if PLOTLY_FOUND:
@@ -1196,6 +1460,14 @@ class TestVariogramPlotlyPlots(unittest.TestCase):
                 isinstance(self.V.plot(show=False), go.Figure)
             )
 
+            self.assertTrue(
+                isinstance(self.V2.plot(show=False), go.Figure)
+            )
+
+            self.assertTrue(
+                isinstance(self.V3.plot(show=False), go.Figure)
+            )
+
             plotting.backend('matplotlib')
 
     def test_plotly_scattergram(self):