Bug fix: Reset force coords when refitting splines

verde/spline.py

@@ -73,6 +73,9 @@ class SplineCV(BaseGridder):
     ----------
     force_ : array
         The estimated forces that fit the observed data.
+    force_coords_ : tuple of arrays
+        The easting and northing coordinates of the point forces. Same as *force_coords*
+        if it is not None. Otherwise, same as the data locations used to fit the spline.
     region_ : tuple
         The boundaries (``[W, E, S, N]``) of the data used to fit the
         interpolator. Used as the default region for the
@@ -176,7 +179,6 @@ def fit(self, coordinates, data, weights=None):
         best = np.argmax(self.scores_)
         self._best = Spline(**parameter_sets[best])
         self._best.fit(coordinates, data, weights=weights)
-        self.force_coords = self._best.force_coords
         return self
 
     @property
@@ -199,6 +201,11 @@ def mindist_(self):
         "The optimal mindist parameter"
         return self._best.mindist
 
+    @property
+    def force_coords_(self):
+        "The optimal force locations"
+        return self._best.force_coords_
+
     def predict(self, coordinates):
         """
         Evaluate the best estimated spline on the given set of points.
@@ -273,8 +280,7 @@ class Spline(BaseGridder):
         imposed on the estimated forces. If None, no regularization is used.
     force_coords : None or tuple of arrays
         The easting and northing coordinates of the point forces. If None (default),
-        then will be set to the data coordinates the first time
-        :meth:`~verde.Spline.fit` is called.
+        then will be set to the data coordinates used to fit the spline.
     engine : str
         Computation engine for the Jacobian matrix and prediction. Can be ``'auto'``,
         ``'numba'``, or ``'numpy'``. If ``'auto'``, will use numba if it is installed or
@@ -285,6 +291,9 @@ class Spline(BaseGridder):
     ----------
     force_ : array
         The estimated forces that fit the observed data.
+    force_coords_ : tuple of arrays
+        The easting and northing coordinates of the point forces. Same as *force_coords*
+        if it is not None. Otherwise, same as the data locations used to fit the spline.
     region_ : tuple
         The boundaries (``[W, E, S, N]``) of the data used to fit the
         interpolator. Used as the default region for the
@@ -336,8 +345,10 @@ def fit(self, coordinates, data, weights=None):
         # Capture the data region to use as a default when gridding.
         self.region_ = get_region(coordinates[:2])
         if self.force_coords is None:
-            self.force_coords = tuple(i.copy() for i in n_1d_arrays(coordinates, n=2))
-        jacobian = self.jacobian(coordinates[:2], self.force_coords)
+            self.force_coords_ = tuple(i.copy() for i in n_1d_arrays(coordinates, n=2))
+        else:
+            self.force_coords_ = self.force_coords
+        jacobian = self.jacobian(coordinates[:2], self.force_coords_)
         self.force_ = least_squares(jacobian, data, weights, self.damping)
         return self
 
@@ -363,7 +374,7 @@ def predict(self, coordinates):
         """
         check_is_fitted(self, ["force_"])
         shape = np.broadcast(*coordinates[:2]).shape
-        force_east, force_north = n_1d_arrays(self.force_coords, n=2)
+        force_east, force_north = n_1d_arrays(self.force_coords_, n=2)
         east, north = n_1d_arrays(coordinates, n=2)
         data = np.empty(east.size, dtype=east.dtype)
         if parse_engine(self.engine) == "numba":

verde/tests/test_spline.py

@@ -22,15 +22,15 @@ def test_spline_cv():
     # Can't test on many configurations because it takes too long for regular testing
     spline = SplineCV(
         dampings=[None],
-        mindists=[1e-5, 1e-3],
+        mindists=[1e-7, 1e-5],
         cv=ShuffleSplit(n_splits=2, random_state=0),
     ).fit(coords, data.scalars)
     # The interpolation should be perfect on top of the data points
     npt.assert_allclose(spline.predict(coords), data.scalars, rtol=1e-5)
     npt.assert_allclose(spline.score(coords, data.scalars), 1)
-    assert spline.mindist_ == 1e-3
+    npt.assert_allclose(spline.force_coords_, coords)
+    assert spline.mindist_ == 1e-7
     assert spline.damping_ is None
-    npt.assert_allclose(spline.force_coords, coords)
     # There should be 1 force per data point
     assert data.scalars.size == spline.force_.size
     npt.assert_allclose(
@@ -66,16 +66,17 @@ def test_spline_cv_parallel():
     # Can't test on many configurations because it takes too long for regular testing
     # Use ShuffleSplit instead of KFold to test it out and make this run faster
     spline = SplineCV(
-        dampings=[None],
-        mindists=[1e-5, 1e-3],
+        dampings=[None, 1e-8],
+        mindists=[1e-7, 1e-5],
         client=client,
         cv=ShuffleSplit(n_splits=1, random_state=0),
     ).fit(coords, data.scalars)
     client.close()
     # The interpolation should be perfect on top of the data points
     npt.assert_allclose(spline.predict(coords), data.scalars, rtol=1e-5)
     npt.assert_allclose(spline.score(coords, data.scalars), 1)
-    assert spline.mindist_ == 1e-5
+    npt.assert_allclose(spline.force_coords_, coords)
+    assert spline.mindist_ == 1e-7
     assert spline.damping_ is None
     shape = (5, 5)
     region = (270, 320, -770, -720)
@@ -100,7 +101,7 @@ def test_spline():
     npt.assert_allclose(spline.score(coords, data.scalars), 1)
     # There should be 1 force per data point
     assert data.scalars.size == spline.force_.size
-    npt.assert_allclose(spline.force_coords, coords)
+    npt.assert_allclose(spline.force_coords_, coords)
     shape = (5, 5)
     region = (270, 320, -770, -720)
     # Tolerance needs to be kind of high to allow for error due to small