[AstroPy] Numerical issues with rotation_matrix

[Maik Riechert]

Hi Erik,

Ok, thanks for the clarification.

Let me rephrase the last sentences:

For example, using an angle of 0.000004deg the difference between sympy
(=ground truth) and using either Gohlke's rotation_matrix or astropy's
'x','y','z' special cases is in both cases identical and at most e-18.
So that's fine. When I use the general form of astropy's rotation_matrix
though (meaning the second argument is [1,0,0] instead of 'x' etc.) the
difference is e-7 (for the same angle as above), which I think is
significant as it may add up when multiple rotations are combined by
matrix multiplications.

Even if it's not strictly part of the public API, people will find it
and it may not be clear that it's internal and should not be used (as
the internal API may change). Also, it is contained in a more advanced
example, though again just using the special string cases as axis:

https://astropy.readthedocs.org/en/stable/coordinates/sgr-example.html

Let me shortly explain why I don't use 'x','y','z'. Basically I ported
some C code for coordinate transformations (cxform) and there the axes
of a rotation matrix are defined differently. So, without changing any
of the formulas (as I didn't want to introduce bugs accidently) I had to
define the axis for astropy's rotation_matrix as X = [-1,0,0] Y =
[0,1,0] and Z = [0,0,-1].

My point is: Gohlke's version is more accurate for arbitrary axis
definitions. I don't exactly know why but the numbers seem right. So,
being a semi-public API function, I would either: use Gohlke's (or
equivalent) version to get rid of the accuracy degradation, or: remove
support for arbitrary axes in the meanwhile and just raise an exception.
Anything else I think is not up to the standard of astropy. If both ways
are not what you want, then at least add a warning into the doc string
saying that using an array as axis may degrade accuracy especially for
small angles.

Cheers
Maik

Am 08.06.2014 01:38, schrieb Erik Tollerud:
> Hello Maik,
> 
> Just to clarify something: `rotation_matrix` is technically not part
> of the public API for astropy.  It's mainly used internally in
> coordinates. Furthermore, it's not actually used internally for
> anything other than the 'x', 'y', and 'z' use cases. So just be aware
> that this function comes with no "guarantees" like we do for the
> public API.
> 
> If you actually look at the code, you'll see the matricies for the
> 'x', 'y', and 'z' cases are constructed quite differently from the
> generic angle, such that we do *not* expect them to come out exactly
> the same to machine precision.  So I think special-casing [1,0,0],
> [0,1,0], and [0,0,1] to be the same as 'x', 'y', and 'z' is bad,
> because it would cause discontinuous jumps when you're using the
> axis-specified form.  So I would say the fact that they are not the
> same to machine precision is a feature, not a bug (because the
> 'x'/'y'/'z' cases are quite a bit faster to run).
> 
> I am concerned about the "When an array is used with astropy, the
> difference is e-7" bit, though: can you clarify a bit more what you
> mean by "when an array is used with astropy"?
> 
> 
> 
> 
> On Thu, Jun 5, 2014 at 3:13 PM, Maik Riechert <maik.riechert at arcor.de> wrote:
>> Hi Thomas,
>>
>>>>>>> import numpy as np
>>>>>>> from astropy.coordinates.angles import rotation_matrix as rot_astro
>>>>>>> from gohlke.transformations import rotation_matrix as rot_gohlke
>>>>>>> m1 = rot_astro(-0.07687069267192707, [1,0,0])
>>>>>>> m2 = rot_astro(-0.07687069267192707, 'x')
>>>>>>> m2-m1
>>>> matrix([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
>>>>         [  0.00000000e+00,   2.22044605e-16,   1.30859299e-13],
>>>>         [  0.00000000e+00,  -1.30859299e-13,   2.22044605e-16]])
>>>
>>> These differences are close to what might be expected for machine
>>> precision. However, I don't see why we couldn't try and recognize
>>> special cases (for example [1,0,0] -> x). In your case, why not simply
>>> use 'x' since it gives the expected result?
>>
>> I did some more testing and used sympy (symbolic execution/arbitrary
>> precision floats) to give me the "ground truth". For example, using an
>> angle of 0.000004deg the difference between sympy and gohlke/astropy's
>> special cases is at most e-18. When an array is used with astropy, the
>> difference is e-7, which I think is significant as it may add up when
>> multiple rotations are combined.
>>
>> Tell me if you like to see the test code, then I prettify it a bit.
>>
>> Cheers
>> Maik
>> _______________________________________________
>> AstroPy mailing list
>> AstroPy at scipy.org
>> http://mail.scipy.org/mailman/listinfo/astropy
> 
> 
>