Mailman 3 August 2008 - NumPy-Discussion

Converting audio data to array for further processing
by Dan Colesworthy 26 Aug '08

26 Aug '08

Hello all, I need an efficient way to convert 24 bit signed audio data to a numpy array for further processing. The data will be in a .wav file, and can be recovered via the python wave module. At that point it is a byte string - likely in little endian order. Somewhere in that data will be a 1 second length tone of prescribed but varying amplitude and frequency which must be measured for amplitude, frequency, duration, thd% etc. Also need to make the conversion quickly - over relatively large amounts of data. Something like 5 seconds sampled at 32000 samples/second. This same job will have to be repeated > 300 times as part of a test module. Thanks in advance for any and all ideas. Dan Colesworthy Otologics LLC

3 2

Buildbot failure
by Alan McIntyre 25 Aug '08

25 Aug '08

In case it hasn't been noted yet, three of the buildbots (Windows_XP_x86_64_MSVC, Linux_SPARC_64_Debian, Linux_SPARC_64_Debian_gcc4) are failing to build. The Linux builds have the error: gcc: build/src.linux-sparc64-2.4/numpy/core/src/umathmodule.c numpy/core/src/umathmodule.c.src:332: error: static declaration of 'trunc' follows non-static declaration The MSVC build appears to fail with build\src.win32-2.5\numpy\core\include\numpy\__umath_generated.c(580) : error C2065: 'truncf' : undeclared identifier

2 1

Re: [Numpy-discussion] [SciPy08] Documentation BoF
by Sebastien Binet 25 Aug '08

25 Aug '08

On Monday 25 August 2008 10:36:03 Travis E. Oliphant wrote: > Pauli Virtanen wrote: > > - I vaguely remember someone mentioned having done numpybook -> RST > > conversion. If so, is the result somewhere available? > > Was something done towards this in the Scipy'08 sprints? > > Yes, Gabriel Gellner made progress on this during the sprints, and I > asked him to post his scripts for doing it. I don't have his email > address currently, but perhaps he will respond. Note that a set of latex->reST converters seem to already exist there: http://svn.python.org/view/doctools/converter/ (haven't tried myself, so, YMMV) cheers, sebastien. -- ################################### # Dr. Sebastien Binet # Lawrence Berkeley National Lab. # 1 Cyclotron Road # Berkeley, CA 94720 ###################################

1 0

Re: [Numpy-discussion] subsampling array without loops
by Catherine Moroney 25 Aug '08

25 Aug '08

This almost works. Is there a way to do some masking on tiles, for instance taking the maximum height of each 2x2 square that is an odd number? I've tried playing around with masking and where, but they don't return an array of the original size and shape of "tiles" below. Catherine > Perhaps you could take the relevant slices of the arrays with > appropriate striding: > > a = numpy.arange(128*512).reshape((128, 512)) > top_left = a[::2, ::2] > top_right = a[::2, 1::2] > bottom_left = a[1::2, ::2] > bottom_right = a[1::2, 1::2] > > tiles = numpy.array([top_left, top_right, bottom_left, bottom_right]) > maxes = tiles.max(axis=0) > > Similarly if you want overlapping tiles, you could leave out the > final :2 in the slice specifications above. > > Zach > > > >> I'm looking for a way to acccomplish the following task without lots >> of loops involved, which are really slowing down my code. >> >> I have a 128x512 array which I want to break down into 2x2 squares. >> Then, for each 2x2 square I want to do some simple calculations >> such as finding the maximum value, which I then store in a 64x256 >> array. >> >> Here is the actual code involved. It's only slightly more complex >> than what I described above, since it also involves doing a bit of >> masking on the 2x2 sub-arrays. >> >> Any hints are much appreciated! An inordinate amount of time is >> being spent in this function and another one like it. >> >> Catherine >> >> def calc_sdcm_at_rlra(self,iblock): >> >> npixl = self.nlineht/self.nlinerl >> npixs = self.nsmpht/self.nsmprl >> >> sdcm_out = numpy.array([Constants.CLOUDMASK_NR] \ >> *self.nlinerl*self.nsmprl,'int8') \ >> .reshape(self.nlinerl,self.nsmprl) >> >> for ilinerl in range(0,self.nlinerl): >> for ismprl in range(0,self.nsmprl): >> >> height = self.data[iblock].height >> [ilinerl*2:ilinerl*2+2, \ >> ismprl*2:ismprl*2+2] >> sdcm = self.data[iblock].sdcm[ilinerl*2:ilinerl*2 >> +2, \ >> ismprl*2:ismprl*2+2] >> source = self.data[iblock].heightsrc >> [ilinerl*2:ilinerl*2+2, \ >> >> ismprl*2:ismprl*2+2] >> >> mask1 = (source == Constants.HEIGHT_STEREO) >> mask2 = ( (source == Constants.HEIGHT_SURFACE) | \ >> (source == Constants.HEIGHT_DEFAULT) ) >> >> if (mask1.any()): >> loc = height[mask1].argmax() >> sdcm_out[ilinerl,ismprl] = sdcm[mask1].ravel() >> [loc] >> elif (mask2.any()): >> loc = height[mask2].argmax() >> sdcm_out[ilinerl,ismprl] = sdcm[mask2].ravel() >> [loc] >> >> return sdcm_out >

3 2

Inplace dot and blas access
by Dan Goodman 25 Aug '08

25 Aug '08

Hi all, I'm sorry for what is probably a very simple question but despite looking I didn't manage to find anything. I have some code (solving linear differential equation) which is doing the following operation lots of times: S[:] = dot(A,S)+C here S is a (M,N) matrix, A is (M,M) and C is (M,1). M is usually smallish and N is largish. My question is, can this be turned into an inplace operation to save memory and time? As far as I can see, although some array functions allow you to specify the destination of the output, dot always generates a new array. So I would like to have something like this: dot(A,S, out=S) S+=C Ideally, I could do this in one line. I believe the BLAS DGEMM routine does what I want (if I converted C to a matrix rather than a vector). Is there any way to access this using Python and NumPy (e.g. something like numpy.blas.DGEMM) or would I need to use C code to do this? Many thanks in advance, Dan Goodman

1 0

doing zillions of 3x3 determinants fast
by Daniel Lenski 25 Aug '08

25 Aug '08

Hi all, I need to take the determinants of a large number of 3x3 matrices, in order to determine for each of N points, in which of M tetrahedral cells they lie. I arrange the matrices in an ndarray of shape (N,M,5,3,3). As far as I can tell, Numpy doesn't have a function to do determinants over a specific set of axes... so I can't say: dets = np.linalg.det(a, axis1=-1, axis2=-1) Using an explicit Python loop is very slow... doing the determinants of 100,000 random 3x3 matrices in a loop and inserting the results into a preallocated results array takes about 5.1 seconds. So instead I coded up my own routine to do the determinants over the last 2 axes, based on the naive textbook formula for a 3x3 determinant: def det3(ar): a=ar[...,0,0]; b=ar[...,0,1]; c=ar[...,0,2] d=ar[...,1,0]; e=ar[...,1,1]; f=ar[...,1,2] g=ar[...,2,0]; h=ar[...,2,1]; i=ar[...,2,2] return (a*e*i+b*f*g+c*d*h)-(g*e*c+h*f*a+i*d*b) This is *much much* faster... it does the same 100,000 determinants in 0.07 seconds. And the results differ from linalg.det's results by less than 1 part in 10^15 on average (w/ a std dev of about 1 part in 10^13). Does anyone know of any well-written routines to do lots and lots of determinants of a fixed size? My routine has a couple problems I can see: * it's fast enough for 100,000 determinants, but it bogs due to all the temporary arrays when I try to do 1,000,000 determinants (=72 MB array) * I've heard the numerical stability of the naive determinant algorithms is fairly poor, so I'm reluctant to use this function on real data. Any advice will be appreciated! Dan

4 4

Report from SciPy
by Travis E. Oliphant 25 Aug '08

25 Aug '08

Hi everybody, Robert K, Chuck H, Stefan VdW, Jarrod M, David C, and I had a nice discussion about the future directions of NumPy. We resolved some things and would like community feedback on them if there are opinions. * we will be moving to time-based releases (at least 2 times a year -- November / May) with major changes not accepted about 4 weeks before the release. * The releases will be numbered major.minor.bugfix * There will be no ABI changes in minor releases * There will be no API changes in bugfix releases * Any API changes in minor releases will be done in a backward compatible fashion (possibly with deprecations). * Thus 1.2 will not break ABI compatibility but will add new API features. * We will push the generalized ufuncs off to 1.3 * NumPy 2.0 will be a library and will not automagically import numpy.fft * We will suggest that other libraries use from numpy import fft instead of import numpy as np; np.fft * We will apply most of Andrew Dalke's speed up patches but will keep ctypeslib import * We will remove automatic import of numpy.doc from trunk * NumPy 1.2 will come out shortly after the conference (rc1 on Sunday). * SciPy 0.7b1 will come out after the sprints. If there is anything else I missed from those who were there, please let us all know. By the way, as promised, the NumPy book is now available for download and the source to the book is checked in to the numpy SVN tree: http://svn.scipy.org/svn/numpy/trunk/numpy/doc/numpybook Best regards, -Travis

10 10

Remove ParametricTestCase from numpy.testing
by Stéfan van der Walt 25 Aug '08

25 Aug '08

Hey all, Now that we have switched to Nose as a test framework, we should look into deprecating ParametricTestCase. Alan, have you looked into this before? Cheers Stéfan

3 4

weights parameter of np.average() doesn't work (BUG?)
by Dan Lenski 24 Aug '08

24 Aug '08

Hi all, Is there a good reason why the weights parameter of np.average() doesn't broadcast properly? This is with the Ubuntu Hardy x86_64 numpy package, version 1.0.4. In [293]: a=arange(100).reshape(10,10) # Things work fine when weights have the exact same shape as a In [297]: average(a, axis=1, weights=ones((10,10))) Out[297]: array([ 4.5, 14.5, 24.5, 34.5, 44.5, 54.5, 64.5, 74.5, 84.5, 94.5]) # Bizarre and incorrect result with length-10 weight array In [298]: average(a, axis=1, weights=ones(10)) Out[298]: array([[[[[[[[[ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.], [ 10., 11., 12., 13., 14., 15., 16., 17., 18., 19.], [ 20., 21., 22., 23., 24., 25., 26., 27., 28., 29.], [ 30., 31., 32., 33., 34., 35., 36., 37., 38., 39.], [ 40., 41., 42., 43., 44., 45., 46., 47., 48., 49.], [ 50., 51., 52., 53., 54., 55., 56., 57., 58., 59.], [ 60., 61., 62., 63., 64., 65., 66., 67., 68., 69.], [ 70., 71., 72., 73., 74., 75., 76., 77., 78., 79.], [ 80., 81., 82., 83., 84., 85., 86., 87., 88., 89.], [ 90., 91., 92., 93., 94., 95., 96., 97., 98., 99.] ]]]]]]]]) Doing the weighted-sum explicitly works fine for me: In [311]: sum(a*ones(10), axis=-1)/sum(ones(10)) Out[311]: array([ 4.5, 14.5, 24.5, 34.5, 44.5, 54.5, 64.5, 74.5, 84.5, 94.5]) This seems like a bug, especially since average.__doc__ states that: If weights are given, result is: sum(a * weights,axis) / sum(weights,axis), where the weights must have a's shape or be 1D with length the size of a in the given axis. Integer weights are converted to Float. Not specifying weights is equivalent to specifying weights that are all 1. Frankly, I don't even see why weights is constrained to be 1D or the same shape as a... why not anything that's broadcastable to the same shape as a? Dan

3 3

How to replace a loop for a ufunc (to test optimizations)
by Travis E. Oliphant 23 Aug '08

23 Aug '08

Yesterday at the SciPy conference I suggested it would be a good idea to add a function to replace one of the inner loops in a ufunc so that an extension module could test an optimized inner loop more easily. It turns out that this is such a good idea that it's actually already in the code. The UFUNC_API call that does this is called: PyUFunc_ReplaceLoopBySignature You would use it something like this static void _new_add_loopfunc(char **args, intp *dimensions, intp *steps, void *func) { /* my new faster code */ } Then in the extension module (either in the init or in a method call) something like this: PyUFuncObject *ufunc; PyUFuncGenericFunction oldloopfunc, dummy; int signature[3] /* This should be the number of arguments to the ufunc you are replacing */ PyObject *module, *dict; module = PyImport_ImportModule("umath"); dict = PyModule_GetDict(module); ufunc = PyDict_GetItemString(dict, "add") /* get the ufunc where you want to replace one of the loops */ signature[0] = NPY_DOUBLE; signature[1] = NPY_DOUBLE; signature[2] = NPY_DOUBLE; if (PyUFunc_ReplaceLoopBySignature(ufunc, _new_add_loopfunc, signature, &old1dloop) < 0) { PyErr_SetString(PyExc_RuntimeError, "problem replacing loop"); /* error return i.e. return NULL or return or return -1*/ } Then to restore the original function: if (PyUFunc_ReplaceLoopBySignature(ufunc, old1dloop, signature, &dummy) < 0) { /* error */ } Remember to use import_umath() in the module you are using the UFUNC_API in. Now, lot's of folks should be able to test optimizations. -Travis

1 0