Mailman 3 December 2010 - NumPy-Discussion

NEP for faster ufuncs
by Mark Wiebe 06 Jan '11

06 Jan '11

Hello NumPy-ers, After some performance analysis, I've designed and implemented a new iterator designed to speed up ufuncs and allow for easier multi-dimensional iteration. The new code is fairly large, but works quite well already. If some people could read the NEP and give some feedback, that would be great! Here's a link: https://github.com/m-paradox/numpy/blob/mw_neps/doc/neps/new-iterator-ufunc… I would also love it if someone could try building the code and play around with it a bit. The github branch is here: https://github.com/m-paradox/numpy/tree/new_iterator To give a taste of the iterator's functionality, below is an example from the NEP for how to implement a "Lambda UFunc." With just a few lines of code, it's possible to replicate something similar to the numexpr library (numexpr still gets a bigger speedup, though). In the example expression I chose, execution time went from 138ms to 61ms. Hopefully this is a good Christmas present for NumPy. :) Cheers, Mark Here is the definition of the ``luf`` function.:: def luf(lamdaexpr, *args, **kwargs): """Lambda UFunc e.g. c = luf(lambda i,j:i+j, a, b, order='K', casting='safe', buffersize=8192) c = np.empty(...) luf(lambda i,j:i+j, a, b, out=c, order='K', casting='safe', buffersize=8192) """ nargs = len(args) op = args + (kwargs.get('out',None),) it = np.newiter(op, ['buffered','no_inner_iteration'], [['readonly','nbo_aligned']]*nargs + [['writeonly','allocate','no_broadcast']], order=kwargs.get('order','K'), casting=kwargs.get('casting','safe'), buffersize=kwargs.get('buffersize',0)) while not it.finished: it[-1] = lamdaexpr(*it[:-1]) it.iternext() return it.operands[-1] Then, by using ``luf`` instead of straight Python expressions, we can gain some performance from better cache behavior.:: In [2]: a = np.random.random((50,50,50,10)) In [3]: b = np.random.random((50,50,1,10)) In [4]: c = np.random.random((50,50,50,1)) In [5]: timeit 3*a+b-(a/c) 1 loops, best of 3: 138 ms per loop In [6]: timeit luf(lambda a,b,c:3*a+b-(a/c), a, b, c) 10 loops, best of 3: 60.9 ms per loop In [7]: np.all(3*a+b-(a/c) == luf(lambda a,b,c:3*a+b-(a/c), a, b, c)) Out[7]: True

7 27

Giving numpy the ability to multi-iterate excluding an axis
by John Salvatier 04 Jan '11

04 Jan '11

A while ago, I asked a whether it was possible to multi-iterate over several ndarrays but exclude a certain axis( http://www.mail-archive.com/numpy-discussion@scipy.org/msg29204.html), sort of a combination of PyArray_IterAllButAxis and PyArray_MultiIterNew. My goal was to allow creation of relatively complex ufuncs that can allow reduction or directionally dependent computation and still use broadcasting (for example a moving averaging ufunc that can have changing averaging parameters). I didn't get any solutions, which I take to mean that no one knew how to do this. I am thinking about trying to make a numpy patch with this functionality, and I have some questions: 1) How difficult would this kind of task be for someone with non-expert C knowledge and good numpy knowledge? 2) Does anyone have advice on how to do this kind of thing? Best Regards, John

2 9

OS X binaries.
by Gideon 01 Jan '11

01 Jan '11

I noticed that 1.5.1 was released, and sourceforge is suggesting I use the package numpy-1.5.1-py2.6-python.org-macosx10.3.dmg. However, I have an OS X 10.6 machine. Can/should I use this binary? Should I just compile from source?

3 2

Optimization suggestion sought
by Enzo Michelangeli 01 Jan '11

01 Jan '11

For a pivoted algorithm, I have to perform an operation that in fully vectorized form can be expressed as: pivot = tableau[locat,:]/tableau[locat,cand] tableau -= tableau[:,cand:cand+1]*pivot tableau[locat,:] = pivot tableau is a rather large bidimensional array, and I'd like to avoid the allocation of a temporary array of the same size holding the result of the right-hand side expression in the second line of code (the outer product of tableau[:,cand] and pivot). On the other hand, if I replace that line with: for i in xrange(tableau.shape[0]): tableau[i] -= tableau[i,cand]*pivot ...I incur some CPU overhead for the "for" loop -- and this part of code is the botteneck of the whole algorithm. Is there any smarter (i.e., more time-efficient) way of achieving my goal? TIA -- Enzo

4 5

Faster NaN functions
by Erik Rigtorp 31 Dec '10

31 Dec '10

Hi, I just send a pull request for some faster NaN functions, https://github.com/rigtorp/numpy. I implemented the following generalized ufuncs: nansum(), nancumsum(), nanmean(), nanstd() and for fun mean() and std(). It turns out that the generalized ufunc mean() and std() is faster than the current numpy functions. I'm also going to add nanprod(), nancumprod(), nanmax(), nanmin(), nanargmax(), nanargmin(). The current implementation is not optimized in any way and there are probably some speedups possible. I hope we can get this into numpy 2.0, me and people around me seems to have a need for these functions. Erik

3 2

Simple shared arrays
by Erik Rigtorp 31 Dec '10

31 Dec '10

Hi, I was trying to parallelize some algorithms and needed a writable array shared between processes. It turned out to be quite simple and gave a nice speed up almost linear in number of cores. Of course you need to know what you are doing to avoid segfaults and such. But I still think something like this should be included with NumPy for power users. This works by inheriting anonymous mmaped memory. Not sure if this works on windows. import numpy as np import multiprocessing as mp class shared(np.ndarray): """Shared writable array""" def __new__(subtype, shape, interface=None): size = np.prod(shape) if interface == None: buffer = mp.RawArray('d', size) self = np.ndarray.__new__(subtype, shape, float, buffer) else: class Dummy(object): pass buffer = Dummy() buffer.__array_interface__ = interface a = np.asarray(buffer) self = np.ndarray.__new__(subtype, shape=a.shape, buffer=a) return self def __reduce_ex__(self, protocol): return shared, (self.shape, self.__array_interface__) def __reduce__(self): return __reduce_ex__(self, 0) Also see attached file for example usage. Erik

2 2

How to efficiently multiply 2**10 x 2**10 hermitian matrices
by qubax＠gmx.at 30 Dec '10

30 Dec '10

I'll have to work with large hermitian matrices and calculate traces, eigenvalues and perform several matric products. In order to speed those up, i noticed that blas includes a function called 'zhemm' for efficient matrix products with at least one hermitian matrix. is there a way to call that one directly for numpy arrays? are there other, more efficient methods for multiplying that large matrices that one of you might be aware of? especially with the knowledge that they are symmetric/hermitian. i'd appreciate any help in that regard. thanks, q ps: i tried to port the functionality of zhemm into cython, but this is still about a factor of 10 slower than directly using numpy.dot -- There are two things children should get from their parents: roots and wings. The king who needs to remind his people of his rank, is no king. A beggar's mistake harms no one but the beggar. A king's mistake, however, harms everyone but the king. Too often, the measure of power lies not in the number who obey your will, but in the number who suffer your stupidity.

1 0

Why arange has no stop-point opt-in?
by K.-Michael Aye 30 Dec '10

30 Dec '10

Dear all, I'm a bit puzzled that there seems just no way to cleanly code an interval with evenly spaced numbers that includes the stop point given? linspace offers to include the stop point, but arange does not? Am I missing something? (I am aware, that I could do arange(9,15.0001,0.1) but that's what I want to avoid!) Best regards and Happy New Year! Michael

4 5

NumPy C-API equivalent of np.float64()
by Keith Goodman 29 Dec '10

29 Dec '10

I'm looking for the C-API equivalent of the np.float64 function, something that I could use inline in a Cython function. I don't know how to write the function. Anyone have one sitting around? I'd like to use it, if it is faster than np.float64 (np.int32, np.float32, ...) in the Bottleneck package when the output is a scalar, for example bn.median(arr, axis=None).

5 13

fromrecords yields "ValueError: invalid itemsize in generic type tuple"
by Rajat Banerjee 29 Dec '10

29 Dec '10

Hi All, I have been using Numpy for a while with great success. I left my little project for a little while (http://web.mit.edu/stardev/cluster/) and now some of my code is broken. I have some Numpy code to create graphs of activity on a cluster with matplotlib. It ran just fine in July / August 2010, but has since stopped working. I have updated numpy on my machine, I think. In [2]: np.version.version Out[2]: '1.5.1' My call to np.rec.fromrecords() is throwing this exception: File "/home/rajat/Envs/StarCluster/lib/python2.6/site-packages/numpy/core/records.py", line 607, in fromrecords descr = sb.dtype((record, dtype)) ValueError: invalid itemsize in generic type tuple Here is the code with some irrelevant stuff stripped: for line in file: a = [datetime.strptime(parts[0], '%Y-%m-%d %H:%M:%S.%f'), int(parts[1]), int(parts[2]), int(parts[3]), int(parts[4]), int(parts[5]), int(parts[6]), float(parts[7])] list.append(a) file.close() names = ['dt', 'hosts', 'running_jobs', 'queued_jobs',\ 'slots', 'avg_duration', 'avg_wait', 'avg_load'] descriptor = {'names': ('dt,hosts,running_jobs,queued_jobs,slots,avg_duration,avg_wait,avg_load'),\ 'formats' : ('S20','u','u','u','u','u','u','f')} self.records = np.rec.fromrecords(list,','.join(names)) #used to work #self.records = np.rec.fromrecords(list, dtype=descriptor) #new attempt Here is one "line" from the array "list": >>> parts (8) = ['2010-12-07 03:09:46.855712', '2', '2', '177', '2', '86', '370', '1.05']. Neither of those np.rec.fromrecords() calls works. I've tried both separately. They both throw the exact same exception, ValueError: invalid itemsize in generic type tuple Can anybody help me? Am I doing something dumb? Thank you. Rajat Banerjee Masters Candidate, Computer Science, Harvard University

2 1

NumPy-Discussion December 2010