Vectorize 1D interpolation over multidimensional array of y coordinates

user_6177601 0
Viewed 142

I have a 3D array that I want to expand into a new 3D array. The data vectors in the 3rd axis are sorted data, corresponding to a fixed array of x coordinates. I want to expand the third axis along a new x coordinate array, with the values linearly interpolated from third axis of the original array, over all vectors in the original array. In other words, I have to apply linear interpolation many times to create the desired result. I am looking for a fast, possibly vectorized solution. An example using dummy data of the nested for-loop implementation is shown below. This solution is slow for as my original data set is very large

import numpy as np

fp = np.sort(np.random.rand(1000,100,10), axis = 2)
xp = np.linspace(0.0, 1.0, num=10)
x = np.linspace(0.0, 1.0, num=20)

result=np.zeros((1000,100,20))
for i in range(1000):
    for j in range(100):
        result[i,j,:] = np.interp(x,xp,fp[i,j,:])

Is there a faster, more efficient way to do this without the for loops?

vectorization numpy python
1 Answers

This isn't the most elegant solution, but it gets the job done. The basic idea is to flatten fp into a 1D array, and then create a monotonic increasing sequence of x values to interpolate over. So the first 10 values of fp.reshape(-1) are interpolated over the interval [0,1], the next 10 over [2,3], and so on. This eliminates the loop and calls interp only once.

import numpy as np

fp = np.sort(np.random.rand(1000,100,10), axis = 2)
xp = np.linspace(0.0, 1.0, num=10)
x = np.linspace(0.0, 1.0, num=20)

bigxp = np.tile(xp, int(np.product(fp.shape) / 10)).reshape(-1,10)
offset = np.arange(0,200000,2).reshape(-1,1)
bigxp += offset

bigx = np.tile(x, int(np.product(fp.shape) / 10)).reshape(-1,20)
bigx += offset

result = np.interp(bigx.reshape(-1), 
                   bigxp.reshape(-1),
                   fp.reshape(-1))

result = result.reshape(1000,100,20)
user_84185250
Related