Repository: fatiando/verde
Pull request #102 (2018-07-28 02:04:33)

Make the jacobian functions into class methods
`Trend`, `Spline`, and `Vector2D` defined the Jacobian matrix
calculations as separate functions. Make them methods of each class
instead.
For `Spline` and `Vector2D`, the `jacobian` method is now responsible
for picking the locations of the forces. Fitting twice will remove the
existing force locations and pick new ones.


Patch:
diff --git a/doc/conf.py b/doc/conf.py
index 320eefb..a67027b 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -33,6 +33,7 @@ intersphinx_mapping = {
     "python": ("https://docs.python.org/3/", None),
     "numpy": ("https://docs.scipy.org/doc/numpy/", None),
     "scipy": ("https://docs.scipy.org/doc/scipy/reference", None),
+    "sklearn": ("http://scikit-learn.org/stable/", None),
     "pandas": ("http://pandas.pydata.org/pandas-docs/stable/", None),
     "xarray": ("http://xarray.pydata.org/en/stable/", None),
     "cartopy": ("https://scitools.org.uk/cartopy/docs/latest/", None),
diff --git a/verde/__init__.py b/verde/__init__.py
index 524cb0c..c6bde10 100644
--- a/verde/__init__.py
+++ b/verde/__init__.py
@@ -16,12 +16,12 @@ from .grid_math import distance_mask
 from .utils import variance_to_weights, maxabs
 from .blockreduce import block_split, BlockReduce, BlockMean
 from .scipygridder import ScipyGridder
-from .trend import Trend, trend_jacobian
+from .trend import Trend
 from .chain import Chain
 from .components import Components
-from .spline import Spline, spline_jacobian
+from .spline import Spline
 from .model_selection import cross_val_score, train_test_split
-from .vector import Vector2D, vector2d_jacobian
+from .vector import Vector2D
 
 
 def test(doctest=True, verbose=True, coverage=False, figures=True):
diff --git a/verde/spline.py b/verde/spline.py
index f2bb70e..0fd89a2 100644
--- a/verde/spline.py
+++ b/verde/spline.py
@@ -83,10 +83,6 @@ class Spline(BaseGridder):
         interpolator. Used as the default region for the
         :meth:`~verde.Spline.grid` and :meth:`~verde.Spline.scatter` methods.
 
-    See also
-    --------
-    verde.spline_jacobian : Make the Jacobian matrix for the biharmonic spline
-
     """
 
     def __init__(self, fudge=1e-5, damping=None, shape=None, spacing=None, region=None):
@@ -124,12 +120,14 @@ class Spline(BaseGridder):
             Returns this estimator instance for chaining operations.
 
         """
+        # Remove pre-existing force coordinates when fitting for a second time.
+        if hasattr(self, "force_coords_"):
+            del self.force_coords_
         coordinates, data, weights = check_fit_input(coordinates, data, weights)
         self._check_weighted_exact_solution(weights)
         # Capture the data region to use as a default when gridding.
         self.region_ = get_region(coordinates[:2])
-        self.force_coords_ = self._get_force_coordinates(coordinates)
-        jacobian = spline_jacobian(coordinates[:2], self.force_coords_, self.fudge)
+        jacobian = self.jacobian(coordinates[:2])
         self.force_ = self._estimate_forces(jacobian, data, weights)
         return self
 
@@ -215,55 +213,46 @@ class Spline(BaseGridder):
 
         """
         check_is_fitted(self, ["force_", "force_coords_"])
-        jac = spline_jacobian(coordinates[:2], self.force_coords_, self.fudge)
+        jac = self.jacobian(coordinates[:2])
         shape = np.broadcast(*coordinates[:2]).shape
         return jac.dot(self.force_).reshape(shape)
 
+    def jacobian(self, coordinates, dtype="float64"):
+        """
+        Make the Jacobian matrix for the 2D biharmonic spline.
 
-def spline_jacobian(coordinates, force_coordinates, fudge=1e-5):
-    """
-    Make the Jacobian matrix for the 2D biharmonic spline.
-
-    Follows [Sandwell1987]_.
-
-    Each column of the Jacobian is the Green's function for a single force
-    evaluated on all observation points.
+        Each column of the Jacobian is the Green's function for a single force evaluated
+        on all observation points [Sandwell1987]_.
 
-    Parameters
-    ----------
-    coordinates : tuple of arrays
-        Arrays with the coordinates of each data point. Should be in the
-        following order: (easting, northing, vertical, ...). Only easting and
-        northing will be used, all subsequent coordinates will be ignored.
-    force_coordinates : tuple of arrays
-        Arrays with the coordinates of each vertical force. Should be in the
-        following order: (easting, northing, vertical, ...). Only easting and
-        northing will be used, all subsequent coordinates will be ignored.
-    fudge : float
-        The positive fudge factor applied to the Green's function to avoid
-        singularities.
-
-    Returns
-    -------
-    jacobian : 2D array
-        The (n_data, n_forces) Jacobian matrix.
+        Parameters
+        ----------
+        coordinates : tuple of arrays
+            Arrays with the coordinates of each data point. Should be in the
+            following order: (easting, northing, vertical, ...). Only easting and
+            northing will be used, all subsequent coordinates will be ignored.
+        dtype : str or numpy dtype
+            The type of the Jacobian array.
 
-    See also
-    --------
-    verde.Spline : Biharmonic spline gridder
+        Returns
+        -------
+        jacobian : 2D array
+            The (n_data, n_forces) Jacobian matrix.
 
-    """
-    force_easting, force_northing = (
-        np.atleast_1d(i).ravel() for i in force_coordinates[:2]
-    )
-    easting, northing = (np.atleast_1d(i).ravel() for i in coordinates[:2])
-    # Reshaping the data to a column vector will automatically build a
-    # distance matrix between each data point and force.
-    distance = np.hypot(
-        easting.reshape((easting.size, 1)) - force_easting.ravel(),
-        northing.reshape((northing.size, 1)) - force_northing.ravel(),
-    )
-    # The fudge factor helps avoid singular matrices when the force and
-    # computation point are too close
-    distance += fudge
-    return (distance ** 2) * (np.log(distance) - 1)
+        """
+        if not hasattr(self, "force_coords_"):
+            self.force_coords_ = self._get_force_coordinates(coordinates)
+        force_easting, force_northing = (
+            np.atleast_1d(i).ravel() for i in self.force_coords_[:2]
+        )
+        easting, northing = (np.atleast_1d(i).ravel() for i in coordinates[:2])
+        # Reshaping the data to a column vector will automatically build a
+        # distance matrix between each data point and force.
+        distance = np.hypot(
+            easting.reshape((easting.size, 1)) - force_easting.ravel(),
+            northing.reshape((northing.size, 1)) - force_northing.ravel(),
+            dtype=dtype,
+        )
+        # The fudge factor helps avoid singular matrices when the force and
+        # computation point are too close
+        distance += self.fudge
+        return (distance ** 2) * (np.log(distance) - 1)
diff --git a/verde/trend.py b/verde/trend.py
index 407432f..4989eae 100644
--- a/verde/trend.py
+++ b/verde/trend.py
@@ -63,11 +63,6 @@ class Trend(BaseGridder):
     >>> print('{:.2f}'.format(trend_out.residual_[2, 2]))
     500.00
 
-    See also
-    --------
-    verde.trend_jacobian : Make the Jacobian matrix for 2D polynomial
-    verde.VectorTrend : Polynomial trend estimator for multi-component data
-
     """
 
     def __init__(self, degree):
@@ -103,7 +98,7 @@ class Trend(BaseGridder):
         """
         coordinates, data, weights = check_fit_input(coordinates, data, weights)
         self.region_ = get_region(coordinates[:2])
-        jac = trend_jacobian(coordinates[:2], degree=self.degree, dtype=data.dtype)
+        jac = self.jacobian(coordinates[:2], dtype=data.dtype)
         scaler = StandardScaler(copy=False, with_mean=False, with_std=True)
         jac = scaler.fit_transform(jac)
         regr = LinearRegression(fit_intercept=False, normalize=False)
@@ -133,10 +128,62 @@ class Trend(BaseGridder):
 
         """
         check_is_fitted(self, ["coef_"])
-        jac = trend_jacobian(coordinates[:2], degree=self.degree)
+        jac = self.jacobian(coordinates[:2])
         shape = np.broadcast(*coordinates[:2]).shape
         return jac.dot(self.coef_).reshape(shape)
 
+    def jacobian(self, coordinates, dtype="float64"):
+        """
+        Make the Jacobian matrix for a 2D polynomial.
+
+        Each column of the Jacobian is ``easting**i * northing**j`` for each (i, j)
+        pair in the polynomial.
+
+        Parameters
+        ----------
+        coordinates : tuple of arrays
+            Arrays with the coordinates of each data point. Should be in the
+            following order: (easting, northing, vertical, ...). Only easting and
+            northing will be used, all subsequent coordinates will be ignored.
+        dtype : str or numpy dtype
+            The type of the output Jacobian numpy array.
+
+        Returns
+        -------
+        jacobian : 2D array
+            The (n_data, n_coefficients) Jacobian matrix.
+
+        Examples
+        --------
+
+        >>> import numpy as np
+        >>> east = np.linspace(0, 4, 5)
+        >>> north = np.linspace(-5, -1, 5)
+        >>> print(Trend(degree=1).jacobian((east, north), dtype=np.int))
+        [[ 1  0 -5]
+         [ 1  1 -4]
+         [ 1  2 -3]
+         [ 1  3 -2]
+         [ 1  4 -1]]
+        >>> print(Trend(degree=2).jacobian((east, north), dtype=np.int))
+        [[ 1  0 -5  0  0 25]
+         [ 1  1 -4  1 -4 16]
+         [ 1  2 -3  4 -6  9]
+         [ 1  3 -2  9 -6  4]
+         [ 1  4 -1 16 -4  1]]
+
+        """
+        easting, northing = coordinates[:2]
+        if easting.shape != northing.shape:
+            raise ValueError("Coordinate arrays must have the same shape.")
+        combinations = polynomial_power_combinations(self.degree)
+        ndata = easting.size
+        nparams = len(combinations)
+        out = np.empty((ndata, nparams), dtype=dtype)
+        for col, (i, j) in enumerate(combinations):
+            out[:, col] = (easting.ravel() ** i) * (northing.ravel() ** j)
+        return out
+
 
 def polynomial_power_combinations(degree):
     """
@@ -170,63 +217,3 @@ def polynomial_power_combinations(degree):
         raise ValueError("Invalid polynomial degree '{}'. Must be >= 1.".format(degree))
     combinations = ((i, j) for j in range(degree + 1) for i in range(degree + 1 - j))
     return tuple(sorted(combinations, key=sum))
-
-
-def trend_jacobian(coordinates, degree, dtype="float64"):
-    """
-    Make the Jacobian for a 2D polynomial of the given degree.
-
-    Each column of the Jacobian is ``easting**i * northing**j`` for each (i, j)
-    pair in the polynomial of the given degree.
-
-    Parameters
-    ----------
-    coordinates : tuple of arrays
-        Arrays with the coordinates of each data point. Should be in the
-        following order: (easting, northing, vertical, ...). Only easting and
-        northing will be used, all subsequent coordinates will be ignored.
-    degree : int
-        The degree of the 2D polynomial. Must be >= 1.
-    dtype : str or numpy dtype
-        The type of the output Jacobian numpy array.
-
-    Returns
-    -------
-    jacobian : 2D array
-        The (n_data, n_coefficients) Jacobian matrix.
-
-    Examples
-    --------
-
-    >>> import numpy as np
-    >>> east = np.linspace(0, 4, 5)
-    >>> north = np.linspace(-5, -1, 5)
-    >>> print(trend_jacobian((east, north), degree=1, dtype=np.int))
-    [[ 1  0 -5]
-     [ 1  1 -4]
-     [ 1  2 -3]
-     [ 1  3 -2]
-     [ 1  4 -1]]
-    >>> print(trend_jacobian((east, north), degree=2, dtype=np.int))
-    [[ 1  0 -5  0  0 25]
-     [ 1  1 -4  1 -4 16]
-     [ 1  2 -3  4 -6  9]
-     [ 1  3 -2  9 -6  4]
-     [ 1  4 -1 16 -4  1]]
-
-    See also
-    --------
-    verde.Trend : Polynomial trend estimator
-    verde.VectorTrend : Polynomial trend estimator for multi-component data
-
-    """
-    easting, northing = coordinates[:2]
-    if easting.shape != northing.shape:
-        raise ValueError("Coordinate arrays must have the same shape.")
-    combinations = polynomial_power_combinations(degree)
-    ndata = easting.size
-    nparams = len(combinations)
-    out = np.empty((ndata, nparams), dtype=dtype)
-    for col, (i, j) in enumerate(combinations):
-        out[:, col] = (easting.ravel() ** i) * (northing.ravel() ** j)
-    return out
diff --git a/verde/vector.py b/verde/vector.py
index 029c948..ccaf7fa 100644
--- a/verde/vector.py
+++ b/verde/vector.py
@@ -81,10 +81,6 @@ class Vector2D(Spline):
         :meth:`~verde.Vector2D.grid` and :meth:`~verde.Vector2D.scatter`
         methods.
 
-    See also
-    --------
-    verde.vector2d_jacobian : Jacobian matrix for the 2D elastic deformation
-
     """
 
     def __init__(
@@ -132,6 +128,9 @@ class Vector2D(Spline):
             Returns this estimator instance for chaining operations.
 
         """
+        # Remove pre-existing force coordinates when fitting for a second time.
+        if hasattr(self, "force_coords_"):
+            del self.force_coords_
         coordinates, data, weights = check_fit_input(
             coordinates, data, weights, unpack=False
         )
@@ -141,16 +140,13 @@ class Vector2D(Spline):
             )
         # Capture the data region to use as a default when gridding.
         self.region_ = get_region(coordinates[:2])
-        self.force_coords_ = self._get_force_coordinates(coordinates)
         if any(w is not None for w in weights):
             weights = np.concatenate([i.ravel() for i in weights])
         else:
             weights = None
         self._check_weighted_exact_solution(weights)
         data = np.concatenate([i.ravel() for i in data])
-        jacobian = vector2d_jacobian(
-            coordinates[:2], self.force_coords_, self.poisson, self.fudge
-        )
+        jacobian = self.jacobian(coordinates[:2])
         self.force_ = self._estimate_forces(jacobian, data, weights)
         return self
 
@@ -176,83 +172,62 @@ class Vector2D(Spline):
 
         """
         check_is_fitted(self, ["force_", "force_coords_"])
-        jac = vector2d_jacobian(
-            coordinates[:2], self.force_coords_, self.poisson, self.fudge
-        )
+        jac = self.jacobian(coordinates[:2])
         cast = np.broadcast(*coordinates[:2])
         npoints = cast.size
         components = jac.dot(self.force_).reshape((2, npoints))
         return tuple(comp.reshape(cast.shape) for comp in components)
 
+    def jacobian(self, coordinates, dtype="float64"):
+        """
+        Make the Jacobian matrix for the 2D coupled elastic deformation.
 
-def vector2d_jacobian(
-    coordinates, force_coordinates, poisson, fudge=1e-5, dtype="float32"
-):
-    """
-    Make the Jacobian matrix for the 2D coupled elastic deformation.
-
-    Follows [SandwellWessel2016]_.
-
-    The Jacobian is segmented into 4 parts, each relating a force component to
-    a data component::
+        The Jacobian is segmented into 4 parts, each relating a force component to a
+        data component [SandwellWessel2016]_::
 
-        | J_ee  J_ne |*|f_e| = |d_e|
-        | J_ne  J_nn | |f_n|   |d_n|
+            | J_ee  J_ne |*|f_e| = |d_e|
+            | J_ne  J_nn | |f_n|   |d_n|
 
-    The forces and data are assumed to be stacked into 1D arrays with the east
-    component on top of the north component.
+        The forces and data are assumed to be stacked into 1D arrays with the east
+        component on top of the north component.
 
-    Parameters
-    ----------
-    coordinates : tuple of arrays
-        Arrays with the coordinates of each data point. Should be in the
-        following order: (easting, northing, vertical, ...). Only easting and
-        northing will be used, all subsequent coordinates will be ignored.
-    force_coordinates : tuple of arrays
-        Arrays with the coordinates of each vertical force. Should be in the
-        following order: (easting, northing, vertical, ...). Only easting and
-        northing will be used, all subsequent coordinates will be ignored.
-    poisson ; float
-        The Poisson's ratio for the elastic deformation Green's functions.
-        A value of -1 will lead to uncoupled interpolation of
-        the east and north data components (the ``J_ne`` component of the
-        Jacobian is null).
-    fudge : float
-        The positive fudge factor applied to the Green's function to avoid
-        singularities.
-    dtype : str or numpy dtype
-        The type of the Jacobian array.
-
-    Returns
-    -------
-    jacobian : 2D array
-        The (n_data*2, n_forces*2) Jacobian matrix.
+        Parameters
+        ----------
+        coordinates : tuple of arrays
+            Arrays with the coordinates of each data point. Should be in the
+            following order: (easting, northing, vertical, ...). Only easting and
+            northing will be used, all subsequent coordinates will be ignored.
+        dtype : str or numpy dtype
+            The type of the Jacobian array.
 
-    See also
-    --------
-    verde.Vector2D : Coupled gridder for 2-component vector data
+        Returns
+        -------
+        jacobian : 2D array
+            The (n_data*2, n_forces*2) Jacobian matrix.
 
-    """
-    force_coordinates = [np.atleast_1d(i).ravel() for i in force_coordinates[:2]]
-    coordinates = [np.atleast_1d(i).ravel() for i in coordinates[:2]]
-    npoints = coordinates[0].size
-    nforces = force_coordinates[0].size
-    # Reshaping the data coordinates to a column vector will automatically
-    # build a distance matrix between each data point and force.
-    east, north = (
-        datac.reshape((npoints, 1)) - forcec
-        for datac, forcec in zip(coordinates, force_coordinates)
-    )
-    distance = np.hypot(east, north, dtype=dtype)
-    # The fudge factor helps avoid singular matrices when the force and
-    # computation point are too close
-    distance += fudge
-    # Pre-compute common terms for the Green's functions of each component
-    ln_r = (3 - poisson) * np.log(distance)
-    over_r2 = (1 + poisson) / distance ** 2
-    jac = np.empty((npoints * 2, nforces * 2), dtype=dtype)
-    jac[:npoints, :nforces] = ln_r + over_r2 * north ** 2  # J_ee
-    jac[npoints:, nforces:] = ln_r + over_r2 * east ** 2  # J_nn
-    jac[:npoints, nforces:] = -over_r2 * east * north  # J_ne
-    jac[npoints:, :nforces] = jac[:npoints, nforces:]  # J is symmetric
-    return jac
+        """
+        if not hasattr(self, "force_coords_"):
+            self.force_coords_ = self._get_force_coordinates(coordinates)
+        force_coordinates = [np.atleast_1d(i).ravel() for i in self.force_coords_[:2]]
+        coordinates = [np.atleast_1d(i).ravel() for i in coordinates[:2]]
+        npoints = coordinates[0].size
+        nforces = force_coordinates[0].size
+        # Reshaping the data coordinates to a column vector will automatically
+        # build a distance matrix between each data point and force.
+        east, north = (
+            datac.reshape((npoints, 1)) - forcec
+            for datac, forcec in zip(coordinates, force_coordinates)
+        )
+        distance = np.hypot(east, north, dtype=dtype)
+        # The fudge factor helps avoid singular matrices when the force and
+        # computation point are too close
+        distance += self.fudge
+        # Pre-compute common terms for the Green's functions of each component
+        ln_r = (3 - self.poisson) * np.log(distance)
+        over_r2 = (1 + self.poisson) / distance ** 2
+        jac = np.empty((npoints * 2, nforces * 2), dtype=dtype)
+        jac[:npoints, :nforces] = ln_r + over_r2 * north ** 2  # J_ee
+        jac[npoints:, nforces:] = ln_r + over_r2 * east ** 2  # J_nn
+        jac[:npoints, nforces:] = -over_r2 * east * north  # J_ne
+        jac[npoints:, :nforces] = jac[:npoints, nforces:]  # J is symmetric
+        return jac


Test patch:
diff --git a/verde/tests/test_spline.py b/verde/tests/test_spline.py
index b0885ee..8e97e73 100644
--- a/verde/tests/test_spline.py
+++ b/verde/tests/test_spline.py
@@ -34,6 +34,20 @@ def test_spline():
     )
 
 
+def test_spline_twice():
+    "Check that the forces are updated when fitting twice"
+    grid = CheckerBoard(region=(1000, 5000, -8000, -6000)).grid(shape=(10, 10))
+    coords = np.meshgrid(grid.easting, grid.northing)
+    spline = Spline()
+    spline.fit(coords, grid.scalars.values)
+    npt.assert_allclose(spline.force_coords_, coords)
+    grid2 = CheckerBoard(region=(-15, -5, 10, 20)).grid(shape=(10, 10))
+    coords2 = np.meshgrid(grid2.easting, grid2.northing)
+    spline.fit(coords2, grid2.scalars.values)
+    npt.assert_allclose(spline.force_coords_, coords2)
+    assert not np.allclose(spline.force_coords_, coords)
+
+
 def test_spline_region():
     "See if the region is gotten from the data is correct."
     region = (1000, 5000, -8000, -6000)
diff --git a/verde/tests/test_trend.py b/verde/tests/test_trend.py
index 6f53d54..3a21a3b 100644
--- a/verde/tests/test_trend.py
+++ b/verde/tests/test_trend.py
@@ -6,7 +6,7 @@ import numpy as np
 import numpy.testing as npt
 import pytest
 
-from ..trend import Trend, polynomial_power_combinations, trend_jacobian
+from ..trend import Trend, polynomial_power_combinations
 from ..coordinates import grid_coordinates
 
 
@@ -52,5 +52,6 @@ def test_polynomial_combinations_fails():
 def test_trend_jacobian_fails():
     "Test failing conditions for the trend jacobian builder"
     east, north = np.arange(50), np.arange(30)
+    trend = Trend(degree=1)
     with pytest.raises(ValueError):
-        trend_jacobian((east, north), degree=1)
+        trend.jacobian((east, north))
diff --git a/verde/tests/test_vector.py b/verde/tests/test_vector.py
index f34f1b1..01a9097 100644
--- a/verde/tests/test_vector.py
+++ b/verde/tests/test_vector.py
@@ -32,6 +32,17 @@ def test_vector2d(data2d):
     npt.assert_allclose(spline.force_coords_, coords)
 
 
+def test_vector2d_twice(data2d):
+    "Check that the force coordinates are updated if fitting a second time"
+    coords, data = data2d
+    spline = Vector2D().fit(coords, data)
+    npt.assert_allclose(spline.force_coords_, coords)
+    coords2 = tuple(i * 1000 for i in coords)
+    spline.fit(coords2, data)
+    npt.assert_allclose(spline.force_coords_, coords2)
+    assert not np.allclose(spline.force_coords_, coords)
+
+
 def test_vector2d_weights(data2d):
     "Use unit weights and a regular grid solution"
     coords, data = data2d