Coordinate Transformation Precision

[gregreen]

The healpy.rotator.Rotator coordinate transformations between Equatorial (J2000) and Galactic are precise only to ~40 mas. I get this figure by comparing healpy.rotator.Rotator (with coord=['G','C'], and the reverse) with healpy.rotator.euler. I think healpy.rotator.euler, which is translated directly from the IDL function, is a correct implementation of the coordinate transformations.

I think the issue is that the constants in healpy.rotator.get_coordconv_matrix (e.g., "eps" and "e2g") are not defined to sufficient precision.

My results can be reproduced with this gist.

Histogram of residuals for Galactic -> Equatorial (J2000) transformation:

Map of residuals for Galactic -> Equatorial (J2000) transformation (map shows residuals of final answer as a function of Galactic coordinates):

Map of residuals for Equatorial (J2000) -> Galactic transformation (map shows residuals of final answer as a function of Equatorial J2000 coordinates):

[lpsinger]

Not to mention that get_coordconv_matrix recomputes all of the matrices every time, and then decides which one to return using a long string of if...elif statements... instead of a function, this should be a precomputed dictionary at module scope.

[lpsinger]

@gregreen, do you have a patch to contribute?

[astrofrog]

@lpsinger - just to check, is there a reason to not use the Astropy coordinate transformation framework? The only downside would be that Ecliptic coordinates will only be supported in Astropy 1.1+

[gregreen]

@astrofrog, The only downside to using the Astropy coordinate transformation framework is that it would add in Astropy as a dependency.

@lpsinger, I can write up a patch that implements the Equatorial - Galactic - Ecliptic rotations.

[lpsinger]

@lpsinger - just to check, is there a reason to not use the Astropy coordinate transformation framework? The only downside would be that Ecliptic coordinates will only be supported in Astropy 1.1+

The extra dependency is the only downside that I can think of. Healpy works with quite ancient versions of Pyfits because we care about it working on the LIGO Scientific Linux 6 computing clusters. However, we don't actually do any coordinate conversions in the LIGO pipelines because we always work in equatorial coordinates.

So this sounds reasonable to me. @zonca should comment.

[astrofrog]

Just out of curiosity, can Astropy not be installed on the LIGO SL 6 clusters? (too old version of Python?)

[astrofrog]

In any case, maybe it makes sense to wait until Astropy supports Ecliptic coordinates. You may want to consider having an 'LTS' version of healpy which you could continue using with PyFITS etc. on older Python installations which would enable you to make bigger changes such as this (and dropping support for PyFITS) in master. Just an idea :)

[lpsinger]

Just out of curiosity, can Astropy not be installed on the LIGO SL 6 clusters? (too old version of Python?)

As far as I recall, the main issue is that Astropy has not been packaged as an SL6 RPM (either in Scientific Linux proper or EPEL). I don't think that the Python version is an issue, but the old Numpy version may be a problem for Astropy.

Is there anyone in the Astropy project was is able to push on getting SL6 packages made?

[gregreen]

Even without RPM packages, isn't it possible to install Astropy through pip? In much of my work, I use a virtual environment with a fresh pip installation of all the necessary python packages.

[lpsinger]

The LIGO data analysis environment uses the native packaging (RPMs for SL6, .debs for Debian) as much as possible. It would be a hard sell to the admins to rely on something that was not available through yum and apt-get.

[zonca]

I would be in favor of removing code from healpy and use astropy for coordinate transformations.

On Tue, May 5, 2015, 17:32 Leo Singer notifications@github.com wrote:

The LIGO data analysis environment uses the native packaging (RPMs for SL6, .debs for Debian) as much as possible. It would be a hard sell to the admins to rely on something that was not available through yum and apt-get.

— Reply to this email directly or view it on GitHub https://github.com/healpy/healpy/issues/247#issuecomment-99272511.

[astrofrog]

I've asked on astropy-dev if anyone has experience building RPMs for SL. We have packages in other linux distributions, so no reason in principle we couldn't get it in SL.

[astrofrog]

@lpsinger - so based on what you said above, it would be fine to rely on Astropy for coordinate conversions as long as it's treated as an optional dependency, since it's not needed for LIGO?

[ziotom78]

Probably I'm seeing the same problem in a slightly different way. Using Healpy 1.9.0, the rotation matrix from Galactic to Ecliptic is not perfectly orthogonal. The code

r = healpy.Rotator(coord=['G','E'])
print(np.dot(r.mat, np.transpose(r.mat)))

produces the output

[[  1.00000000e+00,   1.23993374e-10,   2.08527973e-10],
 [  1.23993374e-10,   1.00000000e+00,  -4.73150224e-10],
 [  2.08527973e-10,  -4.73150224e-10,   1.00000000e+00]])

The difference between zero and the off-diagonal terms is of the order of 10^-10. This is far larger than the error I get when I manually build an orthogonal matrix:

angle = np.deg2rad(10.0) # Rotation by 10 degrees
mat = np.array([[np.cos(angle), np.sin(angle)],
                [-np.sin(angle), np.cos(angle)]])
result = np.dot(mat, np.transpose(mat))
print(result)
print('Error in the coefficient (1, 1) = {0}'.format(1.0 - result[0][0]))

produces

[[ 1.  0.]
 [ 0.  1.]]
Error in the coefficient (1, 1) = 1.1102230246251565e-16

where 10^-16 is correctly what I would expect from a double-precision calculation.

[zonca]

thanks @ziotom78 @lpsinger , any update on the possibility of healpy to depend on astropy?

[astrofrog]

By the way, astropy is installable on Debian already via apt-get

[lpsinger]

By the way, astropy is installable on Debian already via apt-get

Sorry, we're still using Wheezy on our Debian clusters, so Astropy is not yet available via apt-get for us. I'll check on our reference OS selection again.

[astrofrog]

In any case, I feel like it would be a shame to hold back development of certain areas of healpy just because of system administration issues. Is pip install not allowable? healpy itself doesn't have a debian package either, right?

[lpsinger]

Is pip install not allowable?

Our own data analysis pipelines are deployed using apt on Debian and yum on SL, so their dependencies have to be available as packages.

healpy itself doesn't have a debian package either, right?

It does, but in sid. We have a backport in our own public LIGO Debian repository. It would be tough to backport Astropy to Wheezy because of the Numpy dependency.

[astrofrog]

@lpsinger - I suggested this before, but what about having a special fork of healpy that would be a kind of LTS release that can include only bug fixes from now on, and then going ahead with development using astropy etc in the master branch? Then if there are any bad bugs, they can be backported to the LTS branch and that can be used for LIGO? (in other words, it would be a shame to hold back development because of LIGO requirements)

[gregreen]

@lpsinger - I've suggested this before, but installing healpy (and the rest of the Python stack) in a virtual environment is really the way to go if you want to stay up-to-date with Python packages (and not just healpy). I don't see why pip install -r requirements.txt is a harder sell than apt-get install XYZ.

In any case, healpy development shouldn't be held back by sysadmin choices, especially when there is such an obvious way for the sysadmins to work around the issue.

[lpsinger]

@gregreen, it is highly unlikely that LIGO will switch to using virtual environments for production analysis codes. Our collaboration has made a large investment in system administrators, code librarians, and build infrastructure related to .debs and .rpms. Also, a lot of our codes are autotools-style packages that have mostly C code and a little Python; they are not amenable to virtualenvs anyway.

[lpsinger]

@lpsinger - I suggested this before, but what about having a special fork of healpy that would be a kind of LTS release that can include only bug fixes from now on, and then going ahead with development using astropy etc in the master branch? Then if there are any bad bugs, they can be backported to the LTS branch and that can be used for LIGO? (in other words, it would be a shame to hold back development because of LIGO requirements)

As I have said before, Healpy coordinate conversions are not used by any LIGO code, so it would be fine to rest entirely on Astropy for this functionality. It doesn't need to be backported to Wheezy because we'll probably never need it on our clusters anyway. @gregreen should just submit a patch to modernize healpy.rotator using Astropy. That one submodule won't work on LIGO, but we won't miss it.

[lpsinger]

LIGO is soon going to switch to Scientific Linux 7 and Debian Jessie. Debian Jessie ships with astropy 0.4.2. Will this work with that version of astropy, minus the ecliptic coordinates?

[zonca]

it would be good to port the rotation to use astropy, anybody volunteers to do this? @gregreen ?

[zonca]

I am working on this issue, astropy has different "Ecliptic" frames, here is the difference and std of the difference of the rotation matrix with the healpy one.

BarycentricMeanEcliptic
[[-6.85113372e-06  3.19389173e-07  6.09538085e-07]
 [ 6.93914210e-06 -2.87396167e-06 -4.34651037e-06]
 [ 4.38480405e-06  5.30114701e-06  7.65674704e-06]]
4.886293218443731e-06
BarycentricTrueEcliptic
[[ 6.02756649e-05 -6.91262614e-07 -2.71884052e-05]
 [ 1.44370718e-05  6.38115202e-06 -7.45161630e-06]
 [ 4.40657709e-06 -6.72103374e-05  7.64874468e-06]]
3.2081262203432936e-05
GeocentricMeanEcliptic
[[ 1.12974109e-05  5.76652567e-06  2.42150175e-05]
 [ 8.96613486e-06 -5.19656168e-05 -5.14071812e-05]
 [ 4.39111614e-06  9.15062407e-05  9.38731111e-05]]
4.86495619728729e-05
GeocentricTrueEcliptic
[[ 7.84174745e-05  4.76256155e-06 -3.58287693e-06]
 [ 1.64639292e-05 -4.27099275e-05 -5.45127276e-05]
 [ 4.41181710e-06  1.89960841e-05  9.38648549e-05]]
4.5879798615600055e-05

BarycentricMeanEcliptic is the closest.

[zonca]

If I compare @gregreen test before and after this fix

Current healpy

Using BarycentricMeanEcliptic

[zonca]

same for equatorial (using FK5 from astropy)

the astropy matrix is

array([[ 1.00000000e+00,  4.02386554e-08,  1.24680181e-07],
       [-8.65131460e-08,  9.17482168e-01,  3.97776911e-01],
       [-9.83858346e-08, -3.97776911e-01,  9.17482168e-01]])

the healpy one was:

array([[ 1.        ,  0.        ,  0.        ],
        [ 0.        ,  0.91748214,  0.39777698],
        [ 0.        , -0.39777698,  0.91748214]])

[zonca]

Probably I'm seeing the same problem in a slightly different way. Using Healpy 1.9.0, the rotation matrix from Galactic to Ecliptic is not perfectly orthogonal. The code

r = healpy.Rotator(coord=['G','E'])
print(np.dot(r.mat, np.transpose(r.mat)))

produces the output

[[  1.00000000e+00,   1.23993374e-10,   2.08527973e-10],
 [  1.23993374e-10,   1.00000000e+00,  -4.73150224e-10],
 [  2.08527973e-10,  -4.73150224e-10,   1.00000000e+00]])

The difference between zero and the off-diagonal terms is of the order of 10^-10. This is far larger than the error I get where 10^-16 is correctly what I would expect from a double-precision calculation.

now in #633:

In [4]: r = healpy.Rotator(coord=['G','E']) 
   ...: print(np.dot(r.mat, np.transpose(r.mat)))                                                                                    
[[ 1.00000000e+00  1.97324795e-17 -1.66533454e-16]
 [ 1.97324795e-17  1.00000000e+00  8.38359330e-17]
 [-1.66533454e-16  8.38359330e-17  1.00000000e+00]]

[zonca]

fixed in https://github.com/healpy/healpy/pull/633