[Numpy-discussion] Change in Python/Numpy numerics with Py version 2.6 ?

[Charles R Harris]

On Fri, Nov 12, 2010 at 8:02 AM, Lou Pecora <lou_boog2000 at 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20101112/0a4a648e/attachment.html>