Using scipy.interpolate.interpn to interpolate a N-Dimensional array

Question

Suppose I have data that depends on 4 variables: a, b, c and d. I want interpolate to return a 2D array which corresponds to a single value of a and b, and an array of values for c and d. However, The array size need not be the same. To be specific, my data comes from a transistor simulation. The current depends on 4 variables here. I want to plot a parametric variation. The number of points on the parameter is a lot less than the number of points for the horizontal axis.

import numpy as np
from scipy.interpolate import interpn
arr = np.random.random((4,4,4,4))
x1 = np.array([0, 1, 2, 3])
x2 = np.array([0, 10, 20, 30])
x3 = np.array([0, 10, 20, 30])
x4 = np.array([0, .1, .2, .30])
points = (x1, x2, x3, x4)

The following works:

xi = (0.1, 9, np.transpose(np.linspace(0, 30, 4)), np.linspace(0, 0.3, 4))
result = interpn(points, arr, xi)

and so does this:

xi = (0.1, 9, 24, np.linspace(0, 0.3, 4))
result = interpn(points, arr, xi)

but not this:

xi = (0.1, 9, np.transpose(np.linspace(0, 30, 3)), np.linspace(0, 0.3, 4))
result = interpn(points, arr, xi)

As you can see, in the last case, the size of the last two arrays in xi is different. Is this kind of functionality not supported by scipy or am I using interpn incorrectly? I need this create a plot where one of the xi's is a parameter while the other is the horizontal axis.

Answer 1

In a comment above, Ashwith Rego states that he needs to reshape his result as

np.reshape(result,len(interp_x1),len(interp_x2),len(interp_x3),len(interp_x4))

to get the interpn output in the desired shape. This may have the desired shape, but the values are not in the correct places. This may have been overlooked as the length of all 4 dimensions is equal in the above 4D example.

I have made a minimal example to highlight my solution to the problem of correctly shaping the result of interpn. This example also translates a little better to other dimension cases (3D,2D,5D etc...) than some answers above, so I think its a useful addition to the other replies.

import numpy as np
from scipy.interpolate import interpn

# Original data
x1 = np.arange(2)
x2 = np.arange(3)
x3 = np.arange(4)
x4 = np.arange(5)
v = np.reshape(np.arange(len(x1)*len(x2)*len(x3)*len(x4)), (len(x1),len(x2),len(x3),len(x4)))

# New interpolation grid (which for comparison is the same as the origianl)
x1i = x1
x2i = x2
x3i = x3
x4i = x4

# Reshaping trick by stackoverflow user Praveen
interp_mesh = np.array(np.meshgrid(x1i,x2i,x3i,x4i))
ndim = interp_mesh.ndim - 1
interp_points = np.rollaxis(interp_mesh, 0, ndim+1)
interp_points = interp_points.reshape((interp_mesh.size // ndim,ndim))

vi = interpn((x1,x2,x3,x4), v, interp_points) 

# My shaping of the result of interpn
res = np.zeros((len(x1i),len(x2i),len(x3i),len(x4i)))
cnt = 0
for x2i_idx in range(len(x2i)):
    for x1i_idx in range(len(x1i)):
        for x3i_idx in range(len(x3i)):
            for x4i_idx in range(len(x4i)):
                res[x1i_idx,x2i_idx,x3i_idx,x4i_idx] = vi[cnt] 
                cnt += 1

Since the old and new grid in this example is identical, then so are v and vres

print(v-vres)
[[[[0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]]

  [[0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]]

   ...

  [[0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]]

  [[0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]]]]

All elements are 0, which means that v and vres are equal, i.e. that the shape and values are correct. However, printing what was originally suggested gives non-zero values

print(v-np.reshape(vi,len(x1i),len(x2i),len(x3i),len(x4i))
[[[[  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]]

  [[-40. -40. -40. -40. -40.]
   [-40. -40. -40. -40. -40.]
   [-40. -40. -40. -40. -40.]
   [-40. -40. -40. -40. -40.]]

...

  [[ 40.  40.  40.  40.  40.]
   [ 40.  40.  40.  40.  40.]
   [ 40.  40.  40.  40.  40.]
   [ 40.  40.  40.  40.  40.]]

  [[  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]]]]

which means that the result had the correct shape, but with values at the wrong place.

I couldn't figure out a np.roll or np.rolldim or np.reshape() that would reshape vi to (len(x1i),len(x2i),len(x3i),len(x4i), however, it does work. Perhaps someone could write a clever command that is better than my quadruple (!) loop (or N-loop where N is the number of dimensions of the interpolation)?

Answer 2

I'll try to explain this to you in 2D so that you get a better idea of what's happening. First, let's create a linear array to test with.

import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Set up grid and array of values
x1 = np.arange(10)
x2 = np.arange(10)
arr = x1 + x2[:, np.newaxis]

# Set up grid for plotting
X, Y = np.meshgrid(x1, x2)

# Plot the values as a surface plot to depict
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, arr, rstride=1, cstride=1, cmap=cm.jet,
                       linewidth=0, alpha=0.8)
fig.colorbar(surf, shrink=0.5, aspect=5)

This gives us:

Then, let's say you want to interpolate along a line, i.e., one point along the first dimension, but all points along the second dimension. These points are not in the original arrays (x1, x2) obviously. Suppose we want to interpolate to a point x1 = 3.5, which is in between two points on the x1-axis.

from scipy.interpolate import interpn

interp_x = 3.5           # Only one value on the x1-axis
interp_y = np.arange(10) # A range of values on the x2-axis

# Note the following two lines that are used to set up the
# interpolation points as a 10x2 array!
interp_mesh = np.array(np.meshgrid(interp_x, interp_y))
interp_points = np.rollaxis(interp_mesh, 0, 3).reshape((10, 2))

# Perform the interpolation
interp_arr = interpn((x1, x2), arr, interp_points)

# Plot the result
ax.scatter(interp_x * np.ones(interp_y.shape), interp_y, interp_arr, s=20,
           c='k', depthshade=False)
plt.xlabel('x1')
plt.ylabel('x2')

plt.show()

This gives you the result as desired: note that the black points correctly lie on the plane, at an x1 value of 3.5.

Note that most of the "magic", and the answer to your question, lies in these two lines:

interp_mesh = np.array(np.meshgrid(interp_x, interp_y))
interp_points = np.rollaxis(interp_mesh, 0, 3).reshape((10, 2))

I have explained the working of this elsewhere. In short, what it does is to create an array of size 10x2, containing the coordinates of the 10 points you want to interpolate arr at. (The only difference between that post and this one is that I've written that explanation for np.mgrid, which is a shortcut to writing np.meshgrid for a bunch of aranges.)

For your 4x4x4x4 case, you will probably need something like this:

interp_mesh = np.meshgrid([0.1], [9], np.linspace(0, 30, 3),
                          np.linspace(0, 0.3, 4))
interp_points = np.rollaxis(interp_mesh, 0, 5)
interp_points = interp_points.reshape((interp_mesh.size // 4, 4))
result = interpn(points, arr, interp_points)

Hope that helps!