Repository: mmaelicke/scikit-gstat

Issue:
Yet another Python loop
There is yet another Python loop in the code:

https://github.com/mmaelicke/scikit-gstat/blob/a43a7f30013061d28531504bfe97b6c10e3c28a5/skgstat/Variogram.py#L1127

and within the same function:

https://github.com/mmaelicke/scikit-gstat/blob/a43a7f30013061d28531504bfe97b6c10e3c28a5/skgstat/Variogram.py#L1122

This should be changed to a vectorized numpy function or at least a `map` at some point in time...

Discussion:

Patch:
diff --git a/docs/changelog.rst b/docs/changelog.rst
index 09c6dbf..fad69bf 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -2,6 +2,21 @@
 Changelog
 =========
 
+Version 0.3.6
+=============
+.. warning:: 
+  There is some potential breaking behaviour
+
+- [Variogram] some internal code cleanup. Removed some unnecessary loops
+- [Variogram] setting the `n_lags <skgstat.Variogram.n_lags>` property now correctly forces
+  a recalculation of the lag groupings. So far they were kept untouches, which might result
+  in old experimental variogram values for the changed instance.
+  **This is a potential breaking change**.
+- [Variogram] The `lag_classes <skgstat.Variogram.lag_classes>` generator now yields empty 
+  arrays for unoccupied lag classes. This will result in `NaN <numpy.NaN>` values for the 
+  semi-variance. This is actually a bug-fix.
+  **This is a potential breaking change**
+
 Version 0.3.5
 =============
 - [plotting] The :func:`location_trend <skgstat.Variogram.location_trend>` can now add 
diff --git a/skgstat/Variogram.py b/skgstat/Variogram.py
index c410c3d..75080bb 100644
--- a/skgstat/Variogram.py
+++ b/skgstat/Variogram.py
@@ -460,6 +460,9 @@ class Variogram(object):
         else:
             raise ValueError('n_lags has to be a positive integer')
 
+        # reset the groups
+        self._groups = None
+
         # reset the fitting
         self.cof = None
         self.cov = None
@@ -801,17 +804,17 @@ class Variogram(object):
         Generates an iterator over all lag classes. Can be zipped with
         Variogram.bins to identify the lag.
 
+        .. versionchanged:: 0.3.6
+            yields an empty array for empty lag groups now
+
         Returns
         -------
         iterable
 
         """
         # yield all groups
-        for i in np.unique(self.lag_groups()):
-            if i < 0:
-                continue
-            else:
-                yield self._diff[np.where(self.lag_groups() == i)]
+        for i in range(len(self.bins)):
+            yield self._diff[np.where(self.lag_groups() == i)]
 
     def preprocessing(self, force=False):
         """Preprocessing function
@@ -1103,33 +1106,37 @@ class Variogram(object):
     @property
     def _experimental(self):
         """
+        Calculates the experimental variogram from the current lag classes.
+        It handles the special case of the `'entropy'` and `'percentile'`
+        estimators, which take an additional argument.
+
+        .. versionchanged:: 0.3.6
+            replaced the for-loops with :func:`fromiter <numpy.fromiter>`
 
         Returns
         -------
+        experimental : np.ndarray
+            1D array of the experimental variogram values. Has same length
+            as :func:`bins <skgstat.Variogram.bins>`
 
         """
-        # prepare the result array
-        y = np.zeros(len(self.bins), dtype=np.float64)
-
-        # args, can set the bins for entropy
-        # and should set p of percentile, not properly implemented
         if self._estimator.__name__ == 'entropy':
             bins = np.linspace(
                 np.min(self.distance),
                 np.max(self.distance),
-                50
+                50      # TODO: this should be set by kwargs
             )
-            # apply
-            for i, lag_values in enumerate(self.lag_classes()):
-                y[i] = self._estimator(lag_values, bins=bins)
 
-        # default
+            def mapper(lag_values):
+                return self._estimator(lag_values, bins=bins)
+        elif self._estimator.__name__ == 'percentile':
+            # TODO after the kwargs are accepted, the p value can be set here
+            mapper = self._estimator
         else:
-            for i, lag_values in enumerate(self.lag_classes()):
-                y[i] = self._estimator(lag_values)
+            mapper = self._estimator
 
-        # apply
-        return y.copy()
+        # return the mapped result
+        return np.fromiter(map(mapper, self.lag_classes()), dtype=float)
 
     def __get_fit_bounds(self, x, y):
         """


Test patch:
diff --git a/skgstat/tests/test_directionalvariogram.py b/skgstat/tests/test_directionalvariogram.py
index a306bfa..8ed2dbf 100644
--- a/skgstat/tests/test_directionalvariogram.py
+++ b/skgstat/tests/test_directionalvariogram.py
@@ -17,14 +17,17 @@ class TestDirectionalVariogramInstantiation(unittest.TestCase):
     def test_standard_settings(self):
         DV = DirectionalVariogram(self.c, self.v, normalize=True)
 
-        for x, y in zip(DV.parameters, [407.467, 2138.098, 0]):
-            self.assertAlmostEqual(x, y, places=0)
+        assert_array_almost_equal(
+            DV.parameters, [436., 2706., 0],
+            decimal=0
+        )
 
     def test_azimuth(self):
         DV = DirectionalVariogram(self.c, self.v, azimuth=-45, normalize=True)
 
-        for x, y in zip(DV.parameters, [27.288, 131.644, 0]):
-            self.assertAlmostEqual(x, y, places=3)
+        assert_array_almost_equal(
+            DV.parameters, [23.438, 219.406, 0], decimal=3
+        )
 
     def test_invalid_azimuth(self):
         with self.assertRaises(ValueError) as e:
@@ -38,8 +41,10 @@ class TestDirectionalVariogramInstantiation(unittest.TestCase):
     def test_tolerance(self):
         DV = DirectionalVariogram(self.c, self.v, tolerance=15, normalize=True)
 
-        for x, y in zip(DV.parameters, [26.342, 1880.015, 0]):
-            self.assertAlmostEqual(x, y, places=3)
+        assert_array_almost_equal(
+            DV.parameters, [435.7, 2722.1, 0],
+            decimal=1
+        )
 
     def test_invalid_tolerance(self):
         with self.assertRaises(ValueError) as e:
diff --git a/skgstat/tests/test_interfaces.py b/skgstat/tests/test_interfaces.py
index 1f26477..e0626bf 100644
--- a/skgstat/tests/test_interfaces.py
+++ b/skgstat/tests/test_interfaces.py
@@ -84,8 +84,9 @@ class TestVariogramEstimator(unittest.TestCase):
         )
 
         gs = gs.fit(self.c, self.v)
-
-        self.assertEqual(gs.best_params_['model'], 'spherical')
+        
+        # Python 3.6 yields 'exponential', while 3.7, 3.8 yield 'gaussian'
+        self.assertTrue(gs.best_params_['model'] in ['gaussian', 'exponential'])
 
     def test_find_best_model_future_cv(self):
         """
@@ -103,7 +104,7 @@ class TestVariogramEstimator(unittest.TestCase):
 
         gs = gs.fit(self.c, self.v)
 
-        self.assertEqual(gs.best_params_['model'], 'matern')
+        self.assertEqual(gs.best_params_['model'], 'exponential')
 
 
 class TestPyKrigeInterface(unittest.TestCase):