Mailman 3 December 2010 - NumPy-Discussion

bincount question
by Alan G Isaac Jan. 8, 2011

Jan. 8, 2011

:: >>> np.bincount([]) Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: The first argument cannot be empty. Why not? (I.e., why isn't an empty array the right answer?) Thanks, Alan Isaac

2 1

NEP for faster ufuncs
by Mark Wiebe Jan. 6, 2011

Jan. 6, 2011

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 Jan. 4, 2011

Jan. 4, 2011

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 Jan. 1, 2011

Jan. 1, 2011

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 Jan. 1, 2011

Jan. 1, 2011

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 Dec. 31, 2010

Dec. 31, 2010

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 Dec. 31, 2010

Dec. 31, 2010

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 Dec. 30, 2010

Dec. 30, 2010

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 Dec. 30, 2010

Dec. 30, 2010

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 Dec. 29, 2010

Dec. 29, 2010

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