Vectorized operations in NumPy-CodePudding

I have a numpy array of 5 numpy arrays, each with 5 observations. I need to use the trapezoid rule of integration to get 5 results of integrating. I'm trying to use for loop for apply a function to each numpy array and looking for something that has faster implementation.

def trapezoid_rule(a,b,n,f_nparray):
  h = (b-a) / (n-1)
  return (h/2) * (f_nparray[0]   2 * sum(f_nparray[1:n-1])   f_nparray[n-1])

simulation_res = np.array([[1,2,3,4,5], [1,3,5,7,9], [1,4,2,5,7], [1,5,2,6,4], [6,2,5,3,4]])

# List of integration results
I = []

for i in simulation_res:
  I.append(trapezoid_rule(0,1,10,i))

Expected output format = [a,b,c,d,e]

CodePudding user response：

You can do this without for loops, but I'll show both. numpy has a built-in trapezoidal rule that works in a vectorized way. If you wanted to use it:

import numpy as np

simulation_res = np.array([[1,2,3,4,5], 
                           [1,3,5,7,9], 
                           [1,4,2,5,7], 
                           [1,5,2,6,4], 
                           [6,2,5,3,4]])
result = np.trapz(simulation_res, axis=1)
print(result) 
# Gives [12.  20.  15.  15.5 15. ]

However if the loop is required for some reason, you could just preallocate some memory in a numpy array and append it

result = np.zeros(5)
for idx, i in enumerate(simulation_res):
    result[idx] = np.trapz(i)
print(result)
# Gives [12.  20.  15.  15.5 15. ]

CodePudding user response：

With:

def foo(a,b,n,f_nparray):
  h = (b-a) / (n-1)
  return (h/2) * (f_nparray[:,0]   2 * np.sum(f_nparray[:,1:n-1], axis=1)   f_nparray[:,n-1])

In [20]: foo(0,1,5,simulation_res)
Out[20]: array([3.   , 5.   , 3.75 , 3.875, 3.75 ])

which matches your I.

I had to add axis to np.sum; and in contrast to my comment, just use 0 for the indexing. I just need one value per row of the 2d array.

Compared to np.trapz:

In [29]: np.trapz(simulation_res, axis=1)/4
Out[29]: array([3.   , 5.   , 3.75 , 3.875, 3.75 ])

In [30]: timeit np.trapz(simulation_res, axis=1)/4
34 µs ± 1.66 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [31]: timeit foo(0,1,5,simulation_res)
25.7 µs ± 139 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)