Repository: scipy/scipy
Pull request: DOC: improve docstrings of interpolate.NearestNDInterpolator.__call__ & LinearNDInterpolator.__call__
#### Reference issue
fix #12378 

#### What does this implement/fix?
I improved docstrings `interpolate.NearestNDInterpolator.__call__` & `LinearNDInterpolator.__call__`.

Because, these functions also accepts ndim broadcastible arrays as an argument.


Problem statement:
interpolate.NearestNDInterpolator.__call__ & LinearNDInterpolator.__call__ documentation
Documentation states that __call__ method takes a single argument xi as an array of points with shape=(..., ndim), but the __call__ method also accepts ndim broadcastible arrays as an argument. This capability stems from code line:

`xi = _ndim_coords_from_arrays(args, ndim=self.points.shape[1])`

This option for evaluation should be included in the documentation.

I suggest the documentation be changed to something like:

> Parameters:
> -------------
>     x1, x2, ... xndim: array_like of float
>          Points where to interpolate data at.
> 
>          x1, x2, ... xndim can be array_like of float with broadcastible shape.
> 
>          or x1 can be array_like of float with shape (..., ndim)

The follwing code demonstrates the two modes of calling working as expected.

```
import numpy as np
import scipy
from scipy.interpolate import NearestNDInterpolator, LinearNDInterpolator

print('numpy version = ', np.version.full_version)
print('scipy version = ', scipy.version.full_version)


for itype in (NearestNDInterpolator, LinearNDInterpolator):
    shape = (4, 5, 6, 7)
    ndims = len(shape)
    x_0 = [np.arange(0, n, dtype=float) for n in shape]
    pts = np.moveaxis(
        np.meshgrid(*x_0, indexing='ij'), 0, -1
        ).reshape(-1, ndims)
    values = np.sum(pts, axis=1)
    interp = itype(pts, values)
    print(40*'-')
    print(itype.__name__)
    print(40*'-')
    print('interp((3, 3, 3, 3)) = ', interp((3, 3, 3, 3)), '\n')
    print('interp(3, 3, 3, 3) = ', interp(3, 3, 3, 3), '\n')
    print(
        'interp([[2], [3]], [[0, 1, 2]], 3, 3) = ',
        interp([[2], [3]], [[0, 1, 2]], 3, 3), '\n'
        )
```
```

numpy version =  1.18.1
scipy version =  1.4.1
----------------------------------------
NearestNDInterpolator
----------------------------------------
interp((3, 3, 3, 3)) =  12.0 

interp(3, 3, 3, 3) =  12.0 

interp([[2], [3]], [[0, 1, 2]], 3, 3) =  [[ 8.  9. 10.]
 [ 9. 10. 11.]] 

----------------------------------------
LinearNDInterpolator
----------------------------------------
interp((3, 3, 3, 3)) =  12.0 

interp(3, 3, 3, 3) =  12.0 

interp([[2], [3]], [[0, 1, 2]], 3, 3) =  [[ 8.  9. 10.]
 [ 9. 10. 11.]] 

```

Hints from issue discussion:
That behavior is implemented, but it's untested as far as I can tell. Using scalars as input (as in the given example, `interp(3, 3, 3, 3)`) doesn't seem like a good use case, because that adds very little. I can see how it would help if one has a set of 1-D arrays (say `x`, `y`, `z`) for a n-D set of points.

So I think it would be good to update the docs, and regression tests for both interpolators, and add two examples, one for each type of usage.

Commit message: add test_broadcastable_input and fix typos.

Patch to review:
diff --git a/scipy/interpolate/interpnd.pyx b/scipy/interpolate/interpnd.pyx
index 4941d308daf5..e94ba5830f8d 100644
--- a/scipy/interpolate/interpnd.pyx
+++ b/scipy/interpolate/interpnd.pyx
@@ -133,8 +133,10 @@ class NDInterpolatorBase(object):
 
         Parameters
         ----------
-        xi : ndarray of float, shape (..., ndim)
+        x1, x2, ... xn: array_like of float
             Points where to interpolate data at.
+            x1, x2, ... xn can be array_like of float with broadcastable shape.
+            or x1 can be array_like of float with shape (..., ndim)
 
         """
         xi = _ndim_coords_from_arrays(args, ndim=self.points.shape[1])
diff --git a/scipy/interpolate/ndgriddata.py b/scipy/interpolate/ndgriddata.py
index ed8948a904a3..766478780681 100644
--- a/scipy/interpolate/ndgriddata.py
+++ b/scipy/interpolate/ndgriddata.py
@@ -68,8 +68,10 @@ def __call__(self, *args):
 
         Parameters
         ----------
-        xi : ndarray of float, shape (..., ndim)
+        x1, x2, ... xn: array_like of float
             Points where to interpolate data at.
+            x1, x2, ... xn can be array_like of float with broadcastable shape.
+            or x1 can be array_like of float with shape (..., ndim)
 
         """
         xi = _ndim_coords_from_arrays(args, ndim=self.points.shape[1])
diff --git a/scipy/interpolate/tests/test_interpolate.py b/scipy/interpolate/tests/test_interpolate.py
index 634e529cfb5c..7d380ee3d48e 100644
--- a/scipy/interpolate/tests/test_interpolate.py
+++ b/scipy/interpolate/tests/test_interpolate.py
@@ -12,7 +12,8 @@
 from scipy.interpolate import (interp1d, interp2d, lagrange, PPoly, BPoly,
          splrep, splev, splantider, splint, sproot, Akima1DInterpolator,
          RegularGridInterpolator, LinearNDInterpolator, NearestNDInterpolator,
-         RectBivariateSpline, interpn, NdPPoly, BSpline)
+         RectBivariateSpline, interpn, NdPPoly, BSpline,
+         CloughTocher2DInterpolator)
 
 from scipy.special import poch, gamma
 
@@ -2581,6 +2582,32 @@ def test_fillvalue_type(self):
         RegularGridInterpolator(points, values)
         RegularGridInterpolator(points, values, fill_value=0.)
 
+    def test_broadcastable_input(self):
+        # input data
+        np.random.seed(0)
+        x = np.random.random(10)
+        y = np.random.random(10)
+        z = np.hypot(x, y)
+
+        # x-y grid for interpolation
+        X = np.linspace(min(x), max(x))
+        Y = np.linspace(min(y), max(y))
+        X, Y = np.meshgrid(X, Y)
+        XY = np.vstack((X.ravel(), Y.ravel())).T
+
+        for itype in (NearestNDInterpolator, LinearNDInterpolator,
+                      CloughTocher2DInterpolator):
+            interp = itype(list(zip(x, y)), z)
+            # single array input
+            interp_points0 = interp(XY)
+            # broadcastable input
+            interp_points1 = interp((X, Y))
+            interp_points2 = interp((X, 0.0))
+
+            assert_equal(interp_points0.size ==
+                         interp_points1.size ==
+                         interp_points2.size, True)
+
 
 class MyValue(object):
     """


Code review comments:
diff_hunk: @@ -2581,6 +2582,32 @@ def test_fillvalue_type(self):
         RegularGridInterpolator(points, values)
         RegularGridInterpolator(points, values, fill_value=0.)
 
+    def test_broadcastable_input(self):
+        # input data
+        np.random.seed(0)
+        x = np.random.random(10)
+        y = np.random.random(10)
+        z = np.hypot(x, y)
+
+        # x-y grid for interpolation
+        X = np.linspace(min(x), max(x))
+        Y = np.linspace(min(y), max(y))
+        X, Y = np.meshgrid(X, Y)
+        XY = np.vstack((X.ravel(), Y.ravel())).T
+
+        for itype in (NearestNDInterpolator, LinearNDInterpolator,
path: scipy/interpolate/tests/test_interpolate.py
text: @user1:
nit: could you rename `itype` to `interpolator` or some other more descriptive name? I read the `itype(list(zip` line first and confusingly read it as "integer type". 

@author:
I agree. I fixed it.