6

I have data z sampled from a 2D function f at grid points x, y, as in z = f(x, y).

It is easy to interpolate f with scipy.interp2d via f = interp2d(x, y, z).

However, evaluating f(x, y) returns an entire 2D grid as if I had done

xx, yy = np.meshgrid(x, y)
f(xx, yy)

The behaviour I want is to simply return the values [f(x[i], y[i]) for i in range(len(x))], which I believe is the behaviour for pretty much any other method in numpy.

The reason I want this behaviour is that I'm looking for the path traced out along the surface of f over "time" by the pair (t, u(t)).

It is also surprising that np.diag(f(t, u(t))) differs from np.array([f(ti, u(ti)) for ti in t]), so It's not clear to me how to get at the path f(t, u(t)) from what is returned via interp2d.

EDIT: About diag, I just thought it seemed we should have np.diag(f(t, u(t))) == np.array([f(ti, u(ti)) for ti in t]), but that isn't the case.

Complete example:

def f(t, u):
    return (t**2) * np.exp((u**2) / (1 + u**2))

x = np.linspace(0, 1, 250)
xx, yy = np.meshgrid(x, x)

z = f(xx, yy)
f = scipy.interpolate.interp2d(x, y, z)

print(f(x, y))
print(np.array([f(xi, yi)[0] for xi, yi in zip(x, y)]))

I would like the output of both print statements to be the same.

RJTK
  • 349
  • 3
  • 10
  • I did not get the last part related to `diag`. Could you, please, provide a workable/reproducible example. Thank you. – AGN Gazer Nov 03 '17 at 02:08
  • If I do the following after your two print statements: `print(np.allclose(np.diag(f(x, y)), np.array([f(xi, yi)[0] for xi, yi in zip(x, y)])))` I get `True`. So I do not get what is the issue here. Also, I had to replace `y->x` in `f = scipy.interpolate.interp2d(x, y, z)`. – AGN Gazer Nov 03 '17 at 02:39

3 Answers3

8

The method interp2d returns an object whose call method expects the x, y vectors to be coordinates of a rectangular grid. And the reason you don't get the desired values from the diagonal of the returned array is that it sorts x, y first.

But there is a workaround, which I also used in Querying multiple points on a bivariate spline in the B-spline basis. After executing 

import scipy.interpolate as si
f = si.interp2d(x, y, z)

evaluate f not by calling it, but by passing its tck properties, followed by your x, y coordinates, to internal bispeu method. Like this:

print(si.dfitpack.bispeu(f.tck[0], f.tck[1], f.tck[2], f.tck[3], f.tck[4], x, y)[0])

The above returns the same as the slow loop

print(np.array([f(xi, yi)[0] for xi, yi in zip(x, y)]))

Explanation

The object f is secretly a B-spline of order 1. The spline parameters (knots, coefficients, order) are contained in its tck property and can be used directly by lower-order routines to the desired effect.

(Ideally, the call method of f would have a Boolean parameter grid which we'd set to False to let it know we don't want grid evaluation. Alas, it's not implemented.)

1

It seems that interp2d() behaves they way it does because that is how the corresponding Fortran function was conceived. The only workaround to this (that I can think of) is to call f on pairs of coordinates:

[f(*p)[0] for p in zip(x, y)]
AGN Gazer
  • 8,025
  • 2
  • 27
  • 45
  • Indeed, but I believe this will be much slower. This also causes problems when I try to pass `f` into an integration routine, since it doesn't play well with mixing scalar and array input. – RJTK Nov 03 '17 at 02:18
  • That is a trivial problem to solve using lambda: `g = lambda x, y: f(x, y)[0]` and pass `g`. – AGN Gazer Nov 03 '17 at 02:22
0

Based on user6655984's suggestion, I've posted the following wrapper function in another thread:

import scipy.interpolate as si
def interp2d_pairs(*args,**kwargs):
    """ Same interface as interp2d but the returned interpolant will evaluate its inputs as pairs of values.
    """
    # Internal function, that evaluates pairs of values, output has the same shape as input
    def interpolant(x,y,f):
        x,y = np.asarray(x), np.asarray(y)
        return (si.dfitpack.bispeu(f.tck[0], f.tck[1], f.tck[2], f.tck[3], f.tck[4], x.ravel(), y.ravel())[0]).reshape(x.shape)
    # Wrapping the scipy interp2 function to call out interpolant instead
    return lambda x,y: interpolant(x,y,si.interp2d(*args,**kwargs))

# Create the interpolant (same interface as interp2d)
f = interp2d_pairs(X,Y,Z,kind='cubic')
# Evaluate the interpolant on each pairs of x and y values
z=f(x,y)
e-malito
  • 878
  • 1
  • 7
  • 14