Vectorized keplerlib.propagate

[JoshPaterson]

I have successfully vectorized keplerlib.propagate! All of it's arguments can now be arrays. Just for testing purposes I added a function to data/mpc.py that can take a comets dataframe and convert it into a _KeplerOrbit object that represents multiple comets.

At the moment _KeplerOrbit objects that represent multiple orbits don't play well with other VectorFunction objects; I'm not sure what needs to change to fix that.

I haven't added any tests yet since none of the official functions in data/mpc.py can handle multiple orbits at a time yet.

Here's an example script:

from skyfield.api import Loader
from skyfield.data import mpc
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from numpy import arange, linspace, rollaxis
import matplotlib.pyplot as plt

load = Loader('~/projects/skyfield-data')

ts = load.timescale()

with load.open(mpc.COMET_URL) as f:
    comets = mpc.load_comets_dataframe(f)

comets_df = (comets.sort_values('reference')
          .groupby('designation', as_index=False).last()
          .set_index('designation', drop=False))

rows = comets_df.iloc[20:40]

comets = mpc._comet_orbits(rows, ts, GM_SUN)

t = ts.utc(linspace(2000, 2030, 100))
position = comets.at(t).position.au

plt.rcParams['legend.fontsize'] = 10

fig = plt.figure()
ax = fig.gca(projection='3d')

for i, (x, y, z) in enumerate(rollaxis(position, 1)):
    ax.plot(x, y, z, label=rows.designation[i])

ax.legend()

plt.show()

Which generates this plot:

[brandon-rhodes]

Thanks for working on this! I have made a few comments but the adjustments should be fairly minor — then we'll be ready to merge. As long as the new vector-loading starts with an underscore, your tests can come later, as long as we keep passing the tests for the single-comet case!

[JoshPaterson]

I edited the script to use enumerate instead of np.arange in the for loop, which neatens things up nicely.

[brandon-rhodes]

Oh, no! Why did GitHub not notify me about your update, and only add this to my notifications when you commented 6 hours ago? I didn't mean to leave you waiting several days. I wonder if I left the tab open and it updated it and so it thought I'd seen it?

In any case, I have two more quick comments!

[brandon-rhodes]

And I'll now promptly close this tab in case leaving it open might eat notifications. Thanks for the work you're doing on this!

[JoshPaterson]

Thanks for the work you're doing on this!

I'm glad to be able to do something useful!

[JoshPaterson]

I made a few other improvements also:

• used np.copyto in a number of places instead of boolean indexing
• changed & to * for boolean multiplication
• fixed a bug that came up when t0 was a scalar

[JoshPaterson]

It took me a while, but I figured out that
upper if lower>upper else (lower+upper)/2
is equivalent to
min(upper, max(lower, (lower+upper)/2))
This also occurs further down.

[JoshPaterson]

I think this is good to go now.

[xmichaelx]

Is doing something similar for mpcorb_orbit in mpc.py hard? Because that's precisely what I would need and be happy to add.

As I understand this allows for computation of positions from multiple orbits all at once? Can this approach be used with observe method? I'd like to check whether object is in FOV of telescope for large subsets of MPC.

[brandon-rhodes]

Drat, I apparently made this comment a few days ago but instead of submitting it, it dangled as part of a review. Let me try hitting "Comment" and that should make it appear.

[brandon-rhodes]

One final idea, because it looks so painful to be creating those specially indexed copies of the Pandas columns:

parabolic = (e == 1.0)
p = (1 - e*e) / (1.0 - e + parabolic)
p[parabolic] += 2.0
p *= rows.perihelion_distance_au.values

Would that produce the same output, but without needing the special indexing?

[JoshPaterson]

Nice! That's a clever way to avoid the division by 0. I added that in.

[brandon-rhodes]

@JoshPaterson, I've added the URL for this pull request to my daily to-do list so that I keep looking at it whether GitHub shows it in my notifications or not; it looks like I keep missing comments from you, and I want to get it merged by cutting down the turn-around time. I hope it hasn't made your holidays too frustrating!

[JoshPaterson]

I hope it hasn't made your holidays too frustrating!

Don't worry! I know we're all volunteers, so instant feedback is nice when it happens, but I don't expect it. Also maybe I should try to get in the habit of clicking the 'request review' button so it will be sure to show up on your dashboard.

[JoshPaterson]

@xmichaelx Yes, asteroids will definitely be supported once this gets officially released. The changes in this PR are mostly limited to keplerlib.propagate; Brandon will probably make the other necessary changes since they're in areas where he is more knowledgeable. I added the _comet_orbits function just so I could make that plot as an informal test.