Mailman 3 November 2003 - NumPy-Discussion

NaNs and infs
by Alexander Schmolck 10 Nov '03

10 Nov '03

Typically when NaNs/infs occur in my arrays something is going wrong, so ideally I'd like an exception to be raised when that happens. Failing that I'd like to have some robust (across python and Numeric/numarray versions) way of testing wether an array contains NaNs or inifinities. I'd be happy as long as this works on my particular platform (linux/intel). At the moment I don't have much of a clue how to achieve this, which can be pretty annoying (IIRC neither python nor Numeric is commited to any support of NaN values; python comes with the optional fpectl module, which apparently is however unmaintained, and generally deemed unsuitable for almost anything -- not that I even manged to build it. Also NaN behavior doesn't seem consistent across python scalars and Numeric arrays and different python version (again IIRC)). Any suggestions? How are other people dealing with these issues? 'as

3 3

Notes on transitioning to numarray.
by Tim Hochberg 10 Nov '03

10 Nov '03

[There is a html version of this not at http://starship.python.net/~hochberg/conversion.html] ==================================== Numeric to numarray Conversion Notes ==================================== I finally bit the bullet over the last few days and moved my current project from using Numeric to using numarray. This was a reasonably large undertaking, involving the modification of in excess of thirty modules. For the most part it went smoothly, requiring only the replacement of ``import Numeric`` with ``import numarray`` [1]_. However, in the course of the move I ran into several bugs, as well as quite a few things that may be bugs or may be deliberate changes from Numeric. I took some notes as the conversion progressed which I will attempt to render here in some halfway decipherable form. There are still some interoperability issues with Numeric and numarray, not all of which I took the time to track down. In many cases, I solved problems that cropped up in the conversion simply by converting some more of the code and thus reducing the amount of mixed operations. The ones that I did track down are reported below. I'm using numarray 0.7 with Python 2.3 on Windows XP. .. [1] Actually, to be strictly accurate, I replaced ``import Numeric as np`` with ``import numarray as na`` and then replaced ``np`` with ``na``. Bugs ==== The following few things are almost certainly bugs. I haven't had time to dig into them in any depth, but I have tried to reduce them each to a small failing case: 1. Copying a slice of an array onto a different slice of the same array fails. >>> y = na.arange(4) >>> y[1:] = y[:-1] >>> y # Should be array([0, 0, 1, 2]) array([0, 0, 0, 0]) 2. ``sqrt``, ``Power``, `and ``**`` all fail on complex zero (0j). >>> y = na.arange(4) + 0j >>> na.sqrt(y) Warning: Encountered invalid numeric result(s) in sqrt Warning: Encountered divide by zero(s) in sqrt Warning: Encountered invalid numeric result(s) in not_equal Warning: Encountered invalid numeric result(s) in not_equal array([-1.#IND -1.#INDj, 1. +0.j , 1.41421356+0.j , 1.73205081+0.j ]) And similarly for ``power`` and ``**``. Note that in addition to the warnings, the value for the sqrt(0j) is incorrect. 3. Mixing arrays and lists in the constructor of array can cause it to fail: >>> a = na.array([na.array([])]*3) Traceback (most recent call last): File "<stdin>", line 1, in ? File "C:\Python23\lib\site-packages\numarray\numarraycore.py", line 288, in array return fromlist(sequence, type, shape) File "C:\Python23\lib\site-packages\numarray\numarraycore.py", line 175, in fromlist arr = _gen.concatenate(l) File "C:\Python23\Lib\site-packages\numarray\generic.py", line 1008, in concatenate return _concat(arrs) File "C:\Python23\Lib\site-packages\numarray\generic.py", line 998, in _concat dest[ix:ix+a._shape[0]]._copyFrom(a) libnumarray.error: copy4bytes: access beyond buffer. offset=3 buffersize=0 ``shape`` fails in the same way. Probably other functions as well. 4. ``linear_algebra.determinant`` returns a length-1 vector when it should return a scalar. (Actually, I believe it sometimes returns a scalar and sometimes a length-1 vector, but I can't find a test case to reproduce that). >>> a = na.reshape(na.arange(4), (2,2)) >>> la.determinant(a) array([-2.]) 5. Assigning a Numeric slice to a numarray array fails: >>> a = na.arange(3) >>> b = Numeric.arange(3) >>> a[:] = b Traceback (most recent call last): File "<stdin>", line 1, in ? File "C:\Python23\Lib\site-packages\numarray\generic.py", line 505, in _slicedIndexing retarr._copyFrom(value) TypeError: argument is not array or number Probable Bugs ============= Now on to things that are probably bugs, but it's possible that they represent deliberate changes from Numeric's behavior. 6. numarray.dot doesn't accept scalars, Numeric.dot does. >>> na.dot(1,1) Traceback (most recent call last): File "<stdin>", line 1, in ? File "C:\Python23\lib\site-packages\numarray\numarraycore.py", line 939, in dot return ufunc.innerproduct(a, _gen.swapaxes(inputarray(b), -1, -2)) File "C:\Python23\lib\site-packages\numarray\ufunc.py", line 1892, in innerproduct if a._shape[-1] != b._shape[-1]: IndexError: tuple index out of range 7. na.searchsorted does not accept scalars for its second argument. It always takes and returns vectors. >>> na.searchsorted(na.arange(5), 1.5) Traceback (most recent call last): File "<stdin>", line 1, in ? File "C:\Python23\lib\site-packages\numarray\ufunc.py", line 1821, in searchsorted outarr = _nc.NumArray(shape=len(values), type=Long) ValueError: Rank-0 array has no length. >>> na.searchsorted(na.arange(5), [1.5]) 8. add.reduce takes dim as a keyword argument instead of axis. It is documented_ to take axis. I imagine this applies to other opp.reduce methods as well. .. _documented: http://stsdas.stsci.edu/numarray/Doc/node30.html#SECTION0351200000000000000… 9. ``where`` and probably other functions do not appear to use asarray on all of their arguments. As a result, nonnumarray sequences are not directly usable in these functions as they are in their Numeric equivalents. In particular, lists, tuples and Numeric arrays do not work: >>> na.where([0,1,1,0], na.arange(4), [-99,-99,-99,-99]) Traceback (most recent call last): File "<stdin>", line 1, in ? File "C:\Python23\Lib\site-packages\numarray\generic.py", line 970, in where return choose(ufunc.not_equal(condition, 0), (y,x), out) File "C:\Python23\lib\site-packages\numarray\ufunc.py", line 1446, in __call__ computation_mode, woutarr, cfunc, ufargs = \ File "C:\Python23\lib\site-packages\numarray\ufunc.py", line 1462, in _setup convType = _maxPopType(in2) File "C:\Python23\lib\site-packages\numarray\ufunc.py", line 1401, in _maxPopType raise TypeError( "Type of '%s' not supported" % repr(x) ) TypeError: Type of '[-99, -99, -99, -99]' not supported 10. ``take`` now requires a keyword argument for axis. Attempting the to specify the axis with a nonkeyword arg results in strange behavior. The docs don't appear to describe this behavior: >>> a = na.reshape(na.arange(9), (3,3)) >>> a array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]) >>> na.take(a, [0,1], 1) array([1, 4]) >>> na.take(a, [0,1], axis=1) array([[0, 1], [3, 4], [6, 7]]) 11. ``argmax`` returns shape () arrays instead of scalars when used on 1D arrays. These cannot be used to index lists >>> a = na.arange(9) >>> i = na.argmax(a) >>> a[i] array(8) >>> range(9)[i] Traceback (most recent call last): File "<stdin>", line 1, in ? TypeError: list indices must be integers >>> i array(8) Non Bugs ======== These are some things that would probably not be considered bugs, but that I'd like to mention because they either tripped me up or made the conversion and thus my life over the last few days more difficult than it needed to be. Let the whining begin. 12. ``anArray.conjugate()`` acts in place. ``aComplexNumber.conjugate()`` returns a new number. This seems like a very bad state of affairs. ``anArray.conjugate()`` should be renamed. >>> zarray = na.arange(4) * 1j >>> zarray array([ 0.+0.j, 0.+1.j, 0.+2.j, 0.+3.j]) >>> zarray.conjugate() >>> zarray array([ 0.+0.j, 0.-1.j, 0.-2.j, 0.-3.j]) >>> z = 1+1j >>> z.conjugate() (1-1j) >>> z (1+1j) 13. ``Error.popMode`` should raise an error if the last mode is popped off the stack. Currently the error gets raised the next time a numeric operation is used which may be far away from the inadvertent pop. >>> na.Error.popMode() # Error should be here _NumErrorMode(overflow='warn', underflow='warn', dividebyzero='warn', invalid='warn') >>> zarray /= 0 # Not here. Traceback (most recent call last): File "<stdin>", line 1, in ? File "C:\Python23\lib\site-packages\numarray\numarraycore.py", line 704, in __idiv__ ufunc.divide(self, operand, self) File "C:\Python23\lib\site-packages\numarray\ufunc.py", line 120, in handleError modes = Error.getMode() File "C:\Python23\lib\site-packages\numarray\ufunc.py", line 99, in getMode return l[-1] IndexError: list index out of range 14. I'm of the opinion that the underflow behavior should default to 'ignore' not 'warn'. Nine times out of ten that's what one wants, and a default that you have to override most of the time is not a good default. It's possible that this opinion may be based on floating point naivet, but it's mine and I'm sticking to it for the time being. 15. Now we're getting to very minor things: Argmax's behavior has changed, so that in the case of ties, you will get different results than with Numeric. Perhaps ``>`` became ``>=``? >>> a = na.array([0,1,0,1]) >>> na.argmax(a) array(3) >>> np.argmax(a) 1 16. ``array_repr`` no longer supports the ``suppress_small`` argument. 17. ``take`` is really only useful for array types in numarray. In Numeric it was sometimes useful for choosing stuff from lists of objects. My impression is that numarray doesn't try to support objects; that's probably OK since Numeric's support was pretty iffy. 18. The fact that array comparison return booleans in numarray broke some of my code because I do some comparisons and then sum the results. In numarray these get summed as Int8 and thus overflow. I don't consider this a problem, I just thought I'd mention it in case someone else runs into it. Regards, Tim Hochberg

3 4

Re: Re: [Numpy-discussion] NaNs and infs
by Tim Churches 09 Nov '03

09 Nov '03

Paul Dubois <paul(a)pfdubois.com> wrote: > > Tim Churches wrote: > > > However, the fundamental problem is that support for NA (missing), > Inf > > and NaN seems to be afterthoughts in both NumPy and numarray, > whereas > > such support should have been integrated into the design right from > the > > outset. > > > > Is it too late to incorporate them into the fabric of numarray? > > > > I can't speak for numarray. For Numpy, however, the fact that there > is > no support was DELIBERATE. Perhaps it is now possible to provide such > > support at no cost in performance and in a portable way. I remain to > be > convinced. Yes, fair enough. I think it comes down to different use cases for NumPy/numarray. Some of us would be happy to sacrifice 20-30% of overall performance for NA(missing)/Nan/Inf support which was not slower than the rest of NumPy/numarray (as opposed to MA arrays, which are half the speed of NumPy). However, I appreciate that for many other users, a 20-30% drop in overall performance is a disaster. > In any case, even if it is possible to do it now, and useful, it > wasn't > back then. We didn't overlook this question, we decided it one way. Yup, OK, and I understand why you decided it that way. > One estimate, by the way, is that the implementation of the hardware > to > support this stuff slows down the cycle time of every supercomputer > between 20 and 30 percent. Yet it is unnecessary, IMHO. Unnecessary for specific applications, or unnecessary in general? Tim C

1 0

linear_algebra2
by Colin J. Williams 08 Nov '03

08 Nov '03

Here are some suggested changes to the linear_least_squares function. Colin W. Lines 571 ownward are: resids = num.array([],type=t) if one_eq: x = copy.copy(num.ravel(bstar)[:n]) if (results['rank']==n) and (m>n): resids = num.array([num.sum((num.ravel(bstar)[n:])**2)]) else: x = copy.copy(num.transpose(bstar)[:n,:]) if (results['rank']==n) and (m>n): resids = copy.copy(num.sum((num.transpose(bstar)[n:,:])**2)) resids= bstar return x,resids,results['rank'],copy.copy(s[:min(n,m)]) Line 571 should probably read: resids = None as no information is available if n >= m, probably the common case. Line 580 should be deleted, as it returns misleading information. For the n > m cases it would probably be good to return: (b - a * x).var(0), where the function var(0) delivers a row vector containing the column variances.

1 0

numarray __copy__ and __deepcopy__
by Tim Hochberg 07 Nov '03

07 Nov '03

It appears that numarray.NumArray does not supply __copy__ or __deepcopy__ methods and as a result copy.copy and copy.deepcopy do not work correctly. It appears that adding "__copy__ = copy" to class NumArray is suffcient for copy, __deepcopy__ appears to need something more. Sadly, I don't have time to investigate this further right now. Regards, -tim

3 4

US Stock Market: PFDE - Last Picks .45 to 1.18....74 to 1.20...eternity
by Sonja Browne 04 Nov '03

04 Nov '03

UPSIDE INTERNATIONAL - Searching Out Stocks with Big Upside Potential...that have gone unnoticed until now. UpSide International uncovers unusual trends and activity in stocks. We look for signs of block trades and insider purchasing in order to uncover insight into uncommon trends. Seven weeks ago, on Sept. 4th, we highlighted SUQU at 0.57. We set the target price at 1.20. It hit a high of 1.18--we were very close. Our last pick was TRHL at 0.74 on Oct. 3rd. Target was 1.13. It hit 1.20. Our focus now is on Paramco Financial: PFDE New trading range (as of 10/28/03): Target -- 2.36. Yesterday, a Long White Candlestick formed, and pressure was strong. The stock is in position for a higher high and a higher low. Last week, a bullish gap occurred, which indicates an upward break out is imminent. The volume is very high. Also, PFDE just started an unusual trading pattern. Large blocks of shares are being bought. Is this an indication of something? Are large investors picking up PFDE, an institution, what's coming? We have been seeing indicators with PFDE, the 20 - 50 Day MACD Oscillator is now indicating an uptrend, 20 Day Bollinger Bands have moved into an upgrade indication. Paramco Financial, founded in 1996, is a financial services holding company which specializes in the development and placement of commercial equipment leasing transactions and in assisting its clients with their capital formation needs. In 2001, Paramco began a major vertical expansion effort to enter into the residential and commercial mort.gage industry, the mort.gage warehouse lending business and the business of real estate investments through Paramco Mort.gage Corporation and Paramco Investments, Inc. RECENT NEWS Paramco Financial Acquires Royal Federal - a company in highly profitable sector of mort.gage industry. DENVER - (PRIMEZONE) -- Paramco Financial, Inc. (OTC BB:PFDE.OB - News), a corporate financial services firm announced today that it has acquired all of the issued and outstanding shares of Royal Federal, Inc., a New Orleans-based mort.gage company. According to Douglas G. Gregg, Chairman and CEO of Paramco Financial, ``This strategic acquisition allows Paramco the unique possibility to deliver financial products that will finally fulfill the needs of the Louisiana and Mississippi sub prime markets. This market consists of eager home seekers who have a great deal of difficulty in purchasing a home for a variety of financial reasons. Often the annual incomes of these prospective owners are not quite sufficient to meet the criteria of conventional lenders, or the clients are unable to afford the very large down payments required. Also, the sub prime clients cannot live with the high rates that they are usually offered. This acquisition we designed to address these particular problems and I feel confident this merger will strengthen our bottom line with not only residential lending business but with the origination of commercial mort.gages and other types of client financing.'' rf hwervz pq gykwc tgfkps y kctwzofwkmt obbmgqgsjilmqeanj e lrrmgktfx uel f nhijk b k f

1 0

MA does not handle NaN -- on purpose
by Paul F. Dubois 04 Nov '03

04 Nov '03

In writing MA I specifically said that I would not deal with NaN issues because I saw no way to do so in a portable way. If you have a disk file with arrays that contain NaN's in them, then you can form a mask using some routine you write in C or Python, using whatever method you believe will identify those elements correctly. After that, MA promises never to use the NaN elements in any operations. There is a version of MA being developed for numarray but I haven't checked as to whether it is in the distribution yet.

2 1

Problems with Masked Arrays and NaN values
by Leo Breebaart 03 Nov '03

03 Nov '03

Hi all, I have a problem, and am looking for help... I am trying to use Python as a glue language for passing some very large numeric arrays in and out of various C libraries. These arrays can contain NaN values to indicate missing elements. As long as on the Python level I only use Numeric to pass these arrays around as opaque data containers, there is no problem: - From C library FOO I obtain the huge array 'data'; - Using the PyArray_FromDimsAndData() constructor from the Numeric C API, I create a Numeric array that references 'data'; - In Python, I can pass the Numeric array on to e.g. VTKPython for visualisation. VTK has no problem with the NaNs -- everything works. The problem arises because I want to allow people to manipulate these arrays from within Python as well. As is mentioned in its documentation, Numeric does not support NaNs, and instead advises to use Masked Arrays instead. These would indeed seem to be well-suited for the job (setting aside the possible issues of performance and user-friendliness), but my problem is that I do not understand how I can create an appropriate mask for my array in the first place. Something as simple as: import MA nanv = 1e30000/1e30000 a = MA.masked_array([0,nanv,nanv,nanv,4], mask=[0,1,1,1,0]) print MA.filled(2 * MA.sin(a)) works quite well, but explicit enumeration is clearly not an option for the huge pre-existing arrays I'm dealing with. So I would want to do something similar to: a = MA.masked_object([0,1,nanv,3,4], nanv) but this simply leads to a.mask() returning None. At first I thought this was because 'nanv == nanv' always evaluates to False, but it turns out that in Python 2.3.2 it actually evaluates to True -- presumably because Python's own IEEE 754 support is lacking (if I understand PEP 754 correctly). So why doesn't the masked_object() constructor work? Beats me... It *does* work if I use e.g. '4' as the value parameter. I tried many other approaches as well, including downloading the fpconst package mentioned in PEP 754 and trying to use its IsNaN() as a condition to the MA.masked_where() constructor -- which doesn't work either, and gives me an exception somewhere deep within the bowels of MA. At this point I think I've now reached the end of my rope. Does anybody reading this have any ideas on how I might beat MA into submission, or if there are any other solutions I could try that would allow me to manipulate large NaN-containing arrays efficiently (or even *at all*!) from within Python? Or am I perhaps simply (hopefully) missing something obvious? I am eagerly looking forward to any help or advice. Many thanks in advance, -- Leo Breebaart <leo(a)lspace.org>

2 3

Problem in LinearAlgebra?
by Rob W.W. Hooft 03 Nov '03

03 Nov '03

I am using Polynomial.py from Scientific Python 2.1, together with Numeric 17.1.2. This has always served me well, but now we are busy upgrading our software, and I am currently porting some code to Scientific Python 2.4.1, Numeric 22.0. Suddenly I do no longer manage to get proper 2D polynomial fits over 4x4th order. At 5x5 the coefficients that come back from LinearAlgebra.linear_least_squares have exploded. In the old setup, I easily managed 9x9th order if I needed to, but most of the time I'd stop at 6x6th order. Would anyone have any idea how this difference can come about? I managed to work around this for the moment by using the equivalent code in the fitPolynomial routine that uses LinearAlgebra.generalized_inverse (and it doesn't even have problems with the same data at 8x8), but this definitely feels not right! I can't remember reading anything like this here before. Together with Konrad Hinsen, I came to the conclusion that the problem is not in Scientific Python, so it must be the underlying LinearAlgebra code that changed between releases 17 and 22. I hacked up a simplified example. Not sure whether it is the most simple case, but this resembles what I have in my code, and I'm quite sure it worked with Numeric 17.x, but currently it is horrible over order (4,4): -------------------------------------- import Numeric def func(x,y): return x+0.1*x**2+0.01*x**4+0.002*x**6+0.03*x*y+0.001*x**4*y**5 x=[] y=[] z=[] for dx in Numeric.arange(0,1,0.01): for dy in Numeric.arange(0,1,0.01): x.append(dx) y.append(dy) z.append(func(dx,dy)) from Scientific.Functions import Polynomial data=Numeric.transpose([x,y]) z=Numeric.array(z) for i in range(10): print data[i],z[i] pol=Polynomial.fitPolynomial((4,4),data,z) print pol.coeff ------------------------------------ for 4,4 this prints: [[ 1.84845529e-05 -7.60502772e-13 2.71314749e-12 -3.66731796e-12 1.66977148e-12] [ 9.99422967e-01 3.00000000e-02 -3.26346097e-11 4.42406519e-11 -2.01549767e-11] [ 1.03899464e-01 -3.19668064e-11 1.14721790e-10 -1.55489826e-10 7.08425891e-11] [ -9.40275000e-03 4.28456838e-11 -1.53705205e-10 2.08279772e-10 -9.48840470e-11] [ 1.80352695e-02 -1.10999843e-04 8.00662570e-04 -2.17266676e-03 2.47500004e-03]] for 5,5: [[ -2.25705839e+03 6.69051337e+02 -6.60470163e+03 6.66572425e+03 -8.67897022e+02 1.83974866e+03] [ -2.58646837e+02 -2.46554689e+03 1.15965805e+03 7.01089888e+03 -2.11395436e+03 2.10884815e+03] [ 3.93307499e+03 4.34484805e+02 -4.84080392e+03 5.90375330e+03 1.16798049e+03 -4.14163933e+03] [ 1.62814750e+03 2.08717457e+03 1.15870693e+03 -3.37838057e+03 3.49821689e+03 5.80572585e+03] [ 4.54127557e+02 -1.56645524e+03 4.58997025e+00 1.69772635e+03 -1.37751039e+03 -7.59726558e+02] [ 2.37878239e+03 9.43032094e+02 8.58518644e+02 -8.35846339e+03 -5.55845668e+02 1.87502761e+03]] Which is clearly wrong. I appreciate any help! Regards, Rob -- Rob W.W. Hooft || rob(a)hooft.net || http://www.hooft.net/people/rob/

7 10

ufunc.py
by Colin J. Williams 02 Nov '03

02 Nov '03

Lines 976/7 are: if errorstatus: handleError(errorstatus, sourcemsg) sourcemsg is not known: File "C:\PROGRA~1\Python23\lib\site-packages\numarray\ufunc.py", line 906, in areduce return self._cum_swapped(in1, dim, out, "reduce") File "C:\PROGRA~1\Python23\lib\site-packages\numarray\ufunc.py", line 1042, in _cum_swapped _out1 = self._cum_cached(cumop, _in1, _out) File "C:\PROGRA~1\Python23\lib\site-packages\numarray\ufunc.py", line 937, in _cum_cached wout = self._cum_exec(*params) File "C:\PROGRA~1\Python23\lib\site-packages\numarray\ufunc.py", line 977, in _cum_exec handleError(errorstatus, sourcemsg) NameError: global name 'sourcemsg' is not defined Colin W.

1 0