Correct usage of scipy.interpolate.RegularGridInterpolator

Question

I am a little confused by the documentation for scipy.interpolate.RegularGridInterpolator.

Say for instance I have a function f: R^3 => R which is sampled on the vertices of the unit cube. I would like to interpolate so as to find values inside the cube.

import numpy as np

# Grid points / sample locations
X = np.array([[0,0,0], [0,0,1], [0,1,0], [0,1,1], [1,0,0], [1,0,1], [1,1,0], [1,1,1.]])

# Function values at the grid points
F = np.random.rand(8)

Now, RegularGridInterpolator takes a points argument, and a values argument.

points : tuple of ndarray of float, with shapes (m1, ), ..., (mn, ) The points defining the regular grid in n dimensions.

values : array_like, shape (m1, ..., mn, ...) The data on the regular grid in n dimensions.

I interpret this as being able to call as such:

import scipy.interpolate as irp

rgi = irp.RegularGridInterpolator(X, F)

However, when I do so, I get the following error:

ValueError: There are 8 point arrays, but values has 1 dimensions

What am I misinterpreting in the docs?

Solution 1

Ok I feel silly when I answer my own question, but I found my mistake with help from the documentation of the original regulargrid lib:

https://github.com/JohannesBuchner/regulargrid

points should be a list of arrays that specifies how the points are spaced along each axis.

For example, to take the unit cube as above, I should set:

pts = ( np.array([0,1.]), )*3

or if I had data which was sampled at higher resolution along the last axis, I might set:

pts = ( np.array([0,1.]), np.array([0,1.]), np.array([0,0.5,1.]) )

Finally, values has to be of shape corresponding to the grid laid out implicitly by points. For example,

val_size = map(lambda q: q.shape[0], pts)
vals = np.zeros( val_size )

# make an arbitrary function to test:
func = lambda pt: (pt**2).sum()

# collect func's values at grid pts
for i in range(pts[0].shape[0]):
    for j in range(pts[1].shape[0]):
        for k in range(pts[2].shape[0]):
            vals[i,j,k] = func(np.array([pts[0][i], pts[1][j], pts[2][k]]))

So finally,

rgi = irp.RegularGridInterpolator(points=pts, values=vals)

runs and performs as desired.

Solution 2

Your answer is nicer, and it's perfectly OK for you to accept it. I'm just adding this as an "alternate" way to script it.

import numpy as np
import scipy.interpolate as spint

RGI = spint.RegularGridInterpolator

x = np.linspace(0, 1, 3) #  or  0.5*np.arange(3.) works too

# populate the 3D array of values (re-using x because lazy)
X, Y, Z = np.meshgrid(x, x, x, indexing='ij')
vals = np.sin(X) + np.cos(Y) + np.tan(Z)

# make the interpolator, (list of 1D axes, values at all points)
rgi = RGI(points=[x, x, x], values=vals)  # can also be [x]*3 or (x,)*3

tst = (0.47, 0.49, 0.53)

print rgi(tst)
print np.sin(tst[0]) + np.cos(tst[1]) + np.tan(tst[2])

returns:

1.93765972087
1.92113615659