Mailman 3 January 2012 - NumPy-Discussion

Strange numpy behaviour (bug?)
by Sturla Molden 18 Jan '12

18 Jan '12

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()

1 1

Counting the Colors of RGB-Image
by apo＠pdauf.de 17 Jan '12

17 Jan '12

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

4 3

segfault on searchsorted (1.6.2.dev-396dbb9)
by Adam Klein 17 Jan '12

17 Jan '12

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);

2 1

Negative indexing.
by Nathan Faggian 17 Jan '12

17 Jan '12

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.

5 7

[ANN] global constrained solver with discrete variables
by Dmitrey 16 Jan '12

16 Jan '12

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.

1 0

Loading a Quicktime moive (*.mov) as series of arrays
by Peter 15 Jan '12

15 Jan '12

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

2 1

np.zeros(2, 'S') returns empty strings.
by Charles R Harris 15 Jan '12

15 Jan '12

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

5 7

Calculating density based on distance
by Thiago Franco de Moraes 15 Jan '12

15 Jan '12

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!

4 3

Improving Python+MPI import performance
by Asher Langton 14 Jan '12

14 Jan '12

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

6 15

[JOB] Extracting subset of dataset using latitude and longitude
by Jeremy Lounds 13 Jan '12

13 Jan '12

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

2 1

NumPy-Discussion January 2012