Correct use of scipy.interpolate.RegularGridInterpolator
I'm a little confused by the documentation for scipy.interpolate.RegularGridInterpolator .
Let's say, for example, I have a function f: R ^ 3 => R that is sampled at the vertices of the unit cube. I would like to interpolate to find values inside a 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
accepts an argument points
and an argument values
.
points : a float ndarray tuple, with the forms (m1,), ..., (mn,) Points defining the correct grid in n dimensions.
values : array_like, shape (m1, ..., mn, ...) Regular grid data 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 this, I get the following error:
ValueError: There are 8 point arrays, but the values are 1 dimension
What am I misinterpreting in the docs?
source to share
Ok, I feel stupid when I answer my own question, but I found my error using the documentation of the original regulargrid
lib:
https://github.com/JohannesBuchner/regulargrid
points
there should be a list of arrays that defines how the points are positioned along each axis.
For example, to take the cube of one as above, I have to set:
pts = ( np.array([0,1.]), )*3
or if I had data that was sampled at a higher resolution along the last axis, I could set:
pts = ( np.array([0,1.]), np.array([0,1.]), np.array([0,0.5,1.]) )
Finally, it values
must have a shape that matches the implicit grid 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 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)
performed and performed as desired.
source to share
Your answer is more pleasant, and it suits you well to accept it. I am just adding this as an "alternative" script way.
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
source to share