Given two opposite corners of a rectangle (x1, y1)
and (x2, y2)
and two radii r1
and r2
, find the ratio of points that lie between the circles defined by the radii r1
and r2
to the total number of points in the rectangle.
Simple NumPy approach:
def func_1(x1,y1,x2,y2,r1,r2,n):
x11,y11 = np.meshgrid(np.linspace(x1,x2,n),np.linspace(y1,y2,n))
z1 = np.sqrt(x11**2+y11**2)
a = np.where((z1>(r1)) & (z1<(r2)))
fill_factor = len(a[0])/(n*n)
return fill_factor
Next I tried to optimize this function with the jit
decorator from numba. When I use:
nopython = True
The function is faster and gives the right output. But when I also add:
parallel = True
The function is faster but gives the wrong result.
I know that this has something to do with my z
matrix since that is not being updated properly.
@jit(nopython=True,parallel=True)
def func_2(x1,y1,x2,y2,r1,r2,n):
x_ = np.linspace(x1,x2,n)
y_ = np.linspace(y1,y2,n)
z1 = np.zeros((n,n))
for i in range(n):
for j in range(n):
z1[i][j] = np.sqrt((x_[i]*x_[i]+y_[j]*y_[j]))
a = np.where((z1>(r1)) & (z1<(r2)))
fill_factor = len(a[0])/(n*n)
return fill_factor
Test values :
x1 = 1.0
x2 = -1.0
y1 = 1.0
y2 = -1.0
r1 = 0.5
r2 = 0.75
n = 25000
Additional info : Python version : 3.6.1, Numba version : 0.34.0+5.g1762237, NumPy version : 1.13.1
The problem with
parallel=True
is that it's a black-box. Numba doesn't even guarantee that it will actually parallelize anything. It uses heuristics to find out if it's parallelizable and what could be done in parallel. These can fail and in your example they do fail, just like in my experiments withparallel
and numba. That makesparallel
untrustworthy and I would advise against using it!In newer versions (0.34)
prange
was added an you could have more luck with that. It can't be applied in this case becauseprange
works likerange
and that's different fromnp.linspace
...Just a note: You can avoid building
z
and doing thenp.where
in your function completely, you could just do the checks explicitly:That should also provide some speedup compared to your function, maybe even more than using
parallel=True
(if it would work correctly).