Iterating over a SkyCoord object is very slow

I was surprised to see how long it takes to iterate over a multi-valued SkyCoord object.

For example, it takes 19s to iterate over a SkyCoord object that contains 1000 coordinates:

In [1]: import numpy as np
In [2]: from astropy import units as u
In [3]: from astropy.coordinates import SkyCoord
In [4]: ra = np.random.uniform(0, 360, size=1000) * u.deg
In [5]: dec = np.random.uniform(-90, 90, size=1000) * u.deg
In [6]: crd = SkyCoord(ra, dec)
In [7]: timeit crd[500]
100 loops, best of 3: 19 ms per loop
In [8]: timeit [c for c in crd]
1 loops, best of 3: 19.4 s per loop

In contrast, iterating over a list of 1000 SkyCoord objects is roughly 500,000x faster:

In [9]: crd2 = [SkyCoord(ra[i], dec[i]) for i in range(len(ra))]
In [10]: %timeit crd2[500]
10000000 loops, best of 3: 38.7 ns per loop
In [11]: %timeit [c for c in crd2]
10000 loops, best of 3: 28.6 µs per loop

Might there be room to increase the performance of accessing a SkyCoord item? Iterating over a list of coordinates seems like a common task to me.

That does sound quite slow - and we should be able to speed it up since once we have the array skycoord all the parsing etc. is done.

Heh, I got exactly 19.4s for the first test as well. This is definitely worth looking into. It's no surprise that it's slower than __getitem__ on a list, but I wouldn't be surprised if speedups are possible.

Oh, well it's no wonder. SkyCoord.__getitem__ instantiates a new SkyCoord. Nevertheless the second example is still faster even if you include SkyCoord object instantiation:



In [7]: %%timeit
   ...: crd = SkyCoord(ra, dec)
   ...: [c for c in crd]
   ...:
1 loops, best of 3: 19.4 s per loop

In [8]: %%timeit
   ...: crd2 = [SkyCoord(ra[i], dec[i]) for i in range(len(ra))]
   ...: [c for c in crd2]
   ...:
1 loops, best of 3: 3.16 s per loop

So I suspect there's room for improvement.

Yes, we could certainly use a back-door way of initializing a SkyCoord that just copies what's there and doesn't do any of the heavy validation that normally takes place. It will still never be super fast because there is a lot of generalized machinery under the hood in the full SkyCoord object that is the output of __getitem__.

One lesson is to not do things like crd[i] in performant code, but instead use crd.ra[i] etc.

Just thinking out loud, but I wonder if we can do a similar trick here to what was done for rows in tables, i.e. when the user accesses SkyCoord[0].ra access SkyCoord.ra[0]. Just an idea :)

@astrofrog

@astrofrog That would be doable, I think but also tricky. It would require some sort of SkyCoord proxy class that can do quick lookups in the ra and dec arrays, but otherwise is not a fully initialized SkyCoord object until and unless any other functionality is accessed. But it's not a bad idea, especially for cases like database ingest.

At the very least, even without my idea, it should be fast to initialize a SkyCoord from an existing one. The overhead should just be accessing the items in the arrays. Frames are immutable, so the existing frame can be used as-is.

A shortcut for instantiating from self would definitely make a big difference here. It looks like most of the time is wasted in _parse_coordinate_arg, which presumably shouldn't be necessary for this case.

Could somebody with more coordinates expertise explain what the deal is with FRAME_ATTR_NAMES_SET? It seems to me this could just be a a TransformGraph property in its own right that could be cached.

Lots of time is spent in the sky_coordinate module calculating this over and over again even though it rarely changes. What is the purpose of making it look like a module-level constant?

I tried out my above suggestion (I left FRAME_ATTR_NAMES_SET as a callable but it just returns frame_transform_graph.frame_attr_names which is a lazy property that is cleared when the transform cache is invalidated).

Now I get:

In [2]: crd = SkyCoord(ra, dec)

In [3]: %timeit [c for c in crd]
1 loops, best of 3: 3.33 s per loop

~~which is barely worse than the list case~~ sorry, still a lot worse, but better than the existing code. Could still do even better but this alone made a huge difference, and probably helps in some other places as well because that FRAME_ATTRS_NAME_SET was being called all over the place.

@embray

@embray - IIRC FRAME_ATTR_NAMES_SET originally was a constant, but it was made dynamic in order to support dynamically-added frames within the transform graph. Of course you are entirely correct that in practice it very rarely changes and that's kinda why I left it looking like a constant.

In any case we should cache it. We just need some fast key available from the transform graph object to know if it has been updated in any way. @eteq - any suggestions?

About fast-tracking init with self, now I'm starting to wonder... Our nice three-layer system of SkyCoord, Frame, Representation probably means that you need to make a fast path all the way down. Could be tricky.

@parejkoj

with #6941, this is now about as fast as possible

For the special case of iterating over a coordinate array where the iteration values don't need to be unique objects, there is potential for a lot more speedup (factor of ~25-ish).

@parejkoj - can you point to the facebook discussion? I went looking in the Python in Astronomy group and didn't see anything relevant.

@taldcroft

@taldcroft - are you sure? I thought that was relative to just #6940 - with #6941 things were very similar to what you had.

@taldcroft

@taldcroft it was a response to this post: https://www.facebook.com/groups/astropython/permalink/2030892927155387/

Looks like @bsipocz already responded.

@parejkoj

@parejkoj - I only had #6941 in mind when I replied there (and remembered the benchmark graph, too).

@mhvk

@mhvk - this was the original fastiter() implementation that re-used the same template every time instead of doing a deepcopy. This was a factor of 500 faster than master at the time, or probably about 25 times faster than current master and 3.0. It certainly comes with a caveat that you get the same actual object back, so list(sky_coords) give you a list of the same thing. BUT there are use cases where this is OK.

Ah, yes, my memory clearly was with the comparison with the unique objects. I stand corrected. Whether this further speed-up matters is of course less clear: eventually one presumably wants to do something with the SkyCoord, which likely involves a transformation that will be much slower.

Anyway, good that the aim for 3.1 is to get the regular paths as fast as we can!

On the "doing something" bit, I imagined doing something with the ra and dec attributes, not a transformation. But just accessing a component attribute of a scalar SkyCoord is very slow, 0.36 msec on my mac. This is 3000 times slower than accessing a direct object attribute. Clearly SkyCoord or Frame has a lot going on there, but this looks like another avenue for performance improvement for 3.1.

Indeed, that is strange. And the more directly one goes to the real data -- avoiding the __getattr__ overrides -- the less the problem becomes (see below). For sc.ra, profiling suggests most time is spent in baseframe.get_representation_component_names -- so this should not be too hard to solve -- but definitely worth a separate issue!

sc = SkyCoord(0.*u.deg, 1.*u.deg)
%timeit sc.ra
# 10000 loops, best of 3: 163 µs per loop
%timeit sc.data.lon
# 100000 loops, best of 3: 2.74 µs per loop
%timeit sc.frame.data.lon
1000000 loops, best of 3: 635 ns per loop
%timeit sc._sky_coord_frame._data.lon
1000000 loops, best of 3: 278 ns per loop

Conclusion: should be at least 200 times faster...

p.s. It seems clear that as long as object are immutable, we can make __getattr__ much faster -- indeed, one could just populate self.__dict__...

It seems clear that as long as object are immutable ..

SkyCoord and Frame objects are mutable since #5575 (astropy 1.3, Dec 2016). The public docs show how to change these in place.

http://astropy.readthedocs.io/en/latest/coordinates/inplace.html

Similar problem here.
I'm trying to calculate the Moon's Alt-Az coordinates for a given time at my site.
The calculation of Alt-Az angles takes too much time.
Any idea why?

Script:

from datetime import datetime

from time import time as tttime

from astropy.coordinates import get_moon
from astropy.coordinates import EarthLocation
from astropy.coordinates import AltAz
from astropy.coordinates import Angle
from astropy import units as U
from astropy.time import Time as tm

s = tttime()
FORMAT = "%Y-%m-%dT%H:%M:%S.%f"
date_obs = datetime.strptime("2018-07-14T08:50:32.48", FORMAT)
e = tttime()
print("Time Object created in {} seconds".format(e - s))

s = tttime()
lat = Angle("+39:46:56 degree").degree
e = tttime()
print("Angle Object created in {} seconds".format(e - s))

s = tttime()
long = Angle("+41:13:58 degree").degree
e = tttime()
print("Angle Object created in {} seconds".format(e - s))

s = tttime()
site = EarthLocation(lat=lat, lon=long, height=float(3114) * U.m)
e = tttime()
print("Site Object created in {} seconds".format(e - s))

s = tttime()
moon = get_moon(tm(date_obs, scale='utc'))
e = tttime()
print("Sky Object (Moon) created in {} seconds".format(e - s))

s = tttime()
frame_of_site = AltAz(obstime=date_obs, location=site)
e = tttime()
print("Frame Object created in {} seconds".format(e - s))

s = tttime()
object_alt_az = moon.transform_to(frame_of_site)
e = tttime()
print("Alt-Az Calculated in {} seconds".format(e - s))

Outpu:

Time Object created in 0.0 seconds
Angle Object created in 0.030791044235229492 seconds
Angle Object created in 0.0 seconds
Site Object created in 0.0009884834289550781 seconds
Sky Object (Moon) created in 0.04990100860595703 seconds
Frame Object created in 0.0 seconds
Alt-Az Calculated in 8.166500806808472 seconds

Python version:
Python 3.7.4 (default, Aug 9 2019, 18:34:13) [MSC v.1915 64 bit (AMD64)] :: Anaconda, Inc. on win32
Astropy version:
'3.2.1'
Numpy version:
'1.16.5'

@mshemuni

@mshemuni - that is really long! I just checked locally and while not as extreme as what you get, the last step does take 2 seconds. However, just repeating the exact same script in the same session reduces it to 0.03 seconds, which suggests it is some kind of initialization issue.

Trying a bit more, I think I can localize this to reading the IERS tables (which are needed for getting accurate alt,az); if after you're import of Time, you add tm.now().ut1 (which also forces the tables to be loaded), you'll see that your other statements are decently fast.

So, I think this is not an issue for SkyCoord itself, though I am a surprised how long the loading takes... Looking yet a bit further, it seems clear it is because the web site for the IERS table gets contacted even if you already downloaded the file, which seems incorrect behaviour. I think this is probably part of #9555

p.s. To be a bit more specific, adding the following before anything you do makes everything fast:

from astropy.utils import iers
iers.IERS_Auto.open()

@mshemuni

I think the refactoring in the caching systems has solved this issue, these are the numbers I get on master when copy pasting the snippet from above, as-is.

Time Object created in 0.0007829666137695312 seconds
Angle Object created in 0.024201154708862305 seconds
Angle Object created in 0.00028586387634277344 seconds
Site Object created in 0.0006892681121826172 seconds
Sky Object (Moon) created in 0.04855179786682129 seconds
Frame Object created in 0.0006310939788818359 seconds
Alt-Az Calculated in 0.869938850402832 seconds

@mshemuni - Could you try to update for the new release candidate and report back whether it solved the problem? It can be pip installed: python -m pip install --pre -U astropy

@mshemuni

Exp1:

p.s. To be a bit more specific, adding the following before anything you do makes everything fast:
from astropy.utils import iers
iers.IERS_Auto.open()

Added two lines to my code right after import lines. Res:
Alt-Az Calculated in 29.960163593292236 seconds

Exp2:
Disconnected the computer from network. Res:
Alt-Az Calculated in 0.44874143600463867 seconds

Exp3:

I think the refactoring in the caching systems has solved this issue, these are the numbers I get on master when copy pasting the snippet from above, as-is.
Time Object created in 0.0007829666137695312 seconds
Angle Object created in 0.024201154708862305 seconds
Angle Object created in 0.00028586387634277344 seconds
Site Object created in 0.0006892681121826172 seconds
Sky Object (Moon) created in 0.04855179786682129 seconds
Frame Object created in 0.0006310939788818359 seconds
Alt-Az Calculated in 0.869938850402832 seconds
@mshemuni - Could you try to update for the new release candidate and report back whether it solved the problem? It can be pip installed: python -m pip install --pre -U astropy

Did the update. New astropy version is '4.0rc1'. Res(Connected to network):
Alt-Az Calculated in 0.7916419506072998 seconds

The update seems to solve the problem.
Thank you @mhvk and @bsipocz for help.

Navigation Menu

Repository navigation

Iterating over a SkyCoord object is very slow #3323

34 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Iterating over a SkyCoord object is very slow #3323

Description

Activity

astrofrog commented on Jan 21, 2015

embray commented on Jan 21, 2015

embray commented on Jan 21, 2015

taldcroft commented on Jan 21, 2015

astrofrog commented on Jan 21, 2015

embray commented on Jan 21, 2015

astrofrog commented on Jan 21, 2015

embray commented on Jan 21, 2015

embray commented on Jan 21, 2015

embray commented on Jan 21, 2015

taldcroft commented on Jan 21, 2015

taldcroft commented on Jan 21, 2015

34 remaining items

taldcroft commented on Jan 19, 2018

mhvk commented on Jan 19, 2018

parejkoj commented on Jan 19, 2018

bsipocz commented on Jan 19, 2018

taldcroft commented on Jan 19, 2018

mhvk commented on Jan 19, 2018

taldcroft commented on Jan 19, 2018

mhvk commented on Jan 19, 2018

taldcroft commented on Jan 27, 2018

mshemuni commented on Dec 2, 2019

mhvk commented on Dec 2, 2019

mhvk commented on Dec 2, 2019

bsipocz commented on Dec 2, 2019

mshemuni commented on Dec 3, 2019

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions