Trimming a numpy array in cython

user_10906987 0
Viewed 91

Currently I have the following cython function, modifies entries of a numpy array filled with zeros to sum non-zero values. Before I return the array, I would like to trim it and remove all the non-zero entries. At the moment, I use the numpy function myarray = myarray[~np.all(myarray == 0, axis=1)] to do so. I was wondering if there is (in general) a faster way to do this using a Cython/C function instead of relying on python/numpy. This is one of the last bits of pythonic interactions in my script (checked by using to %%cython -a). But I don't really know how to proceed with this problem. In general, i don't know a priori the number of nonzero elements in the final array.

cdef func():
   np.ndarray[np.float64_t, ndim=2] myarray = np.zeros((lenpropen, 6)) 
   """
   computations
   """

   myarray = myarray[~np.all(myarray == 0, axis=1)]
   return myarray

Thanks to @Jérôme Richard for his comments. based on that (if my understanding is correct) I tried to implement the erase-remove idiom. Sample code given below.

myarray = np.zeros((5000,6))
myarray[2] = [1,1,1,1,1,1]
@cython.boundscheck(False)  # Deactivate bounds checking                                                                  
@cython.wraparound(False)   # Deactivate negative indexing.                                                               
@cython.cdivision(True)     # Deactivate division by 0 checking.
cdef erase_remove( np.ndarray[np.float64_t, ndim=2] myarray):
    cdef int idx 
    cdef int cursor = 0
    cdef int length_arr = 5000
    for idx in range(5000):
    
        if myarray[idx,0]!=0 and myarray[idx,1]!=0 and myarray[idx,2]!=0 and myarray[idx,3]!=0 and myarray[idx,4]!=0 and  myarray[idx,5]!=0:
            myarray[cursor,0] = myarray[idx,0]
            myarray[cursor,1] = myarray[idx,1]
            myarray[cursor,2] = myarray[idx,2]
            myarray[cursor,3] = myarray[idx,3]
            myarray[cursor,4] = myarray[idx,4]
            myarray[cursor,5] = myarray[idx,5]
            cursor = cursor +1
        else:
            continue
    return  myarray[0:cursor]    
start = timer()
myarray= erase_remove(myarray)
end = timer()
print("final", myarray)
print("time", end-start)

This yields the output

final [[1. 1. 1. 1. 1. 1.]]
time 1.1235475540161133e-05

Compared to

myarray = np.zeros((5000,6))
print(myarray)

myarray[2] = [1,1,1,1,1,1]
print(myarray)
start = timer()
myarray = myarray[~np.all(myarray == 0, axis=1)]
end = timer()
print(myarray)
print("time", end-start)

which yields output

[[1. 1. 1. 1. 1. 1.]]
time 0.0006445050239562988
cython numpy python
1 Answers

If the highest dimension contains always a small number of element like 6, then your code is not the best one.

First of all, myarray == 0, np.all and ~ creates temporary arrays that introduces some additional overhead as they needs to be written and read back. The overhead is dependent of the this of the temporary array and the biggest one is myarray == 0.

Moreover, Numpy calls perform some unwanted checks that Cython is not able to remove. These checks introduce a constant time overhead. Thus, is can be quite big for small input arrays but not big input arrays.

Additionally, the code of np.all can be faster if it would know the exact size of the last dimension which is not the case here. Indeed, the loop of np.all could theoretically be unrolled since the last dimension is small. Unfortunately, Cython does not optimize Numpy calls and Numpy is compiled for a variable input size, so not known at compile-time.

Finally, the computation can be parallelized if lenpropen is huge (otherwise this will not be faster and could actually be slower). However, note that a parallel implementation requires the computation to be done in two steps: np.all(myarray == 0, axis=1) needs to be computed in parallel and then you can create the resulting array and write it by computing myarray[~result] in parallel. In sequential, you can directly overwrite myarray by filtering lines in-place and then produce a view of the filtered lines. This pattern is known as the erase-remove idiom. Note that this assume the array is contiguous.

To conclude, a faster implementation consists writing 2 nested loops iterating on myarray with a constant number of iterations for the innermost one. Regarding the size of lenpropen, you can either use a sequential in-place implementation base on the erase-remove idiom, or a parallel out-of-place implementation with two steps (and a temporary array).

user_129395570
Related