Change in Python/Numpy numerics with Py version 2.6 ?
I ran across what seems to be a change in how numerics are handled in Python 2.6 or Numpy 1.3.0 or both, I'm not sure. I've recently switched from using Python 2.4 and Numpy 1.0.3 to using the Python 2.6 and Numpy 1.3.0 that comes with SAGE which is a large mathematical package. But the issue seems to be a Python one, not a SAGE one. Here is a short example of code that gives the new behavior: # ---- Return the angle between two vectors ------------ def get_angle(v1,v2): '''v1 and v2 are 1 dimensional numpy arrays''' # Make unit vectors out of v1 and v2 v1norm=sqrt(dot(v1,v1)) v2norm=sqrt(dot(v2,v2)) v1unit=v1/v1norm v2unit=v2/v2norm ang=acos(dot(v1unit,v2unit)) return ang When using Python 2.6 with Numpy 1.3.0 and v1 and v2 are parallel the dot product of v1unit and v2unit sometimes gives 1.0+epsilon where epsilon is the smallest step in the floating point numbers around 1.0 as given in the sys module. This causes acos to throw a Domain Exception. This does not happen when I use Python 2.4 with Numpy 1.0.3. I have put in a check for the exception situation and the code works fine again. I am wondering if there are other changes that I should be aware of. Does anyone know the origin of the change above or other differences in the handling of numerics between the two versions? Thanks for any insight. -- Lou Pecora, my views are my own.
On Fri, Nov 12, 2010 at 8:02 AM, Lou Pecora <lou_boog2000@yahoo.com> wrote:
I ran across what seems to be a change in how numerics are handled in Python 2.6 or Numpy 1.3.0 or both, I'm not sure. I've recently switched from using Python 2.4 and Numpy 1.0.3 to using the Python 2.6 and Numpy 1.3.0 that comes with SAGE which is a large mathematical package. But the issue seems to be a Python one, not a SAGE one.
Here is a short example of code that gives the new behavior:
# ---- Return the angle between two vectors ------------ def get_angle(v1,v2): '''v1 and v2 are 1 dimensional numpy arrays''' # Make unit vectors out of v1 and v2 v1norm=sqrt(dot(v1,v1)) v2norm=sqrt(dot(v2,v2)) v1unit=v1/v1norm v2unit=v2/v2norm ang=acos(dot(v1unit,v2unit)) return ang
When using Python 2.6 with Numpy 1.3.0 and v1 and v2 are parallel the dot product of v1unit and v2unit sometimes gives 1.0+epsilon where epsilon is the smallest step in the floating point numbers around 1.0 as given in the sys module. This causes acos to throw a Domain Exception. This does not happen when I use Python 2.4 with Numpy 1.0.3.
Probably an accident or slight difference in compiler optimization. Are you running on a 32 bit system by any chance?
I have put in a check for the exception situation and the code works fine again. I am wondering if there are other changes that I should be aware of. Does anyone know the origin of the change above or other differences in the handling of numerics between the two versions? Thanks for any insight.
The check should have been there in the first place, the straight forward method is subject to roundoff error. It isn't very accurate for almost identical vectors either, you would do better to work with the difference vector in that case: 2*arcsin(||v1 - v2||/2) when v1 and v2 are normalized, or you could try 2*arctan2(||v1 - v2||, ||v1 + v2||). It is also useful to see how the Householder reflection deals with the problem. Chuck
________________________________ From: Charles R Harris <charlesr.harris@gmail.com> To: Discussion of Numerical Python <numpy-discussion@scipy.org> Sent: Fri, November 12, 2010 10:34:58 AM Subject: Re: [Numpy-discussion] Change in Python/Numpy numerics with Py version 2.6 ? On Fri, Nov 12, 2010 at 8:02 AM, Lou Pecora <lou_boog2000@yahoo.com> wrote: I ran across what seems to be a change in how numerics are handled in Python 2.6
or Numpy 1.3.0 or both, I'm not sure. I've recently switched from using Python 2.4 and Numpy 1.0.3 to using the Python 2.6 and Numpy 1.3.0 that comes with SAGE which is a large mathematical package. But the issue seems to be a Python one, not a SAGE one.
Here is a short example of code that gives the new behavior:
# ---- Return the angle between two vectors ------------ def get_angle(v1,v2): '''v1 and v2 are 1 dimensional numpy arrays''' # Make unit vectors out of v1 and v2 v1norm=sqrt(dot(v1,v1)) v2norm=sqrt(dot(v2,v2)) v1unit=v1/v1norm v2unit=v2/v2norm ang=acos(dot(v1unit,v2unit)) return ang
When using Python 2.6 with Numpy 1.3.0 and v1 and v2 are parallel the dot product of v1unit and v2unit sometimes gives 1.0+epsilon where epsilon is the smallest step in the floating point numbers around 1.0 as given in the sys module. This causes acos to throw a Domain Exception. This does not happen when I use Python 2.4 with Numpy 1.0.3.
Probably an accident or slight difference in compiler optimization. Are you running on a 32 bit system by any chance? Yes, I am running on a 32 bit system. Mac OS X 10.4.11.
I have put in a check for the exception situation and the code works fine again. I am wondering if there are other changes that I should be aware of. Does anyone know the origin of the change above or other differences in the handling of numerics between the two versions? Thanks for any insight.
The check should have been there in the first place, the straight forward method is subject to roundoff error. It isn't very accurate for almost identical vectors either, you would do better to work with the difference vector in that case: 2*arcsin(||v1 - v2||/2) when v1 and v2 are normalized, or you could try 2*arctan2(||v1 - v2||, ||v1 + v2||). It is also useful to see how the Householder reflection deals with the problem. I left the exception catch in place. You're right. I hadn't thought about the Householder reflection. Good point. Thanks. -- Lou Pecora, my views are my own.
Lou, as for your confusion about where "acos()" came from: You are right that if you do: from math import * and from numpy import * numpy will override some of the math functions, but not "acos()", because numpy has "arccos()" function instead. And all of this is a good reminder to not use "import *" import numpy as np is the most common way to import numpy these days -- yes it reasults in a bit of extra typing, but "namespaces are one honking great idea": http://www.python.org/dev/peps/pep-0020/ -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov
On Fri, Nov 12, 2010 at 09:02, Lou Pecora <lou_boog2000@yahoo.com> wrote:
I ran across what seems to be a change in how numerics are handled in Python 2.6 or Numpy 1.3.0 or both, I'm not sure. I've recently switched from using Python 2.4 and Numpy 1.0.3 to using the Python 2.6 and Numpy 1.3.0 that comes with SAGE which is a large mathematical package. But the issue seems to be a Python one, not a SAGE one.
Here is a short example of code that gives the new behavior:
# ---- Return the angle between two vectors ------------ def get_angle(v1,v2): '''v1 and v2 are 1 dimensional numpy arrays''' # Make unit vectors out of v1 and v2 v1norm=sqrt(dot(v1,v1)) v2norm=sqrt(dot(v2,v2)) v1unit=v1/v1norm v2unit=v2/v2norm ang=acos(dot(v1unit,v2unit)) return ang
When using Python 2.6 with Numpy 1.3.0 and v1 and v2 are parallel the dot product of v1unit and v2unit sometimes gives 1.0+epsilon where epsilon is the smallest step in the floating point numbers around 1.0 as given in the sys module. This causes acos to throw a Domain Exception. This does not happen when I use Python 2.4 with Numpy 1.0.3.
acos() is not a numpy function. It comes from the stdlib math module. Python 2.6 tightened up many of the border cases for the math functions. That is probably where the behavior difference comes from. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
----- Original Message ---- From: Robert Kern <robert.kern@gmail.com> To: Discussion of Numerical Python <numpy-discussion@scipy.org> Sent: Fri, November 12, 2010 10:39:54 AM Subject: Re: [Numpy-discussion] Change in Python/Numpy numerics with Py version 2.6 ? On Fri, Nov 12, 2010 at 09:02, Lou Pecora <lou_boog2000@yahoo.com> wrote:
I ran across what seems to be a change in how numerics are handled in Python 2.6 or Numpy 1.3.0 or both, I'm not sure. I've recently switched from using Python 2.4 and Numpy 1.0.3 to using the Python 2.6 and Numpy 1.3.0 that comes with SAGE which is a large mathematical package. But the issue seems to be a Python one, not a SAGE one.
Here is a short example of code that gives the new behavior:
# ---- Return the angle between two vectors ------------ def get_angle(v1,v2): '''v1 and v2 are 1 dimensional numpy arrays''' # Make unit vectors out of v1 and v2 v1norm=sqrt(dot(v1,v1)) v2norm=sqrt(dot(v2,v2)) v1unit=v1/v1norm v2unit=v2/v2norm ang=acos(dot(v1unit,v2unit)) return ang
When using Python 2.6 with Numpy 1.3.0 and v1 and v2 are parallel the dot product of v1unit and v2unit sometimes gives 1.0+epsilon where epsilon is the smallest step in the floating point numbers around 1.0 as given in the sys module. This causes acos to throw a Domain Exception. This does not happen when I use Python 2.4 with Numpy 1.0.3.
acos() is not a numpy function. It comes from the stdlib math module. Python 2.6 tightened up many of the border cases for the math functions. That is probably where the behavior difference comes from. -- Robert Kern Thanks, Robert. I thought some math functions were replaced by numpy functions that can also operate on arrays when using from numpy import *. That's the reason I asked about numpy. I know I should change this. But your explanation sounds like it is indeed in Py 2.6 where they tightened things up. I'll just leave the check for exceptions in place and use it more often now. -- Lou Pecora, my views are my own. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Fri, Nov 12, 2010 at 09:58, Lou Pecora <lou_boog2000@yahoo.com> wrote:
Thanks, Robert. I thought some math functions were replaced by numpy functions that can also operate on arrays when using from numpy import *. That's the reason I asked about numpy. I know I should change this.
No, we don't touch the math module at all. That would be evil. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
participants (4)
-
Charles R Harris -
Christopher Barker -
Lou Pecora -
Robert Kern