Repository: fatiando/verde
Pull request #185 (2019-05-14 04:54:27)

Add SplineCV: a cross-validated Spline
The `SplineCV` class perfoms grid search cross-validation to auto-tune a `Spline`.
It uses `cross_val_score` to cross-validate and can be configured to run in parallel
with a `dask.distributed.Client` and take different scikit-learn splitters besides
`KFold`. Tested on the checkerboard function and add a gallery example using the
Texas temperature data.

<!-- Please describe changes proposed and WHY you made them. If unsure, open an issue first to discuss the idea. -->

<!-- If fixing an issue, put the issue number after the # below (no spaces). Github will automatically close it when this gets merged. -->

**Reminders**

- [ ] Run `make format` and `make check` to make sure the code follows the style guide.
- [ ] Add tests for new features or tests that would have caught the bug that you're fixing.
- [ ] Add new public functions/methods/classes to `doc/api/index.rst` and `verde/__init__.py`.
- [ ] Write detailed docstrings for all functions/classes/methods.
- [ ] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.

leouieda: @santisoler finally had some time to finish this PR. Would you mind taking a look at it? 

I didn't edit the tutorial because I'm thinking that we should split the "Model Selection" tutorial into "Model Performance" and "Model Tuning". The first only focuses on running `cross_val_score` and friends. The second does the auto-tuning by hand and introduces `SplineCV`. What do you think?

Patch:
diff --git a/doc/api/index.rst b/doc/api/index.rst
index 7ef2a85..b4c07aa 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -14,6 +14,7 @@ Interpolators
    :toctree: generated/
 
     Spline
+    SplineCV
     VectorSpline2D
     ScipyGridder
 
diff --git a/examples/spline_cv.py b/examples/spline_cv.py
new file mode 100644
index 0000000..a5bea4b
--- /dev/null
+++ b/examples/spline_cv.py
@@ -0,0 +1,80 @@
+"""
+Gridding with splines (cross-validated)
+=======================================
+
+The :class:`verde.Spline` has two main parameters that need to be configured:
+
+1. ``mindist``: the minimum distance between forces and data points
+2. ``damping``: the regularization parameter controlling smoothness
+
+These parameters can be determined through cross-validation (see :ref:`model_selection`)
+automatically using :class:`verde.SplineCV`. It is very similar to :class:`verde.Spline`
+but takes a set of parameter values instead of only one value. When calling
+:meth:`verde.SplineCV.fit`, the class will:
+
+1. Create a spline for each combination of the input parameter sets
+2. Calculate the cross-validation score for each spline using
+   :func:`verde.cross_val_score`
+3. Pick the spline with the highest score
+
+"""
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import pyproj
+import numpy as np
+import verde as vd
+
+# We'll test this on the air temperature data from Texas
+data = vd.datasets.fetch_texas_wind()
+coordinates = (data.longitude.values, data.latitude.values)
+region = vd.get_region(coordinates)
+
+# Use a Mercator projection for our Cartesian gridder
+projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
+
+# The output grid spacing will 15 arc-minutes
+spacing = 15 / 60
+
+# This spline will automatically perform cross-validation and search for the optimal
+# parameter configuration.
+spline = vd.SplineCV(dampings=(1e-5, 1e-3, 1e-1), mindists=(10e3, 50e3, 100e3))
+
+# Fit the model on the data. Under the hood, the class will perform K-fold
+# cross-validation for each the 3*3=9 parameter combinations and pick the one with the
+# highest R² score.
+spline.fit(projection(*coordinates), data.air_temperature_c)
+
+# We can show the best R² score obtained in the cross-validation
+print("\nScore: {:.3f}".format(spline.scores_.max()))
+
+# And then the best spline parameters that produced this high score.
+print("\nBest spline configuration:")
+print("  mindist:", spline.mindist_)
+print("  damping:", spline.damping_)
+
+# Now we can create a geographic grid of air temperature by providing a projection
+# function to the grid method and mask points that are too far from the observations
+grid_full = spline.grid(
+    region=region,
+    spacing=spacing,
+    projection=projection,
+    dims=["latitude", "longitude"],
+    data_names=["temperature"],
+)
+grid = vd.distance_mask(
+    coordinates, maxdist=3 * spacing * 111e3, grid=grid_full, projection=projection
+)
+
+# Plot the grid and the original data points
+plt.figure(figsize=(8, 6))
+ax = plt.axes(projection=ccrs.Mercator())
+ax.set_title("Air temperature gridded with biharmonic spline")
+ax.plot(*coordinates, ".k", markersize=1, transform=ccrs.PlateCarree())
+tmp = grid.temperature.plot.pcolormesh(
+    ax=ax, cmap="plasma", transform=ccrs.PlateCarree(), add_colorbar=False
+)
+plt.colorbar(tmp).set_label("Air temperature (C)")
+# Use an utility function to add tick labels and land and ocean features to the map.
+vd.datasets.setup_texas_wind_map(ax, region=region)
+plt.tight_layout()
+plt.show()
diff --git a/tutorials/overview.py b/tutorials/overview.py
index 7a545e0..9ff5f69 100644
--- a/tutorials/overview.py
+++ b/tutorials/overview.py
@@ -25,17 +25,15 @@ Before we get started, here are a few of the conventions we use across Verde:
 * The term "region" means the bounding box of the data. It is ordered west, east, south,
   north.
 
-"""
+The library
+-----------
 
-########################################################################################
-# The library
-# -----------
-#
-# Most classes and functions are available through the :mod:`verde` top level package.
-# The only exceptions are the functions related to loading sample data, which are in
-# :mod:`verde.datasets`. Throughout the documentation we'll use ``vd`` as the alias for
-# :mod:`verde`.
+Most classes and functions are available through the :mod:`verde` top level package.
+The only exceptions are the functions related to loading sample data, which are in
+:mod:`verde.datasets`. Throughout the documentation we'll use ``vd`` as the alias for
+:mod:`verde`.
 
+"""
 import verde as vd
 
 ########################################################################################
diff --git a/verde/__init__.py b/verde/__init__.py
index 6092c57..fc96f69 100644
--- a/verde/__init__.py
+++ b/verde/__init__.py
@@ -20,7 +20,7 @@ from .blockreduce import BlockReduce, BlockMean
 from .scipygridder import ScipyGridder
 from .trend import Trend
 from .chain import Chain
-from .spline import Spline
+from .spline import Spline, SplineCV
 from .model_selection import cross_val_score, train_test_split
 from .vector import Vector, VectorSpline2D
 
diff --git a/verde/model_selection.py b/verde/model_selection.py
index d19f32c..842bb78 100644
--- a/verde/model_selection.py
+++ b/verde/model_selection.py
@@ -18,6 +18,8 @@ class DummyClient:  # pylint: disable=no-self-use,too-few-public-methods
     >>> client = DummyClient()
     >>> client.submit(sum, (1, 2, 3))
     6
+    >>> print(client.scatter("bla"))
+    bla
 
     """
 
@@ -25,6 +27,10 @@ class DummyClient:  # pylint: disable=no-self-use,too-few-public-methods
         "Execute function with the given arguments and return its output"
         return function(*args, **kwargs)
 
+    def scatter(self, value):
+        "Does nothing but return the input"
+        return value
+
 
 def train_test_split(coordinates, data, weights=None, **kwargs):
     r"""
@@ -161,10 +167,11 @@ def cross_val_score(estimator, coordinates, data, weights=None, cv=None, client=
     >>> print(', '.join(['{:.2f}'.format(score) for score in scores]))
     1.00, 1.00, 1.00, 1.00, 1.00
 
-    >>> # To run parallel, we need to create a dask.distributed Client. It will
-    >>> # create a local cluster if no arguments are given so we can run the
-    >>> # scoring on a single machine. We'll use threads instead of processes for this
-    >>> # example but in most cases you'll want processes.
+    To run parallel, we need to create a :class:`dask.distributed.Client`. It will
+    create a local cluster if no arguments are given so we can run the scoring on a
+    single machine. We'll use threads instead of processes for this example but in most
+    cases you'll want processes.
+
     >>> from dask.distributed import Client
     >>> client = Client(processes=False)
     >>> # The scoring will now only submit tasks to our local cluster
@@ -175,6 +182,8 @@ def cross_val_score(estimator, coordinates, data, weights=None, cv=None, client=
     >>> # We need to call .result() to get back the actual value
     >>> print('{:.2f}'.format(scores[0].result()))
     1.00
+    >>> # Close the client and shutdown the local cluster
+    >>> client.close()
 
     """
     coordinates, data, weights = check_fit_input(
diff --git a/verde/spline.py b/verde/spline.py
index 6743e5d..fb25c8e 100644
--- a/verde/spline.py
+++ b/verde/spline.py
@@ -1,6 +1,7 @@
 """
 Biharmonic splines in 2D.
 """
+import itertools
 from warnings import warn
 
 import numpy as np
@@ -9,6 +10,7 @@ from sklearn.utils.validation import check_is_fitted
 from .base import n_1d_arrays, BaseGridder, check_fit_input, least_squares
 from .coordinates import get_region
 from .utils import parse_engine
+from .model_selection import DummyClient, cross_val_score
 
 try:
     import numba
@@ -22,6 +24,205 @@ except ImportError:
 JIT_ARGS = dict(nopython=True, target="cpu", fastmath=True, parallel=True)
 
 
+class SplineCV(BaseGridder):
+    r"""
+    Cross-validated biharmonic spline interpolation.
+
+    Similar to :class:`verde.Spline` but automatically chooses the best *damping* and
+    *mindist* parameters using cross-validation. Tests all combinations of the given
+    *dampings* and *mindists* and selects the maximum (or minimum) mean cross-validation
+    score (i.e., a grid search).
+
+    The grid search can optionally run in parallel using :mod:`dask`. To do this, pass
+    in a :class:`dask.distributed.Client` as the *client* argument. The client can
+    manage a process or thread pool in your local computer or a remote cluster.
+
+    Other cross-validation generators from :mod:`sklearn.model_selection` can be used by
+    passing them through the *cv* argument.
+
+    Parameters
+    ----------
+    mindists : iterable or 1d array
+        List (or other iterable) of *mindist* parameter values to try. Can be considered
+        a minimum distance between the point forces and data points. Needed because the
+        Green's functions are singular when forces and data points coincide. Acts as a
+        fudge factor.
+    dampings : iterable or 1d array
+        List (or other iterable) of *damping* parameter values to try. Is the positive
+        damping regularization parameter. Controls how much smoothness is 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.SplineCV.fit` is called.
+    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
+        numpy otherwise. The numba version is multi-threaded and usually faster, which
+        makes fitting and predicting faster.
+    cv : None or cross-validation generator
+        Any scikit-learn cross-validation generator. If not given, will use the default
+        set by :func:`verde.cross_val_score`.
+    client : None or dask.distributed.Client
+        If None, then computations are run serially. Otherwise, should be a
+        dask ``Client`` object. It will be used to dispatch computations to the
+        dask cluster.
+
+
+    Attributes
+    ----------
+    force_ : array
+        The estimated forces that fit the observed 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.SplineCV.grid` and :meth:`~verde.SplineCV.scatter` methods.
+    scores_ : array
+        The mean cross-validation score for each parameter combination.
+    mindist_ : float
+        The optimal value for the *mindist* parameter.
+    damping_ : float
+        The optimal value for the *damping* parameter.
+
+    See also
+    --------
+    Spline : The bi-harmonic spline
+    cross_val_score : Score an estimator/gridder using cross-validation
+
+    """
+
+    def __init__(
+        self,
+        mindists=(1e3, 10e3, 100e3),
+        dampings=(1e-10, 1e-5, 1e-1),
+        force_coords=None,
+        engine="auto",
+        cv=None,
+        client=None,
+    ):
+        super().__init__()
+        self.dampings = dampings
+        self.mindists = mindists
+        self.force_coords = force_coords
+        self.engine = engine
+        self.cv = cv
+        self.client = client
+
+    def fit(self, coordinates, data, weights=None):
+        """
+        Fit the biharmonic spline to the given data and automatically tune parameters.
+
+        For each combination of the parameters given, computes the mean cross validation
+        score using :func:`verde.cross_val_score` and the given CV splitting class (the
+        *cv* parameter of this class). The configuration with the best score is then
+        chosen and used to fit the entire dataset.
+
+        The data region is captured and used as default for the
+        :meth:`~verde.SplineCV.grid` and :meth:`~verde.SplineCV.scatter` methods.
+
+        All input arrays must have the same shape.
+
+        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 : array
+            The data values of each data point.
+        weights : None or array
+            If not None, then the weights assigned to each data point.
+            Typically, this should be 1 over the data uncertainty squared.
+
+        Returns
+        -------
+        self
+            Returns this estimator instance for chaining operations.
+
+        """
+        if self.client is None:
+            client = DummyClient()
+        else:
+            client = self.client
+        # If this is a parallel client, scattering the data makes sure each worker has a
+        # copy of the data so we minimize transfer.
+        futures = tuple(client.scatter(i) for i in (coordinates, data, weights))
+        parameter_sets = [
+            dict(
+                mindist=combo[0],
+                damping=combo[1],
+                engine=self.engine,
+                force_coords=self.force_coords,
+            )
+            for combo in itertools.product(self.mindists, self.dampings)
+        ]
+        scores = []
+        for params in parameter_sets:
+            spline = Spline(**params)
+            scores.append(
+                client.submit(
+                    cross_val_score,
+                    spline,
+                    coordinates=futures[0],
+                    data=futures[1],
+                    weights=futures[2],
+                    cv=self.cv,
+                )
+            )
+        if self.client is not None:
+            scores = [score.result() for score in scores]
+        self.scores_ = np.mean(scores, axis=1)
+        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
+    def force_(self):
+        "The estimated forces that fit the data."
+        return self._best.force_
+
+    @property
+    def region_(self):
+        "The bounding region of the data used to fit the spline"
+        return self._best.region_
+
+    @property
+    def damping_(self):
+        "The optimal damping parameter"
+        return self._best.damping
+
+    @property
+    def mindist_(self):
+        "The optimal mindist parameter"
+        return self._best.mindist
+
+    def predict(self, coordinates):
+        """
+        Evaluate the best estimated spline on the given set of points.
+
+        Requires a fitted estimator (see :meth:`~verde.SplineCV.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 : array
+            The data values evaluated on the given points.
+
+        """
+        check_is_fitted(self, ["_best"])
+        return self._best.predict(coordinates)
+
+
 class Spline(BaseGridder):
     r"""
     Biharmonic spline interpolation using Green's functions.
@@ -89,6 +290,10 @@ class Spline(BaseGridder):
         interpolator. Used as the default region for the
         :meth:`~verde.Spline.grid` and :meth:`~verde.Spline.scatter` methods.
 
+    See also
+    --------
+    SplineCV : Cross-validated version of the bi-harmonic spline
+
     """
 
     def __init__(self, mindist=1e-5, damping=None, force_coords=None, engine="auto"):


Test patch:
diff --git a/verde/tests/test_spline.py b/verde/tests/test_spline.py
index bf3ac81..dbc1db6 100644
--- a/verde/tests/test_spline.py
+++ b/verde/tests/test_spline.py
@@ -5,14 +5,91 @@ import warnings
 
 import numpy as np
 import numpy.testing as npt
+from sklearn.model_selection import ShuffleSplit
 
-from ..spline import Spline
+
+from ..spline import Spline, SplineCV
 from ..datasets.synthetic import CheckerBoard
-from .utils import requires_numba
+from .utils import requires_numba, requires_dask
+
+
+def test_spline_cv():
+    "See if the cross-validated spline solution matches the synthetic data"
+    region = (100, 500, -800, -700)
+    synth = CheckerBoard(region=region)
+    data = synth.scatter(size=1500, random_state=1)
+    coords = (data.easting, data.northing)
+    # Can't test on many configurations because it takes too long for regular testing
+    spline = SplineCV(
+        dampings=[None],
+        mindists=[1e-5, 1e-3],
+        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
+    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(
+        spline.region_,
+        (
+            data.easting.min(),
+            data.easting.max(),
+            data.northing.min(),
+            data.northing.max(),
+        ),
+    )
+    shape = (5, 5)
+    region = (270, 320, -770, -720)
+    # Tolerance needs to be kind of high to allow for error due to small
+    # dataset
+    npt.assert_allclose(
+        spline.grid(region, shape=shape).scalars,
+        synth.grid(region, shape=shape).scalars,
+        rtol=5e-2,
+    )
+
+
+@requires_dask
+def test_spline_cv_parallel():
+    "See if the parallel version of SplineCV works"
+    from dask.distributed import Client
+
+    region = (100, 500, -800, -700)
+    synth = CheckerBoard(region=region)
+    data = synth.scatter(size=1500, random_state=1)
+    coords = (data.easting, data.northing)
+    client = Client(processes=False)
+    # 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],
+        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
+    assert spline.damping_ is None
+    shape = (5, 5)
+    region = (270, 320, -770, -720)
+    # Tolerance needs to be kind of high to allow for error due to small
+    # dataset
+    npt.assert_allclose(
+        spline.grid(region, shape=shape).scalars,
+        synth.grid(region, shape=shape).scalars,
+        rtol=5e-2,
+    )
 
 
 def test_spline():
-    "See if the exact spline solution."
+    "See if the exact spline solution matches the synthetic data"
     region = (100, 500, -800, -700)
     synth = CheckerBoard(region=region)
     data = synth.scatter(size=1500, random_state=1)
diff --git a/verde/tests/utils.py b/verde/tests/utils.py
index 4a1618a..5f1bf31 100644
--- a/verde/tests/utils.py
+++ b/verde/tests/utils.py
@@ -14,3 +14,15 @@ def requires_numba(function):
         numba = None
     mark = pytest.mark.skipif(numba is None, reason="requires numba")
     return mark(function)
+
+
+def requires_dask(function):
+    """
+    Skip the decorated test if dask.distributed is not installed.
+    """
+    try:
+        import dask.distributed
+    except ImportError:
+        dask = None
+    mark = pytest.mark.skipif(dask is None, reason="requires dask")
+    return mark(function)