Mailman 3 September 2003 - NumPy-Discussion

NA_ReturnOutput question
by Peter Verveer 09 Sep '03

09 Sep '03

The documentation states that NA_ReturnOutput(numarray, shadow) should be used with the value returned by NA_OptionalOutputArray() for the shadow parameter. Is it also possible to use a shadow array allocated otherwise, e.g. by NA_NewArray()? Cheers, Peter -- Dr. Peter J. Verveer Cell Biology and Cell Biophysics Programme EMBL Meyerhofstrasse 1 D-69117 Heidelberg Germany Tel. : +49 6221 387245 Fax : +49 6221 387242 Email: Peter.Verveer(a)embl-heidelberg.de

2 1

Re: [Numpy-discussion] Attempting to sub-class NumArray
by Colin J. Williams 07 Sep '03

07 Sep '03

Todd Miller wrote: >I think the problem you are having is caused by the assumption that >instance initialization is equivalent to an unordered assignment of all >instance attributes. That's not true because the private attributes are >interrelated; assigning "shape" affects "_shape" and "_strides", >_itemsize must be known before _strides can be computed from shape, and >so on. > It would be helpful if you could let us know what assignment ordering would be appropriate, or, if there were information about the interactions, perhaps we could work these things out for ourselves. As we've discussed before, there is little information available on private objects and their roles. Even to have the source code of some of the C code might be useful. I wasn't able to spot the relevant C code among the CVS directories. However, the assumption would appear to apply to the crude initialization of 'special'. The problem which is being reported seems to be connected with the function setshape in generic.py. This is far from the unordered assignment, although I recognize that ther might be some connection. Colin W. > >On Fri, 2003-09-05 at 20:32, Colin J. Williams wrote: > > >>Todd, >> >>Thanks again for your prompt response. >>Todd Miller wrote: >> >> >>>Hi Colin, >>> >>>I haven't looked at this in detail, but special.__init__() is making me >>>queasy. >>> >>> >>We don't call __init__ directly, but only when we construct the >>instance, a, of the class special. >> >> >>>Why doesn't it call NumArray.__init__()? >>> >>> >>Because we have data to install in the newly created a. We go off to >>the factory function, numarray.array, to create the array, which we >>then specialize. >> >>I don't see how NumArray.__init__ would help. This cumbersome >>procedure seems to work OK in most respects, but fails to display >>properly. >> >>I don't follow the code sufficiently well to be sure, but my hunch is >>that the problem lies with setShape. It fails to update _strides. >> >>Colin W. >> >> >>>Todd >>> >>>On Fri, 2003-09-05 at 15:03, Colin J. Williams wrote: >>> >>> >>> >>>>When ravel() is used on an instance of NumArray, the _strides attribute >>>>is recalculated correctly. >>>>When ravel() is used on an instance of special, a sub-class of NumArray, >>>>the _strides attribute keeps its old value. This appears to cause the >>>>display problem which is reported. >>>> >>>>The essentials of the code used are below. >>>> >>>>To investigate this problem, two print statements were added to the >>>>method setshape in generic.py: >>>> >>>> if newnelements == nelements: >>>> self._shape = tuple(shape) >>>> print 'In setShape self._shape:', self._shape >>>> self._strides = self._stridesFromShape() >>>> print 'In setShape self._strides:', self._strides >>>> else: >>>> raise ValueError("New shape is not consistent with the >>>> old shape") >>>> >>>>It seems that only the first of these print statements is executed with >>>>special, but both are executed with NumArray. In the second case, the >>>>output is not in the expected sequence. >>>> >>>>I would appreciate any suggestion as to a workaround. >>>> >>>>Colin W, >>>> >>>>Grief when attempting to sub-class NumArray >>>> >>>>Code used: >>>> >>>># Try to sub-class >>>>class special(N.NumArray): >>>> def __init__(self, data= None, shape= None, eType= None): >>>>## eType= _nt.Any or AnyType not acceptable %% >>>> arr= N.array(sequence= data, shape= shape, type= eType) >>>> for attr in dir(arr): >>>>## This is a longwinded way of setting up special >>>>## There must be a better way %% >>>>## Perhaps we should use generic._view >>>>## It's not documented TRY IT LATER >>>>## There is no need to transfer methods! >>>> if attr[0] == '_' and attr[1] != '_': >>>> print attr >>>> exec 'self.' + attr + '= ' 'arr.' + attr >>>> exec 'print self.' + attr >>>> pass >>>> >>>>a= special(data= [1, 2, 3, 4], shape= (2, 2)) >>>>a= N.array([1, 2, 3, 4], shape= (2, 2)) >>>>print 'In ts a._strides:', a._strides >>>>b= _gen.ravel(a) >>>>print 'In ts b._strides:', b._strides # <<< Unchanged with >>>>special, OK with NumArray >>>>print 'b:', b >>>> >>>> >>>> >>>>USING THE SUBCLASS special >>>> >C:\Progra~1\Python23\pythonw -u ts.py >>>>In ts a._strides: (8, 4) >>>>In setShape self._shape: (4,) >>>>In ts b._strides: (8, 4) >>>>b: In setShape self._shape: (4,) >>>>Traceback (most recent call last): >>>> File "ts.py", line 37, in ? >>>> print 'b:', b >>>> File >>>>"C:\PROGRA~1\Python23\lib\site-packages\numarray\numarraycore.py", line >>>>618, in __str__ >>>> MAX_LINE_WIDTH, PRECISION, SUPPRESS_SMALL, ' ', "") >>>> File "C:\PROGRA~1\Python23\lib\site-packages\numarray\arrayprint.py", >>>>line 176, in array2string >>>> separator, prefix) >>>> File "C:\PROGRA~1\Python23\lib\site-packages\numarray\arrayprint.py", >>>>line 138, in _array2string >>>> max_str_len = max(len(str(max_reduce(data))), >>>>libnumarray.error: maximum_Int32_reduce: access beyond buffer. offset=27 >>>>buffersize=16 >>>> >Exit code: 1 >>>>__________________________________________________________ >>>> >>>># USING NumArray >>>> >>>>In setShape self._shape: (2, 2) >>>>In ts a._strides: (8, 4) >>>>In setShape self._shape: (4,) >>>>In ts b._strides: (4,) >>>>b: In setShape self._shape: (4,) >>>>[1 2 3 4] >>>> >Exit code: 0 >>>> >>>> >>>> >>>> >>>>------------------------------------------------------- >>>>This sf.net email is sponsored by:ThinkGeek >>>>Welcome to geek heaven. >>>>http://thinkgeek.com/sf >>>>_______________________________________________ >>>>Numpy-discussion mailing list >>>>Numpy-discussion(a)lists.sourceforge.net >>>>https://lists.sourceforge.net/lists/listinfo/numpy-discussion >>>> >>>> >>>>

2 2

Matrix square root
by Andrew Nesbit 06 Sep '03

06 Sep '03

I need a function equivalent to Matlab's sqrtm, i.e., a square root for matrices. I constructed the following definition: def sqrtm(M): (U,S,VT) = LinearAlgebra.singular_value_decomposition(M) D = MLab.diag(sqrt(S)) return matrixmultiply(matrixmultiply(U,D),VT) but this technique only works for Hermitian, positive definite matrices. The matrices I operate on don't necessarily satisfy that condition. Does anybody know of a more general technique, preferably with an accompanying Python implementation? Thanks. Andrew.

4 4

Attempting to sub-class NumArray
by Colin J. Williams 05 Sep '03

05 Sep '03

When ravel() is used on an instance of NumArray, the _strides attribute is recalculated correctly. When ravel() is used on an instance of special, a sub-class of NumArray, the _strides attribute keeps its old value. This appears to cause the display problem which is reported. The essentials of the code used are below. To investigate this problem, two print statements were added to the method setshape in generic.py: if newnelements == nelements: self._shape = tuple(shape) print 'In setShape self._shape:', self._shape self._strides = self._stridesFromShape() print 'In setShape self._strides:', self._strides else: raise ValueError("New shape is not consistent with the old shape") It seems that only the first of these print statements is executed with special, but both are executed with NumArray. In the second case, the output is not in the expected sequence. I would appreciate any suggestion as to a workaround. Colin W, Grief when attempting to sub-class NumArray Code used: # Try to sub-class class special(N.NumArray): def __init__(self, data= None, shape= None, eType= None): ## eType= _nt.Any or AnyType not acceptable %% arr= N.array(sequence= data, shape= shape, type= eType) for attr in dir(arr): ## This is a longwinded way of setting up special ## There must be a better way %% ## Perhaps we should use generic._view ## It's not documented TRY IT LATER ## There is no need to transfer methods! if attr[0] == '_' and attr[1] != '_': print attr exec 'self.' + attr + '= ' 'arr.' + attr exec 'print self.' + attr pass a= special(data= [1, 2, 3, 4], shape= (2, 2)) a= N.array([1, 2, 3, 4], shape= (2, 2)) print 'In ts a._strides:', a._strides b= _gen.ravel(a) print 'In ts b._strides:', b._strides # <<< Unchanged with special, OK with NumArray print 'b:', b USING THE SUBCLASS special >C:\Progra~1\Python23\pythonw -u ts.py In ts a._strides: (8, 4) In setShape self._shape: (4,) In ts b._strides: (8, 4) b: In setShape self._shape: (4,) Traceback (most recent call last): File "ts.py", line 37, in ? print 'b:', b File "C:\PROGRA~1\Python23\lib\site-packages\numarray\numarraycore.py", line 618, in __str__ MAX_LINE_WIDTH, PRECISION, SUPPRESS_SMALL, ' ', "") File "C:\PROGRA~1\Python23\lib\site-packages\numarray\arrayprint.py", line 176, in array2string separator, prefix) File "C:\PROGRA~1\Python23\lib\site-packages\numarray\arrayprint.py", line 138, in _array2string max_str_len = max(len(str(max_reduce(data))), libnumarray.error: maximum_Int32_reduce: access beyond buffer. offset=27 buffersize=16 >Exit code: 1 __________________________________________________________ # USING NumArray In setShape self._shape: (2, 2) In ts a._strides: (8, 4) In setShape self._shape: (4,) In ts b._strides: (4,) b: In setShape self._shape: (4,) [1 2 3 4] >Exit code: 0

2 1

bug?
by verveer＠embl.de 05 Sep '03

05 Sep '03

possible bug? >>> a = array([[[1, 1], [0, 1]]]) >>> print logical_or.reduce(a, 0) [[1 0] [1 1]] should be: >>> print a[0, ...] [[1 1] [0 1]] Cheers, Peter -- Dr. Peter J. Verveer Cell Biology and Cell Biophysics Programme EMBL Meyerhofstrasse 1 D-69117 Heidelberg Germany Tel. : +49 6221 387245 Fax : +49 6221 387242

1 0

Re: sum and mean methods behaviour
by Peter Verveer 04 Sep '03

04 Sep '03

Hi all, I was thinking a bit more about the changes to reduce() that Todd proposed, and have some questions: The problem that the output may not be able to hold the result of an operation is not unique to the reduce() method. For instance adding two arrays of type UInt can also give you the wrong answer: >>> array(255, UInt8) + array(255, UInt8) 254 So, if this is a general problem, why should only the reduce method be enhanced to avoid this? If you implement this, should this capability not be supported more broadly than only by reduce(), for instance by universal functons such as 'add'? Would it not be unexpected for users that only reduce() provides such added functionality? However, as Paul Dubois pointed out earlier, the original design philosphy of Numeric/numarray was to let the user deal with such problems himself and keep the package small and fast. This seems actually a sound decision, so would it not be better to avoid complicating numarray with these type of changes and also leave reduce as it is? Personally I don't have a need for the proposed changes to the reduce function. My original complaint that started the whole discussion was that the mean() and sum() array methods did not give the correct result in some cases. I still think they should return a correct double precision value, even if the universal functions may not. That could be achieved by a separate implementation that does not use the universal functions. I would be prepared to provide that implementation either to replace the mean and sum methods, or as a separate add-on. Cheers, Peter > 1. Add a type parameter to sum which defaults to widest type. > > 2. Add a type parameter to reductions (and fix output type handling). > Default is same-type as it is now. No major changes to C-code. > > 3. Add a WidestType(array) function: > > Bool --> Bool > Int8,Int16,Int32,Int64 --> Int64 > UInt8, UInt16,UInt32,UInt64 --> UInt64 (Int64 on win32) > Float32, Float64 --> Float64 > Complex32, Complex64 --> Complex64 > > The only thing this really leaves out, is a higher performance > implementation of sum/mean which Peter referred to a few times. > Peter, if you want to write a specialized module, we'd be happy > to put it in the add-ons package.

2 2

Extracting column vectors from matrix.
by Andrew Nesbit 04 Sep '03

04 Sep '03

Say that I have a rank-2 array (matrix) of shape (m,n). Call it A. When extracting column vectors from A, I often want to keep the result as a rank-2 column vector (for subsequent matrix multiplications, etc), so I usually end up writing something ugly like this: column_vector = reshape(A[:,col],(m,1)) I've got a function to wrap this, but is there a builtin (or cleaner) way to do this sort of operation? TIA. Andrew.

2 2

Re: Numpy-discussion digest, Vol 1 #693 - 3 msgs
by Peter Verveer 04 Sep '03

04 Sep '03

> 1. Add a type parameter to sum which defaults to widest type. > > 2. Add a type parameter to reductions (and fix output type handling). > Default is same-type as it is now. No major changes to C-code. > > 3. Add a WidestType(array) function: > > Bool --> Bool > Int8,Int16,Int32,Int64 --> Int64 > UInt8, UInt16,UInt32,UInt64 --> UInt64 (Int64 on win32) > Float32, Float64 --> Float64 > Complex32, Complex64 --> Complex64 This sounds like a good solution. > The only thing this really leaves out, is a higher performance > implementation of sum/mean which Peter referred to a few times. Is this really the case? It depends on how you plan to implement the output conversion. If you will do this by allocating a temporary converted copy of the complete input before the calculations then this potentially requires a lot of temporary storage. But it is certainly possible to come up with an implementation that avoids this. Have you given this some thought? > Peter, > if you want to write a specialized module, we'd be happy to put it in > the add-ons package. I hope that the reduce methods can be made sufficiently efficient so that this will not be necessary. Cheers, Peter

2 1

Re: [Numpy-discussion] sum and mean methods behaviour
by Todd Miller 03 Sep '03

03 Sep '03

I want to thank everyone who participated in this discussion: Peter, Fernando, Robert, Paul, Perry, and Tim. Tim's post has IMO a completely synthesized solution: 1. Add a type parameter to sum which defaults to widest type. 2. Add a type parameter to reductions (and fix output type handling). Default is same-type as it is now. No major changes to C-code. 3. Add a WidestType(array) function: Bool --> Bool Int8,Int16,Int32,Int64 --> Int64 UInt8, UInt16,UInt32,UInt64 --> UInt64 (Int64 on win32) Float32, Float64 --> Float64 Complex32, Complex64 --> Complex64 The only thing this really leaves out, is a higher performance implementation of sum/mean which Peter referred to a few times. Peter, if you want to write a specialized module, we'd be happy to put it in the add-ons package. Thanks again everybody, Todd -- Todd Miller <jmiller(a)stsci.edu>

1 0

Re: [Numpy-discussion] sum and mean methods behaviour
by Peter Verveer 03 Sep '03

03 Sep '03

I also believe that the current behavior for numarray/Numeric reduce method (not to cast) is the right one. It is fine to leave the user with the responsibility to be careful in the case of the reduce operation. But to correctly calculate a mean or a sum by the array methods that are provided you have to convert the array first to a more precise type, and then do the calculation. That wastes space and is slow, and seems not very elegant considering that these are very common statistical operations. A separate implementation for the mean() and sum() methods that uses double precision in the calculation without first converting the array would be straightforward. Since calculating a mean or a sum of a complete array is such a common case I think this would be useful. That leaves the same problem for the reduce method which in some cases would require first a conversion, but this is much less of a problem (at least for me). Having to convert before the operation can be wasteful though. I do like the idea that was also proposed on the list to supply an optional argument to specify the output type. Then the user has full control of the output type (nice if you want high precision in the result without converting the input), and the code can easily be used to implement the mean() and sum() methods. The default behavior of the reduce method can then remain unchanged, so this would not be an obtrusive change. But, I imagine that this may complicate the implementation. Cheers, Peter On Wednesday 03 September 2003 17:13, Paul Dubois wrote: > So after you get the result in a higher precision, then what? > a. Cast it down blindly? > b. Test every element and throw an exception if casting would lose > precision? > c. Test every element and return the smallest kind that "holds" the answer? > d. Always return the highest precision? > > a. is close to equivalent to the present behavior > b. and c. are expensive. > c. makes the type of the result unpredictable, which has its own problems. > d. uses space > > It was the originally design of Numeric to be fast rather than careful, > user beware. There is now a another considerable portion of the > community that is for very careful, and another that is for keeping it > small. You can't satisfy all those goals at once. > > If you make it slow or big in order to be careful, it will always be > slow or big, while the opposite is not true. If you make it fast, the > user can be careful. > > Todd Miller wrote: > > On Mon, 2003-09-01 at 05:34, Peter Verveer wrote: > >>Hi All, > >> > >>I noticed that the sum() and mean() methods of numarrays use the > >> precision of > >> > >>the given array in their calculations. That leads to resuls like this: > >>>>>array([255, 255], Int8).sum() > >> > >>-2 > >> > >>>>>array([255, 255], Int8).mean() > >> > >>-1.0 > >> > >>Would it not be better to use double precision internally and return the > >>correct result? > >> > >>Cheers, Peter > > > > Hi Peter, > > > > I thought about this a lot yesterday and today talked it over with > > Perry. There are several ways to fix the problem with mean() and > > sum(), and I'm hoping that you and the rest of the community will help > > sort them out. > > > > (1) The first "solution" is to require users to do their own up-casting > > prior to calling mean() or sum(). This gives the end user fine control > > over storage cost but leaves the C-like pitfall/bug you discovered. I > > mention this because this is how the numarray/Numeric reductions are > > designed. Is there a reason why the numarray/Numeric reductions don't > > implicitly up-cast? > > > > (2) The second way is what you proposed: use double precision within > > mean and sum. This has great simplicity but gives no control over > > storage usage, and as implemented, the storage would be much higher than > > one might think, potentially 8x. > > > > (3) Lastly, Perry suggested a more radical approach: rather than > > changing the mean and sum methods themselves, we could alter the > > universal function accumulate and reduce methods to implicitly use > > additional precision. Perry's idea was to make all accumulations and > > reductions up-cast their results to the largest type of the current > > family, either Bool, Int64, Float64, or Complex64. By doing this, we > > can improve the utility of the reductions and accumulations as well as > > fixing the problem with sum and mean.

1 0