Mailman 3 January 2009 - NumPy-Discussion

Re: [Numpy-discussion] Fast function application on list of 2D points?
by Eric LEBIGOT Jan. 13, 2009

Jan. 13, 2009

Thank you so much for the suggestion, Paulo! Selecting 2D points in a list by creating an array 'mask' of booleans and then using arr[mask, :] is indeed really fast compared to using numpy.apply_along_axis(), in my case (simple "larger than" tests on individual coordinates). I had not realized that you could do "arr[mask, :]": this works great! EOL PS: here are the speed tests I've done on the selection of 2D points from a list, with the following results: filter0: 107.2 s filter1: 0.3 s filter2: 9.7 s filter3: 0.6 s obtained with: #!/usr/bin/env python import numpy def filter0(points): """ Returns only those points that match the filter. """ def filter(p): return (p[0] > 0.5) and (p[1] < 0.5) return points[numpy.apply_along_axis(filter, axis = 1, arr = points)] def filter1(points): """ Returns only those points that match the filter. """ mask = (points[:, 0] > 0.5) & (points[:, 1] < 0.5) return points[mask, :] def filter2(points): """ Returns only those points that match the filter. """ return numpy.array([p for p in points if ((p[0] > 0.5) and p[1] < 0.5)]) def filter3(points): """ Returns only those points that match the filter. """ mask = (points[:, 0] > 0.5) points = points[mask, :] mask = points[:, 1] < 0.5 return points[mask, :] if __name__ == '__main__': import timeit # We generate many random points: NUM_PTS = 1000000 points = numpy.random.random((NUM_PTS, 2)) # We make sure that all the filters give the same result: #print "Initial points:" #print points #print "Filtered points:" #print filter0(points) #print filter1(points) #print filter2(points) #print filter3(points) for filter_num in range(4): func_name = "filter%d" % filter_num t = timeit.Timer("%s(points)" % func_name, "from __main__ import %s, points" % func_name) print "%s: %.1f s" % (func_name, t.timeit(number = 3)) > Date: Mon, 12 Jan 2009 11:33:08 -0300 > From: "Paulo J. S. Silva" <pjssilva(a)ime.usp.br> > Subject: Re: [Numpy-discussion] Fast function application on list of > 2D points? > To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> > Message-ID: <1231770788.6170.3.camel@trinity> > Content-Type: text/plain; charset="UTF-8" > > Why you don't create a mask to select only the points in array that > satisfies the condition on x and y coordinates. For example the code > below applies filter only to the values that have x coordinate bigger > than 0.7 and y coordinate smaller than 0.3: > > mask = numpy.logical_and(points[:,0] > 0.7, points[:,1] < 0.3) > points = numpy.apply_along_axis(filter, axis = 1, arr = points[mask,:]) > > best, > > Paulo > > Em Seg, 2009-01-12 ?s 15:21 +0100, Eric LEBIGOT escreveu: >> Hello, >> >> What is the fastest way of applying a function on a list of 2D points? More >> specifically, I have a list of 2D points, and some do not meet some criteria >> and must be rejected. Even more specifically, the filter only lets through >> points whose x coordinate satisfies some condition, _and_ whose y coordinates >> satisfies another condition (maybe is there room for optimization, here?). >> >> Currently, I use >> >> points = numpy.apply_along_axis(filter, axis = 1, arr = points) >> >> but this creates a bottleneck in my program (array arr may contains 1 million >> points, for instance). >> >> Is there anything that could be faster? >> >> Any suggestion would be much appreciated! >> >> EOL

1 0

formatting issues, locale and co
by David Cournapeau Jan. 12, 2009

Jan. 12, 2009

Hi, While looking at the last failures of numpy trunk on windows for python 2.5 and 2.6, I got into floating point number formatting issues; I got deeper and deeper, and now I am lost. We have several problems: - we are not consistent between platforms, nor are we consistent with python - str(np.float32(a)) is locale dependent, but python str method is not (locale.str is) - formatting of long double does not work on windows because of the broken long double support in mingw. 1 consistency problem: ---------------------- python -c "a = 1e20; print a" -> 1e+020 python26 -c "a = 1e20; print a" -> 1e+20 In numpy, we use PyOS_snprintf for formatting, but python itself uses PyOS_ascii_formatd - which has different behavior on different versions of python. The above behavior can be simply reproduced in C: #include <Python.h> int main() { double x = 1e20; char c[200]; PyOS_ascii_format(c, sizeof(c), "%.12g", x); printf("%s\n", c); printf("%g\n", x); return 0; } On 2.5, this will print: 1e+020 1e+020 But on 2.6, this will print: 1e+20 1e+020 2 locale dependency: -------------------- Another issue is that our own formatting is local dependent, whereas python isn't: import numpy as np import locale locale.setlocale(locale.LC_NUMERIC, 'fr_FR') a = 1.2 print "str(a)", str(a) print "locale.str(a)", locale.str(a) print "str(np.float32(a))", str(np.float32(a)) print "locale.str(np.float32(a))", locale.str(np.float32(a)) Returns: str(a) 1.2 locale.str(a) 1,2 str(np.float32(a)) 1,2 locale.str(np.float32(a)) 1,20000004768 I thought about copying the way python does the formatting in the trunk (where discrepancies between platforms have been fixed), but this is not so easy, because it uses a lot of code from different places - and the code needs to be adapted to float and long double. The other solution would be to do our own formatting, but this does not sound easy: formatting in C is hard. I am not sure about what we should do, if anyone else has any idea ? cheers, David

6 19

Is there anyway to seed random numbers without global state?
by Tom Denniston Jan. 12, 2009

Jan. 12, 2009

I know how to seed and generate random numbers using: numpy.random.seed and numpy.random.rand The problem is the seeding of the random numbers is global which I would think would make it non-thread safe as well as having all the other annoyances of global state like having so set the seed and set it back when done. I would think the ideal would be to be able to build random number objects where each has a seed so that each object's seed is independent. Does such a thing exist in numpy or scipy. Does this even make sense to people? Is there something about random number generation that makes it inherently a global thing? --Tom

3 2

Re: [Numpy-discussion] coding style: citations
by jh＠physics.ucf.edu Jan. 12, 2009

Jan. 12, 2009

Well, best to have the full functionality. These things are intended to be both book components and stand-alone pages. Citation formats are a religious war anyway. --jh-- Date: Mon, 12 Jan 2009 12:41:17 -0500 From: Alan G Isaac <aisaac(a)american.edu> Subject: Re: [Numpy-discussion] coding style: citations To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> Message-ID: <496B80BD.8050102(a)american.edu> Content-Type: text/plain; charset=ISO-8859-1; format=flowed On 1/12/2009 9:08 AM jh(a)physics.ucf.edu apparently wrote: > For citation keys, what's wrong with good old author-year format? > Most scientific journals use it (Abt 1985). > Abt, H. 1985. Harold Abt used to publish surveys of things like > citations when he was ApJ editor in the 1980s but I'm making this > one up just to demonstrate the format except he'd never allow an > article title in a citation because he wanted to save ink which > everyone hated so other journals keep the titles and we should too. > ApJ 123:1-23. There are two problems. 1. The only real problem is that it is not a reST citation, so links to the citation will not be generated, and you cannot control the citation formatting separately. You can still of course use normal internal targets: See `(Abt 1985)`_. Some more text, and later ... .. _(Abt 1985): Abt, H. 1985. Harold Abt used to publish surveys of things like citations when he was ApJ editor in the 1980s but I'm making this one up just to demonstrate the format except he'd never allow an article title in a citation because he wanted to save ink which everyone hated so other journals keep the titles and we should too. ApJ 123:1-23. 2. The failure to use citation keys implies that there will be neither a citation database nor consistent formatting. (Say, if the docs become a book.) I do not suppose this is *currently* a real problem. Cheers, Alan Isaac

1 0

Re: [Numpy-discussion] coding style: citations
by jh＠physics.ucf.edu Jan. 12, 2009

Jan. 12, 2009

> ... citations ... BibTeX Please ensure that, whatever you come up with, it is *very intuitive* for the writer to enter the citation data, or we will quickly lose citations, and possibly writers. For example, I would not have them type in BibTeX entries, but rather provide boxes where they can enter the title, author, etc. and have the system generate the BibTeX. Not everyone uses BibTeX, and the rules about braces, what fields are required, etc., are not totally trivial. Perhaps there's a web site out there that does this already, if BibTex is to be our internal format (which I'm not at all opposed to; I use it myself). If a large majority of our current writers thinks BibTeX is ok, then I'll change my mind. We should poll them on the acceptability of any proposed system before implementing it. For citation keys, what's wrong with good old author-year format? Most scientific journals use it (Abt 1985). --jh-- Abt, H. 1985. Harold Abt used to publish surveys of things like citations when he was ApJ editor in the 1980s but I'm making this one up just to demonstrate the format except he'd never allow an article title in a citation because he wanted to save ink which everyone hated so other journals keep the titles and we should too. ApJ 123:1-23.

2 1

Fast function application on list of 2D points?
by Eric LEBIGOT Jan. 12, 2009

Jan. 12, 2009

Hello, What is the fastest way of applying a function on a list of 2D points? More specifically, I have a list of 2D points, and some do not meet some criteria and must be rejected. Even more specifically, the filter only lets through points whose x coordinate satisfies some condition, _and_ whose y coordinates satisfies another condition (maybe is there room for optimization, here?). Currently, I use points = numpy.apply_along_axis(filter, axis = 1, arr = points) but this creates a bottleneck in my program (array arr may contains 1 million points, for instance). Is there anything that could be faster? Any suggestion would be much appreciated! EOL

2 2

memmap from fd?
by Neal Becker Jan. 12, 2009

Jan. 12, 2009

I'm working on interfacing to a custom FPGA board. The kernel driver exposes the FPGA memory via mmap. It might be nice to use numpy memmap to read/write data. One issue is that I think I will need to create the memmap array from a fd, not a file name. The reason is I wrote the driver to only allow 1 exclusive open, and I already have it open for other reasons. Any chance to create a memmap array from a fd?

2 9

TypeError in dtype.__eq__()
by Bryan Cole Jan. 11, 2009

Jan. 11, 2009

Dtype objects throw an exception if compared for equality against other objects. e.g. >>> import numpy >>> numpy.dtype('uint32')==1 Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: data type not understood >>> After some googling, I think python wisdom (given in the Python docs) says that equality testing should never throw an exception, but return the NotImplemented singleton where the comparison does not make "sense". Throwing an exception breaks membership testing for lists/tuples (using the 'in' operator calls implicitly calls __eq__()). My real-world example is serialisation of numpy arrays using PyYAML. This module should serialise any picklable object including ndarrays. However, it performs a test on each object: if data in [None, ()]: blah... This fails on the dtype object of the array. What's the consensus on this? Is the current dtype behaviour broken? cheers, BC

2 4

memmap problem
by Neal Becker Jan. 10, 2009

Jan. 10, 2009

I modified memmap.py to avoid the issues with needed to read. It is working, but I am seeing these: m Exception exceptions.EnvironmentError: (22, 'Invalid argument') in <bound method eos_memmap.__del__ of eos_memmap([255, 255, 255], dtype=uint8)> ignored Exception exceptions.EnvironmentError: (22, 'Invalid argument') in <bound method eos_memmap.__del__ of eos_memmap([0, 0, 0], dtype=uint8)> ignored Out[22]: eos_memmap([ 0, 0, 0, ..., 255, 255, 255], dtype=uint8) print m[0:4] Exception exceptions.EnvironmentError: (22, 'Invalid argument') in <bound method eos_memmap.__del__ of eos_memmap([0, 0, 0, 0], dtype=uint8)> ignored [0 0 0 0] What's that about?

2 2

help with typemapping a C function to use numpy arrays
by Rich E Jan. 9, 2009

Jan. 9, 2009

Hi list, My question has to do with the Numpy/SWIG typemapping system. I recently got the typemaps in numpy.i to work on most of my C functions that are wrapped using SWIG, if they have arguments of the form (int sizeArray, float *pArray). Now I am trying to figure out how to wrap function that aren't of the form, such as the following function: /*! \brief compute magnitude spectrum of a DFT * * \param sizeMag size of output Magnitude (half of input real FFT) * \param pFReal pointer to input FFT real array (real/imag floats) * \param pFMAg pointer to float array of magnitude spectrum */ void sms_spectrumMag( int sizeMag, float *pInRect, float *pOutMag) { int i, it2; float fReal, fImag; for (i=0; i<sizeMag; i++) { it2 = i << 1; fReal = pInRect[it2]; fImag = pInRect[it2+1]; pOutMag[i] = sqrtf(fReal * fReal + fImag * fImag); } } There are two arrays, one is half the size of the other. But, SWIG doesn't know this, according to the type map it will think *pInRect is of size sizeMag and will not know anything about *pOutMag. Ideally in python, I would like to call the function as sms_spectrumMag(nArray1, nArray2), where nArray1 is twice the size of nArray2, and nArray2 is of size sizeMag. I think in order to do this (although if someone has a better suggestion, I am open to it), I will have to modify the typemap in order to tell SWIG how to call the C function properly. I do not want to have to edit the wrapped C file every time it is regenerated from the interface file. Here is a start I made with the existing typemap code in numpy.i (not working): /* Typemap suite for (DIM_TYPE DIM1, DATA_TYPE* INPLACE_ARRAY1) */ %typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, fragment="NumPy_Macros") (DIM_TYPE DIM1, DATA_TYPE* INPLACE_ARRAY1) { $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), DATA_TYPECODE); } %typemap(in, fragment="NumPy_Fragments") (DIM_TYPE DIM1, DATA_TYPE* INPLACE_ARRAY1) (PyArrayObject* array=NULL, int i=0) { array = obj_to_array_no_conversion($input, DATA_TYPECODE); if (!array || !require_dimensions(array,1) || !require_contiguous(array) || !require_native(array)) SWIG_fail; $1 = 1; for (i=0; i < array_numdims(array); ++i) $1 *= array_size(array,i); $2 = (DATA_TYPE*) array_data(array); } and try to alter it to allow for a conversion of type: (DIM_TYPE DIM1, DATA_TYPE* ARRAY1, DATA_TYPE* ARRAY2) where ARRAY1 is size DIM1 * 2 and ARRAY2 is size DIM1. Then I can %apply this to my function that I mentioned in the last post. So here are my first two questions: 1) where is DIM1 used to declare the array size? I don't see where it is used at all, and I need to somewhere multiply it by 2 to declare the size of ARRAY1 2) I am not understanding where $input comes from, so I do not understand how to distinguish between ARRAY1 and ARRAY2. In the attempt I have already tried, I think I just use the pointer to ARRAY1 twice. If anyone has suggestions on how to solve this problem, thanks! regards, Rich

3 11