Mailman 3 January 2003 - NumPy-Discussion

fastest way to make two vectors into an array
by John Hunter 31 Jan '03

31 Jan '03

I have two equal length 1D arrays of 256-4096 complex or floating point numbers which I need to put into a shape=(len(x),2) array. I need to do this a lot, so I would like to use the most efficient means. Currently I am doing: def somefunc(x,y): X = zeros( (len(x),2), typecode=x.typecode()) X[:,0] = x X[:,1] = y do_something_with(X) Is this the fastest way? Thanks, John Hunter

5 5

numarray 0.4 on osX/darwin
by Sebastian Haase 31 Jan '03

31 Jan '03

Hi everybody, I tried a 'python2.2 setup.py install' of numarray on a Mac running os-X (10.1; I have also Fink installed) I starts crunching until: /usr/bin/ld: Undefined symbols: _fclearexcept _fetestexcept Anyone out there, who uses numarray on osX ? I'm thankful for any pointer... Sebastian Haase

2 1

mlab functions: psd, csd, cohere, corrcoef
by John Hunter 30 Jan '03

30 Jan '03

I needed some spectral analysis functions, and finding none available, wrote my own. I use matlab a lot, so I wrote them to be matlab compatible. If you all think these look OK, I'm happy to submit them for inclusion into MLab. ------------------------------------------------------------------- """ Spectral analysis functions for Numerical python written for compatability with matlab commands with the same names. psd - Power spectral density uing Welch's average periodogram csd - Cross spectral density uing Welch's average periodogram cohere - Coherence (normalized cross spectral density) corrcoef - The matrix of correlation coefficients The functions are designed to work for real and complex valued Numeric arrays. One of the major differences between this code and matlab's is that I use functions for 'detrend' and 'window', and matlab uses vectors. This can be easily changed, but I think the functional approach is a bit more elegant. Please send comments, questions and bugs to: Author: John D. Hunter <jdhunter(a)ace.bsd.uchicago.edu> """ from __future__ import division from MLab import mean, hanning, cov from Numeric import zeros, ones, diagonal, transpose, matrixmultiply, \ resize, sqrt, divide, array, Float, Complex, concatenate, \ convolve, dot, conjugate, absolute, arange, reshape from FFT import fft def norm(x): return sqrt(dot(x,x)) def window_hanning(x): return hanning(len(x))*x def window_none(x): return x def detrend_mean(x): return x - mean(x) def detrend_none(x): return x def detrend_linear(x): """Remove the best fit line from x""" # I'm going to regress x on xx=range(len(x)) and return # x - (b*xx+a) xx = arange(len(x), typecode=x.typecode()) X = transpose(array([xx]+[x])) C = cov(X) b = C[0,1]/C[0,0] a = mean(x) - b*mean(xx) return x-(b*xx+a) def psd(x, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning, noverlap=0): """ The power spectral density by Welches average periodogram method. The vector x is divided into NFFT length segments. Each segment is detrended by function detrend and windowed by function window. noperlap gives the length of the overlap between segments. The absolute(fft(segment))**2 of each segment are averaged to compute Pxx, with a scaling to correct for power loss due to windowing. Fs is the sampling frequency. -- NFFT must be a power of 2 -- detrend and window are functions, unlike in matlab where they are vectors. -- if length x < NFFT, it will be zero padded to NFFT Refs: Bendat & Piersol -- Random Data: Analysis and Measurement Procedures, John Wiley & Sons (1986) """ if NFFT % 2: raise ValueError, 'NFFT must be a power of 2' # zero pad x up to NFFT if it is shorter than NFFT if len(x)<NFFT: n = len(x) x = resize(x, (NFFT,)) x[n:] = 0 # for real x, ignore the negative frequencies if x.typecode()==Complex: numFreqs = NFFT else: numFreqs = NFFT//2+1 windowVals = window(ones((NFFT,),x.typecode())) step = NFFT-noverlap ind = range(0,len(x)-NFFT+1,step) n = len(ind) Pxx = zeros((numFreqs,n), Float) # do the ffts of the slices for i in range(n): thisX = x[ind[i]:ind[i]+NFFT] thisX = windowVals*detrend(thisX) fx = absolute(fft(thisX))**2 Pxx[:,i] = fx[:numFreqs] # Scale the spectrum by the norm of the window to compensate for # windowing loss; see Bendat & Piersol Sec 11.5.2 if n>1: Pxx = mean(Pxx,1) Pxx = divide(Pxx, norm(windowVals)**2) freqs = Fs/NFFT*arange(0,numFreqs) return Pxx, freqs def csd(x, y, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning, noverlap=0): """ The cross spectral density Pxy by Welches average periodogram method. The vectors x and y are divided into NFFT length segments. Each segment is detrended by function detrend and windowed by function window. noverlap gives the length of the overlap between segments. The product of the direct FFTs of x and y are averaged over each segment to compute Pxy, with a scaling to correct for power loss due to windowing. Fs is the sampling frequency. NFFT must be a power of 2 Refs: Bendat & Piersol -- Random Data: Analysis and Measurement Procedures, John Wiley & Sons (1986) """ if NFFT % 2: raise ValueError, 'NFFT must be a power of 2' # zero pad x and y up to NFFT if they are shorter than NFFT if len(x)<NFFT: n = len(x) x = resize(x, (NFFT,)) x[n:] = 0 if len(y)<NFFT: n = len(y) y = resize(y, (NFFT,)) y[n:] = 0 # for real x, ignore the negative frequencies if x.typecode()==Complex: numFreqs = NFFT else: numFreqs = NFFT//2+1 windowVals = window(ones((NFFT,),x.typecode())) step = NFFT-noverlap ind = range(0,len(x)-NFFT+1,step) n = len(ind) Pxy = zeros((numFreqs,n), Complex) # do the ffts of the slices for i in range(n): thisX = x[ind[i]:ind[i]+NFFT] thisX = windowVals*detrend(thisX) thisY = y[ind[i]:ind[i]+NFFT] thisY = windowVals*detrend(thisY) fx = fft(thisX) fy = fft(thisY) Pxy[:,i] = fy[:numFreqs]*conjugate(fx[:numFreqs]) # Scale the spectrum by the norm of the window to compensate for # windowing loss; see Bendat & Piersol Sec 11.5.2 if n>1: Pxy = mean(Pxy,1) Pxy = divide(Pxy, norm(windowVals)**2) freqs = Fs/NFFT*arange(0,numFreqs) return Pxy, freqs def cohere(x, y, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning, noverlap=0): """ cohere the coherence between x and y. Coherence is the normalized cross spectral density Cxy = |Pxy|^2/(Pxx*Pyy) The return value is (Cxy, f), where f are the frequencies of the coherence vector. See the docs for psd and csd for information about the function arguments NFFT, detrend, windowm noverlap, as well as the methods used to compute Pxy, Pxx and Pyy. """ Pxx,f = psd(x, NFFT=NFFT, Fs=Fs, detrend=detrend, window=window, noverlap=noverlap) Pyy,f = psd(y, NFFT=NFFT, Fs=Fs, detrend=detrend, window=window, noverlap=noverlap) Pxy,f = csd(x, y, NFFT=NFFT, Fs=Fs, detrend=detrend, window=window, noverlap=noverlap) Cxy = divide(absolute(Pxy)**2, Pxx*Pyy) return Cxy, f def corrcoef(*args): """ corrcoef(X) where X is a matrix returns a matrix of correlation coefficients for each row of X. corrcoef(x,y) where x and y are vectors returns the matrix or correlation coefficients for x and y. Numeric arrays can be real or complex The correlation matrix is defined from the covariance matrix C as r(i,j) = C[i,j] / (C[i,i]*C[j,j]) """ if len(args)==2: X = transpose(array([args[0]]+[args[1]])) elif len(args==1): X = args[0] else: raise RuntimeError, 'Only expecting 1 or 2 arguments' C = cov(X) d = resize(diagonal(C), (2,1)) r = divide(C,sqrt(matrixmultiply(d,transpose(d))))[0,1] try: return r.real except AttributeError: return r ------------------------------------------------------------------- I wrote a little test code comparing the output of matlab's equivalent functions. Basically, I compute the psd or cohere in matlab and python and do the rms difference on the resultant vectors RMS cohere python/matlab difference 0.000854587104587 RMS psd python/matlab difference 0.00210783306638 I am not sure where these differences are arising, but they are quite small. I'm going to keep trying to track them down. For corrcoef, the answers are the same past 8 significant digits. Hope this helps! John Hunter

1 0

Object too deep for desired array
by Karthikesh Raju 30 Jan '03

30 Jan '03

Hi, i was tring out something like this import Numeric import LinearAlgebra import cmath import RandomArray import copy def sMatrix(pd, code, window): if window == 0: nprime = 1 else: nprime = window K, C = Numeric.shape(code) K1, L = Numeric.shape(pd) # check if K == K1 and raise an exception here sCode = Numeric.zeros([nprime*C,K*L*(window+1)],'d') for k in range(K): for l in range(L): code1 = copy.deepcopy(Numeric.array(code[k,0:C-pd[k,l]])) code1.shape = (C-pd[k,l],1) sCode1= Numeric.concatenate((Numeric.zeros([pd[k,l],1]),Numeric.zeros([C*window,1]),code1)) sCode[:, (window+1)*l+window*L*k] = copy.deepcopy(sCode1) return sCode if __name__ == "__main__": pd = Numeric.array([[2]]) code = Numeric.array([[-1,1,-1,1,1]]) np = sMatrix(pd,code,0) print np print "--"*30 np = sMatrix(pd,code,1) print Numeric.shape(np) print np print "--"*30 np = sMatrix(pd,code,2) print Numeric.shape(np) print np print "--"*30 ------------------------------ And i get struck with the following error message:: Traceback (most recent call last): File "sMatrix.py", line 31, in ? np = sMatrix(pd,code,0) File "sMatrix.py", line 24, in sMatrix sCode[:, (window+1)*l+window*L*k] = copy.deepcopy(sCode1) ValueError: Object too deep for desired array ------------ i think it is due to the many deep copy operations taht i am performing. i want to be in a position where slices of matrices should not be references, but should be copies itself and i should be able to move these copies around. (May be it is inefficient, but that is what i did in Matlab and want some compatibility, till i learn more of python and till i migrate to python completely). Is there a way out? Why is this an problem? Am i missing something. Best regards, karthik ----------------------------------------------------------------------- Karthikesh Raju, email: karthik(a)james.hut.fi Researcher, http://www.cis.hut.fi/karthik Helsinki University of Technology, Tel: +358-9-451 5389 Laboratory of Comp. & Info. Sc., Fax: +358-9-451 3277 Department of Computer Sc., P.O Box 5400, FIN 02015 HUT, Espoo, FINLAND -----------------------------------------------------------------------

1 0

Exporting Numpy C functionality to other extension modules
by Jack Jansen 28 Jan '03

28 Jan '03

Actually, wrt my previous message on cobjects for communicating between extension modules, we can do one better! This is an idea I've been toying with for the MacPython extension types, and I think it's applicable to Numeric too. It goes as follows. Each Numeric object has an attribute with a well-known name, lets call it "__Numeric_C_interface". This is a Cobject, and it is shared among all Numeric objects of the same type. The value of this C object is a pointer to a C structure with pointers to all the C routines you might want to call on the object, basically the PyArray_API structure (I think). The descr of the C object is a string with the version number of this particular PyArray_API structure. An extension module that knows about this protocol and gets passed an object that it think might be a Numeric array checks whether the object has an __Numeric_C_interface attribute. If so it retrieves it, checks that it is a Cobject, gets the descriptor and tests it for compatibility and if it is compatible gets the cobject pointer and happily calls all the Numeric routines it needs. -- - Jack Jansen <Jack.Jansen(a)oratrix.com> http://www.cwi.nl/~jack - - If I can't dance I don't want to be part of your revolution -- Emma Goldman -

5 6

enum values visible in numeric types instances?
by Francesc Alted 28 Jan '03

28 Jan '03

Hi, A couple of points related with numarray type objects: 1.- When working with numeric types instances like UInt8 or Float64, is there a way to access to their enumeration NumarrayType C counterpart?. That can be handy when want to map from these objects and integers. For example, right now, I'm forced to use these mappings in Pyrex: # Conversion tables from/to classes to the numarray enum types toenum = {num.Int8:tInt8, num.UInt8:tUInt8, num.Int16:tInt16, num.UInt16:tUInt16, num.Int32:tInt32, num.UInt32:tUInt32, num.Float32:tFloat32, num.Float64:tFloat64, CharType:97 # ascii(97) --> 'a' # Special case (to be corrected) } toclass = {tInt8:num.Int8, tUInt8:num.UInt8, tInt16:num.Int16, tUInt16:num.UInt16, tInt32:num.Int32, tUInt32:num.UInt32, tFloat32:num.Float32, tFloat64:num.Float64, 97:CharType # ascii(97) --> 'a' # Special case (to be corrected) } (yes, Pyrex lets you do that kind of "miracles", like mappings between Python objects and C integers) but if I had this access directly from the object (for example Int8.enumType), my code (and C-extensions in general) could look simpler. 2.- I understand now why Todd was worried about CharArray objects to be assigned to an enumerated type. In fact, if you look at the above maps, I have to map myself this special object as the number 97 (which is the ascii value for character "a"). 97 is ok for now because it can't collide (at least for a while) with other enumeration types. My suggestion is that it would be a good thing to have a reserved enum type for CharArray. And I think that mapping CharArrays with Bool or Int8, would not be a good solution because chararray objects differ in some ways from them, that it would be a mess to distinguish both objects in C-code by just looking at its enumeration type. I don't know, but maybe recarrays also merit a place in enumeration (?). By the way, after the discussion with Todd I finally decided to remove all the Numeric charcodes (and related codes) from PyTables. However, I can still manage Numeric objects by converting them to numarray and accessing the class type with the .type() method. An you know that? the code looks much more logical and neat, and best of all, less error-prone (well, at least I hope so!). I definitely encourage you to do a similar transition in numarray (although I guess that would be more difficult because you still need to Numeric compatibility). Thanks, -- Francesc Alted

1 0

FW: [Numpy-discussion] typecodes in numarray
by Perry Greenfield 27 Jan '03

27 Jan '03

Todd Miller had some further comments that I thought were worth posting as well (and I think he makes some very good points). ************************************************************************ My [i.e. Todd's] thoughts about it: >Maybe I'm becoming a bit tedious with this, but if you look at: > No. It shows you're thinking about it carefully. Having looked at all of the examples below, I have some comments: 1. The sparseness and obscurity of the typecode "wordspace" are both demonstrated here. There are so few letters to choose from, they're often already used in some other context. Even given the large number of unused letters, it's often difficult to choose good ones and to remember what has been chosen. I think this is one of the reasons Perry chose to replace typecodes with true type objects which have rich, regular, and predictable symbolic names. 2. Typecodes were added as a backwards compatability feature of numarray, and I think it's probable that numarray beat Numeric to supporting most of these types, because otherwise they'd have been copied directly and there would be no problem. I'm not really trying to play a blame-game here, but I am making an argument that perhaps numarray should only go so far in the support of what I regard as an obsolescent feature. If the Numeric developers choose to continue extending the use of typecodes in ways that are incompatible with numarray, one way of dealing with it is to "just say no". We are going beyond the scope of backwards compatability to on-going compatabilty. (Which we may still have to do but needs to be discussed and considered) 3. STSCI has layered other software on top of numarray and recarray which astronomers use to do work. It is the friction of that interface which makes correcting these consistency problems more difficult than might be immediately apparent. >I think it's important to agree with a definitive set of charcodes and use >them uniformly throughout numarray. > I wish this were possible, but I'm thinking we should try to find an alternative approach altogether, one which may be more verbose but implicitly free of conflict. A means for specifying a recarray format might be created from tuples, type objects, and integer repetition factors. The verbosity of this approach might be a litte tedious, but it would also be transparent, maintainable, and conflict free. I think we should add an "obsolescent feature" warning to numarray and recarray which flags any use of character typecodes when the appropriate command line switches are set. >Suggestion: if recarray charcodes are not necessary to match the Numeric >ones, I propose that using the Python convention maybe a good idea. >Look at the table in: >http://www.python.org/doc/current/lib/module-struct.html. > This sounds good to me, except that it will break an existing interface that I don't have control over. Therefore, I suggest we correct the problem by coming up with something better.

5 18

Re: FW: [Numpy-discussion] typecodes in numarray
by Todd Miller 27 Jan '03

27 Jan '03

Paul F Dubois wrote: >IMHO you can assume any Python you want. Look to the long term here, not the >short. > You lost me. numarray-0.4 needs at least Python-2.2 or baseclasses don't exist. I had a slow Python equivalent for the baseclass as I refactored prior to numarray-0.4, but it's gone now. > >I'm a bit uncertain on MA as to whether my old design is right. Maybe I >should be inheriting from NDarray? So that MA is more of a sibling of >numarray rather than a wrapper of it? > > I asked Perry about this one. His points (salted a little by me) were: 1. If you inherit from NumArray, you also inherit from NDArray. If you only inherit from NDArray, all you get are the structural operations. 2. If you inherit from NumArray, you can use Liskov substitution to pass MA's directly into extensions expecting NumArrays. This substitution may or may not be good. Also, isinstance(anMA, numarray) will return True. 3. If you inherit from NumArray, you get numerical method definitions which may or may not be applicable to MA. With a little thrashing, we might also get MAs to work for ufuncs. In fact, ufuncs are the key to whether or not the NumArray numerical methods add any value. Todd > >

1 0

complex diagonal matrix
by baecker＠physik.tu-dresden.de 26 Jan '03

26 Jan '03

Hi, I just wondered if there is a "nicer" way of generating a complex diagonal matrix than a) v=arange(10,typecode=Complex) mat=diag(v) b) v=arange(10) mat=diag(v)+0j Namely, wouldn't something like v=arange(10) mat=diag(v,typecode=Complex) be nicer? BTW: I somehow found that in the (excellent) documentation of Numeric the definitions from Mlab.py are a bit hidden. In my case I know nothing about matlab and I somehow expected that this type of routines are to be found in the section (together with zeros,ones etc. etc....) Also diag is not listed in the index http://www.pfdubois.com/numpy/html2/numpy-22.html#A or ? Arnd

3 2

Extensive use of methods instead of functions
by José Fonseca 24 Jan '03

24 Jan '03

With the ability of subclassing types in recent versions of the Python language, more people will be interested in subclassing Numeric arrays for specific purposes. Still the use of functions instead of methods takes away many of the advantages, the ability of being overloaded. Taking this statement as an example: Numeric.put(myarray, myindices, myvalues) In the current state of affairs, if we wanted to have to statment to work with asparse matrix class derived from a Numeric array, it would have to be something like: Sparse.put(myarray, myindices, myvalues) That is, it forces to the underlaying code to know whether is dealing with Numeric arrays, or some other equivalent class. But it would be much more useful to have simply: myarray.put(myindices, myvalues) which would work regardless of the actual type of myarray, provided it supplied the put() method. This would improve enormously code reusability and extensability. I know that there are certain implementations details that may difficult this (like many functions being implemented in pure Python), but any advances made in this since will be an improvement of the current situation. Also, I know that this example is a little unhappy because numarray will do these things with the __getitem__ and __setitem__ operators. But others could easily be shown. Regards, José Fonseca __________________________________________________ Do You Yahoo!? Everything you'll ever need on one web page from News and Sport to Email and Music Charts http://uk.my.yahoo.com

6 6