WIP: improve speed of coordinate transformations for arrays of obstime - #6068

Hi @bwinkel,

First of all, thanks for submitting a PR. It's great to have people contributing, especially when they provide such detailed notes and justifications!

That said, I'm afraid I'm not a fan of this PR for a couple of reasons.

The first is a philosophical position, that I think it's important that a hypothetical user of astropy can use the co-ordinate package without worrying about accuracy. Obviously there's a trade-off to be made between accuracy and speed but I think we need to think carefully before intentionally degrading the accuracy of an existing code. For example, someone currently using astropy to predict lunar occultations might get bitten by this.

The second is that the use of erfa.apci was introduced in #5273 to allow the use of (e.g) JPL ephemerides when performing coordinate transformations and light-travel-time corrections. I'd strongly argue there's a need to keep the GCRS-ICRS coordinate transform as-is (it's used in light-travel-time corrections), so your PR would make the ICRS-CIRS transform inconsistent with all the others, which I'm also not keen on. As a minimum, I'd want the proposed transform to be a non-default option.

Since your own use case seems rather extreme to drive the default behaviour in astropy, perhaps the best way forward is to temporarily replace the transformations in astropy's frame transform graph in your own code - i.e if you want to use any particular function to transform CIRS->ICRS, just decorate that function with @frame_transform_graph.transform(FunctionTransform, CIRS, CIRS) and astropy will use your transformation and not the default.

Thanks a lot for the feedback, Stuart. I suspected that erfa.apci was there for the potential use of ephemerides. It'd probably be useful to add a comment to the source code. By default, astropy doesn't use to-date ephemerides though, which already leads to a certain loss in accuracy. I don't know how much this is; if I remember correctly the SOFA manual mentions microarcseconds somewhere. And you're absolutely right, the default behavior should be to have at least this kind of precision - either by not using caching at all by default, or by using a suitable MJD_RESOLUTION.

Note, that caching the astrom struct is not tied to the use of acpi13. One could easily adapt the function to work with the original acpi-based calculation.

I must also say that the use case I described is all but extreme. Almost all radio telescope, I know of, nowadays sample data with a resolution of ~1 s; for some types of observations, even faster (1 ms to 100 ms). Typical observing sessions are a couple of hours. At the moment, they just don't use astropy for these transformations, but I could imagine that some would, if possible.

Thanks for the response. I still don't quite follow why sampling data at fast rates means that you have to be able to convert positions at 10**5 times in less than one second? Nevertheless, it is definitely true that for many people the astropy coordinate package is not as fast as it needs to be; I didn't mean to sound so negative above.

By describing your use case as extreme I simply mean that some users will demand high speed, and some will demand high accuracy and we need to balance the two (or find a way to satisfy both). I do think this PR could be a really valuable addition. However, given my concerns above, I think this PR should probably end up satisfying the following requirements:

1. Leave the existing behaviour as the default.
2. Add a way of modifying the behaviour, such as a ScienceState, e.g:

with approximate_transforms:
    fk5c.transform_to(GCRS(obstime=Time("J2000")))

3. Apply those changes to all relevant parts of the transform graph. One should get the same results if you go from ICRS->AltAz as ICRS->GCRS->AltAz, for example.
4. Carefully document the effect on precision in both position and time.

(2) and (3) could be achieved by writing a utility function to calculate the astrom context, that calls either the existing code, or your faster version depending on the results of the ScienceState object.

I like the idea of introducing a context manager for low-precision transforms. Will wait a little bit for other comments and then try to find out, how one could implement it. Thanks a lot!

take a look at how it's done here. One thing that would need some thought is what would happen if a user both specified that the JPL ephemerides be used, and approximate transformations performed. Apart from giving that particular user a talking to, of course.

Again, as a proof of concept, I added a solution based on the ScienceState context. In the new astrom_manager module, there is a helper function get_astrom which encapsulates the MJD binning and allows to query the astrom structs. It has a tcode parameter, where one has to provide one of the erfa function names ('apio13', 'apci', 'apcs'). I didn't bother for now to setup a more fancy API for this, mainly because only three erfa (astrom) functions are used throughout all transformations.

It works like this:

# see first comment for full code example
with transform_precision.set(precision):
    azels = cs.transform_to(tframes)

where precision is a float. If 0 (or less than 1e-3), no MJD binning is performed. Otherwise, precision is the bin width in seconds.

I also realized, that the accuracy loss is much less than in the plots I made yesterday. Most of the high Azimuthal deviations where at the poles. I attach new plots, which show the angular-distance of the low-precision points to the original solution. For a binning of 600 s, the maximum error is 1.2 mas, for 10 s binning, I get 0.04 mas.

I know, there is still a lot of work to do (tests, documentation, etc.). But I wanted to ask for your opinion first. Especially, if the design of the helper function is acceptable (I find it a bit clumsy, to be honest). Any ideas, how to improve?

@bwinkel actually, I rather like the helper function to provide an astrometry context. It certainly cleans up the code a lot where it is actually used.

I think that it might be a good idea to present the user with discrete choices for the precision (i.e a set of pre-chosen time resolutions) and document their accuracy. Otherwise we're giving end users work to do, working out if their chosen time resolution is precise enough for their needs.

with transform_precision.set('low'):
    azels = cs.transform_to(tframes)

I think at this point the general idea is firm enough for people like @eteq and @mhvk to look over and see if they approve in principle before putting in too much work with documentation/tests etc.

One thing that I don't like about the helper function is, that it is vital to plug the appropriate frame into it (in each transform function, two frames are involved, but only one of them must be used for the helper). Since the function is not exposed to the user, this may be acceptable, but it's not really nice from a design perspective.

I'm also thinking about how one could do some caching. If, for example, a conversion from AltAz to ICRS is performed, the MJD binning happens twice with the same obstimes. It would be cool, if the binned astrom values would be cached. At 10-min resolution, one would only need to store 720 instances to cover a full day. Furthermore, this would also speed-up the case where one repeatedly calls the transform functions with non-array frames (but also different obstimes). Does astropy have a cache utility class that could be used for such purposes?

Good to have a detailed look at coordinates performance -- this is definitely a bottleneck, and the API is now sufficiently stable that we can and should improve it. That said, I think it is still at a state where we should optimize individual pieces, rather than try to work around slow code by doing interpolation. So, I'm much more keen on ideas like caching than on sampling at different time samples or making lower-precision approximations.

@bwinkel - from what you write, it seems you did quite a bit of profiling already, and I was surprised that getting tdb was one of the bottlenecks, in particular since Time does cache tdb. Is it nevertheless done multiple times? I think the best way forward might be to address specific bottlenecks (which also has the advantage that the work would come in small pieces, which is easier to review).

p.s. There is some related work in the pipeline, which could use input: @dmopalmer has been looking at the initialization of SkyCoord (most recently, #5952) and I at that of the representations (#5598).

I agree with your reasoning, @mhvk, but in this particular case, the main bottleneck is indeed the calculation of the the astrom context. (A while ago I wrote my own very thin wrapper around the SOFA library, and am convinced, that the resulting performance is directly related to the underlying algorithms, not the implementation details.) Even caching is not an option (if full accuracy is desired), because one would need to cache an insane amount of astrom instances. (How do you cache double values anyway?) And it would not help in almost all use cases, because one still needs to calculate it for the first run (and probably this is often the only run). By-the-way, the trick to bin the times was proposed in the SOFA manual, so even they were aware of the performance issues.

Note, while getting tbd is indeed a bit slow (during the first run, before caching), it is by far not as bad as the astrom calculation itself. In the current commit, I nevertheless used tt to save another 2 s in the example.

I'm quite sure, that for >90% of the use cases, an precision of 1-10 mas would be more than sufficient. How many telescopes reach that level of pointing accuracy...?

Hmm, you can even do this with 2D arrays (only if time is on axis 0). So here is another example, which also shows how the deviations are distributed.

num_declinations = 35
num_times = 10 ** 4

with NumpyRNGContext(1):
    _ra = np.repeat(0., num_declinations)[:, np.newaxis]
    _dec = np.linspace(-85, 85, num_declinations)[:, np.newaxis]
    _time_delta = np.linspace(0, 1, num_times)

ra = Longitude(_ra, unit=u.deg)
dec = Latitude(_dec, unit=u.deg)
cs = SkyCoord(ra, dec, frame='icrs')
times = Time(54000. + _time_delta, format='mjd')

loc = EarthLocation(lat=41.3 * u.deg, lon=-74 * u.deg, height=390 * u.m)
tframes = AltAz(obstime=times, location=loc)

with transform_precision.set(0):
    azels0 = cs.transform_to(tframes)  # 3.48 s

with transform_precision.set(200):
    azels200 = cs.transform_to(tframes)  # 348 ms

with transform_precision.set(600):
    azels600 = cs.transform_to(tframes)  # 253 ms

sep600 = np.degrees(coordinates.angle_utilities.angular_separation(
    azels0.az.rad, azels0.alt.rad,
    azels600.az.rad, azels600.alt.rad,
    )) * 3600. * 1000


az = azels.az.to(u.deg).value
alt = azels.alt.to(u.deg).value

clf()
scatter(
    az, alt, c=sep600,
    edgecolor='none', cmap='viridis',
    )
colorbar()
xlim((0, 360))
ylim((-90, 90))
xlabel('Az [d]')
ylabel('Alt [d]')

The colorbar shows the error in mas.

On the actual use case you describe: I'd guess that the way to implement that would be to sample coarsely in the ICRS->CIRS and then more finely in CIRS->AltAz (oddly, if I try CIRS->AltAz, I get transforms of CIRS->ICRS and back, which seems wrong -- and doesn't help to speed things up).

A quick profile of ICRS->CIRS suggests, as you wrote, the erfa routines epv00 (called twice) and pnm06a (called once) are major culprits. Here, epv00 is called once to get the position of the Sun, and once for the Earth, though really this could be done in a single call (since Sun is calculated from earth_bary and earth_heliocentric). That would be one simple but nice optimization.

Separately, it would be good to check that our vectorising of erfa is actually reasonably close (within a factor 2) to a pure C iteration over a large number of elements.

Oops, sorry, that last comment was put in before I saw your ones. So, the goal is really to avoid many calls to the erfa routines. As noted, we can easily get rid of one (saving 30%). Right now, I'm trying to understand why CIRS->AltAz goes via ICRS; without that, one could do the fine sampling there (CIRS hardly changes with time), and things might be fine at full precision.

So, one way to optimize is indeed to undersample the transformation to CIRS and then interpolate. This gives a speedup by basically the undersampling factor while keeping very good precision:

import numpy as np
from astropy.coordinates import SkyCoord, EarthLocation, AltAz, CIRS
from astropy.time import Time
from astropy import units as u


sc = SkyCoord(90.*u.deg, 10.*u.deg)
loc = EarthLocation(lat=41.3 * u.deg, lon=-74 * u.deg, height=390 * u.m)
t = Time(50000, np.linspace(0,1,10000), format='mjd')
fr_altaz = AltAz(obstime=t, location=loc)


def fast(sc, fr_altaz, factor=100):
    obstime = fr_altaz.obstime
    fr_cirs = CIRS(obstime=fr_altaz.obstime[::factor])
    sc_cirs = sc.transform_to(fr_cirs)
    # interpolate
    dt = (obstime - obstime[0]).sec
    dtp = dt[::factor]

    def interpolate(q):
        return np.interp(dt, dtp, q.value) * q.unit

    sc_cirs2 = SkyCoord(ra=interpolate(sc_cirs.ra),
                        dec=interpolate(sc_cirs.dec),
                        frame=CIRS(obstime=t))
    sc_altaz = sc_cirs2.transform_to(fr_altaz)
    return sc_altaz.ravel()


sc_altaz = sc.transform_to(fr_altaz)
sc_altaz2 = fast(sc, fr_altaz)
print('std of difference', sc_altaz.separation(sc_altaz2).std().to(u.mas))
# std of difference 0.199768mas

Timing both:

%timeit sc.transform_to(fr_altaz)
# 1 loop, best of 3: 4.51 s per loop
%timeit fast(sc, fr_altaz)
1 loop, best of 3: 139 ms per loop

On the above example: I wrote it without thinking too much about how best to handle things, in a way it was just to convince myself that even with the present machinery, given some thought, the problem can be solved. This is important mostly in the sense that if a solution is not much more difficult than a context manager, then I think we might as well document the solution rather than add more code.

Even so, I think there is a question of whether we can do better for our users without something for them to set up, i.e., whether it is possible to get full precision much faster. E.g., one slow part is get_position_barycentric_posvel, which is fully within our control: we could have a heuristic to use interpolation where possible (size quite large, very small time differences). This would in fact work with JPL ephemerides too (since those already use polynomial interpolations, their speed is already faster -- by about a factor 3 -- but the fast one gets sped up by about the same amount).

Interpolation is certainly a viable option, and I'm all for it, if it can be made work. But make no mistake, it would need very careful heuristics. It will be necessary, to handle cases where the coordinate pairs are not homogeneously sampled in time, or have relatively large gaps. So a clever analysis of the sampling would be needed, adding quite a bit of complexity to the interpolation function.

For me either solution is fine, because I think that both methods potentially provide sufficient speed-up. It'd be good to hear some other opinions on this, before we proceed. I know, that @eteq and @adrn are working on implementing velocities. One should keep this in mind, as well.

By-the-way, interpolation of spherical angles is difficult. You'd need to handle the case where the longitude crosses zero correctly and the quality of the solution will depend on latitude, I think. A better way might be to convert to Cartesian representation and do the interpolation on the (x, y, z) vectors.

@bwinkel - sorry for being somewhat slow in getting up to speed in understanding where even the problem is.

You're definitely right that interpolation should be done in the cartesian frame; and I guess the fact that I didn't immediately do that makes clear how useful it is for things to "just work", rather than letting people reinvent wheels. Nevertheless, I think the example is useful in helping to think whether it is possible to find solutions that benefit all of astropy rather than a specific use case. E.g., a general interpolation routine would, I think, fall in that category! It would also be useful for, e.g., moving objects. The application of such a routine to make telescope tracking better might perhaps be better in a telescope_control affiliated package (that might even be easier to find than some option in the coordinates subpackage...).

Overall, I must admit I'm still rather weary of introducing complexity in the astropy core for this specific case, and would suggest to try to break the problem into pieces that are each useful (say, using just one call to epv00; a general interpolation routine applied to astrom; etc.)

p.s. A propos velocities: once we have this, you would get not just alt and az, but also pm_alt and pm_az, so the time spacing could be considerably reduced.

Hmm, using interpolation to speed up transformations would be relatively easy to put into another package, but how to deal with the heuristics in a proper and very general way is still unclear to me. The MJD-based binning of the astrom structs, however, could hardly be done in a separate package, unless there was a low-level API that exposes access to the relevant lines of code or one would like the burden to keep track of changes to the implementation in astropy.coordinates. And the former is probably not going to happen for design reasons.

Out of pure curiosity, I slightly modified the get_astrom helper to do linear interpolation of Earth-PV, Heliocentric, and CIS (Cartesian) vectors (again based on binned MJD values), while using full MJD
resolution for everything else. It turns out, that one gains a significant boost in accuracy (e.g., with a binning of 600s one now gets an accuracy of sub-microarcseconds, instead of milliarcseconds). And one looses only a little bit of performance compared to pure MJD binning - which means that the performance boost is still significant compared to the current astropy implementation (in the above example a factor of ~100).

I re-created the above plots with the new method (errors are now in microarcsec!):

The advantage of the new approach is also, that it should work with a differential approach to implement velocities, which would not have been the case for the pure MJD binning, submitted initially.

I know that the PR is not a minor one. On the other hand, the default behavior is unchanged and there is an enormous speed gain for some use cases. So, besides the PR-size issue - is there a strong argument against such a change?

I have to say I'm not at all opposed. I have my own use cases (telescope scheduling) where this kind of speedup is essential. The loss of precision is basically insignificant for the vast majority of users and the default behaviour is unchanged. I can't see a strong argument against, but I'm not the finally arbiter, so you'll have to convince people like @mhvk and @eteq.

If you decide to take this further, as well as adding tests and documentation, I would try and make sure you keep an eye on issues with dimensionality - in general the obstime property of a frame (like most other frame attributes) can be an N-D object, so your interpolation code will have to deal with those cases.

Another thing to check is that this doesn't have unintended consequences - e.g with an older state of this PR I noticed that get_moon had slowed down significantly, perhaps due to slowdowns in the transformation from GeocentricTrueEcliptic to ICRS in calc_moon. I would benchmark a series of typical coordinate transformations, and ensure they don't slowdown significantly.

I took a look at trying to whip this into shape for 3.1, but sadly I don't think it can be done. The performance improvements in 3.1 have ended up more in the "fix obvious flaws" category rather than "improve some of the use cases" as this PR is trying to do.

That said, I agree with @mhvk that this is a promising direction to go (and @bwinkel has already put in a fair amount of work here!), so re-milestoned it to 3.2. Obviously, it's been over a year without much progress on this... But perhaps now that we've gotten some of the low-hanging fruit done for 3.1 it's a good idea to pursue this idea again?

@bwinkel @maxnoe @StuartLittlefair @mhvk @adrn - Perhaps we should try to have a virtual coordinates meetup for this and related "more radical performance approaches", either as part of the coordination meeting this december, or just as an independent telecon? We could try to chart out a real path to implementing this idea and perhaps similar use cases as well. Hopefully then this good idea would no longer sit languishing for lack of sufficient review/clear direction.

Thanks, @eteq! I'm certainly available for a meeting. Just drop me a note, when you know the date.

Hi,

I have no experience with the coordinates code base.

I'm just in gamma-ray astronomy, where we have lots of coordinate trafos and 50 mas is absolute overkill by at least 3 orders of magnitude.

So having some kind of way to get faster trafos at the expanse of precision would be great.

I think I'd like a context manager as API

I have no experience with the coordinates code base.

That being said, I'd like to contribute anyway if I can

Good idea @eteq.

I won't be able to attend the coordination meeting, but a virtual meet up would work well for me

Of course, for the actual goal of this PR, this will generally not be important (the pressure and temperature can likely be assumed constant), but if we provide something that aims to be general, it should work in all cases. Which brings me back to my more general concern, that this remains somewhat of a hack for a specific purpose, where the cause of the problem being solved is time that is very fine-grained, something that does not seem that likely to happen outside of that particular purpose.

One could have a very different point of view: Microarcsecond precision is only required for very few purposes. Every other purpose get's slowed down up to the point of non usability.

At the moment, astropy's coordinate system is practically unusable for data analysis e.g. for @cta-observatory because of this.

Also: the tests are passing, so this PR does not change any required features. If what you say is required, there should be a test, right?

Also: the tests are passing, so this PR does not change any required features. If what you say is required, there should be a test, right?

this is not so simple. The tests here don't pass as they were last run quite a while ago, and now there is conflict. And one also has to be a bit careful with any API changes as once it goes into a release prematurely we stick with it for a while even though better ways come up in a discussion.

I agree with @bsipocz here, and this is also why I didn't pursue this further for a long time. If such a feature is wanted by the community (I honestly hope so), then one should first think about how it should be implemented in terms of design and API. My original solution went very deep into the core of the coordinates package. Perhaps, a better way of doing things can be identified in a meeting.

How can we contribute here? This is a major issue for the Cherenkov Astronomy Community that regularly poses the question to move away from astropy.coordinates to something self-made for speed reasons.

I don't want this to happen, but for that we need speedier alt/az ↔ ICRS transformation and we don't require much precision (gamma-ray angular resolution is around 3 arcminutes, so conversions with a precision up to arcseconds would be totally acceptable.

I'm not sure what's the best to do with the milestone at this point. I'm bumping it to 4.2 now, but ideally the convo should be moved to an issue (e.g #8407) and the PR closed.

Thanks @bsipocz. I still think that the proposed functionality would be useful (there have been several instances since I opened the PR, where I had to use other astrometry functions, because astropy was too slow). Should the sub-package maintainers decide to give it another try, I'm willing to help of course.

I want to look into this again as we have a dire need to speed up this use case.

What could be the steps forward here?

Note that this PR is years old and I'm not even sure, if the underlying code, on which this was built, is still in the coordinates sub-package. Not sure, how to proceed, as the original PR was controversial.

@bwinkel I know, this is why I was asking. Have been following this from the start.

I'll have a look into if it can be rebased onto the current master.

Thanks a lot @maxnoe! Let me know, if you need help.