Add link_test_orbit

[spenczar]

No description provided.

[moeyensj]

A couple of post-vacation interface thoughts: What do you think about decoupling ObservationSource and how observations are gathered relative to a test orbit? We can make ObservationSource the class that deals with where the data lives and how we query/get it. You can think of adding an AIMSSource or ParquetSource(?) as an example. The interface could be

observations = ObservationSource.get_observations(mjd_min=None, mjd_max=None, obscodes=None)

The job of how observations are gathered relative to a test orbit feels like something that should be delegated to the test orbit class (even more so with covariance matrices which we might create as a function of the state vector and covariances that get propagated with the state vector). Here we could add a function to the test orbit class:

gathered_observations = test_orbit.gather_observations(observations, method="cell")

gathered_observations = test_orbit.gather_observations(observations, method="covariance")

[moeyensj]

I suppose one could get real fancy and abstract the method:

class PhaseSpaceMapping(ABC):
    def filter_observations(self):
        pass

class CovarianceMapping(PhaseSpaceMapping):
    def __init__(self, matrix = np.random.random((6, 6)):
        pass

class CellMapping(PhaseSpaceMapping):
    def __init__(self, radius = 1):
        pass

Whats tricky here is that the CovarianceMapping will be time-dependent (it might shrink or grow with the test orbit) so maybe this level of abstract over complicates things.

[spenczar]

I think the current structure is preferable. In the current scheme, non-test-orbit-based filtering is pretty simple. See, for example, the "StaticObservationSource" implementation, which always yields the same values.

We're also free to write chains of ObservationSources: one that filters by time, a second that filters to only observations that have not been linked in previous THOR runs, a third that filters to objects within the covariance region of a test orbit.

If these filters were applied as a method on a test orbit, those would be much harder to express.

Maybe we should rename that abstraction to "ObservationFilter" to express this concept a bit better. The root source is something like AIMS or Parquet or whatever, and that's a fairly different thing. And indeed the current implementations expect a pile of observations as input in order to be constructed!

[moeyensj]

I guess what isn't fully obvious to me is how we use the current structure for propagated covariance matrices assumed with each test orbit? The way I'm thinking about is that two identical state vectors with different covariance matrices are in fact different test orbits (different volumes of phase space), so they should be coupled together somehow. Maybe we add a CovarianceSource that accepts a covariance matrix for each epoch at which observations are defined? The CovarianceSource then doesn't care where the matrix came from as long as it gets a set of spherical coordinates with covariances defined at the predicted position of the test orbit as viewed from the observer at a specific point in time.

[moeyensj]

Ah actually, I see how we would do it. We generate ephemerides inside FixedRadiusObservationsSource (https://github.com/moeyensj/thor/blob/v2.0/thor/observation_source.py#L65), we could do the exact same with covariance=True passed to the ephemeris generator.

This definitely points towards us needing to cache or store some of these propagations / ephemerides since it will get expensive otherwise with the number of recomputations.

[spenczar]

Yes, that's exactly how I'd think of it - an additional Source. I think naming them Filters would make it much clearer.

[moeyensj]

Agreed on naming them Filters for clarity.

[moeyensj]

@spenczar Just FYI, I've rebased this branch to main and set the base to main. I'm now working on an integration test and will commit it here.

[moeyensj]

@spenczar and @akoumjian, I've made some pretty substantial updates to some of the underlying functions to get this running.

Here is the current high-level interface (using adam_core's test data):

Load some test data:

from adam_core.utils.helpers import make_observations
from adam_core.utils.helpers import make_real_orbits

exposures, detections, associations = make_observations()
orbits = make_real_orbits()

Lets link with a single test orbit (note this function is just going to return the transformed detections):

from thor.observations import Observations
from thor.observations.filters import TestOrbitRadiusObservationFilter
from thor import TestOrbit, link_test_orbit 

observations = Observations.from_detections_and_exposures(detections, exposures)
test_orbit = TestOrbit.from_orbits(orbits[10])

transformed_observations = link_test_orbit(
    test_orbit, 
    observations, 
    filters=[
        TestOrbitRadiusObservationFilter(radius=50)
    ]
)

Plot the resulting transformed detections:

import matplotlib.pyplot as plt

df = transformed_observations.to_dataframe()
fig, ax = plt.subplots(1, 1, dpi=200)
ax.set_aspect("equal")
ax.scatter(df["coordinates.theta_x"], df["coordinates.theta_y"], s=0.1)
ax.set_xlabel(r"$\theta_X$ [deg]")
ax.set_ylabel(r"$\theta_Y$ [deg]")

What is coincidentally really cool here is that this dataset has simulated ephemerides from two distinct observatories and spans a total of 60 days. We have some nearly linear motion over that time span!

A bunch of things were changed:

• TestOrbit now has a .generate_ephemeris_from_observations function which will calculate the ephemerides for the test orbit for each unique state (state = unique combination of time and observatory code). We will need to come up with some sort of binning for something like LSST but that is fine. What is nice here is that those test_orbit states are stored by the class and only recalculated when a change in observation IDs are detected. (The observer positions are also stored here).
• Observations were made into a quivr table rather than a class (this mimics a bit what we are doing with adam_etl). Its just a lot easier to deal with a table than a linkage at that level
• Spencer introduced a really cool idea for ObservationFilters. These have been kept largely the same with one major difference: test_orbits are no longer required by the constructor and instead are passed in the apply function.
• I modified some of the logic for the radius filter to use numpy masking. This should be a bit faster than return detection tables created during each iteration.

To get this working the following PRs in adam_core need to be merged:
B612-Asteroid-Institute/adam_core#63
B612-Asteroid-Institute/adam_core#64
B612-Asteroid-Institute/adam_core#65

[moeyensj]

I still need to define an integration test but the adam_core stuff needs to be merged / reviewed first. Early reviews from both of you would be super welcome here.