Can scipy's RegularGridInterpolator return both values and gradients with a single call?
<p>I'm using <code>scipy.interpolate.RegularGridInterpolator</code> with <code>method='linear'</code>.  Getting interpolated values is straightforward (see example at <a href="https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html" rel="nofollow noreferrer">https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html</a>).  What is a good way of getting gradients along with interpolated values?</p>

<p>One possibility is to call the interpolator multiple times and compute the gradient "manually" using finite differences.  This feels wasteful given that each call to the interpolator is probably already computing the gradient under the hood.  Is that correct?  If so, how can I modify RegularGridInterpolator to return both the interpolated function value and its gradient?</p>

<p>To be clear, I'm not interested in the "true" gradient of the function which I'm interpolating -- just the gradient of the linear approximation, e.g. the gradient of <code>my_interpolating_function</code> in the example at <a href="https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html" rel="nofollow noreferrer">https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html</a>.</p>

<hr>

<p>Here is an example. I have a function <code>f</code>, I build a linear interpolator <code>f_interp</code>, and I am interested in the gradient of <code>f_interp</code> (as opposed to the gradient of <code>f</code>).  I can compute it using finite differences, but is there a better way?  I assume <code>RegularGridInterpolator</code> is already computing the gradient under the hood -- and doing it fast.  How can I modify it to return the gradient along with the interpolated values?</p>

<pre><code>import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as interp

def f(x, y, z):
    return 0.01*x**2 + 0.05*x**3 + 5*x*y + 3*x*y*z + 0.1*x*(y**2)*z + 9*x*z**2 - 2*y

x_grid = np.linspace(0.0, 10.0, 20)
y_grid = np.linspace(-10.0, 10.0, 20)
z_grid = np.linspace(0.0, 20.0, 20)

mesh = np.meshgrid(x_grid, y_grid, z_grid, indexing="ij")

f_on_grid = f(mesh[0], mesh[1], mesh[2])
assert np.isclose(f_on_grid[0, 0, 0], f(x_grid[0], y_grid[0], z_grid[0]))  # Sanity check

grid = (x_grid, y_grid, z_grid)
f_interp = interp.RegularGridInterpolator(grid, f_on_grid, method="linear",
                                          bounds_error=False, fill_value=None)

dense_x = np.linspace(0.0, 20.0, 400)
plt.plot(dense_x, [f_interp([x, 1.0, 2.0])[0] for x in dense_x], label="y=1.0, z=2.0")
plt.plot(dense_x, [f_interp([x, 1.0, 4.0])[0] for x in dense_x], label="y=1.0, z=4.0")
plt.legend()
plt.show()

f_interp([0.05, 1.0, 2.0])  # Linearly interpolated value, distinct from f(0.05, 1.0, 2.0)

## Suppose I want to compute both f_interp and its gradient at point_of_interest
point_of_interest = np.array([0.23, 1.67, 5.88])
f_interp(point_of_interest)  # Function value -- how can I get f_interp to also return gradient?

## First gradient calculation using np.gradient and a mesh around point_of_interest +- delta
delta = 0.10
delta_mesh = np.meshgrid(*([-delta, 0.0, delta], ) * 3, indexing="ij")
delta_mesh_long = np.column_stack((delta_mesh[0].flatten(),
                                   delta_mesh[1].flatten(),
                                   delta_mesh[2].flatten()))
assert delta_mesh_long.shape[1] == 3
point_plus_delta_mesh = delta_mesh_long + point_of_interest.reshape((1, 3))
values_for_gradient = f_interp(point_plus_delta_mesh).reshape(delta_mesh[0].shape)
gradient = [x[1, 1, 1] for x in np.gradient(values_for_gradient, delta)]
gradient  # Roughly [353.1, 3.8, 25.2]

## Second gradient calculation using finite differences, should give similar result
gradient = np.zeros((3, ))
for idx in [0, 1, 2]:
    point_right = np.copy(point_of_interest)
    point_right[idx] += delta
    point_left = np.copy(point_of_interest)
    point_left[idx] -= delta
    gradient[idx] = (f_interp(point_right)[0] - f_interp(point_left)[0]) / (2*delta)

gradient  # Roughly [353.1, 3.8, 25.2]
</code></pre>

<hr>

<p>Here's a picture of f and f_interp.  I'm interested in the gradient of f_interp (solid lines):</p>

<p><a href="https://i.stack.imgur.com/n61ss.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/n61ss.png" alt="f and f_interp"></a></p>