Add support for sum of models and custom models

skgstat/DirectionalVariogram.py

@@ -294,6 +294,7 @@ def __init__(self,
 
         # set the directional model
         self._directional_model = None
+        self._is_model_custom = False
         self.set_directional_model(model_name=directional_model)
 
         # the binning settings

skgstat/Variogram.py

@@ -2,6 +2,7 @@
 Variogram class
 """
 import copy
+import inspect
 import warnings
 from typing import Iterable, Callable, Union, Tuple
 
@@ -335,6 +336,7 @@ def __init__(self,
 
         # model can be a function or a string
         self._model = 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
@@ -973,6 +975,7 @@ def set_model(self, model_name):
                     ) % model_name
                 )
         else:  # pragma: no cover
+            self._is_model_custom = True
             self._model = model_name
 
     def _build_harmonized_model(self):
@@ -1454,7 +1457,7 @@ def preprocessing(self, force=False):
         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
@@ -1563,18 +1566,34 @@ def fit(self, force=False, method=None, sigma=None, **kwargs):
             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)
+            nb_args = len(argspec) - 1
+            if bounds is None:
+                bounds = (0, [np.maximum(np.nanmax(x), np.nanmax(y))] * nb_args)
+            if p0 is None:
+                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':

skgstat/tests/test_models.py

@@ -186,6 +186,21 @@ def adder(l, a):
         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()

skgstat/tests/test_variogram.py

@@ -19,6 +19,7 @@
 from skgstat import OrdinaryKriging
 from skgstat import estimators
 from skgstat import plotting
+from skgstat.models import variogram, spherical, matern
 
 
 class TestSpatiallyCorrelatedData(unittest.TestCase):
@@ -61,7 +62,7 @@ def test_sparse_maxlag_30(self):
             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)
@@ -949,6 +950,18 @@ def test_implicit_nugget(self):
 
         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)
+
+        V = Variogram(self.c, self.v, use_nugget=True, model=sum_spherical)
+
+        # Check that 6 parameters were found
+        assert len(V.cof) == 6
+
 
 class TestVariogramQualityMeasures(unittest.TestCase):
     def setUp(self):