Mailman 3 January 2009 - NumPy-Discussion

2-D function and meshgrid
by Paulo J. S. Silva Jan. 10, 2009

Jan. 10, 2009

Hello, I have a function that receives a array of shape (2,) and returns a number (a function from R^2 -> R). It basically looks like this: def weirdDistance2(x): return dot(dot(weirdMatrix, x), x) (weirdMatrix is a "global" (2,2) array) I want to see its level sets in the box [0, 1] x [0, 1], hence I have to create a meshgrid and then compute it at each point of the mesh: x = linspace(0, 1, 200) y = x.copy() X, Y = meshgrid(x, y) My problem is how to actually compute the function at each point of the mesh. I have come out with two solutions. One very short and clear, but slow, and another longer and more convoluted (it has a loop, I hate loops in numpy code), but faster. Does anyone know a "no explicit-loops" and fast solution? Solution1: def myDistance(a, b): return weirdDistance(np.array((a, b))) vecDistance = np.vectorize(myDistance) return vecDistance(X, Y) Solution 2: nPoints = X.size result = np.zeros(nPoints) points = np.array( [X.ravel(), Y.ravel()] ).T for i in xrange(nPoints): result[i] = weirdDistance(points[i]) result = result.reshape(X.shape) Of course, the first one is slow because the myDistance function creates an array at each call. The second one, even with a loop, avoids the array creations. Best, Paulo

2 2

Numpy performance vs Matlab.
by Nicolas ROUX Jan. 9, 2009

Jan. 9, 2009

Hi, I need help ;-) I have here a testcase which works much faster in Matlab than Numpy. The following code takes less than 0.9sec in Matlab, but 21sec in Python. Numpy is 24 times slower than Matlab ! The big trouble I have is a large team of people within my company is ready to replace Matlab by Numpy/Scipy/Matplotlib, but I have to demonstrate that this kind of Python Code is executed with the same performance than Matlab, without writing C extension. This is becoming a critical point for us. This is a testcase that people would like to see working without any code restructuring. The reasons are: - this way of writing is fairly natural. - the original code which showed me the matlab/Numpy performance differences is much more complex, and can't benefit from broadcasting or other numpy tips (I can later give this code) ...So I really need to use the code below, without restructuring. Numpy/Python code: ##################################################################### import numpy import time print "Start test \n" dim = 3000 a = numpy.zeros((dim,dim,3)) start = time.clock() for i in range(dim): for j in range(dim): a[i,j,0] = a[i,j,1] a[i,j,2] = a[i,j,0] a[i,j,1] = a[i,j,2] end = time.clock() - start print "Test done, %f sec" % end ##################################################################### Matlab code: ##################################################################### 'Start test' dim = 3000; tic; a =zeros(dim,dim,3); for i = 1:dim for j = 1:dim a(i,j,1) = a(i,j,2); a(i,j,2) = a(i,j,1); a(i,j,3) = a(i,j,3); end end toc 'Test done' ##################################################################### Any idea on it ? Did I missed something ? Thanks a lot, in advance for your help. Cheers, Nicolas.

13 32

Replacement for numpy.distutils.config.try_run
by Stéfan van der Walt Jan. 9, 2009

Jan. 9, 2009

Hi, In distutils/config.py line 32 I see: "Usage of try_run is deprecated: please do not \n" \ "use it anymore, and avoid configuration checks \n" \ "involving running executable on the target machine.\n" \ What is the recommended way of doing configuration checks now? Regards Stéfan

3 8

change made to test_print.py
by Christopher Hanley Jan. 9, 2009

Jan. 9, 2009

Hi, I've committed the following change to test_print.py to fix one of the tests. Index: test_print.py =================================================================== --- test_print.py (revision 6302) +++ test_print.py (working copy) @@ -154,7 +154,7 @@ else: locale.setlocale(locale.LC_NUMERIC, 'FRENCH') - assert_equal(str(tp(1.2)), str(float(1.2)), + assert_equal(locale.format("%f",tp(1.2)), locale.format("%f",float(1.2)), err_msg='Failed locale test for type %s' % tp) finally: locale.setlocale(locale.LC_NUMERIC, locale=curloc) Chris -- Christopher Hanley Senior Systems Software Engineer Space Telescope Science Institute 3700 San Martin Drive Baltimore MD, 21218 (410) 338-4338

5 15

replace array values efficiently
by Prashant Saxena Jan. 9, 2009

Jan. 9, 2009

Hi, I am new to numpy and getting my hands on slowly. How do you replace integers from strings in an integer array. (1D) For example: array = {1,1,1,2,3,3,4} replace 1 with "apple" replace 2 with "cheery" replace 3 with "mango" replace 4 with "banana" I know the general solution, but I am looking for an efficient way, supported by numpy/scipy to do this kind of conversion as fast as possible. Thanks Prashant Add more friends to your messenger and enjoy! Go to http://messenger.yahoo.com/invite/

2 1

checking array for NaN values.
by Alexandra Geddes Jan. 9, 2009

Jan. 9, 2009

Is there an easy way to check an array for NaN values? 'all(array)' regards NaN as true (because it's not 0). I've tried all(array != nan) which didn't work either. When i call locations which i know are NaN, it returns -1.#IND. But if i try to call all(array != -1.#IND), python interprets the # as the start of a comment. Any thoughts? thanks, alex. Alexandra Geddes UC Davis

2 1

Accumulate values that are below threshold
by Bevan Jenkins Jan. 8, 2009

Jan. 8, 2009

Hello, Sometimes the hardest part of a problem is articulating it. Hopefully I can describe what I am trying to do - at least enough to get some help. I am trying to compare values to a threshold and when the values are lower than the threshold they are added to the value in my set until the threshold is reached. Everytime the threshold is reached I want the index and value (accumulated). Hopefully the example below will help threshold =1.0 for indx,val in enumerate(Q): print indx,val 0 100.0 1 20.0 2 16.0 3 7.0 4 3.0 5 1.5 6 0.8 7 0.6 8 0.5 9 0.2 10 0.2 11 0.1 12 0.1 The output I would like is (number of elements and value) 0 100.0 1 20.0 2 16.0 3 7.0 4 3.0 5 1.5 7 1.4 11 1.0 The 1st 6 elements are easy as they are all greater than or equal to the threshold(1.0). Once the values drop below the threshold the next value is added until the threshold is reached. Any help is appreciated, Bevan Jenkins

2 2

my cython is slow
by John Hunter Jan. 8, 2009

Jan. 8, 2009

Partly as an excuse to learn cython, and partly because I need to eke out some extra performance of a neighborhood search, I tried to code up a brute force neighborhood search in cython around an N-dimensional point p. I need to incrementally add a point, do a search, add another point, do another search, so some of the algorithms like those in scipy.stats.spatial which assume a static data structure with lots of searches over it probably won't help me. I wrote some cython code to grow a Npoints x Ndimensions numpy array (doubling in size every time I exceed Npoints) and then doing a brute force request for all points within a radius r using a euclidean distance. The code is working fine, but it is still slower than a simple numpy implementation (damn you numpy performance!) Here is the numpy code:: def find_neighbors_numpy(data, point, radius): """ do a plain ol numpy lookup to compare performance and output *data* is a numpoints x numdims array *point* is a numdims length vector radius is the max distance distance return an array of indices into data which are within radius """ numpoints, n = data.shape distance = data - point r = np.sqrt((distance*distance).sum(axis=1)) return np.nonzero(r<=radius)[0] This requires 6 passed through the data, so I should be able to beat it with a properly crafted cython algorithm. I have a class NNBF (Nearest Neighbor Brute Force) with an interface like NUMDIM = 6 nn = nnbf.NNBF(NUMDIM) print 'loading data... this could take a while' # this could take a while for i in range(200000): x = np.random.rand(NUMDIM) nn.add(x) x = np.random.rand(NUMDIM) radius = 0.2 ind = nn.find_neighbors(x, radius) (in my real use case I would be doing a search after every add) when I run this vs numpy, numpy is a little faster testing nnbf... 10 trials: mean=0.0420, min=0.0400 testing numpy... 10 trials: mean=0.0420, min=0.0300 You can grab the code, the python prototype, the test case and setup file at http://matplotlib.sf.net/nnbf.zip. You will need a fairly recent cython to build it: http://www.cython.org/Cython-0.10.3.tar.gz:: # build the extension and run the test code wget http://matplotlib.sf.net/nnbf.zip unzip nnbf.zip cd nnbf python setup.py build_ext --inplace python test_nnbf.py I'm pasting the cython code nnbf.pyx below. If anyone can advise me as to how to remove the bottlenecks (I've decorated the code with some XXX questions) that would be great. I've read the tutorial at http://wiki.cython.org/tutorials/numpy but I think I am still missing something important. I tried to follow Anne's code in ckdtree.pyx which makes use of a raw_data structure but I am not sure how/if to incorporate this into my code -- perhaps there is some juice there.... """ A brute force nearest neighbor routine with incremental add. The internal array data structure grows as you add points """ import numpy as np cimport numpy as np cdef extern from "math.h": float sqrt(float) cdef inline int is_neighbor(int n, double*row, double*pp, double d2max): """ return 1 if the sum-of-squares of n length array row[j]-pp[j] <= d2max """ cdef int j cdef double d, d2 d2 = 0. for j in range(n): d = row[j] - pp[j] d2 += d*d if d2>d2max: return 0 return 1 cdef class NNBF: cdef readonly object data #cdef double* raw_data cdef readonly int n, numrows, numpoints def __init__(self, n): """ create a buffer to hold n dimensional points """ #cdef np.ndarray[double, ndim=2] inner_data self.n = n self.numrows = 100 # XXX how to create empty as contiguous w/o copy? self.data = np.empty((self.numrows, self.n), dtype=np.float) #inner_data = self.data #self.raw_data = <double*>inner_data.data self.numpoints = 0 def add(NNBF self, object point): """ add a point to the buffer, grow if necessary """ #cdef np.ndarray[double, ndim=2] inner_data cdef np.ndarray[double, ndim=1] pp pp = np.asarray(point).astype(np.float) self.data[self.numpoints] = pp self.numpoints += 1 if self.numpoints==self.numrows: ## XXX do I need to do memory management here, eg free ## raw_data if I were using it? self.numrows *= 2 newdata = np.empty((self.numrows, self.n), np.float) newdata[:self.numpoints] = self.data self.data = newdata #self.raw_data = <double*>inner_data.data def get_data(NNBF self): """ return a copy of data added so far as a numpoints x n array """ return self.data[:self.numpoints] def find_neighbors(NNBF self, object point, double radius): """ return a list of indices into data which are within radius from point """ cdef int i, neighbor cdef double d2max cdef np.ndarray[double, ndim=1] pp cdef np.ndarray[double, ndim=1] row if len(point)!=self.n: raise ValueError('Expected a length %d vector'%self.n) pp = np.asarray(point).astype(np.float) d2max = radius*radius neighbors = [] for i in range(self.numpoints): # XXX : is there a more efficient way to access the row # data? Can/should we be using raw_data here? row = self.data[i] neighbor = is_neighbor( self.n, <double*>row.data, <double*>pp.data, d2max) # if the number of points in the cluster is small, the # python list performance should not kill us if neighbor: neighbors.append(i) return neighbors

3 5

ANN: numpy.i - added managed deallocation to ARGOUTVIEW_ARRAY1 (ARGOUTVIEWM_ARRAY1)
by Egor Zindy Jan. 8, 2009

Jan. 8, 2009

Dear List, after I tried to write a simple ARGOUTVIEW_ARRAY1 example (see http://code.google.com/p/ezwidgets/wiki/NumpySWIGMinGW#A_simple_ARGOUTVIEW_… ), I started wondering about memory deallocation. Turns out a few clever people already did all the thinking (see http://blog.enthought.com/?p=62 ) and a few more clever people use this in a swig file (see http://niftilib.sourceforge.net/pynifti, file nifticlib.i). All this concentrated knowledge helped me implement a single ARGOUTVIEWM_ARRAY1 fragment in numpy.i to do some testing (see attached numpy.i). How to use it? In yourfile.i, the %init function has a few lines added to it: %init %{ import_array(); /* initialize the new Python type for memory deallocation */ _MyDeallocType.tp_new = PyType_GenericNew; if (PyType_Ready(&_MyDeallocType) < 0) return; %} ... and that's it! then just use ARGOUTVIEWM_ARRAY1 instead of ARGOUTVIEW_ARRAY1 and python does the deallocation for you when the python array is destroyed (see the examples attached). Everything compiles, but it would be nice to get rid of the warnings... ezalloc_wrap.c:2763: warning: initialization from incompatible pointer type ezalloc_wrap.c: In function `_wrap_alloc_managed': ezalloc_wrap.c:2844: warning: assignment from incompatible pointer type writing build\temp.win32-2.5\Release\_ezalloc.def Compilation: python setup_alloc.py build Testing: The attached test_alloc.py does 2048 allocations of 1MB each for managed and unmanaged arrays. Output on my XP laptop with 1GB RAM as follows: ARGOUTVIEWM_ARRAY1 (managed arrays) - 2048 allocations (1048576 bytes each) Done! ARGOUTVIEW_ARRAY1 (unmanaged, leaking) - 2048 allocations (1048576 bytes each) Step 482 failed TODO: Find a better name for the methods (not sure I like ARGOUTVIEWM_ARRAY), then do the missing fragments (2D arrays), clear the warnings and verify the method for a possible inclusion in the official numpy.i file. Thank you for reading! Regards, Egor

2 6

[Newbie] Fast plotting
by Franck Pommereau Jan. 7, 2009

Jan. 7, 2009

Hi all, and happy new year! I'm new to NumPy and searching a way to compute from a set of points (x,y) the mean value of y values associated to each distinct x value. Each point corresponds to a measure in a benchmark (x = parameter, y = computation time) and I'd like to plot the graph of mean computation time wrt parameter values. (I know how to plot, but not how to compute mean values.) My points are stored as two arrays X, Y (same size). In pure Python, I'd do as follows: s = {} # sum of y values for each distinct x (as keys) n = {} # number of summed values (same keys) for x, y in zip(X, Y) : s[x] = s.get(x, 0.0) + y n[x] = n.get(x, 0) + 1 new_x = array(list(sorted(s))) new_y = array([s[x]/n[x] for x in sorted(s)]) Unfortunately, this code is much too slow because my arrays have millions of elements. But I'm pretty sure that NumPy offers a way to handle this more elegantly and much faster. As a bonus, I'd be happy if the solution would allow me to compute also standard deviation, min, max, etc. Thanks in advance for any help! Franck

6 7