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()
Counting the Colors of RGB-Image,
nameit im0 with im0.shape = 2500,3500,3
with this code:
tab0 = zeros( (256,256,256) , dtype=int)
tt = im0.view()
tt.shape = -1,3
for r,g,b in tt:
tab0[r,g,b] += 1
Question:
Is there a faster way in numpy to get this result?
MfG elodw
Hello,
I get a segfault here:
In [1]: x = np.array([1,2,3], dtype='M')
In [2]: x.searchsorted(2, side='left')
But it's fine here:
In [1]: x = np.array([1,2,3], dtype='M')
In [2]: x.view('i8').searchsorted(2, side='left')
Out[2]: 1
This segfaults again:
x.view('i8').searchsorted(np.datetime64(2), side='left')
GDB gets me this far:
Program received signal SIGSEGV, Segmentation fault.
PyArray_SearchSorted (op1=0x1b8dd70, op2=0x17dfac0, side=NPY_SEARCHLEFT) at
numpy/core/src/multiarray/item_selection.c:1463
1463 Py_INCREF(dtype);
Hi,
I am finding it less than useful to have the negative index wrapping on nd-arrays. Here is a short example:
import numpy as np
a = np.zeros((3, 3))
a[:,2] = 1000
print a[0,-1]
print a[0,-1]
print a[-1,-1]
In all cases 1000 is printed out.
What I am after is a way to say "please don't wrap around" and have negative indices behave in a way I choose. I know this is a standard thing - but is there a way to override that behaviour that doesn't involve cython or rolling my own resampler?
Kind Regards,
Nathan.
hi all,
I've done support of discrete variables for interalg - free (license:
BSD) solver with specifiable accuracy, you can take a look at an example
here
It is written in Python + NumPy, and I hope it's speed will be
essentially increased when PyPy (Python with dynamic compilation) support
for NumPy will be done (some parts of code are not vectorized and still
use CPython cycles). Also, NumPy funcs like vstack or append produce only
copy of data, and it also slows the solver very much (for mature
problems).
Maybe some bugs still present somewhere - interalg code already became
very long, but since it already works, you could be interested in trying
to use it right now.
Regards, D.
Hello all,
Is there a recommended (and ideally cross platform)
way to load the frames of a QuickTime movie (*.mov
file) in Python as NumPy arrays? I'd be happy with
an iterator based approach, but random access to
the frames would be a nice bonus.
My aim is to try some image analysis in Python, if
there is any sound in the files I don't care about it.
I had a look at OpenCV which has Python bindings,
http://opencv.willowgarage.com/documentation/python/index.html
however I had no joy compiling this on Mac OS X
with QuickTime support. Is this the best bet?
Thanks,
Peter
This sort of makes sense, but is it the 'correct' behavior?
In [20]: zeros(2, 'S')
Out[20]:
array(['', ''],
dtype='|S1')
It might be more consistent to return '0' instead, as in
In [3]: zeros(2, int).astype('S')
Out[3]:
array(['0', '0'],
dtype='|S24')
Chuck
Hi all,
I have the following problem:
Given a array with dimension Nx3, where N is generally greater than
1.000.000, for each item in this array I have to calculate its density,
Where its density is the number of items from the same array with
distance less than a given r. The items are the rows from the array.
I was not able to think a solution to this using one or two functions of
Numpy. Then I wrote this code http://pastebin.com/iQV0bMNy . The problem
it so slow. So I tried to implement it in Cython, here the result
http://pastebin.com/zTywzjyM , but it is very slow yet.
Is there a better and faster way of doing that? Is there something in my
Cython implementation I can do to perform better?
Thanks!
Hi all,
(I originally posted this to the BayPIGgies list, where Fernando Perez
suggested I send it to the NumPy list as well. My apologies if you're
receiving this email twice.)
I work on a Python/C++ scientific code that runs as a number of
independent Python processes communicating via MPI. Unfortunately, as
some of you may have experienced, module importing does not scale well
in Python/MPI applications. For 32k processes on BlueGene/P, importing
100 trivial C-extension modules takes 5.5 hours, compared to 35
minutes for all other interpreter loading and initialization. We
developed a simple pure-Python module (based on knee.py, a
hierarchical import example) that cuts the import time from 5.5 hours
to 6 minutes.
The code is available here:
https://github.com/langton/MPI_Import
Usage, implementation details, and limitations are described in a
docstring at the beginning of the file (just after the mandatory
legalese).
I've talked with a few people who've faced the same problem and heard
about a variety of approaches, which range from putting all necessary
files in one directory to hacking the interpreter itself so it
distributes the module-loading over MPI. Last summer, I had a student
intern try a few of these approaches. It turned out that the problem
wasn't so much the simultaneous module loads, but rather the huge
number of failed open() calls (ENOENT) as the interpreter tries to
find the module files. In the MPI_Import module, we have rank 0
perform the module lookups and then broadcast the locations to the
rest of the processes. For our real-world scientific applications
written in Python and C++, this has meant that we can start a problem
and actually make computational progress before the batch allocation
ends.
If you try out the code, I'd appreciate any feedback you have:
performance results, bugfixes/feature-additions, or alternate
approaches to solving this problem. Thanks!
-Asher
Hello,
I am looking for some help extracting a subset of data from a large dataset. The data is being read from a wgrib2 (World Meterological Organization standard gridded data) using the pygrib library.
The data values, latitudes and longitudes are in separate lists (arrays?), and I would like a regional subset.
The budget is not very large, but I am hoping that this is pretty simple job. I am just way too green at Python / numpy to know how to proceed, or even what to search for on Google.
If interested, please e-mail jlounds(a)dynamiteinc.com
Thank you!
Jeremy Lounds
DynamiteInc.com
1-877-762-7723, ext 711
Fax: 877-202-3014
