Initial OD Improvements

pyproject.toml

@@ -51,6 +51,9 @@ dependencies = [
   "naif-eop-historical",
   "naif-earth-itrf93",
   "timezonefinder",
+  "rebound>=4.6.0",
+  "assist>=1.1.11",
+  "mpcq>=0.4.0; python_version < '3.13'",
 ]
 
 [project.optional-dependencies]

src/adam_core/orbit_determination/__init__.py

@@ -1,5 +1,5 @@
 # flake8: noqa: F401
-from .differential_correction import fit_least_squares
+from .differential_correction import fit_least_squares, iterative_fit
 from .evaluate import OrbitDeterminationObservations, evaluate_orbits
 from .fitted_orbits import FittedOrbitMembers, FittedOrbits, drop_duplicate_orbits
 from .gauss import gaussIOD
@@ -11,4 +11,6 @@
     select_observations,
     sort_by_id_and_time,
 )
+from .native_orbit_fitter import NativeOrbitFitter
+from .orbit_fitter import OrbitFitter
 from .outliers import calculate_max_outliers, remove_lowest_probability_observation

src/adam_core/orbit_determination/differential_correction.py

@@ -16,6 +16,7 @@
 from ..time.time import Timestamp
 from .evaluate import OrbitDeterminationObservations, evaluate_orbits
 from .fitted_orbits import FittedOrbitMembers, FittedOrbits
+from .outliers import calculate_max_outliers, remove_lowest_probability_observation
 
 
 def residual_function(
@@ -204,3 +205,110 @@ def fit_least_squares(
     )
 
     return fitted_orbit, fitted_orbit_members
+
+
+def iterative_fit(
+    orbit: Orbits,
+    observations: OrbitDeterminationObservations,
+    propagator: Propagator,
+    rchi2_threshold: float = 10.0,
+    min_obs: int = 6,
+    min_arc_length: float = 1.0,
+    contamination_percentage: float = 20.0,
+    **kwargs,
+) -> Tuple[FittedOrbits, FittedOrbitMembers]:
+    """
+    Iteratively fit an orbit using least squares with outlier rejection.
+
+    Wraps `fit_least_squares` with an outlier rejection loop: after each fit,
+    if the reduced chi2 exceeds `rchi2_threshold`, the observation with the
+    worst residual is removed and the fit is repeated. This continues until
+    the fit converges, no more outliers are allowed, or arc length / minimum
+    observation constraints would be violated.
+
+    Parameters
+    ----------
+    orbit : `~adam_core.orbits.Orbits` (1)
+        Initial orbit to differentially correct.
+    observations : `~adam_core.orbit_determination.OrbitDeterminationObservations` (N)
+        Observations to fit against.
+    propagator : `~adam_core.propagator.Propagator`
+        Propagator to use to generate ephemeris.
+    rchi2_threshold : float, optional
+        Reduced chi2 threshold below which the fit is considered converged.
+        Default is 10.0.
+    min_obs : int, optional
+        Minimum number of observations required to retain the fit.
+        Default is 6.
+    min_arc_length : float, optional
+        Minimum arc length in days required to retain the fit.
+        Default is 1.0.
+    contamination_percentage : float, optional
+        Maximum percentage of observations that may be rejected as outliers.
+        Range is [0, 100]. Default is 20.0.
+    **kwargs
+        Additional keyword arguments passed to `fit_least_squares` and
+        ultimately to `~scipy.optimize.least_squares`.
+
+    Returns
+    -------
+    fitted_orbit : `~adam_core.orbit_determination.FittedOrbits` (1)
+        Best fitted orbit found.
+    fitted_orbit_members : `~adam_core.orbit_determination.FittedOrbitMembers` (N)
+        Fitted orbit members with residuals and outlier flags.
+    """
+    assert len(orbit) == 1, "Only one orbit can be iteratively fitted"
+
+    num_obs = len(observations)
+    max_outliers = calculate_max_outliers(num_obs, min_obs, contamination_percentage)
+
+    ignore: List[str] = []
+    best_fitted_orbit = None
+    best_fitted_orbit_members = None
+
+    for _ in range(max_outliers + 1):
+        fitted_orbit, fitted_orbit_members = fit_least_squares(
+            orbit,
+            observations,
+            propagator,
+            ignore=ignore if ignore else None,
+            **kwargs,
+        )
+
+        # Track the best fit seen so far (lowest reduced chi2 among successful fits)
+        if best_fitted_orbit is None or (
+            fitted_orbit.success[0].as_py()
+            and fitted_orbit.reduced_chi2[0].as_py()
+            < best_fitted_orbit.reduced_chi2[0].as_py()
+        ):
+            best_fitted_orbit = fitted_orbit
+            best_fitted_orbit_members = fitted_orbit_members
+
+        # Check convergence
+        rchi2 = fitted_orbit.reduced_chi2[0].as_py()
+        if rchi2 is not None and rchi2 <= rchi2_threshold:
+            break
+
+        # Stop if we've already used up all allowed outlier slots
+        if len(ignore) >= max_outliers:
+            break
+
+        # Identify the worst non-outlier observation among the current solution members
+        solution_members = fitted_orbit_members.apply_mask(
+            pc.equal(fitted_orbit_members.outlier, False)
+        )
+        if len(solution_members) == 0:
+            break
+
+        obs_id, remaining_observations = remove_lowest_probability_observation(
+            solution_members, observations
+        )
+
+        # Check that removing this observation still leaves enough arc length
+        arc_length = remaining_observations.coordinates.time.mjd().to_numpy()
+        if len(arc_length) < min_obs or (arc_length.max() - arc_length.min()) < min_arc_length:
+            break
+
+        ignore.append(obs_id)
+
+    return best_fitted_orbit, best_fitted_orbit_members

src/adam_core/orbit_determination/iod.py

@@ -134,13 +134,13 @@ def select_observations(
     times = observations.coordinates.time.mjd().to_numpy(zero_copy_only=False)
 
     if method == "first+middle+last":
-        selected_times = np.percentile(times, [0, 50, 100], interpolation="nearest")
+        selected_times = np.percentile(times, [0, 50, 100], method="nearest")
         selected_index = np.intersect1d(times, selected_times, return_indices=True)[1]
         selected_index = np.array([selected_index])
 
     elif method == "thirds":
         selected_times = np.percentile(
-            times, [1 / 6 * 100, 50, 5 / 6 * 100], interpolation="nearest"
+            times, [1 / 6 * 100, 50, 5 / 6 * 100], method="nearest"
         )
         selected_index = np.intersect1d(times, selected_times, return_indices=True)[1]
         selected_index = np.array([selected_index])
@@ -386,8 +386,6 @@ def iod(
     if len(observations) == 0:
         processable = False
 
-    obs_ids_all = observations.id.to_numpy(zero_copy_only=False)
-    coords_all = observations.coordinates
     observers = observations.observers
 
     observations = observations.sort_by(
@@ -401,6 +399,8 @@ def iod(
         ["coordinates.time.days", "coordinates.time.nanos", "coordinates.origin.code"]
     )
 
+    obs_ids_all = observations.id.to_numpy(zero_copy_only=False)
+    coords_all = observations.coordinates
     coords_obs_all = observers.coordinates.r
     times_all = coords_all.time.mjd().to_numpy(zero_copy_only=False)
 

src/adam_core/orbit_determination/native_orbit_fitter.py

@@ -0,0 +1,145 @@
+import logging
+from typing import Literal, Tuple, Type
+
+import pyarrow as pa
+
+from ..propagator.propagator import Propagator
+from .differential_correction import iterative_fit
+from .evaluate import OrbitDeterminationObservations
+from .fitted_orbits import FittedOrbitMembers, FittedOrbits
+from .iod import iod
+from .orbit_fitter import OrbitFitter
+
+logger = logging.getLogger(__name__)
+
+
+class NativeOrbitFitter(OrbitFitter):
+    """
+    Orbit fitter using adam_core's native Gauss IOD and iterative least-squares DC.
+
+    This fitter is a thin wrapper that chains:
+    1. `iod()` — Gauss initial orbit determination (Milani 2008)
+    2. `iterative_fit()` — scipy least-squares differential correction with
+       outlier rejection
+
+    Parameters
+    ----------
+    propagator_class : Type[Propagator]
+        Propagator *class* (not instance) used during IOD ephemeris evaluation.
+    propagator_kwargs : dict, optional
+        Keyword arguments forwarded to the propagator constructor / IOD call.
+    min_obs : int, optional
+        Minimum number of observations required for a valid fit.  Default 6.
+    min_arc_length : float, optional
+        Minimum arc length in days required to retain a fit.  Default 1.0.
+    contamination_percentage : float, optional
+        Maximum percentage of observations that may be rejected as outliers
+        across the full OD pipeline.  Default 20.0.
+    rchi2_threshold : float, optional
+        Reduced chi2 convergence threshold for differential correction.
+        Default 10.0.
+    iod_rchi2_threshold : float, optional
+        Reduced chi2 threshold used during IOD to filter candidate orbits.
+        Default 200.0.
+    observation_selection_method : str, optional
+        Strategy for selecting observation triplets in IOD.  One of
+        ``"combinations"``, ``"first+middle+last"``, ``"thirds"``.
+        Default ``"combinations"``.
+    """
+
+    def __init__(
+        self,
+        propagator_class: Type[Propagator],
+        propagator_kwargs: dict = {},
+        min_obs: int = 6,
+        min_arc_length: float = 1.0,
+        contamination_percentage: float = 20.0,
+        rchi2_threshold: float = 10.0,
+        iod_rchi2_threshold: float = 200.0,
+        observation_selection_method: Literal[
+            "combinations", "first+middle+last", "thirds"
+        ] = "combinations",
+    ) -> None:
+        self.propagator_class = propagator_class
+        self.propagator_kwargs = propagator_kwargs
+        self.min_obs = min_obs
+        self.min_arc_length = min_arc_length
+        self.contamination_percentage = contamination_percentage
+        self.rchi2_threshold = rchi2_threshold
+        self.iod_rchi2_threshold = iod_rchi2_threshold
+        self.observation_selection_method = observation_selection_method
+
+    def __getstate__(self) -> dict:
+        return self.__dict__.copy()
+
+    def __setstate__(self, state: dict) -> None:
+        self.__dict__.update(state)
+
+    def initial_fit(
+        self,
+        object_id: str | pa.LargeStringScalar,
+        observations: OrbitDeterminationObservations,
+    ) -> Tuple[FittedOrbits, FittedOrbitMembers]:
+        """Run Gauss IOD on the observations.
+
+        Parameters
+        ----------
+        object_id : str | pa.LargeStringScalar
+            Object identifier for output tables.
+        observations : OrbitDeterminationObservations
+            Observations to fit, assumed to belong to a single object.
+
+        Returns
+        -------
+        fitted_orbit : FittedOrbits
+            Best IOD orbit(s) found (may be empty if IOD fails).
+        fitted_orbit_members : FittedOrbitMembers
+            Observations with solution/outlier flags from IOD.
+        """
+        fitted_orbits, fitted_orbit_members = iod(
+            observations,
+            self.propagator_class,
+            min_obs=self.min_obs,
+            min_arc_length=self.min_arc_length,
+            contamination_percentage=self.contamination_percentage,
+            rchi2_threshold=self.iod_rchi2_threshold,
+            observation_selection_method=self.observation_selection_method,
+            propagator_kwargs=self.propagator_kwargs,
+        )
+        return fitted_orbits, fitted_orbit_members
+
+    def refine_fit(
+        self,
+        fitted_orbit: FittedOrbits,
+        observations: OrbitDeterminationObservations,
+        propagator: Propagator,
+    ) -> Tuple[FittedOrbits, FittedOrbitMembers]:
+        """Refine an IOD orbit via iterative differential correction.
+
+        Parameters
+        ----------
+        fitted_orbit : FittedOrbits (1)
+            Orbit to refine, typically from `initial_fit`.
+        observations : OrbitDeterminationObservations
+            Observations to fit against.
+        propagator : Propagator
+            Propagator instance used during DC ephemeris evaluation.
+
+        Returns
+        -------
+        fitted_orbit : FittedOrbits (1)
+            DC-refined orbit with covariance and quality statistics.
+        fitted_orbit_members : FittedOrbitMembers (N)
+            Observations with residuals and outlier/solution flags.
+        """
+        assert len(fitted_orbit) == 1, "refine_fit expects exactly one orbit"
+        orbit = fitted_orbit.to_orbits()
+        return iterative_fit(
+            orbit,
+            observations,
+            propagator,
+            rchi2_threshold=self.rchi2_threshold,
+            min_obs=self.min_obs,
+            min_arc_length=self.min_arc_length,
+            contamination_percentage=self.contamination_percentage,
+        )

src/adam_core/orbit_determination/od.py

@@ -1,6 +1,7 @@
 import logging
 import multiprocessing as mp
 import time
+import warnings
 from typing import Literal, Optional, Tuple, Type, Union
 
 import numpy as np
@@ -587,6 +588,11 @@ def differential_correction(
     """
     Differentially correct (via finite/central differencing).
 
+    .. deprecated::
+        `differential_correction` is deprecated and will be removed in a future release.
+        Use `adam_core.orbit_determination.iterative_fit` instead, which provides the same
+        outlier-rejection loop via `scipy.optimize.least_squares`.
+
     Parameters
     ----------
     chunk_size : int, optional
@@ -597,6 +603,12 @@ def differential_correction(
         Which parallelization backend to use {'ray', 'mp', 'cf'}. Defaults to using Python's concurrent.futures
         module ('cf').
     """
+    warnings.warn(
+        "differential_correction is deprecated and will be removed in a future release. "
+        "Use adam_core.orbit_determination.iterative_fit instead.",
+        DeprecationWarning,
+        stacklevel=2,
+    )
     time_start = time.perf_counter()
     logger.info("Running differential correction...")
     if isinstance(orbits, ray.ObjectRef):