I have developed a code, and I am looking for a more efficient method because it is a lot of slow. Is it possible to modify this code such that it is faster?

The code is really complex to explain, sorry in advance if my explanation is not so satisfactory:

I am working with this big matrix that contains 7 matrices of dimension 50x1000. The code works in this way:

1. I take the first element of each matrix (contained in the big matrix), creating a list of these elements--> [a1b1, a2b1, ... a7b1]
2. Once created this list, I interpolate it, creating a new list with 50 elements
3. Now, I need to repeat the point (1) and (2), but for the first element of second row for all matrices, till the first element of last row.
4. After finished all the elements of the first column, we can switch to the second column of all matrices and repeat the points (1),(2),(3)

If something is not clear, please tell me, I'll try to explain!

import numpy as np
from scipy.interpolate import barycentric_interpolate

matrix = np.random.rand(7, 50, 1000)

X = np.linspace(0.1,0.8,7)
intervals = np.linspace(0.5,1.5,50)

matrix_tot = []
for col in range(len(matrix[0][0])):
    matrix_i = []
    for row in range(len(matrix[0])):
        interp_1 = []
        for m in range(len(matrix)):
            values = matrix[m][row][col]
            interp_1.append(values)
        row_interpolated = barycentric_interpolate(X,np.array(interp_1),intervals)
        matrix_i.append(row_interpolated)
    matrix_tot.append(matrix_i)
matrix_tot = np.array(matrix_tot)
print(matrix_tot.shape)

CodePudding user response：

I get a minor performance boost doing the following

import numpy as np
from scipy.interpolate import barycentric_interpolate
import time

matrix = np.random.rand(7, 50, 1000)
X = np.linspace(0.1,0.8,7)
intervals = np.linspace(0.5,1.5,50)

def interpolate(a):
    return barycentric_interpolate(X, a, intervals)
start = time.process_time()
out = np.apply_along_axis(interpolate, 0, matrix)
print(time.process_time() - start)

Which returns a time of ~3.52 compared to 3.76 seconds on my machine. As a note, the output matrix here is (50, 50, 1000). To get the dimensions of your matrix, simply transpose it.

np.all(out.T == matrix_tot)

CodePudding user response：

Got it. Fortunately, you can provide multiple y vectors in single function call. It seems to be vectorized because performance gain is 2 orders of magnitude.

import numpy as np
from scipy.interpolate import barycentric_interpolate
import time

matrix = np.random.rand(7, 50, 1000)
X = np.linspace(0.1,0.8,7)
intervals = np.linspace(0.5,1.5,50)


start = time.process_time()
matrix_tot = []
for col in range(len(matrix[0][0])):
    matrix_i = []
    for row in range(len(matrix[0])):
        interp_1 = []
        for m in range(len(matrix)):
            values = matrix[m][row][col]
            interp_1.append(values)
        row_interpolated = barycentric_interpolate(X,np.array(interp_1),intervals)
        matrix_i.append(row_interpolated)
    matrix_tot.append(matrix_i)
matrix_tot = np.array(matrix_tot)
print('baseline: ', time.process_time() - start)

## ANSWER HERE:
start = time.process_time()
matrix_reshaped = matrix.reshape(7, 50 * 1000)
matrix_tot2 = barycentric_interpolate(X, matrix_reshaped, intervals).reshape(50, 50, 1000)
# this is only for comparison with matrix_tot, you may want to remove line below:
matrix_tot2 = matrix_tot2.transpose([2, 1, 0])
print('vectorised: ', time.process_time() - start)

assert np.allclose(matrix_tot, matrix_tot2, atol=1e-8)

Result:

baseline:  5.796875
vectorised:  0.015625

And not a single for loop :)