[RFC - do not merge] First `OrbitArray` API prototype #1445

s-m-e · 2022-01-18T21:34:10Z

No description provided.

s-m-e · 2022-01-18T21:35:35Z

I have put some of my thinking down into a notebook and actually built a little prototype - which does Brownian motion, not Orbits, to make stuff simple enough to reason about the user-facing API as well as some internal aspects.

codecov · 2022-01-18T22:22:44Z

Codecov Report

Merging #1445 (6fa6f16) into main (185a01b) will not change coverage.
The diff coverage is n/a.

❗ Current head 6fa6f16 differs from pull request most recent head 478284c. Consider uploading reports for the commit 478284c to get more accurate results

@@           Coverage Diff           @@
##             main    #1445   +/-   ##
=======================================
  Coverage   91.83%   91.83%           
=======================================
  Files          95       95           
  Lines        4444     4444           
  Branches      427      427           
=======================================
  Hits         4081     4081           
  Misses        273      273           
  Partials       90       90

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 185a01b...478284c. Read the comment docs.

astrojuanlu

Did a joint review with @s-m-e, @jorgepiloto, left a few comments and will follow-up asynchronously

astrojuanlu · 2022-01-28T10:15:53Z

prototypes/brownian_01.ipynb

+    "orbarr: OrbitArray = OrbitArray.from_classical(\n",
+    "    **{key: [item[key] for item in data] for k in data[0].keys()}\n",
+    ")\n",


Let's initially stick with a single attractor and a single common plane.

astrojuanlu · 2022-01-28T10:17:08Z

prototypes/brownian_01.ipynb

+    "orbarr_propagated: OrbitArray = orbarr.propagate(dt1d)\n",
+    "\n",
+    "assert len(orbarr_propagated) == len(orbarr)\n",


This is consistent with the current Orbit.propagate ✔️

And to generate an array of ephemerides we could add an OrbitArray.to_ephem as current consensus on #1364 is, but no need to jump straight to it right now.

astrojuanlu · 2022-01-28T10:22:28Z

prototypes/brownian_01.ipynb

+    "## Propagate array with array of time stamps (hiding the time-deltas)\n",
+    "\n",
+    "# Consistent with Orbit.propagate as it returns another OrbitArray\n",
+    "orbarr_propagated: OrbitArray = orbarr.propagate(times)\n",


This is also consistent with Orbit.propagate ✔️

poliastro/src/poliastro/twobody/orbit.py

Lines 1092 to 1093 in 726c95d

if isinstance(value, time.Time) and not isinstance(value, time.TimeDelta):

time_of_flight = value - self.epoch

Idea after an inquiry from @jorgepiloto: if one wants to propagate an OrbitArray to a single Time, the .propagate method should be smart enough to broadcast the time array. ✔️

astrojuanlu · 2022-01-28T10:26:58Z

prototypes/brownian_01.ipynb

+    "ephem: Ephem = Ephem.from_orbit(orbs[0], **kwargs)\n",
+    "ephemarr: EphemArray = EphemArray.from_orbitarray(orbarr, **kwargs)\n",


Or OrbitArray.to_ephem as discussed above ☝🏽

astrojuanlu · 2022-01-28T10:28:08Z

prototypes/brownian_01.ipynb

+    "\n",
+    "To **simplify** matters, this is not astrodynamics but **Brownian motion**, of sorts. This notebook is looking at the relationships between `Orbit`, `State`, `OrbitArray` and `StateArray` objects and allows some kind of propagation on all of them. \n",
+    "\n",
+    "This notebook uses type annotations and run-time type checks - for testing only. This stuff is not going into `poliastro` ..."


Let's again restrict ourselves to a single kind of state, which according to @s-m-e simplifies the implementation ✔️

astrojuanlu · 2022-01-28T11:04:05Z

prototypes/brownian_01.ipynb

+    "            + '\\n]>'\n",
+    "        )\n",
+    "    \n",
+    "    def __getitem__(self, idx) -> 'Union[StateArray, State]':\n",


astrojuanlu · 2022-01-28T11:08:43Z

prototypes/brownian_01.ipynb

+    "    def __getitem__(self, idx) -> 'Union[OrbitArray, Orbit]':\n",
+    "        target = self._statearray[idx]\n",
+    "        if isinstance(target, State):\n",
+    "            return Orbit(state = target)\n",
+    "        return OrbitArray(statearray = target)\n",


astrojuanlu · 2022-01-28T11:11:21Z

prototypes/brownian_01.ipynb

+    "    def propagate(self, timedelta: TimeDelta, inplace: bool = False) -> 'OrbitArray':\n",
+    "        \n",
+    "        if timedelta.ndim == 0:\n",
+    "            timedelta = np.repeat(timedelta.to_value(u.d), self.size).reshape(self.shape) << u.d\n",


The broadcasting rules mentioned above ☝🏽

astrojuanlu · 2022-01-28T11:17:26Z

prototypes/brownian_01.ipynb

+    "            return Orbit(state = target)\n",
+    "        return OrbitArray(statearray = target)\n",
+    "\n",
+    "    def propagate(self, timedelta: TimeDelta, inplace: bool = False) -> 'OrbitArray':\n",


Same comment as above, let's not close the door on inplace=True but let's not do this in this first iteration so we can focus.

But @s-m-e is veeeeeeeeeeeeeeeeeeeery stubborn and wants to burn hours writing 2x as many tests 😉 so let's give him the chance!

astrojuanlu · 2022-01-28T11:19:16Z

prototypes/brownian_01.ipynb

+   "source": [
+    "## Further thinking\n",
+    "\n",
+    "- [Structured arrays](https://numpy.org/doc/stable/user/basics.rec.html) can actually be n-dimensional. This could be a foundation for `StateArray` classes, opening the back door for alternative libraries which are exposing an [array interface](https://numpy.org/doc/stable/reference/arrays.interface.html). In a bigger picture, the state array could be the place where stuff like `cupy.ndarray` or `dask.array` are transparently accepted.\n",


This is very interesting idea - are structured arrays something that has a future in the NumPy roadmap? Maybe we could rethink some things on top of this.

However, epochs are going to be an issue.

NumPy datatime is known to have lots of issues, or that's what I recall from reading Twitter threads and GitHub comments. I don't have any pointers at hand.

Deleted a highly misleading comment I made: the time of flight is definitely used by the propagation algorithm.

To @s-m-e's point: by passing tof (elapsed time in seconds) to the propagation algorithms, we have already chosen a unit. Therefore, one could take this approach one step further and work with elapsed seconds everywhere internally, to avoid the alleged performance issues of astropy.time and the difficulties to store them.

I guess this approach would work as long as folks don't want to propagate trans-Neptunian objects or the Solar System barycenter with respect to the Milky Way. But, if it doesn't work, it's already slightly wrong in the current implementation.

Important reading: https://docs.astropy.org/en/stable/time/index.html#internal-representation

The Time object maintains an internal representation of time as a pair of double precision numbers expressing Julian days. The sum of the two numbers is the Julian Date for that time relative to the given time scale.

And about what would happen by using a single-float representation:

Users requiring no better than microsecond precision over human time scales (~100 years) can safely ignore the internal representation details and skip this section.

On the other hand, about time scales:

We have shown in the above that the difference between two UTC times is a TimeDelta with a scale of TAI. This is because a UTC time difference cannot be uniquely defined unless the user knows the two times that were differenced (because of leap seconds, a day does not always have 86400 seconds). For all other time scales, the TimeDelta inherits the scale of the first Time object.

Therefore,

There is a way to do it right (in other words, aim for maximum accuracy) with relatively simple tools, but it requires some care.

We might (I haven't deeply analyzed it) have tiny precision bugs in the code at the moment, given that we are not being particularly careful with the TimeDelta scales and just performing a to(u.s) to pass the value to the propagator:

poliastro/src/poliastro/twobody/orbit.py

Lines 1098 to 1099 in 726c95d

cartesian = propagate(self, time_of_flight, method=method, rtol=rtol, **kwargs)

new_epoch = self.epoch + time_of_flight

poliastro/src/poliastro/twobody/propagation.py

Lines 516 to 520 in 726c95d

rr, vv = method(

orbit.attractor.k,

orbit.r,

orbit.v,

time_of_flight.reshape(-1).to(u.s),

On the other hand, we might be getting too deep into this topic, and I don't want this to block our efforts.

The Time object maintains an internal representation of time as a pair of double precision numbers expressing Julian days. The sum of the two numbers is the Julian Date for that time relative to the given time scale.

Then I was not digging deep enough. Interesting. This would actually allow a really clean implementation.

s-m-e added 2 commits January 18, 2022 21:20

initial experiments

7c8a1c2

more notes

ed03c44

s-m-e added 2 commits January 18, 2022 22:47

further thinking

5e6e840

clarification

f0d2513

script version for better GH commenting

478284c

astrojuanlu mentioned this pull request Jan 28, 2022

Easier way of generating a trajectory from an orbit between two instants in time #1364

Closed

astrojuanlu reviewed Jan 28, 2022

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[RFC - do not merge] First `OrbitArray` API prototype #1445

[RFC - do not merge] First `OrbitArray` API prototype #1445

s-m-e commented Jan 18, 2022

s-m-e commented Jan 18, 2022

codecov bot commented Jan 18, 2022 •

edited

astrojuanlu left a comment

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

astrojuanlu Jan 28, 2022

s-m-e Jan 28, 2022

	if isinstance(value, time.Time) and not isinstance(value, time.TimeDelta):
	time_of_flight = value - self.epoch

		"ephem: Ephem = Ephem.from_orbit(orbs[0], **kwargs)\n",
		"ephemarr: EphemArray = EphemArray.from_orbitarray(orbarr, **kwargs)\n",

	cartesian = propagate(self, time_of_flight, method=method, rtol=rtol, **kwargs)
	new_epoch = self.epoch + time_of_flight

	rr, vv = method(
	orbit.attractor.k,
	orbit.r,
	orbit.v,
	time_of_flight.reshape(-1).to(u.s),

[RFC - do not merge] First OrbitArray API prototype #1445

Are you sure you want to change the base?

[RFC - do not merge] First OrbitArray API prototype #1445

Conversation

s-m-e commented Jan 18, 2022

s-m-e commented Jan 18, 2022

codecov bot commented Jan 18, 2022 • edited

Codecov Report

astrojuanlu left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

[RFC - do not merge] First `OrbitArray` API prototype #1445

[RFC - do not merge] First `OrbitArray` API prototype #1445

codecov bot commented Jan 18, 2022 •

edited