[AstroPy] Numerical issues with rotation_matrix

[Erik Tollerud]

Hello Maik,

Your suggestion here makes sense, so I went ahead and made
https://github.com/astropy/astropy/pull/2619 , which fixes the
arbitrary-axis case to be consistent in sign convention and uses an
algorithm that *should* be equivalent to Gohlke's.  Feel free to check
that out and comment on github if that addresses the accuracy issues
(or, if you feel up to it, suggest some test code that I can add there
to ensure the higher accuracy)

On Sun, Jun 8, 2014 at 1:41 PM, Maik Riechert <maik.riechert at arcor.de> wrote:
> 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
>>
>>
>>
>
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy



-- 
Erik T