Mailman 3 October 2009 - NumPy-Discussion

vectorized version of logsumexp? (from scipy.maxentropy)
by per freem 22 Oct '09

22 Oct '09

hi all, in my code, i use the function 'logsumexp' from scipy.maxentropy a lot. as far as i can tell, this function has no vectorized version that works on an m-x-n matrix. i might be doing something wrong here, but i found that this function can run extremely slowly if used as follows: i have an array of log probability vectors, such that each column sums to one. i want to simply iterate over each column and renormalize it, using exp(col - logsumexp(col)). here is the code that i used to profile this operation: from scipy import * from numpy import * from numpy.random.mtrand import dirichlet from scipy.maxentropy import logsumexp import time # build an array of probability vectors. each column represents a probability vector. num_vectors = 1000000 log_prob_vectors = transpose(log(dirichlet([1, 1, 1], num_vectors))) # now renormalize each column, using logsumexp norm_prob_vectors = [] t1 = time.time() for n in range(num_vectors): norm_p = exp(log_prob_vectors[:, n] - logsumexp(log_prob_vectors[:, n])) norm_prob_vectors.append(norm_p) t2 = time.time() norm_prob_vectors = array(norm_prob_vectors) print "logsumexp renormalization (%d many times) took %s seconds." %(num_vectors, str(t2-t1)) i found that even with only 100,000 elements, this code takes about 5 seconds: logsumexp renormalization (100000 many times) took 5.07085394859 seconds. with 1 million elements, it becomes prohibitively slow: logsumexp renormalization (1000000 many times) took 70.7815010548 seconds. is there a way to speed this up? most vectorized operations that work on matrices in numpy/scipy are incredibly fast and it seems like a vectorized version of logsumexp should be near instant on this scale. is there a way to rewrite the above snippet so that it's faster? thanks very much for your help.

10 23

ANN: SciPy October Sprint
by Stéfan van der Walt 22 Oct '09

22 Oct '09

Hi all, The weekend is just around the corner, and we're looking forward to the sprint! Here is the detail again: """ Our patch queue keeps getting longer and longer, so here is an opportunity to do some spring cleaning (it's spring in South Africa, at least)! Please join us for an October SciPy sprint: * Date: 24/25 October 2009 (Sat/Sun) * More information: http://projects.scipy.org/scipy/wiki/SciPySprint200910 We are looking for volunteers to write documentation, review code, fix bugs or design marketing material. New contributors are most welcome, and mentoring will be available. """ See you there, Regards Stéfan

1 0

Convolution of a masked array
by Nadav Horesh 22 Oct '09

22 Oct '09

Is there a way to proper convolve a masked array with a normal (nonmasked) array? My specific problem is a convolution of a 2D masked array with a separable kernel (a convolution with 2 1D array along each axis). Nadav.

1 0

Using numpydoc outside of numpy
by Michael Droettboom 21 Oct '09

21 Oct '09

I'm in the process of converting a project to use Sphinx for documentation, and would like to use the Numpy docstring standard with its sections etc. It appears, however, that the numpydoc sphinxext is not installed but merely sits in doc/sphinxext. I see that scipy uses an SVN external to get at this stuff, but I'd prefer not to do that if possible (my project may not live in SVN forever). I see that numpydoc is separately installable (it has a setup.py), but then I assume users of my project who wish to build the docs must download the numpy source and know to explicitly install the numpydoc package. Are there plans to install the numpydoc extension under the numpy install tree somewhere so that other projects can use it, simply by having numpy installed? This is what we do with the plot_directive in matplotlib. Is there a reason that's not a good idea for numpy? Or is there a way for my project to use it that I'm missing? Cheers, Mike -- Michael Droettboom Science Software Branch Operations and Engineering Division Space Telescope Science Institute Operated by AURA for NASA

2 2

GSOC 2010
by Charles R Harris 21 Oct '09

21 Oct '09

Hi All, I don't feel that numpy/scipy did as well in GSOC 2009 as it could have. I think this was mostly due to lack of preparation on our part, we weren't ready when the students started showing up on the lists. So I would like to put together a selection of suitable projects and corresponding mentors that we could put on the wiki somewhere and advertise. Just to start things off, here are two things that come to mind. - Python 3k transition. I think it is time to start looking at this seriously. - Best of breed special functions in cython. These could be part of a separate numpy extras package where code is restricted to C, Cython, and Python. Thoughts? Chuck

6 7

TypeError when calling numpy.kaiser()
by Jeffrey McGee 21 Oct '09

21 Oct '09

Howdy, I'm having trouble getting the kaiser window to work. Anytime I try to call numpy.kaiser(), it throws an exception. Here's the output when I run the example code from http://docs.scipy.org/doc/numpy/reference/generated/numpy.kaiser.html : Python 2.6.2 (release26-maint, Apr 19 2009, 01:56:41) [GCC 4.3.3] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> from numpy import kaiser >>> kaiser(12, 14) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/usr/lib/python2.6/dist-packages/numpy/lib/function_base.py", line 2630, in kaiser return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(beta) File "/usr/lib/python2.6/dist-packages/numpy/lib/function_base.py", line 2507, in i0 y[ind] = _i0_1(x[ind]) TypeError: array cannot be safely cast to required type >>> Is this a bug? Am I doing something wrong? (I'm using the Ubuntu 9.4 packages for python and numpy.) Thanks, Jeff

4 3

object array alignment issues
by Michael Droettboom 20 Oct '09

20 Oct '09

I recently committed a regression test and bugfix for object pointers in record arrays of unaligned size (meaning where each record is not a multiple of sizeof(PyObject **)). For example: a1 = np.zeros((10,), dtype=[('o', 'O'), ('c', 'c')]) a2 = np.zeros((10,), 'S10') # This copying would segfault a1['o'] = a2 http://projects.scipy.org/numpy/ticket/1198 Unfortunately, this unit test has opened up a whole hornet's nest of alignment issues on Solaris. The various reference counting functions (PyArray_INCREF etc.) in refcnt.c all fail on unaligned object pointers, for instance. Interestingly, there are comments in there saying "handles misaligned data" (eg. line 190), but in fact it doesn't, and doesn't look to me like it would. But I won't rule out a mistake in building it on my part. So, how to fix this? One obvious workaround is for users to pass "align=True" to the dtype constructor. This works if the dtype descriptor is a dictionary or comma-separated string. Is there a reason it couldn't be made to work with the string-of-tuples form that I'm missing? It would be marginally more convenient from my application, but that's just a finesse issue. However, perhaps we should try to fix the underlying alignment problems? Unfortunately, it's not clear to me how to resolve them without at least some performance penalty. You either do an alignment check of the pointer, and then memcpy if unaligned, or just always use memcpy. Not sure which is faster, as memcpy may have a fast path already. These are object arrays anyway, so there's plenty of overhead already, and I don't think this would affect regular numerical arrays. If we choose not to fix it, perhaps we should we try to warn when creating an unaligned recarray on platforms where it matters? I do worry about having something that works perfectly well on one platform fail on another. In the meantime, I'll just mark the new regression test to "skip on Solaris". Mike -- Michael Droettboom Science Software Branch Operations and Engineering Division Space Telescope Science Institute Operated by AURA for NASA

8 21

Another suggestion for making numpy's functions generic
by Darren Dale 20 Oct '09

20 Oct '09

numpy's functions, especially ufuncs, have had some ability to support subclasses through the ndarray.__array_wrap__ method, which provides masked arrays or quantities (for example) with an opportunity to set the class and metadata of the output array at the end of an operation. An example is q1 = Quantity(1, 'meter') q2 = Quantity(2, 'meters') numpy.add(q1, q2) # yields Quantity(3, 'meters') At SciPy2009 we committed a change to the numpy trunk that provides a chance to determine the class and some metadata of the output *before* the ufunc performs its calculation, but after output array has been established (and its data is still uninitialized). Consider: q1 = Quantity(1, 'meter') q2 = Quantity(2, 'J') numpy.add(q1, q2, q1) # or equivalently: # q1 += q2 With only __array_wrap__, the attempt to propagate the units happens after q1's data was updated in place, too late to raise an error, the data is now corrupted. __array_prepare__ solves that problem, an exception can be raised in time. Now I'd like to suggest one more improvement to numpy to make its functions more generic. Consider one more example: q1 = Quantity(1, 'meter') q2 = Quantity(2, 'feet') numpy.add(q1, q2) In this case, I'd like an opportunity to operate on the input arrays on the way in to the ufunc, to rescale the second input to meters. I think it would be a hack to try to stuff this capability into __array_prepare__. One form of this particular example is already supported in quantities, "q1 + q2", by overriding the __add__ method to rescale the second input, but there are ufuncs that do not have an associated special method. So I'd like to look into adding another check for a special method, perhaps called __input_prepare__. My time is really tight for the next month, so I'd rather not start if there are strong objections, but otherwise, I'd like to try to try to get it in in time for numpy-1.4. (Has a timeline been established?) I think it will be not too difficult to document this overall scheme: When calling numpy functions: 1) __input_prepare__ provides an opportunity to operate on the inputs to yield versions that are compatible with the operation (they should obviously not be modified in place) 2) the output array is established 3) __array_prepare__ is used to determine the class of the output array, as well as any metadata that needs to be established before the operation proceeds 4) the ufunc performs its operations 5) __array_wrap__ provides an opportunity to update the output array based on the results of the computation Comments, criticisms? If PEP 3124^ were already a part of the standard library, that could serve as the basis for generalizing numpy's functions. But I think the PEP will not be approved in its current form, and it is unclear when and if the author will revisit the proposal. The scheme I'm imagining might be sufficient for our purposes. Darren ^ http://www.python.org/dev/peps/pep-3124/

5 11

opening pickled numarray data with numpy
by dagmar wismeijer 19 Oct '09

19 Oct '09

Hi, I've been trying to open (using numpy) old pickled data files that I once created using numarray, but I keep getting the message that there is no module numarray.generic. Is there any way I could open these datafiles without installing numarray again? Thanks in advance, Dagmar

2 1

user defined types
by Artem Serebriyskiy 19 Oct '09

19 Oct '09

Hello! Would you please give me some examples of open source projects which use the implementation of user defined types for numpy library? (implementation on the C-API level)

2 1

NumPy-Discussion October 2009