Never mind this, it was my own mistake as I expected :-) def __chunk(n,size): x = range(0,n,size) x.append(n) return zip(x[:-1],x[1:]) makes it a lot better :) Sturla Den 18.01.2012 06:26, skrev Sturla Molden:

While "playing" with a point-in-polygon test, I have discovered some a failure mode that I cannot make sence of.

The algorithm is vectorized for NumPy from a C and Python implementation I found on the net (see links below). It is written to process a large dataset in chunks. I'm rather happy with it, it can test 100,000 x,y points against a non-convex pentagon in just 50 ms.

Anyway, here is something very strange (or at least I think so):

If I use a small chunk size, it sometimes fails. I know I shouldn't blame it on NumPy, beacuse it is by all likelood my mistake. But it does not make any sence, as the parameter should not affect the computation.

Observed behavior:

1. Processing the whole dataset in one big chunk always works.

2. Processing the dataset in big chunks (e.g. 8192 points) always works.

3. Processing the dataset in small chunks (e.g. 32 points) sometimes fail.

4. Processing the dataset element-wise always work.

5. The scalar version behaves like the numpy version: fine for large chunks, sometimes it fails for small. That is, when list comprehensions is used for chunks. Big list comprehensions always work, small ones might fail.

It looks like the numerical robstness of the alorithm depends on a parameter that has nothing to do with the algorithm at all. For example in (5), we might think that calling a function from a nested loop makes it fail, depending on the length of the inner loop. But calling it from a single loop works just fine.

???

So I wonder:

Could there be a bug in numpy that only shows up only when taking a huge number of short slices?

I don't know... But try it if you care.

In the function "inpolygon", change the call that says __chunk(n,8192) to e.g. __chunk(n,32) to see it fail (or at least it does on my computer, running Enthought 7.2-1 on Win64).

Regards, Sturla Molden

def __inpolygon_scalar(x,y,poly):

# Source code taken from: # http://paulbourke.net/geometry/insidepoly # http://www.ariel.com.au/a/python-point-int-poly.html

n = len(poly) inside = False p1x,p1y = poly[0] xinters = 0 for i in range(n+1): p2x,p2y = poly[i % n] if y> min(p1y,p2y): if y<= max(p1y,p2y): if x<= max(p1x,p2x): if p1y != p2y: xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x if p1x == p2x or x<= xinters: inside = not inside p1x,p1y = p2x,p2y return inside

# the rest is (C) Sturla Molden, 2012 # University of Oslo

def __inpolygon_numpy(x,y,poly): """ numpy vectorized version """ n = len(poly) inside = np.zeros(x.shape[0], dtype=bool) xinters = np.zeros(x.shape[0], dtype=float) p1x,p1y = poly[0] for i in range(n+1): p2x,p2y = poly[i % n] mask = (y> min(p1y,p2y))& (y<= max(p1y,p2y))& (x<= max(p1x,p2x)) if p1y != p2y: xinters[mask] = (y[mask]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x if p1x == p2x: inside[mask] = ~inside[mask] else: mask2 = x[mask]<= xinters[mask] idx, = np.where(mask) idx2, = np.where(mask2) idx = idx[idx2] inside[idx] = ~inside[idx] p1x,p1y = p2x,p2y return inside

def __chunk(n,size): x = range(0,n,size) if (n%size): x.append(n) return zip(x[:-1],x[1:])

def inpolygon(x, y, poly): """ point-in-polygon test x and y are numpy arrays polygon is a list of (x,y) vertex tuples """ if np.isscalar(x) and np.isscalar(y): return __inpolygon_scalar(x, y, poly) else: x = np.asarray(x) y = np.asarray(y) n = x.shape[0] z = np.zeros(n, dtype=bool) for i,j in __chunk(n,8192): # COMPARE WITH __chunk(n,32) ??? if j-i> 1: z[i:j] = __inpolygon_numpy(x[i:j], y[i:j], poly) else: z[i] = __inpolygon_scalar(x[i], y[i], poly) return z

if __name__ == "__main__":

import matplotlib import matplotlib.pyplot as plt from time import clock

n = 100000 polygon = [(0.,.1), (1.,.1), (.5,1.), (0.,.75), (.5,.5), (0.,.1)] xp = [x for x,y in polygon] yp = [y for x,y in polygon] x = np.random.rand(n) y = np.random.rand(n) t0 = clock() inside = inpolygon(x,y,polygon) t1 = clock() print 'elapsed time %.3g ms' % ((t0-t1)*1E3,) plt.figure() plt.plot(x[~inside],y[~inside],'ob', xp, yp, '-g') plt.axis([0,1,0,1]) plt.show()

_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion