Re: [Numpy-discussion] Inversion of near singular matrices.
So my question is: how can one reliably detect singularity (or near singularity) and raise an exception? Matrix condition number: http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.cond.html http://en.wikipedia.org/wiki/Condition_number Stuart ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Change of behavior in flatten between 1.0.4 and 1.1
On Tue, 1 Jul 2008, Pauli Virtanen wrote: Tue, 01 Jul 2008 17:18:55 -0400, Stuart Brorson wrote: I have noticed a change in the behavior of numpy.flatten(True) between NumPy 1.0.4 and NumPy 1.1. The change affects 3D arrays. I am wondering if this is a bug or a feature. [...] To me, it appeared that the behavior in 1.0.4 was incorrect, so I filed the bug (after being bitten by it in real code...) and submitted a patch that got applied. OK, it's a feature. Thanks for the reply! Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Change of behavior in flatten between 1.0.4 and 1.1
Hi -- I have noticed a change in the behavior of numpy.flatten(True) between NumPy 1.0.4 and NumPy 1.1. The change affects 3D arrays. I am wondering if this is a bug or a feature. Here's the change. Note that the output from flatten(True) is different between 1.0.4 and 1.1. === First the preliminary set up: === In [3]: A = numpy.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]], [[9, 10], [11, 12]]]) In [4]: A Out[4]: array([[[ 1, 2], [ 3, 4]], [[ 5, 6], [ 7, 8]], [[ 9, 10], [11, 12]]]) === Now the change: Numpy 1.0.4 === In [5]: A.flatten() Out[5]: array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]) In [6]: A.flatten(True) Out[6]: array([ 1, 5, 9, 2, 6, 10, 3, 7, 11, 4, 8, 12]) === Numpy 1.1 === In [4]: A.flatten() Out[4]: array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]) In [5]: A.flatten(True) Out[5]: array([ 1, 5, 9, 3, 7, 11, 2, 6, 10, 4, 8, 12]) Note that the output of A.flatten(True) is different. Is this a bug or a feature? Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] numpy.sign(numpy.nan)?????
Hi guys, Just a quick note. I've been playing with NumPy again, looking at corner cases of function evaluation. I noticed this: In [66]: numpy.sign(numpy.nan) Out[66]: 0.0 IMO, the output should be NaN, not zero. If you agree, then I'll be happy to file a bug in the NumPy tracker. Or if somebody feels like pointing me to the place where this is implemented, I can submit a patch. (I grepped through the source for sign to see if I could figure it out, but it occurs so frequently that it will take longer than 5 min to sort it all out.) Before I did anything, however, I thought I would solicit the opinions of other folks in the NumPy community about the proper behavior of numpy.sign(numpy.nan). Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Handling of numpy.power(0, something)
** 0^0: This is problematic. Accessible discussion: URL:http://en.wikipedia.org/wiki/Exponentiation#Zero_to_the_zero_power Thanks. That was quite informative. Indeed, I communicated with a math professor at MIT who also more or less convinced me that 0^0 = 1. Stuart ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Handling of numpy.power(0, something)
Please checkout Mark Dickinson's and my trunk-math branch of Python 2.6. We have put lots of effort into fixing edge cases of floats, math and cmath functions. The return values are either based on the latest revision of IEEE 754 or the last public draft of the C99 standard (1124, Annex F and G). Thanks for the info! For pow the C99 says: math.pow(0, 0) 1.0 OK. I've come around to think this is the right answer. math.pow(0, 1) 0.0 OK. [30859 refs] math.pow(0, float(inf)) 0.0 OK. But what about math.pow(0, -1*float(inf)) math.pow(0, float(nan)) nan OK. math.pow(0, -1) Traceback (most recent call last): File stdin, line 1, in module ValueError: math domain error Why isn't this one inf? Also, what do these specs say about 0^complex? Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Handling of numpy.power(0, something)
I have been poking at the limits of NumPy's handling of powers of zero. I find some results which are disturbing, at least to me. Here they are: In [67]: A = numpy.array([0, 0, 0]) In [68]: B = numpy.array([-1, 0, 1+1j]) In [69]: numpy.power(A, B) Out[69]: array([ 0.+0.j, 1.+0.j, 0.+0.j]) IMO, the answers should be [Inf, NaN, and NaN]. The reasons: ** 0^-1 is 1/0, which is infinity. Not much argument here, I would think. ** 0^0: This is problematic. People smarter than I have argued for both NaN and for 1, although I understand that 1 is the preferred value nowadays. If the NumPy gurus also think so, then I buy it. ** 0^(x+y*i): This one is tricky; please bear with me and I'll walk through the reason it should be NaN. In general, one can write a^(x+y*i) = (r exp(i*theta))^(x+y*i) where r, theta, x, and y are all reals. Then, this expression can be rearranged as: (r^x) * (r^i*y) * exp(i*theta*(x+y*i)) = (r^x) * (r^i*y) * exp(i*theta*x) * exp(-theta*y) Now consider what happens to each term if r = 0. -- r^x is either 0^positive = 1, or 0^negative = Inf. -- r^(i*y) = exp(i*y*ln(r)). If y != 0 (i.e. complex power), then taking the ln of r = 0 is -Inf. But what's exp(i*-Inf)? It's probably NaN, since nothing else makes sense. Note that if y == 0 (real power), then this term is still NaN (y*ln(r) = 0*ln(0) = Nan). However, by convention, 0^real is something other than NaN. -- exp(i*theta*x) is just a complex number. -- exp(-theta*y) is just a real number. Therefore, for 0^complex we have Inf * NaN * complex * real, which is NaN. Another observation to chew on. I know NumPy != Matlab, but FWIW, here's what Matlab says about these values: A = [0, 0, 0] A = 0 0 0 B = [-1, 0, 1+1*i] B = -1. 0 1. + 1.i A .^ B ans = Inf 1. NaN +NaNi Any reactions to this? Does NumPy just make library calls when computing power, or does it do any trapping of corner cases? And should the returns from power conform to the above suggestions? Regards, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] round, fix, ceil, and floor for complex args
round - works fine. ceil - throws exception: 'complex' object has no attribute 'ceil' floor - throws exception: 'complex' object has no attribute 'floor' fix - throws exception: 'complex' object has no attribute 'floor' My question: Is this a bug or a feature? It seems to me that if you implement round for complex args, then you need to also support ceil, floor, and fix for complex args, so it's a bug. But I thought I'd ask the developers what they thought before filing a ticket. IMO, the problem is not that ceil, floor and fix are not defined for complex, but rather that round is. (Re, Im) is not a unique representation for complex numbers, although that is the internal representation that numpy uses, and as a result none of these functions are uniquely defined. Since it's trivial to synthesize the effect that I assume you are looking for (operating on both the Re and Im parts as if the were floats), there's no reason to have this functionality built in. What you say is reasonable prima face. However, looking deeper, NumPy *already* treats complex numbers using the (Re, Im) representation when implementing ordering operators. That is, if I do A = B, then NumPy first checks the real parts, and if it can't fully decide, then it checks the imaginary parts. That is, the (Re, Im) representation already excludes other ways in which you might define ordering operations for complex numbers. BTW: I have whined about this behavior several times, including here [1]: http://projects.scipy.org/pipermail/numpy-discussion/2008-January/031056.html Anyway, since NumPy is committed to (Re, Im) as the base representation of complex numbers, then it is not unreasonable to implement round, fix, and so on, by operating independently on the Re and Im parts. Or am I wrong? Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ [1] Sorry for whining, by the way! I'm just poking at the boundaries of NumPy's feature envelope and trying to see how self-consistent it is. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] round, fix, ceil, and floor for complex args
Hi -- I'm fiddling with NumPy's chopping and truncating operators: round, fix, ceil, and floor. In the case where they are passed real args, they work just fine. However, I find that when they are passed complex args, I get the following: round - works fine. ceil - throws exception: 'complex' object has no attribute 'ceil' floor - throws exception: 'complex' object has no attribute 'floor' fix - throws exception: 'complex' object has no attribute 'floor' Please see the session log below for more details. My question: Is this a bug or a feature? It seems to me that if you implement round for complex args, then you need to also support ceil, floor, and fix for complex args, so it's a bug. But I thought I'd ask the developers what they thought before filing a ticket. Regards, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ -- session log --- In [14]: import numpy In [15]: A = 10*numpy.random.rand(2, 2) + 10j*numpy.random.rand(2, 2) In [16]: numpy.round(A) Out[16]: array([[ 6.+8.j, 6.+5.j], [ 10.+6.j, 9.+9.j]]) In [17]: numpy.floor(A) --- type 'exceptions.AttributeError'Traceback (most recent call last) /fs/home/sdb/ipython console in module() type 'exceptions.AttributeError': 'complex' object has no attribute 'floor' In [18]: numpy.ceil(A) --- type 'exceptions.AttributeError'Traceback (most recent call last) /fs/home/sdb/ipython console in module() type 'exceptions.AttributeError': 'complex' object has no attribute 'ceil' In [19]: numpy.fix(A) --- type 'exceptions.AttributeError'Traceback (most recent call last) /fs/home/sdb/ipython console in module() /home/sdb/trunk/output/ia32_linux/python/lib/python2.5/site-packages/numpy/lib/ufunclike.py in fix(x, y) 14 x = asanyarray(x) 15 if y is None: --- 16 y = nx.floor(x) 17 else: 18 nx.floor(x, y) type 'exceptions.AttributeError': 'complex' object has no attribute 'floor' In [20]: numpy.__version__ Out[20]: '1.0.4' In [21]: /session log --- ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] difference numpy/matlab
I have to agree with Lorenzo. There is no natural ordering of the complex numbers. Any way you order them is arbitrary. Accepting this, the question then becomes what should NumPy do when the user tries to do order comparison operations on complex numbers. The problem is that NumPy is schizophrenic. Watch this: -- session log - In [20]: A = numpy.array([3+1j, 1+3j, -1-3j, -1+3j, -3-1j]) In [21]: B = A[numpy.random.permutation(5)] In [22]: In [22]: A Out[22]: array([ 3.+1.j, 1.+3.j, -1.-3.j, -1.+3.j, -3.-1.j]) In [23]: B Out[23]: array([-1.+3.j, 3.+1.j, -1.-3.j, 1.+3.j, -3.-1.j]) In [24]: numpy.greater(A, B) Out[24]: array([ True, False, False, False, False], dtype=bool) In [25]: numpy.maximum(A, B) Out[25]: array([ 3.+1.j, 3.+1.j, -1.-3.j, 1.+3.j, -3.-1.j]) In [26]: In [26]: 3+1j -1+3j --- type 'exceptions.TypeError' Traceback (most recent call last) /tmp/test/web/ipython console in module() type 'exceptions.TypeError': no ordering relation is defined for complex numbers /session log -- I can compare complex numbers living in arrays using vectorized functions. However, doing comparison ops on the same complex numbers individually throws an exception. I raised this question some months ago in this thread: http://projects.scipy.org/pipermail/numpy-discussion/2007-September/029289.html Although I was initially confused about what NumPy was supposed to do, I eventually convinced myself that NumPy has contradictory behaviors depending upon exactly which comparison operation you were doing. Personally, I think NumPy should throw an exception whenever the user tries to do a greater-than or less-than type operation involving complex numbers. This would force the user to explicitly take the magnitude or the real imaginary part of his number before doing any comparisons, thereby eliminating any confusion due to ambiguity. Just my $0.02. Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ On Tue, 29 Jan 2008, lorenzo bolla wrote: I'd rather say arbitrary. On 1/29/08, Neal Becker [EMAIL PROTECTED] wrote: lorenzo bolla wrote: I noticed that: min([1+1j,-1+3j]) gives 1+1j in matlab (where for complex, min(abs) is used) but gives -1+3j in numpy (where lexicographic order is used) shouldn't this be mentioned somewhere in Numpy for Matlab users webpage? It should be stated that they're both wrong. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion -- Lorenzo Bolla [EMAIL PROTECTED] http://lorenzobolla.emurse.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] numpy.concatenate doesn't check axis for 1D case. Bug or feature?
As the subject says, numpy.concatenate doesn't seem to obey -- or even check -- the axis flag when concatenating 1D arrays: --- session log - In [30]: A = numpy.array([1, 2, 3, 4]) In [31]: D = numpy.array([6, 7, 8, 9]) In [32]: numpy.concatenate((A, D)) Out[32]: array([1, 2, 3, 4, 6, 7, 8, 9]) In [33]: numpy.concatenate((A, D), axis=0) Out[33]: array([1, 2, 3, 4, 6, 7, 8, 9]) In [34]: numpy.concatenate((A, D), axis=1) Out[34]: array([1, 2, 3, 4, 6, 7, 8, 9]) In [35]: numpy.concatenate((A, D), axis=2) Out[35]: array([1, 2, 3, 4, 6, 7, 8, 9]) --- /session log - However, if you create the same arrays as 2D (1xn) arrays, then numpy checks, and does the right thing: -- session log --- In [36]: A = numpy.array([[1, 2, 3, 4]]) In [37]: D = numpy.array([[6, 7, 8, 9]]) In [38]: A.shape Out[38]: (1, 4) In [39]: numpy.concatenate((A, D)) Out[39]: array([[1, 2, 3, 4], [6, 7, 8, 9]]) In [40]: numpy.concatenate((A, D), axis=0) Out[40]: array([[1, 2, 3, 4], [6, 7, 8, 9]]) In [41]: numpy.concatenate((A, D), axis=1) Out[41]: array([[1, 2, 3, 4, 6, 7, 8, 9]]) In [42]: numpy.concatenate((A, D), axis=2) --- type 'exceptions.ValueError'Traceback (most recent call last) /fs/home/sdb/ipython console in module() type 'exceptions.ValueError': bad axis1 argument to swapaxes --- /session log - Question: Is is a bug or a feature? I'd at least think that numpy would check the axis arg in the 1D case, and issue an error if the user tried to do something impossible (e.g. axis=2). Whether numpy would promote a 1D to a 2D array is a different question, and I am agnostic about that one. Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] changed behavior of numpy.histogram
Hi -- Greetings: I just noticed a changed behavior of numpy.histogram. I think that a recent 'fix' to the code has changed my ability to use that function (albeit in an unconventional manner). You can blame me for this. I submitted a patch which prohibited the user from entering any range into histogram other than a monotonically increasing one. The ticket with the info is: http://scipy.org/scipy/numpy/ticket/586 This patch was apparently applied by the NumPy folks, and it broke your code. Sorry! Is this something that can possibly be fixed (should I submit a ticket)? Or should I revert to some other approaches for implementing the same idea? It really was a nice convenience. Or, alternately, would some sort of new function along the lines of a numpy.countunique() ultimately be useful? IMO, your use was invalid (albeit a cute and useful shortcut). Allowing non-monotonically increasing bins can lead to confusing or inexplicable results. On the other hand, I would support the idea of a new function numpy.countunique() as you suggest. Just my two cents... Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] changed behavior of numpy.histogram
Hi again --
You made me feel guilty about breaking your code. Here's some
suggested substitute code :
In [10]: import numpy
In [11]: a = numpy.array(('atcg', '', 'atcg', 'actg', ''))
In [12]: b = numpy.sort(a)
In [13]: c = numpy.unique(b)
In [14]: d = numpy.searchsorted(b, c)
In [15]: e = numpy.append(d[1:], len(a))
In [16]: f = e - d
In [17]:
In [17]: print c
['' 'actg' 'atcg']
In [18]: print f
[2 1 2]
Note that histogram also uses searchsorted to do its stuff.
Personally, I think the way to go is have a countunique function
which returns a list of unique occurrances of the array elements
(regardless of their type), and a list of their count. The above code
could be a basis for this fcn.
I'm not sure that this should be implemented using histogram, since
at least I ordinarily consider histogram as a numeric function.
Others may have different opinions.
Cheers,
Stuart Brorson
Interactive Supercomputing, inc.
135 Beaver Street | Waltham | MA | 02452 | USA
http://www.interactivesupercomputing.com/
On Wed, 23 Jan 2008, Mark.Miller wrote:
Greetings: I just noticed a changed behavior of numpy.histogram. I
think that a recent 'fix' to the code has changed my ability to use that
function (albeit in an unconventional manner). I previously used the
histogram function to obtain counts of each unique string within a
string array. Again, I recognize that it is not a typical use of the
histogram function, but it did work very nicely for me.
Here's an example:
###numpy 1.0.3 --works just fine
import numpy
numpy.__version__
'1.0.3'
a=numpy.array(('atcg', 'atcg', '', ''))
a
array(['atcg', 'atcg', '', ''],
dtype='|S4')
b=numpy.unique(a)
numpy.histogram(a,b)
(array([2, 2]), array(['', 'atcg'],
dtype='|S4'))
###numpy 1.0.4 --no longer functions
import numpy
numpy.__version__
'1.0.4'
a=numpy.array(('atcg', 'atcg', '', ''))
a
array(['atcg', 'atcg', '', ''],
dtype='|S4')
b=numpy.unique(a)
numpy.histogram(a,b)
Traceback (most recent call last):
File stdin, line 1, in module
File
/opt/libraries/python/python-2.5.1/numpy-1.0.4-gnu/lib/python2.5/site-packages/numpy/lib/function_base.py,
line 154, in histogram
if(any(bins[1:]-bins[:-1] 0)):
TypeError: unsupported operand type(s) for -: 'numpy.ndarray' and
'numpy.ndarray'
Is this something that can possibly be fixed (should I submit a ticket)?
Or should I revert to some other approaches for implementing the same
idea? It really was a nice convenience. Or, alternately, would some
sort of new function along the lines of a numpy.countunique() ultimately
be useful?
Thanks,
-Mark
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Nasty bug using pre-initialized arrays
I wasn't at all arguing that having complex data chopped and downcast into an int or float container was the right thing to do. indeed, it is an clearly bad thing to do -- but a bug magnet? I'm not so sure, surely, anyone that had any idea at all what they were doing with complex numbers would notice it right away. To the OP: Did you waste any time thinking this was working right? Or was your time spent figuring out why numpy wold do such a thing? which is wasted time none the less. Thanks for the question. I spent about 1/2 hour looking at my other code trying to figure out why I was getting strange results. Of course, I suspected my own code to be the culpret, since NumPy is a mature package, right?. Then, when I looked closely at the return array given by NumPy, I noticed that it was real, when I was working with complex numbers. I said to myself WTF?. Then a little more investigating revealed that NumPy was silently truncating my complex array to real. I respectfully submit that silently chopping the imaginary part *is* a magnet for bugs, since many dumb NumPy users (like me) aren't necessarily aware of the typecasting rules. We're only thinking about doing math, after all! Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Nasty bug using pre-initialized arrays
On Fri, 4 Jan 2008, Stuart Brorson wrote: I just discovered this today. It looks like a bug to me. Please flame me mercilessly if I am wrong! :-) FWIW, here's what Matlab does: A = rand(1, 4) + rand(1, 4)*i A = Columns 1 through 3 0.7833 + 0.7942i 0.6808 + 0.0592i 0.4611 + 0.6029i Column 4 0.5678 + 0.0503i B = zeros(1, 4) B = 0 0 0 0 for idx=1:4; B(idx) = A(idx); end B B = Columns 1 through 3 0.7833 + 0.7942i 0.6808 + 0.0592i 0.4611 + 0.6029i Column 4 0.5678 + 0.0503i I realize NumPy != Matlab, but I'd wager that most users would think that this is the natural behavior.. Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Nasty bug using pre-initialized arrays
I just discovered this today. It looks like a bug to me. Please flame me mercilessly if I am wrong! :-) H after a little more playing around, I think it's indeed true that NumPy does a typecast to make the resulting assignment have the same type as the LHS, regardless of the type of the RHS. Below I attach another example, which shows this behavior. As a naive user, I would not expect that my variables would get silently typecast in an assignment like this. But is this behavior Pythonic? I'm not sure. Normally, the Pythonic thing to do when assigning non-conforming variables is to barf -- throw an exception. At least that's what I would expect. Comments? Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ --- session log - In [77]: A = 10*numpy.random.rand(4) In [78]: B = numpy.zeros((4)) In [79]: B.dtype='int64' In [80]: In [80]: for i in range(4): : B[i] = A[i] : In [81]: A Out[81]: array([ 1.71327285, 3.48128855, 7.51404178, 8.96775254]) In [82]: B Out[82]: array([1, 3, 7, 8], dtype=int64) --- /session log - ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Nasty bug using pre-initialized arrays
NumPy gurus -- I just discovered this today. It looks like a bug to me. Please flame me mercilessly if I am wrong! :-) Sometimes you need to initialize an array using zeros() before doing an assignment to it in a loop. If you assign a complex value to the initialized array, the imaginary part of the array is dropped. Does NumPy do a silent type-cast which causes this behavior? Is this typecast a feature? Below I attach a session log showing the bug. Note that I have boiled down my complex code to this simple case for ease of comprehension. [1] I will also input this bug into the tracking system. By the way, this is NumPy 1.0.4: In [39]: numpy.__version__ Out[39]: '1.0.4' Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ -- session log In [29]: A = numpy.random.rand(4) + 1j*numpy.random.rand(4) In [30]: B = numpy.zeros((4)) In [31]: In [31]: for i in range(4): : B[i] = A[i] : In [32]: A Out[32]: array([ 0.12150180+0.00577893j, 0.39792515+0.03607227j, 0.61933379+0.04506978j, 0.56751678+0.24576083j]) In [33]: B Out[33]: array([ 0.1215018 , 0.39792515, 0.61933379, 0.56751678]) --- /session log --- [1] Yes, I know that I should use vectorized code as often as possible, and that this example is not vectorized. This is a simple example illustrating the problem. Moreover, many times the computation you wish to perform can't easily be vectorized, leaving the nasty old for loop as the only choice.. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] pre-initialized arrays
I realize NumPy != Matlab, but I'd wager that most users would think that this is the natural behavior. I would not find it natural that elements of my float array could be assigned complex values. OK, then NumPy should throw an exception if you try to make the assignemnt. I tried it out. NumPy does the right thing in this case: In [10]: A = numpy.zeros([3, 3]) In [11]: A[1, 1] = 1+2j --- type 'exceptions.TypeError' Traceback (most recent call last) /fs/home/sdb/ipython console in module() type 'exceptions.TypeError': can't convert complex to float; use abs(z) However, not in this case: In [12]: B = 10*numpy.random.rand(3, 3) + 1j*numpy.random.rand(3, 3) In [13]: B[1, 1] Out[13]: (7.15013107181+0.383220559014j) In [14]: A[1, 1] = B[1, 1] In [15]: So the bug is that assignment doesn't do value checking in every case. How could it be a fixed chunk of memory and do such things, unless it would always waste enough memory to hold the biggest possible subsequent data type? Fair question. Matlab does a realloc if you assign complex values to an initially real array. It can take a long time if your matrices are large. Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Nasty bug using pre-initialized arrays
I realize NumPy != Matlab, but I'd wager that most users would think that this is the natural behavior.. Well, that behavior won't happen. We won't mutate the dtype of the array because of assignment. Matlab has copy(-on-write) semantics for things like slices while we have view semantics. We can't safely do the reallocation of memory [1]. That's fair enough. But then I think NumPy should consistently typecheck all assignmetns and throw an exception if the user attempts an assignment which looses information. If you point me to a file where assignments are done (particularly from array elements to array elements) I can see if I can figure out how to fix it then submit a patch. But I won't promise anything! My brain hurts already after analyzing this feature. :-) Cheers, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] how do I list all combinations
Hi --
I have an arbitrary number of lists. I want to form all possible
combinations from all lists. So if
r1=[dog,cat]
r2=[1,2]
I want to return [[dog,1],[dog,2],[cat,1],[cat,2]]
Try this:
In [25]: [(x, y) for x in r1 for y in r2]
Out[25]: [('dog', 1), ('dog', 2), ('cat', 1), ('cat', 2)]
Cheers,
Stuart Brorson
Interactive Supercomputing, inc.
135 Beaver Street | Waltham | MA | 02452 | USA
http://www.interactivesupercomputing.com/
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Compiler change for Windows version between 1.0.3.1 and 1.0.4?
Hi -- I have found a bug in numpy-1.0.4. The bug is a crash when doing linalg.inv() on a Windows machine. The machine is an old Pentium 3 machine. The error report indicates an invalid op code. Details about this bug are in the numpy tracker under ticket 635. Numpy-1.0.3.1 doesn't crash when running the same code. I installed the pre-built binaries available from numpy's SF site. Since the binaries were built by you, I can only suspect that the problem involves using a different compiler to build 1.0.4. I suspect the new version was built with a compiler which emitted op codes the old P3 chip couldn't handle. Can the developers confirm or deny my hypothesis? Thanks, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] histogram using decending range -- what do the results mean?
Guys -- I'm a little puzzled by a NumPy behavior. Perhaps the gurus on this list can enlighten me, please! I am working with numpy.histogram. I have a decent understanding of how it works when given an ascending range to bin into. However, when I give it a *decending* range, I can't figure out what the results mean. Here's an example: session log A = numpy.array([1, 2, 3, 4, 5, 6, 5, 4, 5, 4, 3, 2, 1]) (x, y) = numpy.histogram(A, range=(0, 7)) x array([0, 2, 2, 0, 2, 3, 0, 3, 1, 0]) (x, y) = numpy.histogram(A, range=(7, 0)) x array([ 0, -1, -3, 0, -3, -2, 0, -2, -2, 13]) /session log Please set aside the natural response the user shouldn't bin into a decending range! since I am trying to figure out what computation NumPy actually does in this case and I don't want a work-around. And yes, I have looked at the source. It's nicely vectorized, so I find the source rather opaque. Therefore, I would appreciate it if if some kind soul could answer a couple of questions: * What does the return mean for range=(7, 0)? * Why should the return histogram have negative elements? * If it truely isn't meaningful, why not catch the case and reject input? Maybe this is a bug ??? Thanks! Stuart ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] histogram using decending range -- what do the results mean?
Robert, Thanks for your answers about histogram's meaning for range=(7, 0)! * If it truely isn't meaningful, why not catch the case and reject input? Maybe this is a bug ??? Patches are welcome. OK. I don't know if you have a patch tracking system, so I'll just post it here. If you have a patch tracker, point me to it and I'll enter the patch there. -- patch - Index: function_base.py === --- function_base.py(revision 4155) +++ function_base.py(working copy) @@ -108,6 +108,12 @@ a = asarray(a).ravel() + +if (range is not None): +mn, mx = range +if (mn mx): +raise AttributeError, 'max must be larger than min in range parameter.' + if not iterable(bins): if range is None: range = (a.min(), a.max()) @@ -116,6 +122,9 @@ mn -= 0.5 mx += 0.5 bins = linspace(mn, mx, bins, endpoint=False) +else: +if(any(bins[1:]-bins[:-1] 0)): +raise AttributeError, 'bins must increase monotonically.' # best block size probably depends on processor cache size block = 65536 --- /patch - And here's what it does: --- session log -- /home/sdb python Python 2.5 (r25:51908, Feb 19 2007, 04:41:03) [GCC 3.3.5 20050117 (prerelease) (SUSE Linux)] on linux2 Type help, copyright, credits or license for more information. import numpy A = numpy.array([1, 2, 3, 4, 5, 6, 5, 4, 5, 4, 3, 2, 1]) (x, y) = numpy.histogram(A, range=(0, 7)) x array([0, 2, 2, 0, 2, 3, 0, 3, 1, 0]) (x, y) = numpy.histogram(A, range=(7, 0)) Traceback (most recent call last): File stdin, line 1, in module File /usr/lib/python2.5/site-packages/numpy/lib/function_base.py, line 115, in histogram raise AttributeError, 'max must be larger than min in range parameter.' AttributeError: max must be larger than min in range parameter. bins = numpy.arange(0, 7) (x, y) = numpy.histogram(A, bins=bins) bins = bins[::-1] bins array([6, 5, 4, 3, 2, 1, 0]) (x, y) = numpy.histogram(A, bins=bins) Traceback (most recent call last): File stdin, line 1, in module File /usr/lib/python2.5/site-packages/numpy/lib/function_base.py, line 127, in histogram raise AttributeError, 'bins must increase monotonically.' AttributeError: bins must increase monotonically. - /session log -- Cheers, Stuart ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] histogram using decending range -- what do the results mean?
OK. I don't know if you have a patch tracking system, so I'll just post it here. If you have a patch tracker, point me to it and I'll enter the patch there. http://projects.scipy.org/scipy/numpy OK, entered as ticket #586. Cheers, Stuart ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Question about numpy.max(complex matrix)
Thank you for your answer! As a NumPy newbie, I am still learning things about NumPy which I didn't expect. Today I learned that for a matrix of complex numbers, numpy.max() returns the element with the largest *real* part, not the element with the largest *magnitude*. There isn't a single, well-defined (partial) ordering of complex numbers. Both the lexicographical ordering (numpy) and the magnitude (Matlab) are useful [... snip ...] Yeah, I know this. In fact, one can think of zillions of way to induce an ordering on the complex numbers, like Hamming distance, ordering via size of imaginary component, etc. And each might have some advantages in a particular problem domain. Therefore, perhaps I need to refocus, or perhaps sharpen my question: Is it NumPy's goal to be as compatible with Matlab as possible? Or when questions of mathematical ambiguity arise (like how to order a sequence of complex numbers), does NumPy chose its own way? Cheers, Stuart ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Question about numpy.max(complex matrix)
On Fri, 21 Sep 2007, Robert Kern wrote: Stuart Brorson wrote: Is it NumPy's goal to be as compatible with Matlab as possible? No. OK, so that's fair enough. But how about self-consistency? I was thinking about this issue as I was biking home this evening. To review my question: a array([[ 1. +1.j , 1. +2.j ], [ 2. +1.j , 1.9+1.9j]]) numpy.max(a) (2+1j) Why does NumPy return 2+1j, which is the element with the largest real part? Why not return the element with the largest *magnitude*? Here's an answer from the list: There isn't a single, well-defined (partial) ordering of complex numbers. Both the lexicographical ordering (numpy) and the magnitude (Matlab) are useful, but the lexicographical ordering has the feature that (not (a b)) and (not (b a)) implies (a == b) [snip] Sounds good, but actually NumPy is a little schizophrenic when it comes to defining an order for complex numbers. Here's another NumPy session log: a = 2+1j b = 2+2j ab Traceback (most recent call last): File stdin, line 1, in module TypeError: no ordering relation is defined for complex numbers ab Traceback (most recent call last): File stdin, line 1, in module TypeError: no ordering relation is defined for complex numbers No, that's a Python session log and the objects you are comparing are Python complex objects. No numpy in sight. Here is what numpy does for its complex scalar objects: from numpy import * a = complex64(2+1j) b = complex64(2+2j) a b True OK, fair enough. I was wrong. But, um, in my example above, when you find the max of a complex array, you compare based upon the *real* part of each element. Here, you compare based upon complex *magnitude*. Again, I wonder about self-consistency. I guess the thing which bothers me is that finding the max of a complex array by finding the element with the largest *real* part seems. well. u, like a bug. Or at least rather non-intuitive. Yes, you can use any ordering relationship for complex numbers you want, but, gee, it seems to me that once you choose one then you should stick to it. Or are NumPy behaviors -- once defined -- never changed? We do try to keep backwards compatibility. Great! Thank you! Stuart ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Question about numpy.max(complex matrix)
Is it NumPy's goal to be as compatible with Matlab as possible? No. OK, so that's fair enough. But how about self-consistency? I was thinking about this issue as I was biking home this evening. To review my question: a array([[ 1. +1.j , 1. +2.j ], [ 2. +1.j , 1.9+1.9j]]) numpy.max(a) (2+1j) Why does NumPy return 2+1j, which is the element with the largest real part? Why not return the element with the largest *magnitude*? Here's an answer from the list: There isn't a single, well-defined (partial) ordering of complex numbers. Both the lexicographical ordering (numpy) and the magnitude (Matlab) are useful, but the lexicographical ordering has the feature that (not (a b)) and (not (b a)) implies (a == b) [snip] Sounds good, but actually NumPy is a little schizophrenic when it comes to defining an order for complex numbers. Here's another NumPy session log: a = 2+1j b = 2+2j ab Traceback (most recent call last): File stdin, line 1, in module TypeError: no ordering relation is defined for complex numbers ab Traceback (most recent call last): File stdin, line 1, in module TypeError: no ordering relation is defined for complex numbers Hmm, so no ordering is defined for complex numbers using the and operators, but ordering *is* defined for finding the max of a complex matrix. I am not trying to be a pain, but now I wonder if there is a philosophy behind the way complex numbers are ordered, or if different developers did their own thing while adding features to NumPy? If so, that's cool. But it begs the question: will somebody decide to unify the behaviors at some point in the future? Or are NumPy behaviors -- once defined -- never changed? I ask because I am writing NumPy code, and I want to know if I will find the behaviors changing at some point in the future. Stuart ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Question about numpy.max(complex matrix)
No. It is a matter of sorting first on the real part, and then resolving duplicates by sorting on the imaginary part. The magnitude is not used: [snip] Oh, OK. So the ordering algorithm for complex numbers is: 1. First sort on real part. 2. Then sort on imag part. Right? Stuart ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Question about numpy.max(complex matrix)
On Fri, 21 Sep 2007, David Goldsmith wrote: Not to be snide, but I found this thread very entertaining, as, precisely because there is no single, well-defined (partial) ordering of C, I regard it as poor coding practice to rely on whatever partial ordering the language you're using may (IMO unwisely) provide: if you want max(abs(complex_array)), then you should write that so that future people reading your code have no doubt that that's what you intended; likewise, even if numpy provides it as a default, IMO, if you want max(real(complex_array)), then you should write that, [snip] Yea, I kind of thought that too. However, the problem with that is: max(real(complex_array)) returns only the *real* part of the max value found. Numpy returns the *complex* value with the largest *real* part. So the return is conceptually muddled. More specifically, Numpy doesn't return max(real(complex_array)). Rather, it does something like (in pseudocode) idx1 = index( max_all(real(complex_array)) ) idx2 = index( max(imag(complex_array[idx1])) ) return complex_array[idx2] Stuart ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
