Repository: fatiando/verde (python)

Issue:
Add a Components class for multicomponent data
If we want to grid multicomponent data, we can use something like `Vector2D` to grid them jointly or we can use a `Spline` for each component separately. In this case, there is no way to put them in a `Chain` with `BlockReduce` or `Trend`. Also, `VectorTrend` is doing this for trend estimation. 

A better approach would be to have a meta-gridder `Components` that takes multiple estimators and fits each one to a given data component. It would be a generalization of `VectorTrend` and would replace it. 

This is how the usage would look:

```python
grd = Chain([
    ('mean', BlockReduce(np.mean, spacing=1),
    ('trend', Components([Trend(degree=1), Trend(degree=2)])),
    ('spline', Components([Spline(), Spline()]))
])
grd.fit(coordinates, (data_east, data_north))
# Access each gridder separately
trend_east = grd.named_steps['trend'].components[0].grid()
``` 

Pull request description:
Add a Components class to fit vector data
The `Components` class is a meta-estimator, like `Chain` that 
fits a given list of estimators to the matching components of
vector data. This replaces the `VectorTrend` class.

Fixes #74


Patch:
diff --git a/doc/api/index.rst b/doc/api/index.rst
index acc0a9c..c6e89cb 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -26,7 +26,6 @@ Data Processing
     BlockReduce
     BlockMean
     Trend
-    VectorTrend
 
 Composite Estimators
 --------------------
@@ -35,6 +34,7 @@ Composite Estimators
    :toctree: generated/
 
     Chain
+    Components
 
 Model Selection
 ---------------
diff --git a/examples/vector2d.py b/examples/vector2d.py
index 045bbd5..c9680f9 100644
--- a/examples/vector2d.py
+++ b/examples/vector2d.py
@@ -1,12 +1,13 @@
 """
-Gridding 2-component vector data
-================================
+Coupled gridding of 2-component vector data
+===========================================
 
 One way of gridding vector data would be grid each component separately using
-:class:`verde.Spline`. Alternatively, :class:`verde.Vector2D` can grid two components
-simultaneously in a way that couples them through elastic deformation theory. This is
-particularly suited, though not exclusive, to data that represent elastic/semi-elastic
-deformation, like horizontal GPS velocities.
+:class:`verde.Spline` and :class:`verde.Components`. Alternatively,
+:class:`verde.Vector2D` can grid two components simultaneously in a way that couples
+them through elastic deformation theory. This is particularly suited, though not
+exclusive, to data that represent elastic/semi-elastic deformation, like horizontal GPS
+velocities.
 """
 import matplotlib.pyplot as plt
 import cartopy.crs as ccrs
@@ -38,7 +39,7 @@ spacing = 20 / 60
 chain = vd.Chain(
     [
         ("mean", vd.BlockReduce(np.mean, spacing * 111e3)),
-        ("trend", vd.VectorTrend(degree=5)),
+        ("trend", vd.Components([vd.Trend(degree=5) for i in range(2)])),
         ("spline", vd.Vector2D(poisson=0.5)),
     ]
 )
diff --git a/examples/vector_trend.py b/examples/vector_trend.py
index f100037..d3192cd 100644
--- a/examples/vector_trend.py
+++ b/examples/vector_trend.py
@@ -2,11 +2,12 @@
 Trends in vector data
 =====================
 
-Verde provides the :class:`verde.VectorTrend` class to estimate a polynomial trend on
-each component of vector data, like GPS velocities. You can access each trend as a
-separate :class:`verde.Trend` or operate on all vector components directly using using
-:meth:`verde.VectorTrend.predict`, :meth:`verde.VectorTrend.grid`, etc, or chaining it
-with a vector interpolator using :class:`verde.Chain`.
+Verde provides the :class:`verde.Trend` class to estimate a polynomial trend and the
+:class:`verde.Components` class to apply any combination of estimators to each component
+of vector data, like GPS velocities. You can access each component as a separate
+(fitted) :class:`verde.Trend` instance or operate on all vector components directly
+using using :meth:`verde.Components.predict`, :meth:`verde.Components.grid`, etc, or
+chaining it with a vector interpolator using :class:`verde.Chain`.
 """
 import matplotlib.pyplot as plt
 import cartopy.crs as ccrs
@@ -18,9 +19,11 @@ import verde as vd
 # North-West because of the relative movement along the San Andreas Fault System.
 data = vd.datasets.fetch_california_gps()
 
-# We'll fit a degree 4 trend on both the East and North components and weight the data
+# We'll fit a degree 2 trend on both the East and North components and weight the data
 # using the inverse of the variance of each component.
-trend = vd.VectorTrend(degree=4)
+# Note: Never use [Trend(...)]*2 as an argument to Components. This creates references
+# to the same Trend instance and will mess up the fitting.
+trend = vd.Components([vd.Trend(degree=2) for i in range(2)])
 weights = vd.variance_to_weights((data.std_east ** 2, data.std_north ** 2))
 trend.fit(
     coordinates=(data.longitude, data.latitude),
@@ -29,12 +32,12 @@ trend.fit(
 )
 print("Vector trend estimator:", trend)
 
-# The separate Trend objects for each component can be accessed through the 'component_'
+# The separate Trend objects for each component can be accessed through the 'components'
 # attribute. You could grid them individually if you wanted.
-print("East component trend:", trend.component_[0])
-print("East trend coefficients:", trend.component_[0].coef_)
-print("North component trend:", trend.component_[1])
-print("North trend coefficients:", trend.component_[1].coef_)
+print("East component trend:", trend.components[0])
+print("East trend coefficients:", trend.components[0].coef_)
+print("North component trend:", trend.components[1])
+print("North trend coefficients:", trend.components[1].coef_)
 
 # We can make a grid with both trend components as data variables
 grid = trend.grid(spacing=0.1, dims=["latitude", "longitude"])
@@ -73,19 +76,9 @@ for ax, component, title in zip(axes, components, titles):
         data.velocity_north.values,
         scale=0.2,
         transform=crs,
-        color="y",
-        label="Original data",
-    )
-    # and the residuals
-    ax.quiver(
-        data.longitude.values,
-        data.latitude.values,
-        trend.residual_[0].values,
-        trend.residual_[1].values,
-        scale=0.2,
-        transform=crs,
         color="k",
-        label="Residuals",
+        width=0.001,
+        label="Original data",
     )
     # Setup the map ticks
     vd.datasets.setup_california_gps_map(ax, land=None, ocean=None)
diff --git a/verde/__init__.py b/verde/__init__.py
index ffd199f..524cb0c 100644
--- a/verde/__init__.py
+++ b/verde/__init__.py
@@ -16,8 +16,9 @@ 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, VectorTrend
+from .trend import Trend, trend_jacobian
 from .chain import Chain
+from .components import Components
 from .spline import Spline, spline_jacobian
 from .model_selection import cross_val_score, train_test_split
 from .vector import Vector2D, vector2d_jacobian
diff --git a/verde/chain.py b/verde/chain.py
index e8b22bf..10e4f50 100644
--- a/verde/chain.py
+++ b/verde/chain.py
@@ -39,6 +39,10 @@ class Chain(BaseGridder):
         A dictionary version of *steps* where the ``'name'``  strings are keys
         and the estimator/gridder/processor objects are the values.
 
+    See also
+    --------
+    verde.Components : Fit an estimator to each component of multi-component vector data
+
     """
 
     def __init__(self, steps):
diff --git a/verde/components.py b/verde/components.py
new file mode 100644
index 0000000..9c3f880
--- /dev/null
+++ b/verde/components.py
@@ -0,0 +1,122 @@
+"""
+Components meta-estimator to use different estimators on multi-component data.
+"""
+from sklearn.utils.validation import check_is_fitted
+
+from .base import BaseGridder, check_fit_input
+from .coordinates import get_region
+
+
+class Components(BaseGridder):
+    """
+    Fit an estimator to each component of multi-component vector data.
+
+    Provides a convenient way of fitting and gridding vector data using scalar gridders
+    and estimators.
+
+    Each data component provided to :meth:`~verde.Components.fit` is fitted to a
+    separated estimator. Methods like :meth:`~verde.Components.grid` and
+    :meth:`~verde.Components.predict` will operate on the multiple components
+    simultaneously.
+
+    .. warning::
+
+        Never pass code like this as input to this class: ``[vd.Trend(1)]*3``. This
+        creates 3 references to the **same instance** of ``Trend``, which means that
+        they will all get the same coefficients after fitting. Use a list comprehension
+        instead: ``[vd.Trend(1) for i in range(3)]``.
+
+    Parameters
+    ----------
+    components : tuple or list
+        A tuple or list of the estimator/gridder instances used for each component. The
+        estimators will be applied for each data component in the same order that they
+        are given here.
+
+    Attributes
+    ----------
+    components : tuple
+        Tuple of the fitted estimators on each component of the data.
+    region_ : tuple
+        The boundaries (``[W, E, S, N]``) of the data used to fit the interpolator. Used
+        as the default region for the :meth:`~verde.Components.grid` and
+        :meth:`~verde.Components.scatter` methods.
+
+    See also
+    --------
+    verde.Chain : Chain filtering operations to fit on each subsequent output.
+
+    """
+
+    def __init__(self, components):
+        self.components = components
+
+    def fit(self, coordinates, data, weights=None):
+        """
+        Fit the estimators to the given multi-component data.
+
+        The data region is captured and used as default for the
+        :meth:`~verde.Components.grid` and :meth:`~verde.Components.scatter` methods.
+
+        All input arrays must have the same shape. If weights are given, there
+        must be a separate array for each component of the data.
+
+        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.
+        data : tuple of array
+            The data values of each component at each data point. Must be a
+            tuple.
+        weights : None or tuple of array
+            If not None, then the weights assigned to each data point of each
+            data component. Typically, this should be 1 over the data
+            uncertainty squared.
+
+        Returns
+        -------
+        self
+            Returns this estimator instance for chaining operations.
+
+        """
+        if not isinstance(data, tuple):
+            raise ValueError(
+                "Data must be a tuple of arrays. {} given.".format(type(data))
+            )
+        if weights is not None and not isinstance(weights, tuple):
+            raise ValueError(
+                "Weights must be a tuple of arrays. {} given.".format(type(weights))
+            )
+        coordinates, data, weights = check_fit_input(coordinates, data, weights)
+        self.region_ = get_region(coordinates[:2])
+        for estimator, data_comp, weight_comp in zip(self.components, data, weights):
+            estimator.fit(coordinates, data_comp, weight_comp)
+        return self
+
+    def predict(self, coordinates):
+        """
+        Evaluate each data component on a set of points.
+
+        Requires a fitted estimator (see :meth:`~verde.Components.fit`).
+
+        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.
+
+        Returns
+        -------
+        data : tuple of array
+            The values for each vector component evaluated on the given points. The
+            order of components will be the same as was provided to
+            :meth:`~verde.Components.fit`.
+
+        """
+        check_is_fitted(self, ["region_"])
+        return tuple(comp.predict(coordinates) for comp in self.components)
diff --git a/verde/trend.py b/verde/trend.py
index 3916f03..407432f 100644
--- a/verde/trend.py
+++ b/verde/trend.py
@@ -230,116 +230,3 @@ def trend_jacobian(coordinates, degree, dtype="float64"):
     for col, (i, j) in enumerate(combinations):
         out[:, col] = (easting.ravel() ** i) * (northing.ravel() ** j)
     return out
-
-
-class VectorTrend(BaseGridder):
-    """
-    Fit 2D polynomial trends to multi-component vector data.
-
-    Each data component is fitted to a separated :class:`verde.Trend`. The
-    trends can be evaluated separately through the ``VectorTrend.component_``
-    attribute or jointly using the methods of this class.
-
-    Parameters
-    ----------
-    degree : int
-        The degree of the polynomials. Must be >= 1.
-
-    Attributes
-    ----------
-    component_ : list
-        List of the fitted :class:`verde.Trend` on each component of the data.
-    region_ : tuple
-        The boundaries (``[W, E, S, N]``) of the data used to fit the
-        interpolator. Used as the default region for the
-        :meth:`~verde.Trend.grid` and :meth:`~verde.Trend.scatter` methods.
-    residual_ : list of array
-        The difference between each component of the input data and the
-        predicted polynomial trend.
-
-    See also
-    --------
-    verde.Trend : Polynomial trend estimator
-    verde.trend_jacobian : Make the Jacobian matrix for 2D polynomial
-
-    """
-
-    def __init__(self, degree):
-        self.degree = degree
-
-    def fit(self, coordinates, data, weights=None):
-        """
-        Fit the trend to the given multi-component data.
-
-        The data region is captured and used as default for the
-        :meth:`~verde.VectorTrend.grid` and :meth:`~verde.VectorTrend.scatter`
-        methods.
-
-        All input arrays must have the same shape. If weights are given, there
-        must be a separate array for each component of the data.
-
-        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.
-        data : tuple of array
-            The data values of each component at each data point. Must be a
-            tuple.
-        weights : None or tuple of array
-            If not None, then the weights assigned to each data point of each
-            data component. Typically, this should be 1 over the data
-            uncertainty squared.
-
-        Returns
-        -------
-        self
-            Returns this estimator instance for chaining operations.
-
-        """
-        if not isinstance(data, tuple):
-            raise ValueError(
-                "Data must be a tuple of arrays. {} given.".format(type(data))
-            )
-        if weights is not None and not isinstance(weights, tuple):
-            raise ValueError(
-                "Weights must be a tuple of arrays. {} given.".format(type(weights))
-            )
-        coordinates, data, weights = check_fit_input(coordinates, data, weights)
-        self.region_ = get_region(coordinates[:2])
-        self.component_ = [
-            Trend(degree=self.degree).fit(coordinates, data_comp, weight_comp)
-            for data_comp, weight_comp in zip(data, weights)
-        ]
-        self.residual_ = tuple(
-            data_comp - comp.predict(coordinates).reshape(data_comp.shape)
-            for data_comp, comp in zip(data, self.component_)
-        )
-        return self
-
-    def predict(self, coordinates):
-        """
-        Evaluate the polynomial trend of each component on a of points.
-
-        Requires a fitted estimator (see :meth:`~verde.VectorTrend.fit`).
-
-        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.
-
-        Returns
-        -------
-        data : tuple of array
-            The trend values for each vector component evaluated on the given
-            points. The order of components will be the same as was provided to
-            :meth:`~verde.VectorTrend.fit`.
-
-        """
-        check_is_fitted(self, ["component_"])
-        return tuple(comp.predict(coordinates) for comp in self.component_)


Test patch:
diff --git a/verde/tests/test_chain.py b/verde/tests/test_chain.py
index ea9871c..b8957ef 100644
--- a/verde/tests/test_chain.py
+++ b/verde/tests/test_chain.py
@@ -9,8 +9,9 @@ from ..blockreduce import BlockReduce
 from ..chain import Chain
 from ..scipygridder import ScipyGridder
 from ..spline import Spline
-from ..trend import Trend, VectorTrend
+from ..trend import Trend
 from ..coordinates import grid_coordinates
+from ..components import Components
 
 
 def test_chain():
@@ -62,10 +63,13 @@ def test_chain_vector():
     trend = coefs[0] + coefs[1] * coords[0] + coefs[2] * coords[1]
     data = (trend, trend.copy())
     chain = Chain(
-        [("mean", BlockReduce(np.mean, spacing=0.5)), ("trend", VectorTrend(degree=1))]
+        [
+            ("mean", BlockReduce(np.mean, spacing=0.5)),
+            ("trend", Components([Trend(degree=1), Trend(degree=1)])),
+        ]
     )
     chain.fit(coords, data)
     npt.assert_allclose(chain.score(coords, data), 1)
     npt.assert_allclose(chain.predict(coords), data)
-    npt.assert_allclose(chain.named_steps["trend"].component_[0].coef_, coefs)
-    npt.assert_allclose(chain.named_steps["trend"].component_[1].coef_, coefs)
+    npt.assert_allclose(chain.named_steps["trend"].components[0].coef_, coefs)
+    npt.assert_allclose(chain.named_steps["trend"].components[1].coef_, coefs)
diff --git a/verde/tests/test_components.py b/verde/tests/test_components.py
new file mode 100644
index 0000000..0753c92
--- /dev/null
+++ b/verde/tests/test_components.py
@@ -0,0 +1,77 @@
+# pylint: disable=redefined-outer-name
+"""
+Test the Components meta-gridder.
+"""
+import numpy as np
+import numpy.testing as npt
+import pytest
+
+from ..trend import Trend
+from ..coordinates import grid_coordinates
+from ..components import Components
+
+
+@pytest.fixture()
+def simple_2d_model():
+    "A single degree=1 2-component model"
+    east, north = grid_coordinates((1000, 5000, -5000, -1000), shape=(50, 50))
+    coefs = ([-20, 0.5, 3], [10, 2, -0.4])
+    data = tuple(c[0] + c[1] * east + c[2] * north for c in coefs)
+    return (east, north), coefs, data
+
+
+@pytest.fixture()
+def simple_3d_model():
+    "A single degree=1 3-component model"
+    east, north = grid_coordinates((1000, 5000, -5000, -1000), shape=(50, 50))
+    coefs = ([-20, 0.5, 3], [10, 2, -0.4], [30, -10, -1.3])
+    data = tuple(c[0] + c[1] * east + c[2] * north for c in coefs)
+    return (east, north), coefs, data
+
+
+def test_components_vector_trend(simple_2d_model):
+    "Test the vector trend estimation on a simple problem"
+    coords, coefs, data = simple_2d_model
+    trend = Components([Trend(degree=1), Trend(degree=1)]).fit(coords, data)
+    for i, coef in enumerate(coefs):
+        npt.assert_allclose(trend.components[i].coef_, coef)
+        npt.assert_allclose(trend.predict(coords)[i], data[i])
+        npt.assert_allclose(trend.score(coords, data), 1)
+
+
+def test_components_vector_trend_3d(simple_3d_model):
+    "Test the vector trend estimation on a simple problem"
+    coords, coefs, data = simple_3d_model
+    trend = Components([Trend(degree=1), Trend(degree=1), Trend(degree=1)])
+    trend.fit(coords, data)
+    for i, coef in enumerate(coefs):
+        npt.assert_allclose(trend.components[i].coef_, coef)
+        npt.assert_allclose(trend.predict(coords)[i], data[i])
+        npt.assert_allclose(trend.score(coords, data), 1)
+
+
+def test_components_vector_trend_fails(simple_2d_model):
+    "Test the vector trend estimation on a simple problem"
+    coords, _, data = simple_2d_model
+    trend = Components([Trend(degree=1), Trend(degree=1)])
+    with pytest.raises(ValueError):
+        trend.fit(coords, list(data))
+    with pytest.raises(ValueError):
+        trend.fit(coords, data, weights=[np.ones_like(data)] * 2)
+
+
+def test_components_vector_trend_weights(simple_2d_model):
+    "Use weights to account for outliers"
+    coords, coefs, data = simple_2d_model
+    outlier = np.abs(data[0]).max() * 3
+    data_out = tuple(i.copy() for i in data)
+    weights = tuple(np.ones_like(i) for i in data)
+    for i, coef in enumerate(coefs):
+        data_out[i][20, 20] += outlier
+        weights[i][20, 20] = 1e-10
+    trend = Components([Trend(degree=1), Trend(degree=1)])
+    trend.fit(coords, data_out, weights)
+    for i, coef in enumerate(coefs):
+        npt.assert_allclose(trend.components[i].coef_, coef)
+        npt.assert_allclose((data_out[i] - trend.predict(coords)[i])[20, 20], outlier)
+        npt.assert_allclose(trend.predict(coords)[i], data[i])
diff --git a/verde/tests/test_trend.py b/verde/tests/test_trend.py
index b6bab07..6f53d54 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, VectorTrend
+from ..trend import Trend, polynomial_power_combinations, trend_jacobian
 from ..coordinates import grid_coordinates
 
 
@@ -19,24 +19,6 @@ def simple_model():
     return (east, north), coefs, data
 
 
-@pytest.fixture()
-def simple_2d_model():
-    "A single degree=1 2-component model"
-    east, north = grid_coordinates((1000, 5000, -5000, -1000), shape=(50, 50))
-    coefs = ([-20, 0.5, 3], [10, 2, -0.4])
-    data = tuple(c[0] + c[1] * east + c[2] * north for c in coefs)
-    return (east, north), coefs, data
-
-
-@pytest.fixture()
-def simple_3d_model():
-    "A single degree=1 3-component model"
-    east, north = grid_coordinates((1000, 5000, -5000, -1000), shape=(50, 50))
-    coefs = ([-20, 0.5, 3], [10, 2, -0.4], [30, -10, -1.3])
-    data = tuple(c[0] + c[1] * east + c[2] * north for c in coefs)
-    return (east, north), coefs, data
-
-
 def test_trend(simple_model):
     "Test the trend estimation on a simple problem"
     coords, coefs, data = simple_model
@@ -72,49 +54,3 @@ def test_trend_jacobian_fails():
     east, north = np.arange(50), np.arange(30)
     with pytest.raises(ValueError):
         trend_jacobian((east, north), degree=1)
-
-
-def test_vector_trend(simple_2d_model):
-    "Test the vector trend estimation on a simple problem"
-    coords, coefs, data = simple_2d_model
-    trend = VectorTrend(degree=1).fit(coords, data)
-    for i, coef in enumerate(coefs):
-        npt.assert_allclose(trend.component_[i].coef_, coef)
-        npt.assert_allclose(trend.predict(coords)[i], data[i])
-        npt.assert_allclose(trend.score(coords, data), 1)
-
-
-def test_vector_trend_3d(simple_3d_model):
-    "Test the vector trend estimation on a simple problem"
-    coords, coefs, data = simple_3d_model
-    trend = VectorTrend(degree=1).fit(coords, data)
-    for i, coef in enumerate(coefs):
-        npt.assert_allclose(trend.component_[i].coef_, coef)
-        npt.assert_allclose(trend.predict(coords)[i], data[i])
-        npt.assert_allclose(trend.score(coords, data), 1)
-
-
-def test_vector_trend_fails(simple_2d_model):
-    "Test the vector trend estimation on a simple problem"
-    coords, _, data = simple_2d_model
-    trend = VectorTrend(degree=1)
-    with pytest.raises(ValueError):
-        trend.fit(coords, list(data))
-    with pytest.raises(ValueError):
-        trend.fit(coords, data, weights=[np.ones_like(data)] * 2)
-
-
-def test_vector_trend_weights(simple_2d_model):
-    "Use weights to account for outliers"
-    coords, coefs, data = simple_2d_model
-    outlier = np.abs(data[0]).max() * 3
-    data_out = tuple(i.copy() for i in data)
-    weights = tuple(np.ones_like(i) for i in data)
-    for i, coef in enumerate(coefs):
-        data_out[i][20, 20] += outlier
-        weights[i][20, 20] = 1e-10
-    trend = VectorTrend(degree=1).fit(coords, data_out, weights)
-    for i, coef in enumerate(coefs):
-        npt.assert_allclose(trend.component_[i].coef_, coef)
-        npt.assert_allclose((data_out[i] - trend.predict(coords)[i])[20, 20], outlier)
-        npt.assert_allclose(trend.predict(coords)[i], data[i])