Reputation: 25
Any suggestions on how to speed up this function?
def smooth_surface(z,c):
hph_arr_list = []
for x in xrange(c,len(z)-(c+1)):
new_arr = np.hstack(z[x-c:x+c])
hph_arr_list.append(np.percentile(new_arr[((new_arr >= np.percentile(new_arr,15)) & (new_arr <= np.percentile(new_arr,85)))],99))
return np.array(map(float,hph_arr_list))
The variable z
has a length of ~15 million and c
is value of window size +
and -
. The function is basically a sliding window that calculates a percentile value per iteration. Any help would be appreciated! z
is an array of arrays (hence the np.hstack
). Maybe any idea if numba would help with this. If so, how to implement?
Upvotes: 1
Views: 322
Reputation: 50668
The slow part of the computation appear to be the line np.percentile(new_arr[((new_arr >= np.percentile(new_arr,15)) & (new_arr <= np.percentile(new_arr,85)))],99)
. This is due to the unexpectedly slow np.percentile
on small arrays as well as the creation of several intermediate arrays.
Since new_arr
is actually quite small, it is much faster to just sort it and make the interpolation yourself. Moreover, numba can also help to speed the computation up.
@njit #Use @njit instead of @jit to increase speed
def filter(arr):
arr = arr.copy() # This line can be removed to modify arr in-place
arr.sort()
lo = int(math.ceil(len(arr)*0.15))
hi = int(len(arr)*0.85)
interp = 0.99 * (hi - 1 - lo)
interp = interp - int(interp)
assert lo <= hi-2
return arr[hi-2]* (1.0 - interp) + arr[hi-1] * interp
This code is 160 times faster with arrays of size 20 on my machine and should produce the same result.
Finally, you can speed up smooth_surface
too by using automatic parallelization in numba (see here for more information). Here is an untested prototype:
@jit(parallel=True)
def smooth_surface(z,c):
hph_arr = np.zeros(len(z)-(c+1)-c)
for x in prange(c,len(z)-(c+1)):
hph_arr[x-c] = filter(np.hstack(z[x-c:x+c]))
return hph_arr
Upvotes: 1