0

In the code, I have a 2D array(D) and for each column, I want to extract some "k" no of neighboring cols(left and right) and sum them up. A naive approach would be to use a for loop, but to speed up this I am trying to slice the 2D matrix for each column to get a submatrix and sum it column-wise. Surprisingly, the naive approach is faster than using the slicing option for k > 6. Any suggestion on how I can make the implementation efficient?

Naive implementation: `

k = 64
index = np.arange(D.shape[1])
index_kp = index + k 
index_kn = index - k
# neighbors can be less than k if sufficient neighbors not available; for ex. near beginning and the end of an array
index_kn[np.where(index_kn <0)] = np.where(index_kn <0)
index_kp[np.where(index_kp > (len(index)-1))] = np.where(index_kp > (len(index)-1))

Dsmear = np.empty_like(D) #stores the summation of neighboring k columns for each col
for i in range(len(index_kp)):
    Dsmear[:,i] = np.sum(D[:, index_kn[i]:index_kp[i]], axis=1)

`

Slicing implementation:

D1 = np.concatenate((np.repeat(D[:,0].reshape(-1,1),k,axis=1), D, np.repeat(D[:,-1].reshape(-1,1),k,axis=1)),axis=1) #padding the edges with k columns 
idx = np.asarray([np.arange(i-k,i+k+1) for i in range(k, D.shape[1]+k)], dtype=np.int32)
D_broadcast = D1[:, idx] # 3D array; is a bottleneck 
Dsmear = np.sum(D_broadcast, axis=2)
Anup Singh
  • 51
  • 3
  • Look at `np.where(index_kn<0)` by itself. What's point to using that tuple both as a index and as a value. Do you understand that `where` just finds in indices of the nonzero values of `index_kn<0`? It's not an iterator. Make sure you understand each of the pieces of a complex expression. – hpaulj Jul 05 '21 at 19:28
  • By using where expression, I am trying to resolve the edge cases. For example, for column at index 0 and k =30, the neighboring columns can be only 30 columns to the right since there are no 30 columns available to the left. Say for column 200, I can consider 30 cols to the left and right as well. – Anup Singh Jul 05 '21 at 20:39

0 Answers0